pkg/generic_advdiff/gad_dst2u1_impl_r.F

381eb13d88 Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: GAD_DST2U1_IMPL_R
                0005 C     !INTERFACE:
ec0db5c1b3 Jean*0006       SUBROUTINE GAD_DST2U1_IMPL_R(
                0007      I           bi,bj,k, iMin,iMax,jMin,jMax,
                0008      I           advectionScheme, deltaTarg, rTrans, recip_hFac,
381eb13d88 Jean*0009      O           a3d, b3d, c3d,
                0010      I           myThid )
                0011
                0012 C     !DESCRIPTION:
                0013 C     Compute matrix element to solve vertical advection implicitly
                0014 C     using DST 2nd.Order (=Lax-Wendroff) or 1rst Order Upwind 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.
                0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021
                0022 C     == Global variables ===
                0023 #include "SIZE.h"
                0024 #include "GRID.h"
a7ec469280 Jean*0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
381eb13d88 Jean*0027 #include "GAD.h"
                0028
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine Arguments ==
                0031 C     bi,bj        :: tile indices
                0032 C     k            :: vertical level
                0033 C  advectionScheme :: advection scheme to use: either 2nd Order DST
                0034 C                                                or 1rst Order Upwind
                0035 C     iMin,iMax    :: computation domain
                0036 C     jMin,jMax    :: computation domain
                0037 C     deltaTarg    :: time step
                0038 C     rTrans       :: vertical volume transport
ec0db5c1b3 Jean*0039 C     recip_hFac   :: inverse of cell open-depth factor
                0040 C     a3d          :: lower diagonal of the tridiagonal matrix
                0041 C     b3d          :: main  diagonal of the tridiagonal matrix
                0042 C     c3d          :: upper diagonal of the tridiagonal matrix
381eb13d88 Jean*0043 C     myThid       :: thread number
                0044       INTEGER bi,bj,k
                0045       INTEGER iMin,iMax,jMin,jMax
                0046       INTEGER advectionScheme
ec0db5c1b3 Jean*0047       _RL     deltaTarg(Nr)
                0048       _RL     rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0050       _RL     a3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0051       _RL     b3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0052       _RL     c3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
381eb13d88 Jean*0053       INTEGER myThid
                0054
                0055 C     == Local Variables ==
                0056 C     i,j          :: loop indices
                0057 C     w_CFL        :: Courant-Friedrich-Levy number
                0058 C     rLimit       :: centered (vs upwind) fraction
                0059 C     rCenter      :: centered contribution
                0060 C     rUpwind      :: upwind   contribution
                0061       INTEGER i,j
a7ec469280 Jean*0062 c     _RL w_CFL
                0063       _RL rLimit
381eb13d88 Jean*0064       _RL rCenter, rUpwind
                0065       _RL deltaTcfl
                0066
                0067 CEOP
                0068
                0069       rLimit = 0. _d 0
                0070       IF ( advectionScheme.EQ.ENUM_DST2 ) rLimit = 1. _d 0
                0071
                0072 C--   process interior interface only:
                0073       IF ( k.GT.1 .AND. k.LE.Nr ) THEN
                0074
                0075 C--    Add centered & upwind contributions
                0076        deltaTcfl = deltaTarg(k)
                0077        DO j=jMin,jMax
                0078          DO i=iMin,iMax
a7ec469280 Jean*0079 c          w_CFL = deltaTcfl*ABS(rTrans(i,j))
                0080 c    &            *recip_rA(i,j,bi,bj)*recip_drC(k)
                0081 c    &            *recip_deepFac2F(k)*recip_rhoFacF(k)
381eb13d88 Jean*0082            rCenter = 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
                0083            rUpwind = ABS(rCenter)
                0084      &             * ( 1. _d 0 - rLimit )
                0085 c    &             * ( 1. _d 0 - rLimit*( 1. _d 0 + w_CFL ) )
                0086            a3d(i,j,k)   = a3d(i,j,k)
                0087      &                  - (rCenter+rUpwind)*deltaTarg(k)
ec0db5c1b3 Jean*0088      &                   *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0089      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
381eb13d88 Jean*0090            b3d(i,j,k)   = b3d(i,j,k)
                0091      &                  - (rCenter-rUpwind)*deltaTarg(k)
ec0db5c1b3 Jean*0092      &                    *recip_hFac(i,j,k)*recip_drF(k)
a7ec469280 Jean*0093      &                   *recip_deepFac2C(k)*recip_rhoFacC(k)
381eb13d88 Jean*0094            b3d(i,j,k-1) = b3d(i,j,k-1)
                0095      &                  + (rCenter+rUpwind)*deltaTarg(k-1)
ec0db5c1b3 Jean*0096      &                    *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0097      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
381eb13d88 Jean*0098            c3d(i,j,k-1) = c3d(i,j,k-1)
                0099      &                  + (rCenter-rUpwind)*deltaTarg(k-1)
ec0db5c1b3 Jean*0100      &                    *recip_hFac(i,j,k-1)*recip_drF(k-1)
a7ec469280 Jean*0101      &                   *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
381eb13d88 Jean*0102          ENDDO
                0103        ENDDO
                0104
                0105 C--   process interior interface only: end
                0106       ENDIF
                0107
                0108       RETURN
                0109       END