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!-----------------------------------------------------------------------