pkg/generic_advdiff/gad_fluxlimit_impl_r.F

7b4413ef94 Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: GAD_FLUXLIMIT_IMPL_R
                0005 C     !INTERFACE:
ec0db5c1b3 Jean*0006       SUBROUTINE GAD_FLUXLIMIT_IMPL_R(
                0007      I           bi,bj,k, iMin,iMax,jMin,jMax,
                0008      I           deltaTarg, rTrans, recip_hFac, tFld,
7b4413ef94 Jean*0009      O           a3d, b3d, c3d,
                0010      I           myThid )
                0011
1b5fb69d21 Ed H*0012 C     !DESCRIPTION:
                0013 C     Compute matrix element to solve vertical advection implicitly
381eb13d88 Jean*0014 C      using flux--limiter advection scheme.
                0015 C     Method:
ec0db5c1b3 Jean*0016 C      contribution of vertical transport at interface k is added
381eb13d88 Jean*0017 C      to matrix lines k and k-1.
7b4413ef94 Jean*0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021
                0022 C     == Global variables ===
                0023 #include "SIZE.h"
                0024 #include "GRID.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027
                0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     == Routine Arguments ==
2132daf4a7 Jean*0030 C     bi, bj       :: tile indices
1b5fb69d21 Ed H*0031 C     k            :: vertical level
2132daf4a7 Jean*0032 C     iMin, iMax   :: computation domain
                0033 C     jMin, jMax   :: computation domain
1b5fb69d21 Ed H*0034 C     deltaTarg    :: time step
                0035 C     rTrans       :: vertical volume transport
ec0db5c1b3 Jean*0036 C     recip_hFac   :: inverse of cell open-depth factor
1b5fb69d21 Ed H*0037 C     tFld         :: tracer field
ec0db5c1b3 Jean*0038 C     a3d          :: lower diagonal of the tridiagonal matrix
                0039 C     b3d          :: main  diagonal of the tridiagonal matrix
                0040 C     c3d          :: upper diagonal of the tridiagonal matrix
1b5fb69d21 Ed H*0041 C     myThid       :: thread number
7b4413ef94 Jean*0042       INTEGER bi,bj,k
                0043       INTEGER iMin,iMax,jMin,jMax
ec0db5c1b3 Jean*0044       _RL     deltaTarg(Nr)
                0045       _RL     rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0047       _RL     tFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0048       _RL     a3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0049       _RL     b3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0050       _RL     c3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
7b4413ef94 Jean*0051       INTEGER myThid
                0052
                0053 C     == Local Variables ==
2132daf4a7 Jean*0054 C     i, j         :: loop indices
1b5fb69d21 Ed H*0055 C     kp1          :: =min( k+1 , Nr )
                0056 C     km1          :: =max( k-1 , 1 )
                0057 C     km2          :: =max( k-2 , 1 )
                0058 C     Cr           :: slope ratio
2132daf4a7 Jean*0059 C     Rjm, Rj, Rjp :: differences at k-1,k,k+1
1b5fb69d21 Ed H*0060 C     w_CFL        :: Courant-Friedrich-Levy number
                0061 C     upwindFac    :: upwind factor
                0062 C     rCenter      :: centered contribution
                0063 C     rUpwind      :: upwind   contribution
2132daf4a7 Jean*0064       INTEGER i, j, kp1, km1, km2
                0065       _RL Cr, Rjm, Rj, Rjp, w_CFL
7b4413ef94 Jean*0066       _RL upwindFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL rCenter, rUpwind
9de7a55d87 Jean*0068       _RL deltaTcfl
2132daf4a7 Jean*0069       _RL CrMax
                0070       PARAMETER( CrMax = 1.D+6 )
7b4413ef94 Jean*0071
                0072 C Statement function provides Limiter(Cr)
                0073 #include "GAD_FLUX_LIMITER.h"
                0074 CEOP
                0075
2132daf4a7 Jean*0076       km2 = MAX(1,k-2)
                0077       km1 = MAX(1,k-1)
                0078       kp1 = MIN(Nr,k+1)
7b4413ef94 Jean*0079
9de7a55d87 Jean*0080 C--   process interior interface only:
                0081       IF ( k.GT.1 .AND. k.LE.Nr ) THEN
7b4413ef94 Jean*0082
                0083 C--   Compute the upwind fraction:
9de7a55d87 Jean*0084        deltaTcfl = deltaTarg(k)
7b4413ef94 Jean*0085        DO j=jMin,jMax
                0086         DO i=iMin,iMax
9de7a55d87 Jean*0087          w_CFL = deltaTcfl*rTrans(i,j)*recip_rA(i,j,bi,bj)*recip_drC(k)
a7ec469280 Jean*0088      &                    *recip_deepFac2F(k)*recip_rhoFacF(k)
2132daf4a7 Jean*0089          Rjp = (tFld(i,j,kp1)-tFld(i,j,k)  )*maskC(i,j,kp1,bi,bj)
                0090          Rj  = (tFld(i,j,k)  -tFld(i,j,km1))
                0091          Rjm = (tFld(i,j,km1)-tFld(i,j,km2))*maskC(i,j,km2,bi,bj)
7b4413ef94 Jean*0092
2132daf4a7 Jean*0093          IF ( rTrans(i,j).LT.zeroRL ) THEN
                0094            Cr = Rjm
7b4413ef94 Jean*0095          ELSE
2132daf4a7 Jean*0096            Cr = Rjp
7b4413ef94 Jean*0097          ENDIF
2132daf4a7 Jean*0098          IF ( ABS(Rj)*CrMax .LE. ABS(Cr) ) THEN
                0099            Cr = SIGN( CrMax, Cr )*SIGN( oneRL, Rj )
                0100          ELSE
                0101            Cr = Cr/Rj
                0102          ENDIF
                0103
                0104 C        calculate upwind fraction using Limiter Function:
                0105          upwindFac(i,j) = 1. _d 0
                0106      &                  - Limiter(Cr) * ( 1. _d 0 + ABS(w_CFL) )
                0107          upwindFac(i,j) = MAX( -1. _d 0, upwindFac(i,j) )
7b4413ef94 Jean*0108         ENDDO
                0109        ENDDO
                0110
                0111 C--    Add centered & upwind contributions
                0112        DO j=jMin,jMax
                0113          DO i=iMin,iMax
bb6c554092 Jean*0114            rCenter = 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
                0115            rUpwind = ABS(rCenter)*upwindFac(i,j)
7b4413ef94 Jean*0116            a3d(i,j,k)   = a3d(i,j,k)
bb6c554092 Jean*0117      &                  - (rCenter+rUpwind)*deltaTarg(k)
ec0db5c1b3 Jean*0118      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0119      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
7b4413ef94 Jean*0120            b3d(i,j,k)   = b3d(i,j,k)
bb6c554092 Jean*0121      &                  - (rCenter-rUpwind)*deltaTarg(k)
ec0db5c1b3 Jean*0122      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0123      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
7b4413ef94 Jean*0124            b3d(i,j,k-1) = b3d(i,j,k-1)
bb6c554092 Jean*0125      &                  + (rCenter+rUpwind)*deltaTarg(k-1)
ec0db5c1b3 Jean*0126      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0127      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
7b4413ef94 Jean*0128            c3d(i,j,k-1) = c3d(i,j,k-1)
bb6c554092 Jean*0129      &                  + (rCenter-rUpwind)*deltaTarg(k-1)
ec0db5c1b3 Jean*0130      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0131      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
7b4413ef94 Jean*0132          ENDDO
                0133        ENDDO
                0134
9de7a55d87 Jean*0135 C--   process interior interface only: end
                0136       ENDIF
                0137
7b4413ef94 Jean*0138       RETURN
                0139       END