!=========================================================================== ! subroutine y_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. pjR-pjL+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) * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend do i = iiL, iiR df(i,j,k) = ae *( f(i,j+1,k)-f(i,j-1,k) ) end do end do end do ! ! Boundary nodes: ! Left side: ! if ( nproLy .eq. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd1A(1)*f(i,ybegin-1,k) + & bnd1A(2)*f(i,ybegin+0,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .eq. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd1A(1)*f(i,yend+1,k) + & bnd1A(2)*f(i,yend+0,k))*rdy 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 * rdy be = -1.d0/12.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend do i = iiL, iiR df(i,j,k) = ae *( f(i,j+1,k)-f(i,j-1,k) ) & + be *( f(i,j+2,k)-f(i,j-2,k) ) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd3B(1)*f(i,ybegin-2,k) + & bnd3B(2)*f(i,ybegin-1,k) + & bnd3B(3)*f(i,ybegin+0,k) + & bnd3B(4)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .eq. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd3A(1)*f(i,ybegin-2,k) + & bnd3A(2)*f(i,ybegin-1,k) + & bnd3A(3)*f(i,ybegin+0,k) + & bnd3A(4)*f(i,ybegin+1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd3B(1)*f(i,yend+2,k) + & bnd3B(2)*f(i,yend+1,k) + & bnd3B(3)*f(i,yend+0,k) + & bnd3B(4)*f(i,yend-1,k))*rdy end do end do end if if ( nproRy .eq. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd3A(1)*f(i,yend+2,k) + & bnd3A(2)*f(i,yend+1,k) + & bnd3A(3)*f(i,yend+0,k) + & bnd3A(4)*f(i,yend-1,k))*rdy 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 * rdy be = -3.d0/20.d0 * rdy ce = 1.d0/60.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend do i = iiL, iiR df(i,j,k) = ae *( f(i,j+1,k)-f(i,j-1,k) ) & + be *( f(i,j+2,k)-f(i,j-2,k) ) & + ce *( f(i,j+3,k)-f(i,j-3,k) ) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd4C(1)*f(i,ybegin-3,k) + & bnd4C(2)*f(i,ybegin-2,k) + & bnd4C(3)*f(i,ybegin-1,k) + & bnd4C(4)*f(i,ybegin+0,k) + & bnd4C(5)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd3B(1)*f(i,ybegin-3,k) + & bnd3B(2)*f(i,ybegin-2,k) + & bnd3B(3)*f(i,ybegin-1,k) + & bnd3B(4)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .eq. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3A(1)*f(i,ybegin-3,k) + & bnd3A(2)*f(i,ybegin-2,k) + & bnd3A(3)*f(i,ybegin-1,k) + & bnd3A(4)*f(i,ybegin+0,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd4C(1)*f(i,yend+3,k) + & bnd4C(2)*f(i,yend+2,k) + & bnd4C(3)*f(i,yend+1,k) + & bnd4C(4)*f(i,yend+0,k) + & bnd4C(5)*f(i,yend-1,k))*rdy end do end do end if if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd3B(1)*f(i,yend+3,k) + & bnd3B(2)*f(i,yend+2,k) + & bnd3B(3)*f(i,yend+1,k) + & bnd3B(4)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .eq. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3A(1)*f(i,yend+3,k) + & bnd3A(2)*f(i,yend+2,k) + & bnd3A(3)*f(i,yend+1,k) + & bnd3A(4)*f(i,yend+0,k))*rdy 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 * rdy be = -1.d0/ 5.d0 * rdy ce = 4.d0/105.d0 * rdy de = -1.d0/280.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend do i = iiL, iiR df(i,j,k) = ae *( f(i,j+1,k)-f(i,j-1,k) ) & + be *( f(i,j+2,k)-f(i,j-2,k) ) & + ce *( f(i,j+3,k)-f(i,j-3,k) ) & + de *( f(i,j+4,k)-f(i,j-4,k) ) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd6D(1)*f(i,ybegin-4,k) + & bnd6D(2)*f(i,ybegin-3,k) + & bnd6D(3)*f(i,ybegin-2,k) + & bnd6D(4)*f(i,ybegin-1,k) + & bnd6D(5)*f(i,ybegin+0,k) + & bnd6D(6)*f(i,ybegin+1,k) + & bnd6D(7)*f(i,ybegin+2,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd4C(1)*f(i,ybegin-4,k) + & bnd4C(2)*f(i,ybegin-3,k) + & bnd4C(3)*f(i,ybegin-2,k) + & bnd4C(4)*f(i,ybegin-1,k) + & bnd4C(5)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3B(1)*f(i,ybegin-4,k) + & bnd3B(2)*f(i,ybegin-3,k) + & bnd3B(3)*f(i,ybegin-2,k) + & bnd3B(4)*f(i,ybegin-1,k))*rdy end do end do end if ! if ( nproLy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-4,k) = (bnd3A(1)*f(i,ybegin-4,k) + & bnd3A(2)*f(i,ybegin-3,k) + & bnd3A(3)*f(i,ybegin-2,k) + & bnd3A(4)*f(i,ybegin-1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd6D(1)*f(i,yend+4,k) + & bnd6D(2)*f(i,yend+3,k) + & bnd6D(3)*f(i,yend+2,k) + & bnd6D(4)*f(i,yend+1,k) + & bnd6D(5)*f(i,yend+0,k) + & bnd6D(6)*f(i,yend-1,k) + & bnd6D(7)*f(i,yend-2,k))*rdy end do end do end if ! if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd4C(1)*f(i,yend+4,k) + & bnd4C(2)*f(i,yend+3,k) + & bnd4C(3)*f(i,yend+2,k) + & bnd4C(4)*f(i,yend+1,k) + & bnd4C(5)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3B(1)*f(i,yend+4,k) + & bnd3B(2)*f(i,yend+3,k) + & bnd3B(3)*f(i,yend+2,k) + & bnd3B(4)*f(i,yend+1,k))*rdy end do end do end if if ( nproRy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+4,k) = -(bnd3A(1)*f(i,yend+4,k) + & bnd3A(2)*f(i,yend+3,k) + & bnd3A(3)*f(i,yend+2,k) + & bnd3A(4)*f(i,yend+1,k))*rdy 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 * rdy Up = -1.d0/ 2.d0 * rdy aR = 1.d0/ 1.d0 * rdy bR = -1.d0/ 6.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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-sign*1,k) & + Up * f(i,j ,k) & + aR * f(i,j+sign*1,k) & + bR * f(i,j+sign*2,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd3B(1)*f(i,ybegin-2,k) + & bnd3B(2)*f(i,ybegin-1,k) + & bnd3B(3)*f(i,ybegin+0,k) + & bnd3B(4)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .eq. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd3A(1)*f(i,ybegin-2,k) + & bnd3A(2)*f(i,ybegin-1,k) + & bnd3A(3)*f(i,ybegin+0,k) + & bnd3A(4)*f(i,ybegin+1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd3B(1)*f(i,yend+2,k) + & bnd3B(2)*f(i,yend+1,k) + & bnd3B(3)*f(i,yend+0,k) + & bnd3B(4)*f(i,yend-1,k))*rdy end do end do end if if ( nproRy .eq. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd3A(1)*f(i,yend+2,k) + & bnd3A(2)*f(i,yend+1,k) + & bnd3A(3)*f(i,yend+0,k) + & bnd3A(4)*f(i,yend-1,k))*rdy 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 * rdy Up = -5.d0/ 6.d0 * rdy aR = 3.d0/ 2.d0 * rdy bR = -1.d0/ 2.d0 * rdy cR = 1.d0/12.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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-sign*1,k) & + Up * f(i,j ,k) & + aR * f(i,j+sign*1,k) & + bR * f(i,j+sign*2,k) & + cR * f(i,j+sign*3,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd4C(1)*f(i,ybegin-3,k) + & bnd4C(2)*f(i,ybegin-2,k) + & bnd4C(3)*f(i,ybegin-1,k) + & bnd4C(4)*f(i,ybegin+0,k) + & bnd4C(5)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd3B(1)*f(i,ybegin-3,k) + & bnd3B(2)*f(i,ybegin-2,k) + & bnd3B(3)*f(i,ybegin-1,k) + & bnd3B(4)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .eq. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3A(1)*f(i,ybegin-3,k) + & bnd3A(2)*f(i,ybegin-2,k) + & bnd3A(3)*f(i,ybegin-1,k) + & bnd3A(4)*f(i,ybegin+0,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd4C(1)*f(i,yend+3,k) + & bnd4C(2)*f(i,yend+2,k) + & bnd4C(3)*f(i,yend+1,k) + & bnd4C(4)*f(i,yend+0,k) + & bnd4C(5)*f(i,yend-1,k))*rdy end do end do end if if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd3B(1)*f(i,yend+3,k) + & bnd3B(2)*f(i,yend+2,k) + & bnd3B(3)*f(i,yend+1,k) + & bnd3B(4)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .eq. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3A(1)*f(i,yend+3,k) + & bnd3A(2)*f(i,yend+2,k) + & bnd3A(3)*f(i,yend+1,k) + & bnd3A(4)*f(i,yend+0,k))*rdy 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 * rdy aL = -1.d0/ 2.d0 * rdy Up = -1.d0/ 3.d0 * rdy aR = 1.d0/ 1.d0 * rdy bR = -1.d0/ 4.d0 * rdy cR = 1.d0/30.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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-sign*2,k) & + aL * f(i,j-sign*1,k) & + Up * f(i,j ,k) & + aR * f(i,j+sign*1,k) & + bR * f(i,j+sign*2,k) & + cR * f(i,j+sign*3,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd4C(1)*f(i,ybegin-3,k) + & bnd4C(2)*f(i,ybegin-2,k) + & bnd4C(3)*f(i,ybegin-1,k) + & bnd4C(4)*f(i,ybegin+0,k) + & bnd4C(5)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd3B(1)*f(i,ybegin-3,k) + & bnd3B(2)*f(i,ybegin-2,k) + & bnd3B(3)*f(i,ybegin-1,k) + & bnd3B(4)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .eq. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3A(1)*f(i,ybegin-3,k) + & bnd3A(2)*f(i,ybegin-2,k) + & bnd3A(3)*f(i,ybegin-1,k) + & bnd3A(4)*f(i,ybegin+0,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd4C(1)*f(i,yend+3,k) + & bnd4C(2)*f(i,yend+2,k) + & bnd4C(3)*f(i,yend+1,k) + & bnd4C(4)*f(i,yend+0,k) + & bnd4C(5)*f(i,yend-1,k))*rdy end do end do end if if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd3B(1)*f(i,yend+3,k) + & bnd3B(2)*f(i,yend+2,k) + & bnd3B(3)*f(i,yend+1,k) + & bnd3B(4)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .eq. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3A(1)*f(i,yend+3,k) + & bnd3A(2)*f(i,yend+2,k) + & bnd3A(3)*f(i,yend+1,k) + & bnd3A(4)*f(i,yend+0,k))*rdy 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 * rdy aL = - 2.d0/ 5.d0 * rdy Up = - 7.d0/12.d0 * rdy aR = 4.d0/ 3.d0 * rdy bR = -1.d0/ 2.d0 * rdy cR = 2.d0/15.d0 * rdy dR = -1.d0/60.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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-sign*2,k) & + aL * f(i,j-sign*1,k) & + Up * f(i,j ,k) & + aR * f(i,j+sign*1,k) & + bR * f(i,j+sign*2,k) & + cR * f(i,j+sign*3,k) & + dR * f(i,j+sign*4,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd6D(1)*f(i,ybegin-4,k) + & bnd6D(2)*f(i,ybegin-3,k) + & bnd6D(3)*f(i,ybegin-2,k) + & bnd6D(4)*f(i,ybegin-1,k) + & bnd6D(5)*f(i,ybegin+0,k) + & bnd6D(6)*f(i,ybegin+1,k) + & bnd6D(7)*f(i,ybegin+2,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd4C(1)*f(i,ybegin-4,k) + & bnd4C(2)*f(i,ybegin-3,k) + & bnd4C(3)*f(i,ybegin-2,k) + & bnd4C(4)*f(i,ybegin-1,k) + & bnd4C(5)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3B(1)*f(i,ybegin-4,k) + & bnd3B(2)*f(i,ybegin-3,k) + & bnd3B(3)*f(i,ybegin-2,k) + & bnd3B(4)*f(i,ybegin-1,k))*rdy end do end do end if ! if ( nproLy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-4,k) = (bnd3A(1)*f(i,ybegin-4,k) + & bnd3A(2)*f(i,ybegin-3,k) + & bnd3A(3)*f(i,ybegin-2,k) + & bnd3A(4)*f(i,ybegin-1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd6D(1)*f(i,yend+4,k) + & bnd6D(2)*f(i,yend+3,k) + & bnd6D(3)*f(i,yend+2,k) + & bnd6D(4)*f(i,yend+1,k) + & bnd6D(5)*f(i,yend+0,k) + & bnd6D(6)*f(i,yend-1,k) + & bnd6D(7)*f(i,yend-2,k))*rdy end do end do end if ! if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd4C(1)*f(i,yend+4,k) + & bnd4C(2)*f(i,yend+3,k) + & bnd4C(3)*f(i,yend+2,k) + & bnd4C(4)*f(i,yend+1,k) + & bnd4C(5)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3B(1)*f(i,yend+4,k) + & bnd3B(2)*f(i,yend+3,k) + & bnd3B(3)*f(i,yend+2,k) + & bnd3B(4)*f(i,yend+1,k))*rdy end do end do end if if ( nproRy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+4,k) = -(bnd3A(1)*f(i,yend+4,k) + & bnd3A(2)*f(i,yend+3,k) + & bnd3A(3)*f(i,yend+2,k) + & bnd3A(4)*f(i,yend+1,k))*rdy 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 * rdy bL = 1.d0/ 10.d0 * rdy aL = -3.d0/ 5.d0 * rdy Up = -1.d0/ 4.d0 * rdy aR = 1.d0/ 1.d0 * rdy bR = -3.d0/ 10.d0 * rdy cR = 1.d0/ 15.d0 * rdy dR = -1.d0/140.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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-sign*3,k) & + bL * f(i,j-sign*2,k) & + aL * f(i,j-sign*1,k) & + Up * f(i,j ,k) & + aR * f(i,j+sign*1,k) & + bR * f(i,j+sign*2,k) & + cR * f(i,j+sign*3,k) & + dR * f(i,j+sign*4,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd6D(1)*f(i,ybegin-4,k) + & bnd6D(2)*f(i,ybegin-3,k) + & bnd6D(3)*f(i,ybegin-2,k) + & bnd6D(4)*f(i,ybegin-1,k) + & bnd6D(5)*f(i,ybegin+0,k) + & bnd6D(6)*f(i,ybegin+1,k) + & bnd6D(7)*f(i,ybegin+2,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd4C(1)*f(i,ybegin-4,k) + & bnd4C(2)*f(i,ybegin-3,k) + & bnd4C(3)*f(i,ybegin-2,k) + & bnd4C(4)*f(i,ybegin-1,k) + & bnd4C(5)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3B(1)*f(i,ybegin-4,k) + & bnd3B(2)*f(i,ybegin-3,k) + & bnd3B(3)*f(i,ybegin-2,k) + & bnd3B(4)*f(i,ybegin-1,k))*rdy end do end do end if ! if ( nproLy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-4,k) = (bnd3A(1)*f(i,ybegin-4,k) + & bnd3A(2)*f(i,ybegin-3,k) + & bnd3A(3)*f(i,ybegin-2,k) + & bnd3A(4)*f(i,ybegin-1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd6D(1)*f(i,yend+4,k) + & bnd6D(2)*f(i,yend+3,k) + & bnd6D(3)*f(i,yend+2,k) + & bnd6D(4)*f(i,yend+1,k) + & bnd6D(5)*f(i,yend+0,k) + & bnd6D(6)*f(i,yend-1,k) + & bnd6D(7)*f(i,yend-2,k))*rdy end do end do end if ! if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd4C(1)*f(i,yend+4,k) + & bnd4C(2)*f(i,yend+3,k) + & bnd4C(3)*f(i,yend+2,k) + & bnd4C(4)*f(i,yend+1,k) + & bnd4C(5)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3B(1)*f(i,yend+4,k) + & bnd3B(2)*f(i,yend+3,k) + & bnd3B(3)*f(i,yend+2,k) + & bnd3B(4)*f(i,yend+1,k))*rdy end do end do end if if ( nproRy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+4,k) = -(bnd3A(1)*f(i,yend+4,k) + & bnd3A(2)*f(i,yend+3,k) + & bnd3A(3)*f(i,yend+2,k) + & bnd3A(4)*f(i,yend+1,k))*rdy 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 * rdy aR = 4.d0/ 1.d0 * rdy bR = -3.d0/ 1.d0 * rdy cR = 4.d0/ 3.d0 * rdy dR = -1.d0/ 4.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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+sign*1,k) & + bR * f(i,j+sign*2,k) & + cR * f(i,j+sign*3,k) & + dR * f(i,j+sign*4,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd4C(1)*f(i,ybegin-3,k) + & bnd4C(2)*f(i,ybegin-2,k) + & bnd4C(3)*f(i,ybegin-1,k) + & bnd4C(4)*f(i,ybegin-0,k) + & bnd4C(5)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd4C(1)*f(i,ybegin-4,k) + & bnd4C(2)*f(i,ybegin-3,k) + & bnd4C(3)*f(i,ybegin-2,k) + & bnd4C(4)*f(i,ybegin-1,k) + & bnd4C(5)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3B(1)*f(i,ybegin-4,k) + & bnd3B(2)*f(i,ybegin-3,k) + & bnd3B(3)*f(i,ybegin-2,k) + & bnd3B(4)*f(i,ybegin-1,k))*rdy end do end do end if ! if ( nproLy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-4,k) = (bnd3A(1)*f(i,ybegin-4,k) + & bnd3A(2)*f(i,ybegin-3,k) + & bnd3A(3)*f(i,ybegin-2,k) + & bnd3A(4)*f(i,ybegin-1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd4C(1)*f(i,yend+3,k) + & bnd4C(2)*f(i,yend+2,k) + & bnd4C(3)*f(i,yend+1,k) + & bnd4C(4)*f(i,yend+0,k) + & bnd4C(5)*f(i,yend-1,k))*rdy end do end do end if ! if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd4C(1)*f(i,yend+4,k) + & bnd4C(2)*f(i,yend+3,k) + & bnd4C(3)*f(i,yend+2,k) + & bnd4C(4)*f(i,yend+1,k) + & bnd4C(5)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3B(1)*f(i,yend+4,k) + & bnd3B(2)*f(i,yend+3,k) + & bnd3B(3)*f(i,yend+2,k) + & bnd3B(4)*f(i,yend+1,k))*rdy end do end do end if if ( nproRy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+4,k) = -(bnd3A(1)*f(i,yend+4,k) + & bnd3A(2)*f(i,yend+3,k) + & bnd3A(3)*f(i,yend+2,k) + & bnd3A(4)*f(i,yend+1,k))*rdy 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 * rdy Up = -13.d0/12.d0 * rdy aR = 2.d0/ 1.d0 * rdy bR = -1.d0/ 1.d0 * rdy cR = 1.d0/ 3.d0 * rdy dR = -1.d0/20.d0 * rdy ! ! Internal nodes: ! do k = ikL, ikR do j = ybegin, yend 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-sign*1,k) & + Up * f(i,j ,k) & + aR * f(i,j+sign*1,k) & + bR * f(i,j+sign*2,k) & + cR * f(i,j+sign*3,k) & + dR * f(i,j+sign*4,k)) end do end do end do ! ! Boundary nodes: ! ! Left side: ! if ( nproLy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-1,k) = (bnd4C(1)*f(i,ybegin-3,k) + & bnd4C(2)*f(i,ybegin-2,k) + & bnd4C(3)*f(i,ybegin-1,k) + & bnd4C(4)*f(i,ybegin+0,k) + & bnd4C(5)*f(i,ybegin+1,k))*rdy end do end do end if ! if ( nproLy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-2,k) = (bnd4C(1)*f(i,ybegin-4,k) + & bnd4C(2)*f(i,ybegin-3,k) + & bnd4C(3)*f(i,ybegin-2,k) + & bnd4C(4)*f(i,ybegin-1,k) + & bnd4C(5)*f(i,ybegin+0,k))*rdy end do end do end if ! if ( nproLy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-3,k) = (bnd3B(1)*f(i,ybegin-4,k) + & bnd3B(2)*f(i,ybegin-3,k) + & bnd3B(3)*f(i,ybegin-2,k) + & bnd3B(4)*f(i,ybegin-1,k))*rdy end do end do end if ! if ( nproLy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,ybegin-4,k) = (bnd3A(1)*f(i,ybegin-4,k) + & bnd3A(2)*f(i,ybegin-3,k) + & bnd3A(3)*f(i,ybegin-2,k) + & bnd3A(4)*f(i,ybegin-1,k))*rdy end do end do end if ! ! Right side: ! if ( nproRy .ge. 1 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+1,k) = -(bnd4C(1)*f(i,yend+3,k) + & bnd4C(2)*f(i,yend+2,k) + & bnd4C(3)*f(i,yend+1,k) + & bnd4C(4)*f(i,yend+0,k) + & bnd4C(5)*f(i,yend-1,k))*rdy end do end do end if ! if ( nproRy .ge. 2 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+2,k) = -(bnd4C(1)*f(i,yend+4,k) + & bnd4C(2)*f(i,yend+3,k) + & bnd4C(3)*f(i,yend+2,k) + & bnd4C(4)*f(i,yend+1,k) + & bnd4C(5)*f(i,yend+0,k))*rdy end do end do end if if ( nproRy .ge. 3 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+3,k) = -(bnd3B(1)*f(i,yend+4,k) + & bnd3B(2)*f(i,yend+3,k) + & bnd3B(3)*f(i,yend+2,k) + & bnd3B(4)*f(i,yend+1,k))*rdy end do end do end if if ( nproRy .eq. 4 ) then do k = ikL, ikR do i = iiL, iiR df(i,yend+4,k) = -(bnd3A(1)*f(i,yend+4,k) + & bnd3A(2)*f(i,yend+3,k) + & bnd3A(3)*f(i,yend+2,k) + & bnd3A(4)*f(i,yend+1,k))*rdy end do end do end if ! end if ! end if ! ! End Super Upwind Difference Operators ! return end ! !===========================================================================