tp_core.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!BOC
!
! !LOCAL VARIABLES:
 real (r8) r3, r23
 parameter ( r3 = 1./3., r23 = 2./3. )

 integer i, lmt
 integer iu, itmp
 real (r8) ru
 logical steep

  if( iord == 6 ) then
      steep = .true.
  else
      steep = .false.
  endif

  do i=1,im
     al(i) = 0.5*(p(i-1)+p(i)) + (dm(i-1) - dm(i))*r3
  enddo

  if( steep ) call steepx( im, p, al(1), dm )

     do i=1,im-1
        ar(i) = al(i+1)
     enddo
        ar(im) = al(1)

  if(iord == 7) then
     call huynh(im, ar(1), al(1), p(1), a6(1), dm(1))
  else
     if(iord .eq. 3 .or. iord .eq. 5) then
         do i=1,im
            a6(i) = 3.*(p(i)+p(i)  - (al(i)+ar(i)))
         enddo
     endif
     lmt = iord - 3
     call lmppm( dm(1), a6(1), ar(1), al(1), p(1), im, lmt )
  endif

  if( ffsl ) then

      do i=iuw, 0
         al(i) = al(im+i)
         ar(i) = ar(im+i)
         a6(i) = a6(im+i)
      enddo

      do i=im+1, iue
         al(i) = al(i-im)
         ar(i) = ar(i-im)
         a6(i) = a6(i-im)
      enddo

      do i=1,im
            iu = c(i)
            ru = c(i) - iu
         if(c(i) .gt. 0.) then
            itmp = i - iu - 1
            isave(i) = itmp + 1
            fx(i) = ru*(ar(itmp)+0.5*ru*(al(itmp)-ar(itmp) +     &
                        a6(itmp)*(1.-r23*ru)) )
         else
            itmp = i - iu
            isave(i) = itmp - 1
            fx(i) = ru*(al(itmp)-0.5*ru*(ar(itmp)-al(itmp) +     &
                        a6(itmp)*(1.+r23*ru)) )
         endif
      enddo

  else
         al(0) = al(im)
         ar(0) = ar(im)
         a6(0) = a6(im)
      do i=1,im
         if(c(i) .gt. 0.) then
            fx(i) = ar(i-1) + 0.5*c(i)*(al(i-1) - ar(i-1) +   &
                    a6(i-1)*(1.-r23*c(i)) )
      else
            fx(i) = al(i) - 0.5*c(i)*(ar(i) - al(i) +         &
                    a6(i)*(1.+r23*c(i)))
      endif
            fx(i) = mfx(i) * fx(i)
      enddo
  endif
  return
!EOC
 end subroutine fxppm
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: steepx
!
! !INTERFACE: 
 subroutine  steepx(im, p, al, dm)
!-----------------------------------------------------------------------

! !USES:
 use precision
 implicit none

! !INPUT PARAMETERS:
 integer im
 real (r8)  p(-im/3:im+im/3)
 real (r8) dm(-im/3:im+im/3)

! !INPUT/OUTPUT PARAMETERS:
 real (r8) al(im)

! !DESCRIPTION:
!   
!
! !REVISION HISTORY:
!   99.01.01 Lin      Creation
!   01.03.27 Sawyer   Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
 integer i
 real (r8) r3
 parameter ( r3 = 1./3. )

 real (r8) dh(0:im)
 real (r8) d2(0:im+1)
 real (r8) eta(0:im)
 real (r8) xxx, bbb, ccc

   do i=0,im
      dh(i) = p(i+1) - p(i)
   enddo

! Needs dh(0:im)
   do i=1,im
      d2(i) = dh(i) - dh(i-1)
   enddo
   d2(0) = d2(im)
   d2(im+1) = d2(1)

! needs p(-1:im+2), d2(0:im+1)
   do i=1,im
      if( d2(i+1)*d2(i-1).lt.0. .and. p(i+1).ne.p(i-1) ) then
          xxx    = 1. - 0.5 * ( p(i+2) - p(i-2) ) / ( p(i+1) - p(i-1) )
          eta(i) = max(0., min(xxx, 0.5) )
      else
          eta(i) = 0.
      endif
    enddo

    eta(0) = eta(im)

! needs eta(0:im), dh(0:im-1), dm(0:im)
   do i=1,im
      bbb = ( 2.*eta(i  ) - eta(i-1) ) * dm(i-1) 
      ccc = ( 2.*eta(i-1) - eta(i  ) ) * dm(i  ) 
      al(i) = al(i) + 0.5*( eta(i-1) - eta(i)) * dh(i-1) + (bbb - ccc) * r3
   enddo
   return
!EOC
 end subroutine steepx
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: lmppm
!
! !INTERFACE: 
 subroutine lmppm(dm, a6, ar, al, p, im, lmt)
!-----------------------------------------------------------------------

 use precision
 implicit none

! !INPUT PARAMETERS:
 integer im   ! Total longitudes
 integer lmt  ! LMT = 0: full monotonicity
              ! LMT = 1: Improved and simplified full monotonic constraint
              ! LMT = 2: positive-definite constraint
              ! LMT = 3: Quasi-monotone constraint
 real(r8) p(im)
 real(r8) dm(im)

! !OUTPUT PARAMETERS:
 real(r8) a6(im)
 real(r8) ar(im)
 real(r8) al(im)

! !DESCRIPTION:
!
!
! !REVISION HISTORY:
!   99.01.01 Lin      Creation
!   01.03.27 Sawyer   Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
 real (r8) r12
 parameter ( r12 = 1./12. )

 real (r8) da1, da2, fmin, a6da
 real (r8) dr, dl

 integer i

! LMT = 0: full monotonicity
! LMT = 1: Improved and simplified full monotonic constraint
! LMT = 2: positive-definite constraint
! LMT = 3: Quasi-monotone constraint

  if( lmt == 0 ) then

! Full constraint
  do i=1,im
     if(dm(i) .eq. 0.) then
         ar(i) = p(i)
         al(i) = p(i)
         a6(i) = 0.
     else
         da1  = ar(i) - al(i)
         da2  = da1**2
         a6da = a6(i)*da1
         if(a6da .lt. -da2) then
            a6(i) = 3.*(al(i)-p(i))
            ar(i) = al(i) - a6(i)
         elseif(a6da .gt. da2) then
            a6(i) = 3.*(ar(i)-p(i))
            al(i) = ar(i) - a6(i)
         endif
     endif
  enddo

  elseif( lmt == 1 ) then

! Improved (Lin 2001?) full constraint
      do i=1,im
           da1 = dm(i) + dm(i)
            dl = sign(min(abs(da1),abs(al(i)-p(i))), da1)
            dr = sign(min(abs(da1),abs(ar(i)-p(i))), da1)
         ar(i) = p(i) + dr
         al(i) = p(i) - dl
         a6(i) = 3.*(dl-dr)
      enddo

  elseif( lmt == 2 ) then
! Positive definite constraint
      do 250 i=1,im
      if(abs(ar(i)-al(i)) .ge. -a6(i)) go to 250
      fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12
      if(fmin.ge.0.) go to 250
      if(p(i).lt.ar(i) .and. p(i).lt.al(i)) then
            ar(i) = p(i)
            al(i) = p(i)
            a6(i) = 0.
      elseif(ar(i) .gt. al(i)) then
            a6(i) = 3.*(al(i)-p(i))
            ar(i) = al(i) - a6(i)
      else
            a6(i) = 3.*(ar(i)-p(i))
            al(i) = ar(i) - a6(i)
      endif
250   continue

  elseif(lmt .eq. 3) then
! Quasi-monotone constraint
      do i=1,im
         da1 = 4.*dm(i)
          dl = sign(min(abs(da1),abs(al(i)-p(i))), da1)
          dr = sign(min(abs(da1),abs(ar(i)-p(i))), da1)
         ar(i) = p(i) + dr
         al(i) = p(i) - dl
         a6(i) = 3.*(dl-dr)
      enddo
  endif
  return
!EOC
 end subroutine lmppm
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
!BOP
! !IROUTINE: huynh --- Enforce Huynh's 2nd constraint in 1D periodic domain
!
! !INTERFACE: 
 subroutine huynh(im, ar, al, p, d2, d1)
!-----------------------------------------------------------------------

! !USES:

 use precision
 implicit none

! !INPUT PARAMETERS:
 integer im
 real(r8)  p(im)

! !OUTPUT PARAMETERS:
 real(r8) ar(im)
 real(r8) al(im)
 real(r8) d2(im)
 real(r8) d1(im)

! !DESCRIPTION:
!
!   Enforce Huynh's 2nd constraint in 1D periodic domain
!
! !REVISION HISTORY:
!   99.01.01 Lin      Creation
!   01.03.27 Sawyer   Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
 integer  i
 real(r8) pmp
 real(r8) lac
 real(r8) pmin
 real(r8) pmax

! Compute d1 and d2
      d1(1) = p(1) - p(im)
      do i=2,im
         d1(i) = p(i) - p(i-1)
      enddo

      do i=1,im-1
         d2(i) = d1(i+1) - d1(i)
      enddo
      d2(im) = d1(1) - d1(im)

! Constraint for AR
!            i = 1
         pmp   = p(1) + 2.0 * d1(1)
         lac   = p(1) + 0.5 * (d1(1)+d2(im)) + d2(im) 
         pmin  = min(p(1), pmp, lac)
         pmax  = max(p(1), pmp, lac)
         ar(1) = min(pmax, max(ar(1), pmin))

      do i=2, im
         pmp   = p(i) + 2.0*d1(i)
         lac   = p(i) + 0.5*(d1(i)+d2(i-1)) + d2(i-1)
         pmin  = min(p(i), pmp, lac)
         pmax  = max(p(i), pmp, lac)
         ar(i) = min(pmax, max(ar(i), pmin))
      enddo

! Constraint for AL
      do i=1, im-1
         pmp   = p(i) - 2.0*d1(i+1)
         lac   = p(i) + 0.5*(d2(i+1)-d1(i+1)) + d2(i+1)
         pmin  = min(p(i), pmp, lac)
         pmax  = max(p(i), pmp, lac)
         al(i) = min(pmax, max(al(i), pmin))
      enddo

! i=im
         i = im
         pmp    = p(im) - 2.0*d1(1)
         lac    = p(im) + 0.5*(d2(1)-d1(1)) + d2(1)
         pmin   = min(p(im), pmp, lac)
         pmax   = max(p(im), pmp, lac)
         al(im) = min(pmax, max(al(im), pmin))

! compute A6 (d2)
      do i=1, im
         d2(i) = 3.*(p(i)+p(i)  - (al(i)+ar(i)))
      enddo
    return
!EOC
 end subroutine huynh
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
!BOP
! !IROUTINE: ytp
!
! !INTERFACE: 
 subroutine ytp(im, jm, fy, q, c, yfx, 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                        !  order of subgrid dist
 integer iv                          !  Scalar=0, Vector=1
 real (r8) q(im,jfirst-ng:jlast+ng)       !  advected scalar N*jord S*jord
 real (r8) c(im,jfirst:jlast+1)           !  Courant   N (like FY)
 real (r8) yfx(im,jfirst:jlast+1)         !  Backgrond mass flux

! !OUTPUT PARAMETERS:
 real (r8) fy(im,jfirst:jlast+1)          !  Flux      N (see tp2c)

! !DESCRIPTION:
!     This routine calculates the flux FX.  The method chosen
!     depends on the order of the calculation JORD (currently
!     1, 2 or 3).  
!
! !CALLED FROM:
!     cd_core
!     tp2d
!
! !REVISION HISTORY:
!
!  SJL 99.04.13:  Delivery
!  WS  99.04.13:  Added jfirst:jlast concept
!  WS  99.04.21:  Removed j1 and j2 (j1=2, j2=jm-1 consistently)
!                 removed a6,ar,al from argument list
!  WS  99.04.27:  DM made local to this routine
!  WS  99.09.09:  Documentation; indentation; cleaning
!  WS  99.10.22:  Added NG as argument; pruned arrays
!  SJL 99.12.24:  Revised documentation; optimized for better cache usage
!  WS  00.05.14:  Renamed ghost indices as per Kevin's definitions
!
!EOP
!---------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
 integer i, j, jt
 integer js2g0, jn1g1

! work arrays (should pass in eventually for performance enhancement):
 real (r8) dm(im,jfirst-ng:jlast+ng)

!     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

   js2g0  = max(2,jfirst)       ! No ghosting
   jn1g1  = min(jm,jlast+1)     ! Ghost N*1
     
   if(jord == 1) then
!!!        do j=2,jm
        do j=js2g0,jn1g1
          do i=1,im
            jt = float(j) - c(i,j)
            fy(i,j) = q(i,jt)
          enddo
        enddo
   else

!
! YMIST requires q on NS;  Only call to YMIST here
!
        call ymist(im, jm, q, dm, ng, jord, iv, jfirst, jlast)

        if( abs(jord) .ge. 3 ) then
 
          call fyppm(c,q,dm,fy,im,jm,ng,jord,iv,jfirst,jlast)

        else
!
! JORD can either have the value 2 or -2 at this point
!
!!!          do j=2,jm
          do j=js2g0,jn1g1
            do i=1,im
              jt = float(j) - c(i,j)
              fy(i,j) = q(i,jt) + (sign(1.,c(i,j))-c(i,j))*dm(i,jt)
            enddo
          enddo
        endif
   endif

!!!      do j=2,jm
      do j=js2g0,jn1g1
        do i=1,im
          fy(i,j) = fy(i,j)*yfx(i,j)
        enddo
      enddo
    return
!EOC
 end subroutine ytp
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: ymist
!
! !INTERFACE: 
 subroutine ymist(im, jm, q, dm, ng, jord, iv, jfirst, jlast)
!-----------------------------------------------------------------------

! !USES:
 use precision
 implicit none

! !INPUT PARAMETERS:
 integer im, jm                      !  Dimensions
 integer jfirst, jlast               !  Latitude strip
 integer ng                          !  NS dependencies
 integer jord                        !  order of subgrid distribution
 integer iv                          !  Scalar (==0) Vector (==1)
 real (r8) q(im,jfirst-ng:jlast+ng)  !  transported scalar  N*ng S*ng

! !OUTPUT PARAMETERS:
 real (r8) dm(im,jfirst-ng:jlast+ng)      !  Slope only N*(ng-1) S*(ng-1) used

! !DESCRIPTION:
!     Calculate the slope of the pressure.  The number of ghost
!     latitudes (NG) depends on what method (JORD) will be used
!     subsequentally.    NG is equal to MIN(ABS(JORD),3).
!
! !CALLED FROM:
!     ytp
!
! !REVISION HISTORY:
!
!  SJL 99.04.13:  Delivery
!  WS  99.04.13:  Added jfirst:jlast concept
!  WS  99.09.09:  Documentation; indentation; cleaning
!  SJL 00.01.06:  Documentation
!  WS  00.05.14:  Renamed ghost indices as per Kevin's definitions
!
!EOP
!---------------------------------------------------------------------
!BOC

! Local variables

 integer i, j, jm1, im2, js2gng1, jn2gng1
 real (r8) qmax, qmin, tmp

    js2gng1 = max(2,   jfirst-ng+1)     !  Number needed on S
    jn2gng1 = min(jm-1,jlast+ng-1)      !  Number needed on N

    jm1 = jm - 1
    im2 = im / 2

!!!      do j=2,jm1
      do j=js2gng1,jn2gng1
        do i=1,im
           dm(i,j) = 0.25*(q(i,j+1) - q(i,j-1))
        enddo
      enddo

   if( iv == 0 ) then

        if ( jfirst == 1 ) then
! S pole
          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)