pkg/generic_advdiff/gad_dst3fl_impl_r.F

381eb13d88 Jean*0001 #include "GAD_OPTIONS.h"
113b21cd45 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
381eb13d88 Jean*0005
                0006 CBOP
                0007 C     !ROUTINE: GAD_DST3FL_IMPL_R
                0008 C     !INTERFACE:
                0009       SUBROUTINE GAD_DST3FL_IMPL_R(
                0010      I           bi,bj,k, iMin,iMax,jMin,jMax,
ec0db5c1b3 Jean*0011      I           deltaTarg, rTrans, recip_hFac, tFld,
381eb13d88 Jean*0012      O           a5d, b5d, c5d, d5d, e5d,
                0013      I           myThid )
                0014
                0015 C     !DESCRIPTION:
                0016
                0017 C     Compute matrix element to solve vertical advection implicitly
                0018 C     using 3rd order Direct Space and Time (DST) advection scheme
                0019 C           with Flux-Limiter.
                0020 C     Method:
                0021 C      contribution of vertical transport at interface k is added
                0022 C      to matrix lines k and k-1
                0023
                0024 C     !USES:
                0025       IMPLICIT NONE
                0026
                0027 C     == Global variables ===
                0028 #include "SIZE.h"
                0029 #include "GRID.h"
                0030 #include "EEPARAMS.h"
                0031 #include "PARAMS.h"
                0032 #include "GAD.h"
                0033
                0034 C     !INPUT/OUTPUT PARAMETERS:
                0035 C     == Routine Arguments ==
ec0db5c1b3 Jean*0036 C     bi,bj        :: tile indices
                0037 C     k            :: vertical level
                0038 C     iMin,iMax    :: computation domain
                0039 C     jMin,jMax    :: computation domain
                0040 C     deltaTarg    :: time step
                0041 C     rTrans       :: vertical volume transport
                0042 C     recip_hFac   :: inverse of cell open-depth factor
                0043 C     tFld         :: tracer field
                0044 C     a5d          :: 2nd  lower diag of pentadiagonal matrix
                0045 C     b5d          :: 1rst lower diag of pentadiagonal matrix
                0046 C     c5d          :: main diag       of pentadiagonal matrix
                0047 C     d5d          :: 1rst upper diag of pentadiagonal matrix
                0048 C     e5d          :: 2nd  upper diag of pentadiagonal matrix
                0049 C     myThid       :: thread number
381eb13d88 Jean*0050       INTEGER bi,bj,k
                0051       INTEGER iMin,iMax,jMin,jMax
ec0db5c1b3 Jean*0052       _RL     deltaTarg(Nr)
                0053       _RL     rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0054       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0055       _RL     tFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0056       _RL     a5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0057       _RL     b5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0058       _RL     c5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0059       _RL     d5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0060       _RL     e5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
381eb13d88 Jean*0061       INTEGER myThid
                0062
                0063 C     == Local Variables ==
ec0db5c1b3 Jean*0064 C     i,j          :: loop indices
                0065 C     kp1          :: =min( k+1 , Nr )
                0066 C     km2          :: =max( k-2 , 1 )
                0067 C     wCFL         :: Courant-Friedrich-Levy number
                0068 C     lowFac       :: low  order term factor
                0069 C     highFac      :: high order term factor
                0070 C     rCenter      :: centered contribution
                0071 C     rUpwind      :: upwind   contribution
                0072 C     rC4km, rC4kp :: high order contributions
381eb13d88 Jean*0073       INTEGER i,j,kp1,km2
                0074       _RL wCFL, rCenter, rUpwind
                0075       _RL lowFac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076       _RL highFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077       _RL rC4km, rC4kp
                0078       _RL mskM, mskP, maskM2, maskP1
                0079       _RL Rj, Rjh, cL1, cH3, cM2, th1, th2
                0080       _RL deltaTcfl
7448700841 Mart*0081 #if ( defined ALLOW_AUTODIFF_TAMC && defined INCLUDE_IMPLVERTADV_CODE )
113b21cd45 Mart*0082       INTEGER ikey
                0083 #endif
381eb13d88 Jean*0084 CEOP
                0085
7448700841 Mart*0086 #if ( defined ALLOW_AUTODIFF_TAMC && defined INCLUDE_IMPLVERTADV_CODE )
113b21cd45 Mart*0087 CADJ INIT loctape_gad = COMMON, (sNx+2*OLx)*(sNy+2*OLy)
                0088 #endif
381eb13d88 Jean*0089 C--   process interior interface only:
                0090       IF ( k.GT.1 .AND. k.LE.Nr ) THEN
                0091
                0092        km2=MAX(1,k-2)
                0093        kp1=MIN(Nr,k+1)
                0094        maskP1 = 1. _d 0
                0095        maskM2 = 1. _d 0
                0096        IF ( k.LE.2 ) maskM2 = 0. _d 0
                0097        IF ( k.GE.Nr) maskP1 = 0. _d 0
                0098
                0099 C--   Compute the low-order term & high-order term fractions :
                0100        deltaTcfl = deltaTarg(k)
                0101 C     DST-3 Flux-Limiter Advection Scheme:
                0102 C-    Limiter: Psi=max(0,min(1,cL1+theta*cH1,theta*(1-cfl)/cfl) )
                0103 C              with theta=Rjh/Rj ;
                0104 C       is linearize arround the current value of theta(tFld) & cfl:
                0105 C       lowFac & highFac are set such as Psi*Rj = lowFac*Rj + highFac*Rjh
                0106        DO j=jMin,jMax
                0107          DO i=iMin,iMax
7448700841 Mart*0108 #if ( defined ALLOW_AUTODIFF_TAMC && defined INCLUDE_IMPLVERTADV_CODE )
113b21cd45 Mart*0109            ikey = i + (sNx+2*OLx)*(j-1)
                0110 #endif
381eb13d88 Jean*0111            wCFL = deltaTcfl*ABS(rTrans(i,j))
                0112      &           *recip_rA(i,j,bi,bj)*recip_drC(k)
a7ec469280 Jean*0113      &           *recip_deepFac2F(k)*recip_rhoFacF(k)
381eb13d88 Jean*0114            cL1 = (2. _d 0 -wCFL)*(1. _d 0 -wCFL)*oneSixth
                0115            cH3 = (1. _d 0 -wCFL*wCFL)*oneSixth
                0116 c          cM2 = (1. _d 0 - wCFL)/( wCFL +1. _d -20)
                0117            cM2 = (1. _d 0 + wCFL)/( wCFL +1. _d -20)
                0118
                0119            Rj =(tFld(i,j,k)  -tFld(i,j,k-1))
                0120            IF ( rTrans(i,j).GT.0. _d 0 ) THEN
                0121              Rjh = (tFld(i,j,k-1)-tFld(i,j,km2))*maskC(i,j,km2,bi,bj)
                0122            ELSE
                0123              Rjh = (tFld(i,j,kp1)-tFld(i,j,k)  )*maskC(i,j,kp1,bi,bj)
                0124            ENDIF
                0125            IF ( Rj*Rjh.LE.0. _d 0 ) THEN
                0126 C-         1rst case: theta < 0 (Rj & Rjh opposite sign) => Psi = 0
                0127              lowFac(i,j) = 0. _d 0
                0128              highFac(i,j)= 0. _d 0
                0129            ELSE
7448700841 Mart*0130 #if ( defined ALLOW_AUTODIFF_TAMC && defined INCLUDE_IMPLVERTADV_CODE )
113b21cd45 Mart*0131 C     These stores to local tape avoid recomputations and associated warnings
                0132 CADJ STORE Rjh, Rj = loctape_gad, key = ikey
                0133 #endif
381eb13d88 Jean*0134              Rj  = ABS(Rj)
                0135              Rjh = ABS(Rjh)
                0136              th1 = cL1*Rj+cH3*Rjh
                0137              th2 = cM2*Rjh
                0138             IF     ( th1.LE.th2 .AND. th1.LE.Rj ) THEN
                0139 C-          2nd case: cL1+theta*cH3 = min of the three = Psi
                0140              lowFac(i,j) = cL1
                0141              highFac(i,j)= cH3
                0142             ELSEIF ( th2.LT.th1 .AND. th2.LE.Rj ) THEN
                0143 C-          3rd case: theta*cM2 = min of the three = Psi
                0144              lowFac(i,j) = 0. _d 0
                0145              highFac(i,j)= cM2
                0146             ELSE
                0147 C-          4th case (Rj < th1 & Rj < th2) : 1 = min of the three = Psi
                0148              lowFac(i,j) = 1. _d 0
                0149              highFac(i,j)= 0. _d 0
                0150             ENDIF
                0151            ENDIF
                0152          ENDDO
                0153        ENDDO
                0154
                0155 C--    Add centered & upwind contributions
                0156        DO j=jMin,jMax
                0157          DO i=iMin,iMax
                0158            rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
                0159            mskM   = maskC(i,j,km2,bi,bj)*maskM2
                0160            mskP   = maskC(i,j,kp1,bi,bj)*maskP1
                0161            rUpwind= (0.5 _d 0 -lowFac(i,j))*ABS(rCenter)*2. _d 0
                0162            rC4km  = highFac(i,j)*(rCenter+ABS(rCenter))*mskM
                0163            rC4kp  = highFac(i,j)*(rCenter-ABS(rCenter))*mskP
                0164
                0165            a5d(i,j,k)   = a5d(i,j,k)
                0166      &                  + rC4km
                0167      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0168      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0169      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
381eb13d88 Jean*0170            b5d(i,j,k)   = b5d(i,j,k)
                0171      &                  - ( (rCenter+rUpwind) + rC4km )
                0172      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0173      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0174      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
381eb13d88 Jean*0175            c5d(i,j,k)   = c5d(i,j,k)
                0176      &                  - ( (rCenter-rUpwind) + rC4kp )
                0177      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0178      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0179      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
381eb13d88 Jean*0180            d5d(i,j,k)   = d5d(i,j,k)
                0181      &                  + rC4kp
                0182      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0183      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0184      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
381eb13d88 Jean*0185            b5d(i,j,k-1) = b5d(i,j,k-1)
                0186      &                  - rC4km
                0187      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0188      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0189      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
381eb13d88 Jean*0190            c5d(i,j,k-1) = c5d(i,j,k-1)
                0191      &                  + ( (rCenter+rUpwind) + rC4km )
                0192      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0193      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0194      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
381eb13d88 Jean*0195            d5d(i,j,k-1) = d5d(i,j,k-1)
                0196      &                  + ( (rCenter-rUpwind) + rC4kp )
                0197      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0198      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0199      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
381eb13d88 Jean*0200            e5d(i,j,k-1) = e5d(i,j,k-1)
                0201      &                  - rC4kp
                0202      &                   *deltaTarg(k-1)
a7ec469280 Jean*0203      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
                0204      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
381eb13d88 Jean*0205          ENDDO
                0206        ENDDO
                0207
                0208 C--   process interior interface only: end
                0209       ENDIF
                0210
                0211       RETURN
                0212       END