!=========================================================================== ! subroutine z_der1_co & (f,df,dx,dy,dz,vel,upwind,orderi,orderb, & biL, biR, bjL, bjR, bkL, bkR, & fiL, fiR, fjL, fjR, fkL, fkR, & iiL, iiR, ijL, ijR, ikL, ikR, & viL, viR, vjL, vjR, vkL, vkR, & dfiL,dfiR,dfjL,dfjR,dfkL,dfkR, & iperx,ipery,iperz,ierror,error) ! include 'der1_include.h' ! call der1_pre( ierror, error, upwind, & xbegin,xend,ybegin,yend,zbegin,zend, & bnd1A, bnd3A, bnd3B, bnd4C, bnd6D, & iperx,ipery,iperz, orderi, orderb, & biL, biR, bjL, bjR, bkL, bkR, & fiL, fiR, fjL, fjR, fkL, fkR, & iiL, iiR, ijL, ijR, ikL, ikR, & piL, piR, pjL, pjR, pkL, pkR, & viL, viR, vjL, vjR, vkL, vkR, & dfiL,dfiR,dfjL,dfjR,dfkL,dfkR, & nproLx,nproRx,nproLy,nproRy,nproLz,nproRz, & proL,proR,promax,dx,dy,dz,rdx,rdy,rdz) ! if( 2*promax+1 .gt. pkR-pkL+1) then ierror = -100 error = 'DER_1: Stencil is wider than the f domain' return endif ! ! Zero-out first derivative ! do k = ikL, ikR do j = ijL, ijR do i = iiL, iiR df(i,j,k) = 0.0 end do end do end do ! !======================================================================== ! CENTERED DIFFERENCE OPERATORS ! if ( upwind .eq. 0 ) then !----------------------------- ! 2nd-order explicit: (1-2E-1) ! if ( orderi .eq. 2 ) then ! ae = (1.d0/ 2.d0) * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR df(i,j,k) = ae *( f(i,j,k+1)-f(i,j,k-1) ) end do end do end do ! ! Boundary nodes: ! Left side: ! if ( nproLz .eq. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd1A(1)*f(i,j,zbegin-1) + & bnd1A(2)*f(i,j,zbegin+0))*rdz end do end do end if ! ! Right side: ! if ( nproRz .eq. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd1A(1)*f(i,j,zend+1) + & bnd1A(2)*f(i,j,zend+0))*rdz end do end do end if ! end if ! !----------------------------- ! 4th order explicit: (3,3-4E-3,3) ! if ( orderi .eq. 4 ) then ! ae = 2.d0/ 3.d0 * rdz be = -1.d0/12.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR df(i,j,k) = ae *( f(i,j,k+1)-f(i,j,k-1) ) & + be *( f(i,j,k+2)-f(i,j,k-2) ) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd3B(1)*f(i,j,zbegin-2) + & bnd3B(2)*f(i,j,zbegin-1) + & bnd3B(3)*f(i,j,zbegin+0) + & bnd3B(4)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .eq. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd3A(1)*f(i,j,zbegin-2) + & bnd3A(2)*f(i,j,zbegin-1) + & bnd3A(3)*f(i,j,zbegin+0) + & bnd3A(4)*f(i,j,zbegin+1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd3B(1)*f(i,j,zend+2) + & bnd3B(2)*f(i,j,zend+1) + & bnd3B(3)*f(i,j,zend+0) + & bnd3B(4)*f(i,j,zend-1))*rdz end do end do end if if ( nproRz .eq. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd3A(1)*f(i,j,zend+2) + & bnd3A(2)*f(i,j,zend+1) + & bnd3A(3)*f(i,j,zend+0) + & bnd3A(4)*f(i,j,zend-1))*rdz end do end do end if ! end if ! !----------------------------- ! 6th order explicit: (3,3,4-6E-4,3,3) ! ! if ( orderi .eq. 6 ) then ! ae = 3.d0/ 4.d0 * rdz be = -3.d0/20.d0 * rdz ce = 1.d0/60.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR df(i,j,k) = ae *( f(i,j,k+1)-f(i,j,k-1) ) & + be *( f(i,j,k+2)-f(i,j,k-2) ) & + ce *( f(i,j,k+3)-f(i,j,k-3) ) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd4C(1)*f(i,j,zbegin-3) + & bnd4C(2)*f(i,j,zbegin-2) + & bnd4C(3)*f(i,j,zbegin-1) + & bnd4C(4)*f(i,j,zbegin+0) + & bnd4C(5)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd3B(1)*f(i,j,zbegin-3) + & bnd3B(2)*f(i,j,zbegin-2) + & bnd3B(3)*f(i,j,zbegin-1) + & bnd3B(4)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .eq. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3A(1)*f(i,j,zbegin-3) + & bnd3A(2)*f(i,j,zbegin-2) + & bnd3A(3)*f(i,j,zbegin-1) + & bnd3A(4)*f(i,j,zbegin+0))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd4C(1)*f(i,j,zend+3) + & bnd4C(2)*f(i,j,zend+2) + & bnd4C(3)*f(i,j,zend+1) + & bnd4C(4)*f(i,j,zend+0) + & bnd4C(5)*f(i,j,zend-1))*rdz end do end do end if if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd3B(1)*f(i,j,zend+3) + & bnd3B(2)*f(i,j,zend+2) + & bnd3B(3)*f(i,j,zend+1) + & bnd3B(4)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .eq. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3A(1)*f(i,j,zend+3) + & bnd3A(2)*f(i,j,zend+2) + & bnd3A(3)*f(i,j,zend+1) + & bnd3A(4)*f(i,j,zend+0))*rdz end do end do end if ! end if !----------------------------- ! 8th order explicit: (3,3,4,6-8E-6,4,3,3) ! ! if ( orderi .eq. 8 ) then ! ae = 4.d0/ 5.d0 * rdz be = -1.d0/ 5.d0 * rdz ce = 4.d0/105.d0 * rdz de = -1.d0/280.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR df(i,j,k) = ae *( f(i,j,k+1)-f(i,j,k-1) ) & + be *( f(i,j,k+2)-f(i,j,k-2) ) & + ce *( f(i,j,k+3)-f(i,j,k-3) ) & + de *( f(i,j,k+4)-f(i,j,k-4) ) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd6D(1)*f(i,j,zbegin-4) + & bnd6D(2)*f(i,j,zbegin-3) + & bnd6D(3)*f(i,j,zbegin-2) + & bnd6D(4)*f(i,j,zbegin-1) + & bnd6D(5)*f(i,j,zbegin+0) + & bnd6D(6)*f(i,j,zbegin+1) + & bnd6D(7)*f(i,j,zbegin+2))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd4C(1)*f(i,j,zbegin-4) + & bnd4C(2)*f(i,j,zbegin-3) + & bnd4C(3)*f(i,j,zbegin-2) + & bnd4C(4)*f(i,j,zbegin-1) + & bnd4C(5)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3B(1)*f(i,j,zbegin-4) + & bnd3B(2)*f(i,j,zbegin-3) + & bnd3B(3)*f(i,j,zbegin-2) + & bnd3B(4)*f(i,j,zbegin-1))*rdz end do end do end if ! if ( nproLz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-4) = (bnd3A(1)*f(i,j,zbegin-4) + & bnd3A(2)*f(i,j,zbegin-3) + & bnd3A(3)*f(i,j,zbegin-2) + & bnd3A(4)*f(i,j,zbegin-1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd6D(1)*f(i,j,zend+4) + & bnd6D(2)*f(i,j,zend+3) + & bnd6D(3)*f(i,j,zend+2) + & bnd6D(4)*f(i,j,zend+1) + & bnd6D(5)*f(i,j,zend+0) + & bnd6D(6)*f(i,j,zend-1) + & bnd6D(7)*f(i,j,zend-2))*rdz end do end do end if ! if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd4C(1)*f(i,j,zend+4) + & bnd4C(2)*f(i,j,zend+3) + & bnd4C(3)*f(i,j,zend+2) + & bnd4C(4)*f(i,j,zend+1) + & bnd4C(5)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3B(1)*f(i,j,zend+4) + & bnd3B(2)*f(i,j,zend+3) + & bnd3B(3)*f(i,j,zend+2) + & bnd3B(4)*f(i,j,zend+1))*rdz end do end do end if if ( nproRz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+4) = -(bnd3A(1)*f(i,j,zend+4) + & bnd3A(2)*f(i,j,zend+3) + & bnd3A(3)*f(i,j,zend+2) + & bnd3A(4)*f(i,j,zend+1))*rdz end do end do end if ! end if end if ! ! End Centered Difference Operators ! !======================================================================== ! UPWIND/DOWNWIND DIFFERENCE OPERATORS ! if ( abs(upwind) .eq. 1 ) then !----------------------------- ! 3rd order upwind: (3,3-3U-3,3) if ( orderi .eq. 3 ) then ! aL = -1.d0/ 3.d0 * rdz Up = -1.d0/ 2.d0 * rdz aR = 1.d0/ 1.d0 * rdz bR = -1.d0/ 6.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (aL * f(i,j,k-sign*1) & + Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd3B(1)*f(i,j,zbegin-2) + & bnd3B(2)*f(i,j,zbegin-1) + & bnd3B(3)*f(i,j,zbegin+0) + & bnd3B(4)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .eq. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd3A(1)*f(i,j,zbegin-2) + & bnd3A(2)*f(i,j,zbegin-1) + & bnd3A(3)*f(i,j,zbegin+0) + & bnd3A(4)*f(i,j,zbegin+1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd3B(1)*f(i,j,zend+2) + & bnd3B(2)*f(i,j,zend+1) + & bnd3B(3)*f(i,j,zend+0) + & bnd3B(4)*f(i,j,zend-1))*rdz end do end do end if if ( nproRz .eq. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd3A(1)*f(i,j,zend+2) + & bnd3A(2)*f(i,j,zend+1) + & bnd3A(3)*f(i,j,zend+0) + & bnd3A(4)*f(i,j,zend-1))*rdz end do end do end if ! end if ! !---------------------------------------------------------- ! 4th order upwind: (3,3,4-4U-4,3,3) ! if ( orderi .eq. 4 ) then ! aL = -1.d0/ 4.d0 * rdz Up = -5.d0/ 6.d0 * rdz aR = 3.d0/ 2.d0 * rdz bR = -1.d0/ 2.d0 * rdz cR = 1.d0/12.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (aL * f(i,j,k-sign*1) & + Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2) & + cR * f(i,j,k+sign*3)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd4C(1)*f(i,j,zbegin-3) + & bnd4C(2)*f(i,j,zbegin-2) + & bnd4C(3)*f(i,j,zbegin-1) + & bnd4C(4)*f(i,j,zbegin+0) + & bnd4C(5)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd3B(1)*f(i,j,zbegin-3) + & bnd3B(2)*f(i,j,zbegin-2) + & bnd3B(3)*f(i,j,zbegin-1) + & bnd3B(4)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .eq. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3A(1)*f(i,j,zbegin-3) + & bnd3A(2)*f(i,j,zbegin-2) + & bnd3A(3)*f(i,j,zbegin-1) + & bnd3A(4)*f(i,j,zbegin+0))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd4C(1)*f(i,j,zend+3) + & bnd4C(2)*f(i,j,zend+2) + & bnd4C(3)*f(i,j,zend+1) + & bnd4C(4)*f(i,j,zend+0) + & bnd4C(5)*f(i,j,zend-1))*rdz end do end do end if if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd3B(1)*f(i,j,zend+3) + & bnd3B(2)*f(i,j,zend+2) + & bnd3B(3)*f(i,j,zend+1) + & bnd3B(4)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .eq. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3A(1)*f(i,j,zend+3) + & bnd3A(2)*f(i,j,zend+2) + & bnd3A(3)*f(i,j,zend+1) + & bnd3A(4)*f(i,j,zend+0))*rdz end do end do end if ! end if ! !---------------------------------------------------------- ! 5th-order upwind: (3,3,4-5U-4,3,3) ! if ( orderi .eq. 5 ) then ! bL = 1.d0/20.d0 * rdz aL = -1.d0/ 2.d0 * rdz Up = -1.d0/ 3.d0 * rdz aR = 1.d0/ 1.d0 * rdz bR = -1.d0/ 4.d0 * rdz cR = 1.d0/30.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (bL * f(i,j,k-sign*2) & + aL * f(i,j,k-sign*1) & + Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2) & + cR * f(i,j,k+sign*3)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd4C(1)*f(i,j,zbegin-3) + & bnd4C(2)*f(i,j,zbegin-2) + & bnd4C(3)*f(i,j,zbegin-1) + & bnd4C(4)*f(i,j,zbegin+0) + & bnd4C(5)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd3B(1)*f(i,j,zbegin-3) + & bnd3B(2)*f(i,j,zbegin-2) + & bnd3B(3)*f(i,j,zbegin-1) + & bnd3B(4)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .eq. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3A(1)*f(i,j,zbegin-3) + & bnd3A(2)*f(i,j,zbegin-2) + & bnd3A(3)*f(i,j,zbegin-1) + & bnd3A(4)*f(i,j,zbegin+0))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd4C(1)*f(i,j,zend+3) + & bnd4C(2)*f(i,j,zend+2) + & bnd4C(3)*f(i,j,zend+1) + & bnd4C(4)*f(i,j,zend+0) + & bnd4C(5)*f(i,j,zend-1))*rdz end do end do end if if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd3B(1)*f(i,j,zend+3) + & bnd3B(2)*f(i,j,zend+2) + & bnd3B(3)*f(i,j,zend+1) + & bnd3B(4)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .eq. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3A(1)*f(i,j,zend+3) + & bnd3A(2)*f(i,j,zend+2) + & bnd3A(3)*f(i,j,zend+1) + & bnd3A(4)*f(i,j,zend+0))*rdz end do end do end if ! end if ! !---------------------------------------------------------- ! 6th-order upwind: (3,3,4,6-6U-6,4,3,3) ! if ( orderi .eq. 6 ) then ! bL = 1.d0/30.d0 * rdz aL = - 2.d0/ 5.d0 * rdz Up = - 7.d0/12.d0 * rdz aR = 4.d0/ 3.d0 * rdz bR = -1.d0/ 2.d0 * rdz cR = 2.d0/15.d0 * rdz dR = -1.d0/60.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (bL * f(i,j,k-sign*2) & + aL * f(i,j,k-sign*1) & + Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2) & + cR * f(i,j,k+sign*3) & + dR * f(i,j,k+sign*4)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd6D(1)*f(i,j,zbegin-4) + & bnd6D(2)*f(i,j,zbegin-3) + & bnd6D(3)*f(i,j,zbegin-2) + & bnd6D(4)*f(i,j,zbegin-1) + & bnd6D(5)*f(i,j,zbegin+0) + & bnd6D(6)*f(i,j,zbegin+1) + & bnd6D(7)*f(i,j,zbegin+2))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd4C(1)*f(i,j,zbegin-4) + & bnd4C(2)*f(i,j,zbegin-3) + & bnd4C(3)*f(i,j,zbegin-2) + & bnd4C(4)*f(i,j,zbegin-1) + & bnd4C(5)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3B(1)*f(i,j,zbegin-4) + & bnd3B(2)*f(i,j,zbegin-3) + & bnd3B(3)*f(i,j,zbegin-2) + & bnd3B(4)*f(i,j,zbegin-1))*rdz end do end do end if ! if ( nproLz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-4) = (bnd3A(1)*f(i,j,zbegin-4) + & bnd3A(2)*f(i,j,zbegin-3) + & bnd3A(3)*f(i,j,zbegin-2) + & bnd3A(4)*f(i,j,zbegin-1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd6D(1)*f(i,j,zend+4) + & bnd6D(2)*f(i,j,zend+3) + & bnd6D(3)*f(i,j,zend+2) + & bnd6D(4)*f(i,j,zend+1) + & bnd6D(5)*f(i,j,zend+0) + & bnd6D(6)*f(i,j,zend-1) + & bnd6D(7)*f(i,j,zend-2))*rdz end do end do end if ! if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd4C(1)*f(i,j,zend+4) + & bnd4C(2)*f(i,j,zend+3) + & bnd4C(3)*f(i,j,zend+2) + & bnd4C(4)*f(i,j,zend+1) + & bnd4C(5)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3B(1)*f(i,j,zend+4) + & bnd3B(2)*f(i,j,zend+3) + & bnd3B(3)*f(i,j,zend+2) + & bnd3B(4)*f(i,j,zend+1))*rdz end do end do end if if ( nproRz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+4) = -(bnd3A(1)*f(i,j,zend+4) + & bnd3A(2)*f(i,j,zend+3) + & bnd3A(3)*f(i,j,zend+2) + & bnd3A(4)*f(i,j,zend+1))*rdz end do end do end if ! end if !---------------------------------------------------------- ! 7th-order upwind: (3,3,4,6-7U-6,4,3,3) ! if ( orderi .eq. 7 ) then ! cL = -1.d0/105.d0 * rdz bL = 1.d0/ 10.d0 * rdz aL = -3.d0/ 5.d0 * rdz Up = -1.d0/ 4.d0 * rdz aR = 1.d0/ 1.d0 * rdz bR = -3.d0/ 10.d0 * rdz cR = 1.d0/ 15.d0 * rdz dR = -1.d0/140.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (cL * f(i,j,k-sign*3) & + bL * f(i,j,k-sign*2) & + aL * f(i,j,k-sign*1) & + Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2) & + cR * f(i,j,k+sign*3) & + dR * f(i,j,k+sign*4)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd6D(1)*f(i,j,zbegin-4) + & bnd6D(2)*f(i,j,zbegin-3) + & bnd6D(3)*f(i,j,zbegin-2) + & bnd6D(4)*f(i,j,zbegin-1) + & bnd6D(5)*f(i,j,zbegin+0) + & bnd6D(6)*f(i,j,zbegin+1) + & bnd6D(7)*f(i,j,zbegin+2))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd4C(1)*f(i,j,zbegin-4) + & bnd4C(2)*f(i,j,zbegin-3) + & bnd4C(3)*f(i,j,zbegin-2) + & bnd4C(4)*f(i,j,zbegin-1) + & bnd4C(5)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3B(1)*f(i,j,zbegin-4) + & bnd3B(2)*f(i,j,zbegin-3) + & bnd3B(3)*f(i,j,zbegin-2) + & bnd3B(4)*f(i,j,zbegin-1))*rdz end do end do end if ! if ( nproLz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-4) = (bnd3A(1)*f(i,j,zbegin-4) + & bnd3A(2)*f(i,j,zbegin-3) + & bnd3A(3)*f(i,j,zbegin-2) + & bnd3A(4)*f(i,j,zbegin-1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd6D(1)*f(i,j,zend+4) + & bnd6D(2)*f(i,j,zend+3) + & bnd6D(3)*f(i,j,zend+2) + & bnd6D(4)*f(i,j,zend+1) + & bnd6D(5)*f(i,j,zend+0) + & bnd6D(6)*f(i,j,zend-1) + & bnd6D(7)*f(i,j,zend-2))*rdz end do end do end if ! if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd4C(1)*f(i,j,zend+4) + & bnd4C(2)*f(i,j,zend+3) + & bnd4C(3)*f(i,j,zend+2) + & bnd4C(4)*f(i,j,zend+1) + & bnd4C(5)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3B(1)*f(i,j,zend+4) + & bnd3B(2)*f(i,j,zend+3) + & bnd3B(3)*f(i,j,zend+2) + & bnd3B(4)*f(i,j,zend+1))*rdz end do end do end if if ( nproRz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+4) = -(bnd3A(1)*f(i,j,zend+4) + & bnd3A(2)*f(i,j,zend+3) + & bnd3A(3)*f(i,j,zend+2) + & bnd3A(4)*f(i,j,zend+1))*rdz end do end do end if ! end if ! end if ! ! End Simple Upwind Difference Operators ! !----------------------------- if ( abs(upwind) .eq. 2 ) then ! ! 4th-order super-upwind: (3,3,4,4-4UU-4,4,3,3) ! if ( orderi .eq. 4 ) then ! Up = -25.d0/12.d0 * rdz aR = 4.d0/ 1.d0 * rdz bR = -3.d0/ 1.d0 * rdz cR = 4.d0/ 3.d0 * rdz dR = -1.d0/ 4.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2) & + cR * f(i,j,k+sign*3) & + dR * f(i,j,k+sign*4)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd4C(1)*f(i,j,zbegin-3) + & bnd4C(2)*f(i,j,zbegin-2) + & bnd4C(3)*f(i,j,zbegin-1) + & bnd4C(4)*f(i,j,zbegin-0) + & bnd4C(5)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd4C(1)*f(i,j,zbegin-4) + & bnd4C(2)*f(i,j,zbegin-3) + & bnd4C(3)*f(i,j,zbegin-2) + & bnd4C(4)*f(i,j,zbegin-1) + & bnd4C(5)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3B(1)*f(i,j,zbegin-4) + & bnd3B(2)*f(i,j,zbegin-3) + & bnd3B(3)*f(i,j,zbegin-2) + & bnd3B(4)*f(i,j,zbegin-1))*rdz end do end do end if ! if ( nproLz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-4) = (bnd3A(1)*f(i,j,zbegin-4) + & bnd3A(2)*f(i,j,zbegin-3) + & bnd3A(3)*f(i,j,zbegin-2) + & bnd3A(4)*f(i,j,zbegin-1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd4C(1)*f(i,j,zend+3) + & bnd4C(2)*f(i,j,zend+2) + & bnd4C(3)*f(i,j,zend+1) + & bnd4C(4)*f(i,j,zend+0) + & bnd4C(5)*f(i,j,zend-1))*rdz end do end do end if ! if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd4C(1)*f(i,j,zend+4) + & bnd4C(2)*f(i,j,zend+3) + & bnd4C(3)*f(i,j,zend+2) + & bnd4C(4)*f(i,j,zend+1) + & bnd4C(5)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3B(1)*f(i,j,zend+4) + & bnd3B(2)*f(i,j,zend+3) + & bnd3B(3)*f(i,j,zend+2) + & bnd3B(4)*f(i,j,zend+1))*rdz end do end do end if if ( nproRz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+4) = -(bnd3A(1)*f(i,j,zend+4) + & bnd3A(2)*f(i,j,zend+3) + & bnd3A(3)*f(i,j,zend+2) + & bnd3A(4)*f(i,j,zend+1))*rdz end do end do end if ! end if ! ! !----------------------------- ! 5th-order super-upwind: (3,3,4,4-5UU-4,4,3,3) ! if ( orderi .eq. 5 ) then ! aL = - 1.d0/ 5.d0 * rdz Up = -13.d0/12.d0 * rdz aR = 2.d0/ 1.d0 * rdz bR = -1.d0/ 1.d0 * rdz cR = 1.d0/ 3.d0 * rdz dR = -1.d0/20.d0 * rdz ! ! Internal nodes: ! do k = zbegin, zend do j = ijL, ijR do i = iiL, iiR ! sign = vel(i,j,k)*upwind/( abs(vel(i,j,k) ) ) if( vel(i,j,k) .eq. 0.d0) sign = 1 ! df(i,j,k) = sign * (aL * f(i,j,k-sign*1) & + Up * f(i,j,k ) & + aR * f(i,j,k+sign*1) & + bR * f(i,j,k+sign*2) & + cR * f(i,j,k+sign*3) & + dR * f(i,j,k+sign*4)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-1) = (bnd4C(1)*f(i,j,zbegin-3) + & bnd4C(2)*f(i,j,zbegin-2) + & bnd4C(3)*f(i,j,zbegin-1) + & bnd4C(4)*f(i,j,zbegin+0) + & bnd4C(5)*f(i,j,zbegin+1))*rdz end do end do end if ! if ( nproLz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-2) = (bnd4C(1)*f(i,j,zbegin-4) + & bnd4C(2)*f(i,j,zbegin-3) + & bnd4C(3)*f(i,j,zbegin-2) + & bnd4C(4)*f(i,j,zbegin-1) + & bnd4C(5)*f(i,j,zbegin+0))*rdz end do end do end if ! if ( nproLz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-3) = (bnd3B(1)*f(i,j,zbegin-4) + & bnd3B(2)*f(i,j,zbegin-3) + & bnd3B(3)*f(i,j,zbegin-2) + & bnd3B(4)*f(i,j,zbegin-1))*rdz end do end do end if ! if ( nproLz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zbegin-4) = (bnd3A(1)*f(i,j,zbegin-4) + & bnd3A(2)*f(i,j,zbegin-3) + & bnd3A(3)*f(i,j,zbegin-2) + & bnd3A(4)*f(i,j,zbegin-1))*rdz end do end do end if ! ! Right side: ! if ( nproRz .ge. 1 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+1) = -(bnd4C(1)*f(i,j,zend+3) + & bnd4C(2)*f(i,j,zend+2) + & bnd4C(3)*f(i,j,zend+1) + & bnd4C(4)*f(i,j,zend+0) + & bnd4C(5)*f(i,j,zend-1))*rdz end do end do end if ! if ( nproRz .ge. 2 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+2) = -(bnd4C(1)*f(i,j,zend+4) + & bnd4C(2)*f(i,j,zend+3) + & bnd4C(3)*f(i,j,zend+2) + & bnd4C(4)*f(i,j,zend+1) + & bnd4C(5)*f(i,j,zend+0))*rdz end do end do end if if ( nproRz .ge. 3 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+3) = -(bnd3B(1)*f(i,j,zend+4) + & bnd3B(2)*f(i,j,zend+3) + & bnd3B(3)*f(i,j,zend+2) + & bnd3B(4)*f(i,j,zend+1))*rdz end do end do end if if ( nproRz .eq. 4 ) then do j = ijL, ijR do i = iiL, iiR df(i,j,zend+4) = -(bnd3A(1)*f(i,j,zend+4) + & bnd3A(2)*f(i,j,zend+3) + & bnd3A(3)*f(i,j,zend+2) + & bnd3A(4)*f(i,j,zend+1))*rdz end do end do end if ! end if ! end if ! ! End Super Upwind Difference Operators ! return end ! !===========================================================================