pkg/generic_advdiff/gad_u3c4_impl_r.F

7b4413ef94 Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
d285d86583 Jean*0004 C     !ROUTINE: GAD_U3C4_IMPL_R
7b4413ef94 Jean*0005 C     !INTERFACE:
381eb13d88 Jean*0006       SUBROUTINE GAD_U3C4_IMPL_R(
                0007      I           bi,bj,k, iMin,iMax,jMin,jMax,
ec0db5c1b3 Jean*0008      I           advectionScheme, deltaTarg, rTrans, recip_hFac,
7b4413ef94 Jean*0009      O           a5d, b5d, c5d, d5d, e5d,
                0010      I           myThid )
                0011
1b5fb69d21 Ed H*0012 C     !DESCRIPTION:
d285d86583 Jean*0013 C     Compute matrix element to solve vertical advection implicitly
                0014 C      using 3rd order upwind advection scheme,
381eb13d88 Jean*0015 C         or 3rd order Direct Space and Time advection scheme,
d285d86583 Jean*0016 C         or 4th order Centered advection scheme.
                0017 C     Method:
                0018 C      contribution of vertical transport at interface k is added
                0019 C      to matrix lines k and k-1
7b4413ef94 Jean*0020
                0021 C     !USES:
                0022       IMPLICIT NONE
                0023
                0024 C     == Global variables ===
                0025 #include "SIZE.h"
                0026 #include "GRID.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "GAD.h"
                0030
d285d86583 Jean*0031 C     !INPUT/OUTPUT PARAMETERS:
                0032 C     == Routine Arguments ==
ec0db5c1b3 Jean*0033 C     bi,bj        :: tile indices
                0034 C     k            :: vertical level
                0035 C     iMin,iMax    :: computation domain
                0036 C     jMin,jMax    :: computation domain
                0037 C  advectionScheme :: advection scheme to use
                0038 C     deltaTarg    :: time step
                0039 C     rTrans       :: vertical volume transport
                0040 C     recip_hFac   :: inverse of cell open-depth factor
                0041 C     a5d          :: 2nd  lower diag of pentadiagonal matrix
                0042 C     b5d          :: 1rst lower diag of pentadiagonal matrix
                0043 C     c5d          :: main diag       of pentadiagonal matrix
                0044 C     d5d          :: 1rst upper diag of pentadiagonal matrix
                0045 C     e5d          :: 2nd  upper diag of pentadiagonal matrix
                0046 C     myThid       :: thread number
7b4413ef94 Jean*0047       INTEGER bi,bj,k
                0048       INTEGER iMin,iMax,jMin,jMax
                0049       INTEGER advectionScheme
ec0db5c1b3 Jean*0050       _RL     deltaTarg(Nr)
                0051       _RL     rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0052       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0053       _RL     a5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0054       _RL     b5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0055       _RL     c5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0056       _RL     d5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0057       _RL     e5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
7b4413ef94 Jean*0058       INTEGER myThid
                0059
d285d86583 Jean*0060 C     == Local Variables ==
ec0db5c1b3 Jean*0061 C     i,j          :: loop indices
                0062 C     kp1          :: =min( k+1 , Nr )
                0063 C     km2          :: =max( k-2 , 1 )
                0064 C     rCenter      :: centered contribution
                0065 C     rUpwind      :: upwind   contribution
                0066 C     rC4km, rC4kp :: high order contribution
                0067 C     rHigh        :: high order term factor
7b4413ef94 Jean*0068       LOGICAL flagC4
                0069       INTEGER i,j,kp1,km2
64442af1fe Jean*0070 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
ec0db5c1b3 Jean*0071       _RL rC4km2D  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
757b06e91e Mart*0072       _RL rC4kp2D  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073       _RL rCenter2D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0074       _RL rUpwind2D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075 #endif
381eb13d88 Jean*0076       _RL wCFL, rCenter, rUpwind
                0077       _RL rC4km, rC4kp, rHigh
7b4413ef94 Jean*0078       _RL mskM, mskP, maskM2, maskP1
381eb13d88 Jean*0079       _RL deltaTcfl
7b4413ef94 Jean*0080 CEOP
                0081
9de7a55d87 Jean*0082 C--   process interior interface only:
                0083       IF ( k.GT.1 .AND. k.LE.Nr ) THEN
7b4413ef94 Jean*0084
                0085       km2=MAX(1,k-2)
                0086       kp1=MIN(Nr,k+1)
381eb13d88 Jean*0087       maskP1 = 1. _d 0
                0088       maskM2 = 1. _d 0
7b4413ef94 Jean*0089       IF ( k.LE.2 ) maskM2 = 0. _d 0
                0090       IF ( k.GE.Nr) maskP1 = 0. _d 0
381eb13d88 Jean*0091       flagC4 = advectionScheme.EQ.ENUM_CENTERED_4TH
7b4413ef94 Jean*0092      &         .AND. k.GT.2 .AND. k.LT.Nr
                0093
381eb13d88 Jean*0094 C--    Add centered, upwind and high-order contributions
                0095        deltaTcfl = deltaTarg(k)
64442af1fe Jean*0096 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
757b06e91e Mart*0097        DO j=jMin,jMax
                0098         DO i=iMin,iMax
ec0db5c1b3 Jean*0099          rCenter2D(i,j) =
757b06e91e Mart*0100      &        0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
                0101          mskM   = maskC(i,j,km2,bi,bj)*maskM2
                0102          mskP   = maskC(i,j,kp1,bi,bj)*maskP1
                0103          IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN
                0104           rUpwind2D(i,j) = 0. _d 0
                0105           rC4km2D  (i,j) = oneSixth*rCenter*mskM
                0106           rC4kp2D  (i,j) = oneSixth*rCenter*mskP
                0107          ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
                0108           wCFL = deltaTcfl*ABS(rTrans(i,j))
a7ec469280 Jean*0109      &            *recip_rA(i,j,bi,bj)*recip_drC(k)
                0110      &            *recip_deepFac2F(k)*recip_rhoFacF(k)
757b06e91e Mart*0111           rHigh  = (1. _d 0 -wCFL*wCFL)*oneSixth
                0112 c         rUpwind2D(i,j) = (2. _d 0*rHigh - wCFL)*ABS(rCenter)
                0113           rUpwind2D(i,j) = (2. _d 0*rHigh )*ABS(rCenter)
                0114           rC4km2D  (i,j) =  rHigh * (rCenter+ABS(rCenter))*mskM
                0115           rC4kp2D  (i,j) =  rHigh * (rCenter-ABS(rCenter))*mskP
                0116          ELSE
                0117           rUpwind2D(i,j) =  2. _d 0*oneSixth*ABS(rCenter)
                0118           rC4km2D  (i,j) = oneSixth*(rCenter+ABS(rCenter))*mskM
                0119           rC4kp2D  (i,j) = oneSixth*(rCenter-ABS(rCenter))*mskP
                0120          ENDIF
                0121         ENDDO
                0122        ENDDO
64442af1fe Jean*0123 #endif /* ALLOW_AUTODIFF and TARGET_NEC_SX */
7b4413ef94 Jean*0124        DO j=jMin,jMax
                0125          DO i=iMin,iMax
64442af1fe Jean*0126 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
757b06e91e Mart*0127            rC4km   = rC4km2D  (i,j)
                0128            rC4kp   = rC4kp2D  (i,j)
                0129            rCenter = rCenter2D(i,j)
                0130            rUpwind = rUpwind2D(i,j)
                0131 #else
bb6c554092 Jean*0132            rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
7b4413ef94 Jean*0133            mskM   = maskC(i,j,km2,bi,bj)*maskM2
                0134            mskP   = maskC(i,j,kp1,bi,bj)*maskP1
                0135            IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN
                0136             rUpwind= 0. _d 0
381eb13d88 Jean*0137             rC4km  = oneSixth*rCenter*mskM
                0138             rC4kp  = oneSixth*rCenter*mskP
                0139            ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
                0140             wCFL = deltaTcfl*ABS(rTrans(i,j))
                0141      &            *recip_rA(i,j,bi,bj)*recip_drC(k)
a7ec469280 Jean*0142      &            *recip_deepFac2F(k)*recip_rhoFacF(k)
381eb13d88 Jean*0143             rHigh  = (1. _d 0 -wCFL*wCFL)*oneSixth
                0144 c           rUpwind= (2. _d 0*rHigh - wCFL)*ABS(rCenter)
                0145             rUpwind= (2. _d 0*rHigh )*ABS(rCenter)
                0146             rC4km  =  rHigh * (rCenter+ABS(rCenter))*mskM
                0147             rC4kp  =  rHigh * (rCenter-ABS(rCenter))*mskP
7b4413ef94 Jean*0148            ELSE
381eb13d88 Jean*0149             rUpwind=  2. _d 0*oneSixth*ABS(rCenter)
                0150             rC4km  = oneSixth*(rCenter+ABS(rCenter))*mskM
                0151             rC4kp  = oneSixth*(rCenter-ABS(rCenter))*mskP
7b4413ef94 Jean*0152            ENDIF
64442af1fe Jean*0153 #endif /* ALLOW_AUTODIFF and TARGET_NEC_SX */
7b4413ef94 Jean*0154            a5d(i,j,k)   = a5d(i,j,k)
381eb13d88 Jean*0155      &                  + rC4km
9de7a55d87 Jean*0156      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0157      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0158      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
7b4413ef94 Jean*0159            b5d(i,j,k)   = b5d(i,j,k)
381eb13d88 Jean*0160      &                  - ( (rCenter+rUpwind) + rC4km )
9de7a55d87 Jean*0161      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0162      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0163      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
7b4413ef94 Jean*0164            c5d(i,j,k)   = c5d(i,j,k)
381eb13d88 Jean*0165      &                  - ( (rCenter-rUpwind) + rC4kp )
9de7a55d87 Jean*0166      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0167      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0168      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
7b4413ef94 Jean*0169            d5d(i,j,k)   = d5d(i,j,k)
381eb13d88 Jean*0170      &                  + rC4kp
9de7a55d87 Jean*0171      &                   *deltaTarg(k)
ec0db5c1b3 Jean*0172      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0173      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
7b4413ef94 Jean*0174            b5d(i,j,k-1) = b5d(i,j,k-1)
381eb13d88 Jean*0175      &                  - rC4km
9de7a55d87 Jean*0176      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0177      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0178      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
7b4413ef94 Jean*0179            c5d(i,j,k-1) = c5d(i,j,k-1)
381eb13d88 Jean*0180      &                  + ( (rCenter+rUpwind) + rC4km )
9de7a55d87 Jean*0181      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0182      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0183      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
7b4413ef94 Jean*0184            d5d(i,j,k-1) = d5d(i,j,k-1)
381eb13d88 Jean*0185      &                  + ( (rCenter-rUpwind) + rC4kp )
9de7a55d87 Jean*0186      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0187      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0188      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
7b4413ef94 Jean*0189            e5d(i,j,k-1) = e5d(i,j,k-1)
381eb13d88 Jean*0190      &                  - rC4kp
9de7a55d87 Jean*0191      &                   *deltaTarg(k-1)
ec0db5c1b3 Jean*0192      &                   *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0193      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
7b4413ef94 Jean*0194          ENDDO
                0195        ENDDO
                0196
9de7a55d87 Jean*0197 C--   process interior interface only: end
                0198       ENDIF
                0199
7b4413ef94 Jean*0200       RETURN
                0201       END