qmin = q(i,1) - min(q(i,2),q(i,1), q(i+im2,2))            dm(i,1) = sign(min(abs(tmp),qmax,qmin),tmp)          enddo          do i=im2+1,im            dm(i, 1) =  - dm(i-im2, 1)          enddo        endif        if ( jlast == jm ) then! N pole          do i=1,im2            tmp = 0.25*(q(i+im2,jm1)-q(i,jm1))            qmax = max(q(i+im2,jm1),q(i,jm), q(i,jm1)) - q(i,jm)            qmin = q(i,jm) - min(q(i+im2,jm1),q(i,jm), q(i,jm1))            dm(i,jm) = sign(min(abs(tmp),qmax,qmin),tmp)          enddo          do i=im2+1,im            dm(i,jm) =  - dm(i-im2,jm)          enddo        endif   else        if ( jfirst == 1 ) then! South          do i=1,im2            tmp  = 0.25*(q(i,2)+q(i+im2,2))            qmax = max(q(i,2),q(i,1), -q(i+im2,2)) - q(i,1)            qmin = q(i,1) - min(q(i,2),q(i,1),-q(i+im2,2))            dm(i,1) = sign(min(abs(tmp),qmax,qmin),tmp)          enddo          do i=im2+1,im            dm(i, 1) = dm(i-im2, 1)          enddo        endif        if ( jlast == jm ) then! North          do i=1,im2            tmp  = -0.25*(q(i+im2,jm1)+q(i,jm1))            qmax = max(-q(i+im2,jm1),q(i,jm), q(i,jm1)) - q(i,jm)            qmin = q(i,jm) - min(-q(i+im2,jm1),q(i,jm), q(i,jm1))            dm(i,jm) = sign(min(abs(tmp),qmax,qmin),tmp)          enddo          do i=im2+1,im            dm(i,jm) = dm(i-im2,jm)          enddo        endif   endif   if( jord > 0 ) then!! Applies monotonic slope constraint (off if jord less than zero)!!!!        do j=2,jm1        do j=js2gng1,jn2gng1          do i=1,im            qmax = max(q(i,j-1),q(i,j),q(i,j+1)) - q(i,j)            qmin = q(i,j) - min(q(i,j-1),q(i,j),q(i,j+1))            dm(i,j) = sign(min(abs(dm(i,j)),qmin,qmax),dm(i,j))          enddo        enddo   endif    return!EOC end subroutine ymist!-----------------------------------------------------------------------!-----------------------------------------------------------------------!BOP! !IROUTINE: fyppm!! !INTERFACE:  subroutine fyppm(c,  q,  dm, flux, im, jm, ng, jord, iv, jfirst, jlast)!-----------------------------------------------------------------------! !USES: use precision implicit none! !INPUT PARAMETERS: integer im, jm                      !  Dimensions integer jfirst, jlast               !  Latitude strip integer ng                          !  Max. NS dependencies integer jord                        !  Approximation order integer iv                          !  Scalar=0, Vector=1 real (r8)  q(im,jfirst-ng:jlast+ng) !  mean value needed only N*2 S*2 real (r8) dm(im,jfirst-ng:jlast+ng) !  Slope     needed only N*2 S*2 real (r8)  c(im,jfirst:jlast+1)     !  Courant   N (like FLUX)! !INPUT/OUTPUT PARAMETERS: real (r8) ar(im,jfirst-1:jlast+1)   ! AR needs to be ghosted on NS real (r8) al(im,jfirst-1:jlast+2)   ! AL needs to be ghosted on N2S real (r8) a6(im,jfirst-1:jlast+1)   ! A6 needs to be ghosted on NS! !OUTPUT PARAMETERS: real (r8) flux(im,jfirst:jlast+1)   !  Flux      N (see tp2c)! !DESCRIPTION:!!   NG is passed from YTP for convenience -- it is actually 1 more in NS!   than the actual number of latitudes needed here.  But in the shared-memory !   case it becomes 0, which is much cleaner.!! !CALLED FROM:!      ytp!! !REVISION HISTORY:!!  SJL 99.04.13:  Delivery!  WS  99.04.19:  Added jfirst:jlast concept; FYPPM only called from YTP!  WS  99.04.21:  Removed j1, j2  (j1=2, j2=jm-1 consistently)!                 removed a6,ar,al from argument list!  WS  99.09.09:  Documentation; indentation; cleaning!  WS  99.10.22:  Added ng as argument; Pruned arrays!  WS  00.05.14:  Renamed ghost indices as per Kevin's definitions!!EOP!---------------------------------------------------------------------!BOC real (r8)   r3, r23 parameter ( r3 = 1./3., r23 = 2./3. ) integer i, j, imh, jm1, lmt integer js1g1, js2g0, js2g1, jn1g2, jn1g1, jn2g1!     logical steep!     if(jord .eq. 6) then!        steep = .true.!     else!        steep = .false.!     endif      imh = im / 2      jm1 = jm - 1      js1g1  = max(1,jfirst-1)         ! Ghost S*1      js2g0  = max(2,jfirst)           ! No ghosting      js2g1  = max(2,jfirst-1)         ! Ghost S*1      jn1g1  = min(jm,jlast+1)         ! Ghost N*1      jn1g2  = min(jm,jlast+2)         ! Ghost N*2      jn2g1  = min(jm-1,jlast+1)       ! Ghost N*1!!!      do j=2,jm      do j=js2g1,jn1g2                 ! AL needed N2S        do i=1,im                      ! P, dm ghosted N2S2 (at least)          al(i,j) = 0.5*(q(i,j-1)+q(i,j)) + r3*(dm(i,j-1) - dm(i,j))        enddo      enddo! Yeh's steepening procedure; to be implemented!     if(steep) call steepy(im,   jm,   jfirst,   jlast,       &!                           ng,    q,       al,   dm )!!!      do j=1,jm1      do j=js1g1,jn2g1                 ! AR needed NS        do i=1,im          ar(i,j) = al(i,j+1)          ! AL ghosted N2S        enddo      enddo! WS 990726 :  Added condition to decide if poles are on this processor! Poles:   if( iv == 0 ) then        if ( jfirst .eq. 1 ) then          do i=1,imh            al(i,    1) = al(i+imh,2)            al(i+imh,1) = al(i,    2)          enddo        endif        if ( jlast .eq. jm ) then          do i=1,imh            ar(i,    jm) = ar(i+imh,jm1)            ar(i+imh,jm) = ar(i,    jm1)          enddo        endif   else        if ( jfirst .eq. 1 ) then          do i=1,imh            al(i,    1) = -al(i+imh,2)            al(i+imh,1) = -al(i,    2)          enddo        endif        if ( jlast .eq. jm ) then          do i=1,imh            ar(i,    jm) = -ar(i+imh,jm1)            ar(i+imh,jm) = -ar(i,    jm1)          enddo        endif   endif   if( jord == 3 .or. jord == 5 ) then!!!      do j=1,jm      do j=js1g1,jn1g1               ! A6 needed NS        do i=1,im          a6(i,j) = 3.*(q(i,j)+q(i,j) - (al(i,j)+ar(i,j)))        enddo      enddo   endif      lmt = jord - 3!!!        do j=1,jm!       do j=js1g1,jn1g1             !  A6, AR, AL needed NS!         call lmppm(dm(1,j),a6(1,j),ar(1,j),al(1,j),q(1,j),im,lmt)!       enddo        call lmppm(dm(1,js1g1), a6(1,js1g1), ar(1,js1g1),               &                   al(1,js1g1),  q(1,js1g1), im*(jn1g1-js1g1+1), lmt)!!!      do j=2,jm      do j=js2g0,jn1g1                 ! flux needed N        do i=1,im          if(c(i,j).gt.0.) then            flux(i,j) = ar(i,j-1) + 0.5*c(i,j)*(al(i,j-1) - ar(i,j-1) +  &                        a6(i,j-1)*(1.-r23*c(i,j)) )          else            flux(i,j) = al(i,j) - 0.5*c(i,j)*(ar(i,j) - al(i,j) +        &                        a6(i,j)*(1.+r23*c(i,j)))          endif        enddo      enddo    return!EOC end subroutine fyppm !-----------------------------------------------------------------------!-----------------------------------------------------------------------!BOP! !IROUTINE: tpcc!! !INTERFACE:  subroutine tpcc(va,   ymass,  q,   crx,  cry,  im,   jm,      ng,     &                 iord, jord,   fx,  fy,   ffsl, cose, jfirst, jlast,   &                 dm,   qtmp,   al,  ar,   a6 )       !-----------------------------------------------------------------------! !USES: use precision implicit none! !INPUT PARAMETERS:  integer im, jm                    ! Dimensions  integer ng                        ! ng_d (needed to properly declare q  integer jfirst, jlast             ! Latitude strip  integer iord, jord                ! Interpolation order in x,y  logical ffsl(jm)                  ! Flux-form semi-Lagrangian transport?  real (r8) cose(jm)                ! Critical cosine  (replicated)  real (r8) va(im,jfirst:jlast)     ! Courant (unghosted like FX)  real (r8) q(im,jfirst-ng:jlast+ng)! C-grid abs vort: N*north S*(north-1)  real (r8) crx(im,jfirst-1:jlast+2)  real (r8) cry(im,jfirst:jlast)    ! Courant # (ghosted like FY)  real (r8) ymass(im,jfirst:jlast)  ! Background y-mass-flux (ghosted like FY)! Input 1D work arrays:  real (r8)   dm(-im/3:im+im/3)  real (r8) qtmp(-im/3:im+im/3)  real (r8)   al(-im/3:im+im/3)  real (r8)   ar(-im/3:im+im/3)  real (r8)   a6(-im/3:im+im/3)! !OUTPUT PARAMETERS:  real (r8) fx(im,jfirst:jlast)     ! Flux in x (unghosted)  real (r8) fy(im,jfirst:jlast)     ! Flux in y (unghosted since iv==0)! !DESCRIPTION:!     In this routine the number !     of north ghosted latitude min(abs(jord),2), and south ghosted!     latitudes is XXXX!! !CALLED FROM:!     cd_core!! !REVISION HISTORY:!   SJL 99.04.13:  Delivery!   WS  99.04.13:  Added jfirst:jlast concept!   WS  99.05.10:  Replaced JNP with JM, JMR with JM-1, IMR with IM!   WS  99.05.10:  Removed fvcore.h and JNP, IMH, IML definitions!   WS  99.10.20:  Pruned arrays!   WS  00.05.14:  Renamed ghost indices as per Kevin's definitions!!EOP!-----------------------------------------------------------------------!BOC!! !LOCAL VARIABLES:  real (r8) adx(im,jfirst-1:jlast+2)  integer north, south  integer i, j, jp, im2, js2g0, js2gs, jn2g0, jn1g0, jn1gn  real (r8) wk1(im)  real (r8) fx1(im)    im2 = im/2    north = min(2,abs(jord))         ! north == 1 or 2    south = north-1                  ! south == 0 or 1    js2g0 = max(2,jfirst)    js2gs = max(2,jfirst-south)    jn2g0 = min(jm-1,jlast)    jn1gn = min(jm,jlast+north)    jn1g0 = min(jm,jlast)! This loop must be ghosted N*NG, S*NG!!!      do j=2,jm      do j=js2gs,jn1gn        call xtp( im, ffsl(j), wk1, q(1,j),           &                  crx(1,j), 1, crx(1,j), cose(j), 0,  &                  dm, qtmp, al, ar, a6 )        do i=1,im-1          adx(i,j) = q(i,j) + 0.5 *                   &                     (wk1(i)-wk1(i+1) + q(i,j)*(crx(i+1,j)-crx(i,j)))        enddo        adx(im,j) = q(im,j) + 0.5 *                   &                    (wk1(im)-wk1(1) + q(im,j)*(crx(1,j)-crx(im,j)))      enddo      call ycc(im, jm, fy, adx, cry, ymass, jord, 0,jfirst,jlast)! For Scalar only!!!      if ( jfirst == 1 ) then        do i=1,im2          q(i,1) = q(i+im2,  2)        enddo        do i=im2+1,im           q(i,1) = q(i-im2,  2)        enddo      endif      if ( jlast == jm ) then        do i=1,im2          fx1(i) = q(i+im2,jm)        enddo        do i=im2+1,im           fx1(i) = q(i-im2,jm)        enddo        do i=1,im          if(va(i,jm) .gt. 0.) then            adx(i,jm) = q(i,jm) + 0.5*va(i,jm)*(q(i,jm-1)-q(i,jm))          else            adx(i,jm) = q(i,jm) + 0.5*va(i,jm)*(q(i,jm)-fx1(i))          endif        enddo      endif!!!      do j=2,jm-1      do j=js2g0,jn2g0        do i=1,im          jp = j-va(i,j)! jp = j     if va < 0! jp = j -1  if va < 0! q needed max(1, jfirst-1)          adx(i,j) = q(i,j) + 0.5*va(i,j)*(q(i,jp)-q(i,jp+1))        enddo      enddo!!!      do j=2,jm      do j=js2g0,jn1g0        call xtp( im, ffsl(j), fx(1,j), adx(1,j),                 &                  crx(1,j), iord, crx(1,j), cose(j), 0,           &                  dm, qtmp, al, ar, a6 )      enddo    return!EOC end subroutine tpcc!-----------------------------------------------------------------------!-----------------------------------------------------------------------!BOP! !IROUTINE: ycc!! !INTERFACE:  subroutine ycc(im, jm, fy, q, vc, ymass, jord, iv, jfirst, jlast)!-----------------------------------------------------------------------! !USES: use precision implicit none! !INPUT PARAMETERS: integer im, jm                      !  Dimensions integer jfirst, jlast               !  Latitude strip integer jord                        !  Approximation order integer iv                          !  Scalar=0, Vector=1 real (r8) q(im,jfirst-1-iv:jlast+2)      !  Field (N*2 S*(iv+1)) real (r8) vc(im,jfirst-iv:jlast)         !  Courant  (like FY) real (r8) ymass(im,jfirst-iv:jlast)      !  background mass flux! !OUTPUT PARAMETERS: real (r8) fy(im,jfirst-iv:jlast)         !  Flux (S if iv=1)! !DESCRIPTION:!     Will Sawyer's note: In this routine the number !     of ghosted latitudes NG is min(abs(jord),2).  The scalar/vector!     flag determines whether the flux FY needs to be ghosted on the!     south.  If called from CD\_CORE (iv==1) then it does, if called!     from TPCC (iv==0) it does not.  !! !CALLED FROM:!     cd_core!     tpcc!! !REVISION HISTORY:!!   SJL 99.04.13:  Delivery!   WS  99.04.19:  Added jfirst:jlast concept!   WS  99.04.27:  DC removed as argument (local to this routine); DC on N!   WS  99.05.10:  Replaced JNP with JM, JMR with JM-1, IMR with IM!   WS  99.05.10:  Removed fvcore.h!   WS  99.07.27:  Built in tests for SP or NP!   WS  99.09.09:  Documentation; indentation; cleaning; pole treatment!   WS  99.09.14:  Loop limits!   WS  00.05.14:  Renamed ghost indices as per Kevin's definitions!!EOP!---------------------------------------------------------------------!BOC! !LOCAL VARIABLES:  real (r8) dc(im,jfirst-iv:jlast+1)  real (r8) qmax, qmin  integer i, j, jt, im2, js2giv, js3giv, jn2g1, jn2g0   im2 = im/2   js2giv = max(2,jfirst-iv)   js3giv = max(3,jfirst-iv)   jn2g1  = min(jm-1,jlast+1)   jn2g0  = min(jm-1,jlast)         if(jord == 1) then!!!        do j=2,jm-1        do j=js2giv,jn2g0                      ! FY needed on S*iv          do i=1,im! jt=j if vc > 0; jt=j+1 if vc <=0            jt = float(j+1)  - vc(i,j)         ! VC ghosted like fy            fy(i,j) = q(i,jt)*ymass(i,j)       ! ymass ghosted like fy          enddo                                ! q ghosted N*1, S*iv        enddo   else!!!        do j=3,jm-1        do j=js3giv,jn2g1                      ! dc needed N*1, S*iv          do i=1,im            dc(i,j) = 0.25*(q(i,j+1)-q(i,j-1)) ! q ghosted N*2, S*(iv+1)          enddo        enddo        if(iv.eq.0) then! Scalar.! WS 99.07.27 : Split loops in SP and NP regions, added SP/NP tests          if ( jfirst == 1 ) then            do i=1,im2              dc(i, 2) = 0.25 * ( q(i,3) - q(i+im2,2) )            enddo            do i=im2+1,im              dc(i, 2) = 0.25 * ( q(i,3) - q(i-im2,2) )            enddo          endif          if ( jlast == jm ) then            do i=1,im2              dc(i,jm) = 0.25 * ( q(i+im2,jm) - q(i,jm-1) )            enddo            do i=im2+1,im              dc(i,jm) = 0.25 * ( q(i-im2,jm) - q(i,jm-1) )            enddo          endif        else! Vector winds! WS 99.07.27 : Split loops in SP and NP regions, added SP/NP tests          if ( jfirst == 1 ) then            do i=1,im2              dc(i, 2) =  0.25 * ( q(i,3) + q(i+im2,2) )            enddo            do i=im2+1,im              dc(i, 2) =  0.25 * ( q(i,3) + q(i-im2,2) )            enddo          endif          if ( jlast == jm ) then            do i=1,im2              dc(i,jm) = -0.25 * ( q(i,jm-1) + q(i+im2,jm) )            enddo            do i=im2+1,im              dc(i,jm) = -0.25 * ( q(i,jm-1) + q(i-im2,jm) )            enddo          endif        endif        if( jord > 0 ) then! Monotonic constraint!!!          do j=3,jm-1          do j=js3giv,jn2g1            ! DC needed N*1, S*iv            do i=1,im                  ! P ghosted N*2, S*(iv+1)              qmax = max(q(i,j-1),q(i,j),q(i,j+1)) - q(i,j)              qmin = q(i,j) - min(q(i,j-1),q(i,j),q(i,j+1))              dc(i,j) = sign(min(abs(dc(i,j)),qmin,qmax),dc(i,j))            enddo          enddo!! WS 99.08.03 : Following loop split into SP and NP part!          if ( jfirst == 1 ) then            do i=1,im              dc(i, 2) = 0.            enddo          endif          if ( jlast == jm ) then            do i=1,im              dc(i,jm) = 0.            enddo          endif        endif!!!        do j=2,jm-1       do j=js2giv,jn2g0                   ! fy needed S*iv         do i=1,im                                  jt = float(j+1)  - vc(i,j)      ! vc, ymass ghosted like fy           fy(i,j) = (q(i,jt)+(sign(1.,vc(i,j))-vc(i,j))*dc(i,jt))*ymass(i,j)         enddo       enddo    endif    return!EOC end subroutine ycc!-----------------------------------------------------------------------end module tp_core