tp_core.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

module tp_core
!BOP
!
! !MODULE: tp_core --- Utilities for the transport core
!
!
! !PUBLIC MEMBER FUNCTIONS:
      public tp2c, tp2d, xtp, fxppm, xmist, steepx, lmppm
      public huynh, ytp, ymist, fyppm, tpcc, ycc
!
! !DESCRIPTION:
!
!      This module provides 
!
!      \begin{tabular}{|l|l|} \hline \hline
!       tp2c  &   \\ \hline
!       tp2d  &   \\ \hline 
!       xtp  &   \\ \hline 
!       fxppm  &   \\ \hline 
!       xmist  &   \\ \hline 
!       steepx  &   \\ \hline 
!       lmppm  &   \\ \hline 
!       huynh  &   \\ \hline 
!       ytp  &   \\ \hline 
!       ymist  &   \\ \hline 
!       fyppm  &   \\ \hline 
!       tpcc  &   \\ \hline 
!       ycc  &   \\ \hline
!                                \hline
!      \end{tabular}
!
! !REVISION HISTORY:
!   01.01.15   Lin        Routines coalesced into this module
!   01.03.26   Sawyer     Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------

CONTAINS

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: tp2c --- Perform transport on a C grid
!
! !INTERFACE: 
 subroutine tp2c(dh, va, h, crx, cry, im, jm,                            &
                 iord, jord, ng, fx, fy, ffsl,                           &
                 rcap, acosp, xfx, yfx, cosp, id, jfirst, jlast)
!-----------------------------------------------------------------------

 use precision
 implicit none

! !INPUT PARAMETERS:
   integer im, jm                  ! Dimensions
   integer jfirst, jlast           ! Latitude strip
   integer iord, jord              ! Interpolation order in x,y
   integer ng                      ! Max. NS dependencies
   integer id                      ! density (0)  (mfx = C)
   real (r8) rcap                  ! Ask S.-J. (polar constant?)
   real (r8) acosp(jm)             ! Ask S.-J. (difference to cosp??)
   logical ffsl(jm)                ! Use flux-form semi-Lagrangian trans.?
                                      ! (N*NG S*NG)
   real (r8) cosp(jm)                   ! Critical angle
   real (r8) va(im,jfirst:jlast)        ! Courant  (unghosted)
   real (r8) h(im,jfirst-ng:jlast+ng)   ! Pressure ( N*NG S*NG )
   real (r8) crx(im,jfirst-ng:jlast+ng) ! Ask S.-J. ( N*NG S*NG )
   real (r8) cry(im,jfirst:jlast+1)     ! Ask S.-J. ( N like FY )
   real (r8) xfx(im,jfirst:jlast)       ! Ask S.-J. ( unghosted like FX )
   real (r8) yfx(im,jfirst:jlast+1)     ! Ask S.-J. ( N like FY )

! !OUTPUT PARAMETERS:
   real (r8) dh(im,jfirst:jlast)        ! Ask S.-J. ( unghosted )
   real (r8) fx(im,jfirst:jlast)        ! Flux in x ( unghosted )
   real (r8) fy(im,jfirst:jlast+1)      ! Flux in y ( N, see tp2c )

! !DESCRIPTION:
!     Perform transport on a C grid.   The number of ghost
!     latitudes (NG) depends on what method (JORD) will be used
!     subsequentally.    NG is equal to MIN(ABS(JORD),3).
!     Ask S.-J. how exactly this differs from TP2C.
!
! !REVISION HISTORY:
!
!EOP
!-----------------------------------------------------------------------
!BOC
   integer i, j, js2g0, jn2g0
   real (r8)  sum1

   js2g0  = max(2,jfirst)          !  No ghosting
   jn2g0  = min(jm-1,jlast)        !  No ghosting

   call tp2d(va, h, crx, cry, im, jm, iord, jord, ng,fx, fy, ffsl,    &
             xfx, yfx, cosp, id, jfirst, jlast)

   do j=js2g0,jn2g0
      do i=1,im-1
         dh(i,j) = fx(i,j) - fx(i+1,j) + (fy(i,j)-fy(i,j+1))*acosp(j)
      enddo
   enddo

   do j=js2g0,jn2g0
      dh(im,j) = fx(im,j) - fx(1,j) + (fy(im,j)-fy(im,j+1))*acosp(j)
   enddo

! Poles
   if ( jfirst ==  1 ) then
!       sum1 = - SUM( fy(1:im, 2) ) * rcap
        sum1 = 0.
        do i=1,im
          sum1 = sum1 + fy(i,2)
        enddo
          sum1 = -sum1*rcap
        do i=1,im
          dh(i, 1) = sum1
        enddo
   endif
   
   if ( jlast == jm ) then
!       sum1 = SUM( fy(1:im,jm) ) * rcap
        sum1 = 0.
        do i=1,im
          sum1 = sum1 + fy(i,jm)
        enddo
          sum1 = sum1*rcap
        do i=1,im
          dh(i,jm) = sum1
        enddo
   endif
   return
!EOC
 end subroutine tp2c
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: tp2d --- Perform transport on a D grid
!
! !INTERFACE: 
 subroutine tp2d(va, q, crx, cry, im, jm, iord, jord, ng, fx, fy,        &
                 ffsl, xfx, yfx, cosp, id, jfirst, jlast)
!-----------------------------------------------------------------------
! !USES:

 use precision
 implicit none

! !INPUT PARAMETERS:
   integer im, jm                    ! Dimensions
   integer jfirst, jlast             ! Latitude strip
   integer iord, jord                ! Interpolation order in x,y
   integer ng                        ! Max. NS dependencies
   integer id                        ! density (0)  (mfx = C)
                                     ! mixing ratio (1) (mfx = mass flux)
   logical ffsl(jm)                  ! Use flux-form semi-Lagrangian trans.?
                                     ! ghosted N*ng S*ng
   real (r8) cosp(jm)                     ! Critical angle
   real (r8) va(im,jfirst:jlast)          ! Courant  (unghosted)
   real (r8) q(im,jfirst-ng:jlast+ng)     ! transported scalar ( N*NG S*NG )
   real (r8) crx(im,jfirst-ng:jlast+ng)   ! Ask S.-J. ( N*NG S*NG )
   real (r8) cry(im,jfirst:jlast+1)       ! Ask S.-J. ( N like FY )
   real (r8) xfx(im,jfirst:jlast)         ! Ask S.-J. ( unghosted like FX )
   real (r8) yfx(im,jfirst:jlast+1)       ! Ask S.-J. ( N like FY )

! !OUTPUT PARAMETERS:
   real (r8) fx(im,jfirst:jlast)          ! Flux in x ( unghosted )
   real (r8) fy(im,jfirst:jlast+1)        ! Flux in y ( N, see tp2c )

! !DESCRIPTION:
!     Perform transport on a D grid.   The number of ghost
!     latitudes (NG) depends on what method (JORD) will be used
!     subsequentally.    NG is equal to MIN(ABS(JORD),3).
!
!
! !REVISION HISTORY:
!   WS  99.04.13:  Added jfirst:jlast concept
!       99.04.21:  Removed j1 and j2 (j1=2, j2=jm-1 consistently)
!       99.04.27:  Removed dc, wk2 as arguments (local to YTP)
!       99.04.27:  Removed adx as arguments (local here)
!   SJL 99.07.26:  ffsl flag added
!   WS  99.09.07:  Restructuring, cleaning, documentation
!   WS  99.10.22:  NG now argument; arrays pruned
!   WS  00.05.14:  Renamed ghost indices as per Kevin's definitions
!
!EOP
!---------------------------------------------------------------------
!BOC

! Local:
   integer i, j, iad, jp, js2g0, js2gng, jn2g0, jn2gng
   real (r8) adx(im,jfirst-ng:jlast+ng)
   real (r8) wk1(im)
   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)

! Number of ghost latitudes
    js2g0  = max(2,jfirst)          !  No ghosting
    js2gng = max(2,jfirst-ng)       !  Number needed on S
    jn2g0  = min(jm-1,jlast)        !  No ghosting
    jn2gng = min(jm-1,jlast+ng)     !  Number needed on N
    iad = 1

!!!      do j=2,jm-1
    do j=js2gng,jn2gng               !  adx needed on N*ng S*ng

       call xtp(im,  ffsl(j), wk1, q(1,j),                &
                crx(1,j), iad, crx(1,j), cosp(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

! WS 99.09.07 : Split up north and south pole

     if ( jfirst == 1 ) then
        do i=1,im 
          adx(i, 1) = q(i,1)
        enddo
     endif 
     if ( jlast == jm ) then
        do i=1,im 
          adx(i,jm) = q(i,jm)
        enddo
     endif

     call ytp(im,jm,fy, adx,cry,yfx,ng,jord,0,jfirst,jlast)

!!!   do j=2,jm-1
      do j=js2g0,jn2g0
        do i=1,im
           jp = j-va(i,j)
           wk1(i) = q(i,j) +0.5*va(i,j)*(q(i,jp)-q(i,jp+1))
        enddo

        call xtp(im,  ffsl(j), fx(1,j), wk1,                  &
                 crx(1,j), iord, xfx(1,j), cosp(j), id,       &
                 dm, qtmp, al, ar, a6)
      enddo
    return
!EOC
 end subroutine tp2d
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: xtp
!
! !INTERFACE: 
 subroutine xtp(im, ffsl,  fx,  q,  c,  iord,  mfx,            &
                cosa, id, dm, qtmp, al, ar, a6)
!-----------------------------------------------------------------------

 use precision
 implicit none
 
! !INPUT PARAMETERS:
   integer id               ! ID = 0: density (mfx = C)
                            ! ID = 1: mixing ratio (mfx is mass flux)

   integer im               ! Total longitudes
   real (r8) c(im)          ! Courant numbers
   real (r8) q(im)
   real (r8) mfx(im)
   logical ffsl
   integer iord
   real (r8) cosa

! !INPUT/OUTPUT PARAMETERS:
   real (r8) qtmp(-im/3:im+im/3)   ! Input work arrays:
   real (r8)   dm(-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)

! !DESCRIPTION:
!   
!
! !REVISION HISTORY:
!   99.01.01 Lin      Creation
!   01.03.27 Sawyer   Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC

! Local:
   real (r8)       cos_upw               !critical cosine for upwind
   real (r8)       cos_van               !critical cosine for van Leer
   real (r8)       cos_ppm               !critical cosine for ppm

   parameter (cos_upw = 0.05)       !roughly at 87 deg.
   parameter (cos_van = 0.25)       !roughly at 75 deg.
   parameter (cos_ppm = 0.25)

   integer i, imp
   real (r8) qmax, qmin
   real (r8) rut, tmp
   integer iu, itmp, ist
   integer isave(im)
   integer iuw, iue

   imp = im + 1

   do i=1,im
      qtmp(i) = q(i)
   enddo

   if( ffsl ) then
! Flux-Form Semi-Lagrangian transport

! Figure out ghost zone for the western edge:
      iuw =  -c(1)
      iuw = min(0, iuw)
 
      do i=iuw, 0
         qtmp(i) = q(im+i)
      enddo 

! Figure out ghost zone for the eastern edge:
      iue = im - c(im)
      iue = max(imp, iue)

      do i=imp, iue
         qtmp(i) = q(i-im)
      enddo

      if( iord == 1 .or. cosa < cos_upw) then
      do i=1,im
        iu = c(i)
      if(c(i) .le. 0.) then
        itmp = i - iu
        isave(i) = itmp - 1
      else
        itmp = i - iu - 1
        isave(i) = itmp + 1
      endif
        fx(i) = (c(i)-iu) * qtmp(itmp)
      enddo
      else

      do i=1,im
! 2nd order slope
         tmp = 0.25*(qtmp(i+1) - qtmp(i-1))
         qmax = max(qtmp(i-1), qtmp(i), qtmp(i+1)) - qtmp(i)
         qmin = qtmp(i) - min(qtmp(i-1), qtmp(i), qtmp(i+1))
         dm(i) = sign(min(abs(tmp),qmax,qmin), tmp)
      enddo

 
      do i=iuw, 0
         dm(i) = dm(im+i)
      enddo 

      do i=imp, iue
         dm(i) = dm(i-im)
      enddo

      if(iord .ge. 3 .and. cosa .gt. cos_ppm) then
         call fxppm(im, c, mfx, qtmp, dm, fx, iord, al, ar, a6,         &
                    iuw, iue, ffsl, isave)
      else
      do i=1,im
            iu  = c(i)
            rut = c(i) - iu
         if(c(i) .le. 0.) then
            itmp = i - iu
            isave(i) = itmp - 1
            fx(i) = rut*(qtmp(itmp)-dm(itmp)*(1.+rut))
         else
            itmp = i - iu - 1
            isave(i) = itmp + 1
            fx(i) = rut*(qtmp(itmp)+dm(itmp)*(1.-rut))
         endif
      enddo
      endif

      endif

      do i=1,im
      if(c(i) .ge. 1.) then
        do ist = isave(i),i-1
           fx(i) = fx(i) + qtmp(ist)
        enddo
      elseif(c(i) .le. -1.) then
        do ist = i,isave(i)
           fx(i) = fx(i) - qtmp(ist)
        enddo
      endif
      enddo

      if(id .ne. 0) then
         do i=1,im
            fx(i) =  fx(i)*mfx(i)
         enddo
      endif

   else
! Regular PPM (Eulerian without FFSL extension)

      qtmp(imp) = q(1)
      qtmp(  0) = q(im)

      if(iord == 1 .or. cosa < cos_upw) then
         do i=1,im
            iu = float(i) - c(i)
            fx(i) = mfx(i)*qtmp(iu)
         enddo
      else

         qtmp(-1)    = q(im-1)
         qtmp(imp+1) = q(2)

         if(iord > 0 .or. cosa < cos_van) then
            call xmist(im, qtmp, dm, 2)
         else
            call xmist(im, qtmp, dm, iord)
         endif

         dm(0) = dm(im)

         if( abs(iord).eq.2 .or. cosa .lt. cos_van ) then
            do i=1,im
               iu = float(i) - c(i)
               fx(i) =  mfx(i)*(qtmp(iu)+dm(iu)*(sign(1.,c(i))-c(i)))

!              if(c(i) .le. 0.) then
!                 fx(i) = qtmp(i) - dm(i)*(1.+c(i))
!              else
!                 fx(i) = qtmp(i-1) + dm(i-1)*(1.-c(i))
!              endif
!                 fx(i) = fx(i)*mfx(i)

            enddo
         else
            call fxppm(im, c, mfx, qtmp, dm, fx, iord, al, ar, a6,       &
                       iuw, iue, ffsl, isave)
         endif
      endif

   endif
   return
!EOC
 end subroutine xtp
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: xmist
!
! !INTERFACE: 
 subroutine xmist(im,  q,  dm,  id)
!-----------------------------------------------------------------------

 use precision
 implicit none

! !INPUT PARAMETERS:
 integer im                   ! Total number of longitudes
 integer id                   ! ID = 0: density (mfx = C)
                              ! ID = 1: mixing ratio (mfx is mass flux)
 real(r8)  q(-im/3:im+im/3)   ! Input latitude 

! !OUTPUT PARAMETERS:
 real(r8) dm(-im/3:im+im/3)   ! 

! !DESCRIPTION:
!   
!
! !REVISION HISTORY:
!   99.01.01 Lin      Creation
!   01.03.27 Sawyer   Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC

 real(r8) r24
 parameter( r24 = 1./24.)

 integer i
 real(r8) qmin, qmax

    if(id .le. 2) then
       do i=1,im
          dm(i) = r24*(8.*(q(i+1) - q(i-1)) + q(i-2) - q(i+2))
       enddo
    else
       do i=1,im
          dm(i) = 0.25*(q(i+1) - q(i-1))
       enddo
    endif

    if( id < 0 ) return

! Apply monotonicity constraint (Lin et al. 1994, MWR)
      do i=1,im
         qmax = max( q(i-1), q(i), q(i+1) ) - q(i)
         qmin = q(i) - min( q(i-1), q(i), q(i+1) )
         dm(i) = sign( min(abs(dm(i)), qmax, qmin), dm(i) )
      enddo
  return
!EOC
 end subroutine xmist
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: fxppm
!
! !INTERFACE: 
 subroutine fxppm(im, c, mfx,  p, dm, fx, iord, al, ar, a6,        &
                  iuw, iue, ffsl, isave)
!-----------------------------------------------------------------------
!
! !USES:
 use precision
 implicit none

! !INPUT PARAMETERS:
 integer im, iord
 real (r8)  c(im)
 real (r8) p(-im/3:im+im/3)
 real (r8) dm(-im/3:im+im/3)
 real (r8) mfx(im)
 integer iuw, iue
 logical ffsl
 integer isave(im)

! !INPUT/OUTPUT PARAMETERS:
 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)

! !DESCRIPTION:
!   
!
! !REVISION HISTORY:
!   99.01.01 Lin      Creation
!   01.03.27 Sawyer   Additional ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------