tp_core.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

            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