pkg/generic_advdiff/gad_dst3_adv_r.F

7baaf69241 Alis*0001 #include "GAD_OPTIONS.h"
                0002
8a9f54a9ca Jean*0003 CBOP
                0004 C !ROUTINE: GAD_DST3_ADV_R
                0005
                0006 C !INTERFACE: ==========================================================
0af3073e4e Jean*0007       SUBROUTINE GAD_DST3_ADV_R(
8a9f54a9ca Jean*0008      I           bi,bj,k,dTarg,
0af3073e4e Jean*0009      I           rTrans, wFld,
7baaf69241 Alis*0010      I           tracer,
                0011      O           wT,
                0012      I           myThid )
8a9f54a9ca Jean*0013
                0014 C !DESCRIPTION:
                0015 C  Calculates the area integrated vertical flux due to advection of a tracer
                0016 C  using 3rd-order Direct Space and Time (DST-3) Advection Scheme
                0017
                0018 C !USES: ===============================================================
7baaf69241 Alis*0019       IMPLICIT NONE
                0020
                0021 C     == GLobal variables ==
                0022 #include "SIZE.h"
cc94647d10 Jean*0023 #ifdef OLD_DST3_FORMULATION
7baaf69241 Alis*0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
cc94647d10 Jean*0026 #endif
                0027 #include "GRID.h"
7baaf69241 Alis*0028 #include "GAD.h"
                0029
                0030 C     == Routine arguments ==
8a9f54a9ca Jean*0031 C !INPUT PARAMETERS: ===================================================
                0032 C  bi,bj             :: tile indices
                0033 C  k                 :: vertical level
                0034 C  deltaTloc         :: local time-step (s)
                0035 C  rTrans            :: vertical volume transport
0af3073e4e Jean*0036 C  wFld              :: vertical flow
8a9f54a9ca Jean*0037 C  tracer            :: tracer field
                0038 C  myThid            :: thread number
                0039       INTEGER bi,bj,k
7baaf69241 Alis*0040       _RL dTarg
                0041       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0042       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8a9f54a9ca Jean*0043       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
7baaf69241 Alis*0044       INTEGER myThid
                0045
8a9f54a9ca Jean*0046 C !OUTPUT PARAMETERS: ==================================================
                0047 C  wT                :: vertical advective flux
                0048       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049
7baaf69241 Alis*0050 C     == Local variables ==
8a9f54a9ca Jean*0051 C !LOCAL VARIABLES: ====================================================
                0052 C  i,j               :: loop indices
                0053 C  km1               :: =max( k-1 , 1 )
0af3073e4e Jean*0054 C  wLoc              :: velocity, vertical component
8a9f54a9ca Jean*0055 C  wCFL              :: Courant-Friedrich-Levy number
                0056       INTEGER i,j,kp1,km1,km2
cf79b8bd99 Jean*0057       _RL wLoc
7baaf69241 Alis*0058       _RL Rjm,Rj,Rjp,cfl,d0,d1
cf79b8bd99 Jean*0059 #ifdef OLD_DST3_FORMULATION
7baaf69241 Alis*0060       _RL psiP,psiM,thetaP,thetaM
b79a37688e Patr*0061       _RL smallNo
7baaf69241 Alis*0062
0d75a51072 Mart*0063 c     IF (inAdMode .AND. useApproxAdvectionInAdMode) THEN
370e5a47fd Jean*0064 c      smallNo = 1.0D-20
                0065 c     ELSE
b79a37688e Patr*0066        smallNo = 1.0D-20
370e5a47fd Jean*0067 c     ENDIF
cf79b8bd99 Jean*0068 #endif
b79a37688e Patr*0069
7baaf69241 Alis*0070       km2=MAX(1,k-2)
                0071       km1=MAX(1,k-1)
                0072       kp1=MIN(Nr,k+1)
                0073
370e5a47fd Jean*0074       DO j=1-OLy,sNy+OLy
                0075        DO i=1-OLx,sNx+OLx
8a9f54a9ca Jean*0076         Rjp=(tracer(i,j,k)-tracer(i,j,kp1))
                0077      &         *maskC(i,j,kp1,bi,bj)
                0078         Rj =(tracer(i,j,km1)-tracer(i,j,k))
                0079      &         *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
                0080         Rjm=(tracer(i,j,km2)-tracer(i,j,km1))
                0081      &         *maskC(i,j,km1,bi,bj)
7baaf69241 Alis*0082
99c9058df1 Jean*0083         wLoc = wFld(i,j)
                0084 c       wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
cf79b8bd99 Jean*0085         cfl=ABS(wLoc*dTarg*recip_drC(k))
7baaf69241 Alis*0086         d0=(2.-cfl)*(1.-cfl)*oneSixth
                0087         d1=(1.-cfl*cfl)*oneSixth
cf79b8bd99 Jean*0088 #ifdef OLD_DST3_FORMULATION
b79a37688e Patr*0089         IF ( ABS(Rj).LT.smallNo .OR.
                0090      &       ABS(Rjm).LT.smallNo ) THEN
                0091          thetaP=0.
                0092          psiP=0.
                0093         ELSE
                0094          thetaP=(Rjm+smallNo)/(smallNo+Rj)
                0095          psiP=d0+d1*thetaP
                0096         ENDIF
                0097         IF ( ABS(Rj).LT.smallNo .OR.
                0098      &       ABS(Rjp).LT.smallNo ) THEN
                0099          thetaM=0.
                0100          psiM=0.
                0101         ELSE
                0102          thetaM=(Rjp+smallNo)/(smallNo+Rj)
                0103          psiM=d0+d1*thetaM
                0104         ENDIF
                0105          wT(i,j)=
cf79b8bd99 Jean*0106      &    0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
8a9f54a9ca Jean*0107      &       *( tracer(i,j, k ) + psiM*Rj )
cf79b8bd99 Jean*0108      &   +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
8a9f54a9ca Jean*0109      &       *( tracer(i,j,km1) - psiP*Rj )
cf79b8bd99 Jean*0110 #else /* OLD_DST3_FORMULATION */
                0111         wT(i,j)=
                0112      &    0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
                0113      &       *( tracer(i,j, k ) + (d0*Rj+d1*Rjp) )
                0114      &   +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
                0115      &       *( tracer(i,j,km1) - (d0*Rj+d1*Rjm) )
                0116 #endif /* OLD_DST3_FORMULATION */
983c7d32b1 Jean*0117
7baaf69241 Alis*0118        ENDDO
                0119       ENDDO
                0120
                0121       RETURN
                0122       END