pkg/generic_advdiff/gad_dst3fl_adv_r.F

7baaf69241 Alis*0001 #include "GAD_OPTIONS.h"
                0002
8a9f54a9ca Jean*0003 CBOP
                0004 C !ROUTINE: GAD_DST3FL_ADV_R
                0005
                0006 C !INTERFACE: ==========================================================
0af3073e4e Jean*0007       SUBROUTINE GAD_DST3FL_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 DST Scheme with flux limiting
                0017
                0018 C !USES: ===============================================================
7baaf69241 Alis*0019       IMPLICIT NONE
                0020
                0021 C     == GLobal variables ==
                0022 #include "SIZE.h"
                0023 #include "GRID.h"
                0024 #include "GAD.h"
                0025
                0026 C     == Routine arguments ==
8a9f54a9ca Jean*0027 C !INPUT PARAMETERS: ===================================================
                0028 C  bi,bj             :: tile indices
                0029 C  k                 :: vertical level
                0030 C  deltaTloc         :: local time-step (s)
                0031 C  rTrans            :: vertical volume transport
0af3073e4e Jean*0032 C  wFld              :: vertical flow
8a9f54a9ca Jean*0033 C  tracer            :: tracer field
                0034 C  myThid            :: thread number
                0035       INTEGER bi,bj,k
7baaf69241 Alis*0036       _RL dTarg
                0037       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0038       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8a9f54a9ca Jean*0039       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
7baaf69241 Alis*0040       INTEGER myThid
                0041
8a9f54a9ca Jean*0042 C !OUTPUT PARAMETERS: ==================================================
                0043 C  wT                :: vertical advective flux
                0044       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0045
7baaf69241 Alis*0046 C     == Local variables ==
8a9f54a9ca Jean*0047 C !LOCAL VARIABLES: ====================================================
                0048 C  i,j               :: loop indices
                0049 C  km1               :: =max( k-1 , 1 )
0af3073e4e Jean*0050 C  wLoc              :: velocity, vertical component
8a9f54a9ca Jean*0051 C  wCFL              :: Courant-Friedrich-Levy number
                0052       INTEGER i,j,kp1,km1,km2
c6bc48bdb5 Jean*0053       _RL Rjm,Rj,Rjp,wCFL,d0,d1
7baaf69241 Alis*0054       _RL psiP,psiM,thetaP,thetaM
0af3073e4e Jean*0055       _RL wLoc
20eb3a9e6b Jean*0056       _RL thetaMax
                0057       PARAMETER( thetaMax = 1.D+20 )
7baaf69241 Alis*0058
                0059       km2=MAX(1,k-2)
                0060       km1=MAX(1,k-1)
                0061       kp1=MIN(Nr,k+1)
                0062
64442af1fe Jean*0063       DO j=1-OLy,sNy+OLy
                0064        DO i=1-OLx,sNx+OLx
                0065 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
bf0222e9e9 Mart*0066 C     These lines make TAF create vectorizable code
                0067         thetaP = 0. _d 0
                0068         thetaM = 0. _d 0
                0069 #endif
8a9f54a9ca Jean*0070         Rjp=(tracer(i,j,k)-tracer(i,j,kp1))
                0071      &         *maskC(i,j,kp1,bi,bj)
                0072         Rj =(tracer(i,j,km1)-tracer(i,j,k))
                0073      &         *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
                0074         Rjm=(tracer(i,j,km2)-tracer(i,j,km1))
                0075      &         *maskC(i,j,km1,bi,bj)
                0076
99c9058df1 Jean*0077         wLoc = wFld(i,j)
c6bc48bdb5 Jean*0078         wCFL = ABS(wLoc*dTarg*recip_drC(k))
                0079         d0=(2. _d 0 -wCFL)*(1. _d 0 -wCFL)*oneSixth
                0080         d1=(1. _d 0 -wCFL*wCFL)*oneSixth
20eb3a9e6b Jean*0081
                0082 C-      the old version: can produce overflow, division by zero,
                0083 C       and is wrong for tracer with low concentration:
                0084 c       thetaP=Rjm/(1.D-20+Rj)
                0085 c       thetaM=Rjp/(1.D-20+Rj)
                0086 C-      the right expression, but not bounded:
51862bbf1d Patr*0087 c       thetaP=0.D0
                0088 c       thetaM=0.D0
20eb3a9e6b Jean*0089 c       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
51862bbf1d Patr*0090 c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
20eb3a9e6b Jean*0091 C-      prevent |thetaP,M| to reach too big value:
                0092         IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
                0093           thetaP=SIGN(thetaMax,Rjm*Rj)
                0094         ELSE
                0095           thetaP=Rjm/Rj
                0096         ENDIF
                0097         IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
                0098           thetaM=SIGN(thetaMax,Rjp*Rj)
                0099         ELSE
                0100           thetaM=Rjp/Rj
                0101         ENDIF
                0102
                0103         psiP=d0+d1*thetaP
                0104         psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),
c6bc48bdb5 Jean*0105      &                       thetaP*(1. _d 0 -wCFL)/(wCFL+1. _d -20) ))
7baaf69241 Alis*0106         psiM=d0+d1*thetaM
20eb3a9e6b Jean*0107         psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),
c6bc48bdb5 Jean*0108      &                       thetaM*(1. _d 0 -wCFL)/(wCFL+1. _d -20) ))
20eb3a9e6b Jean*0109
7baaf69241 Alis*0110         wT(i,j)=
c6bc48bdb5 Jean*0111      &   0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
8a9f54a9ca Jean*0112      &      *( tracer(i,j, k ) + psiM*Rj )
c6bc48bdb5 Jean*0113      &  +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
8a9f54a9ca Jean*0114      &      *( tracer(i,j,km1) - psiP*Rj )
7baaf69241 Alis*0115
                0116        ENDDO
                0117       ENDDO
                0118
                0119       RETURN
                0120       END