pkg/generic_advdiff/gad_fluxlimit_adv_r.F

31566b6684 Alis*0001 #include "GAD_OPTIONS.h"
                0002
527a84022c Alis*0003 CBOP
                0004 C !ROUTINE: GAD_FLUXLIMIT_ADV_R
                0005
                0006 C !INTERFACE: ==========================================================
0af3073e4e Jean*0007       SUBROUTINE GAD_FLUXLIMIT_ADV_R(
8f93275e63 Jean*0008      I           bi, bj, k, dTarg,
0af3073e4e Jean*0009      I           rTrans, wFld,
31566b6684 Alis*0010      I           tracer,
                0011      O           wT,
                0012      I           myThid )
                0013
527a84022c Alis*0014 C !DESCRIPTION:
                0015 C Calculates the area integrated vertical flux due to advection of a tracer
                0016 C using second-order interpolation with a flux limiter:
                0017 C \begin{equation*}
                0018 C F^x_{adv} = W \overline{ \theta }^k
                0019 C - \frac{1}{2} \left(
                0020 C     [ 1 - \psi(C_r) ] |W|
                0021 C    + W \frac{w \Delta t}{\Delta r_c} \psi(C_r)
                0022 C              \right) \delta_k \theta
                0023 C \end{equation*}
                0024 C where the $\psi(C_r)$ is the limiter function and $C_r$ is
                0025 C the slope ratio.
                0026
                0027 C !USES: ===============================================================
                0028       IMPLICIT NONE
31566b6684 Alis*0029 #include "SIZE.h"
                0030 #include "GRID.h"
24016b3859 Alis*0031 #include "EEPARAMS.h"
                0032 #include "PARAMS.h"
31566b6684 Alis*0033
527a84022c Alis*0034 C !INPUT PARAMETERS: ===================================================
8f93275e63 Jean*0035 C  bi, bj             :: tile indices
1b5fb69d21 Ed H*0036 C  k                  :: vertical level
                0037 C  rTrans             :: vertical volume transport
0af3073e4e Jean*0038 C  wFld               :: vertical flow
1b5fb69d21 Ed H*0039 C  tracer             :: tracer field
                0040 C  myThid             :: thread number
8f93275e63 Jean*0041       INTEGER bi, bj, k
24016b3859 Alis*0042       _RL dTarg
31566b6684 Alis*0043       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0044       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8a9f54a9ca Jean*0045       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31566b6684 Alis*0046       INTEGER myThid
                0047
527a84022c Alis*0048 C !OUTPUT PARAMETERS: ==================================================
1b5fb69d21 Ed H*0049 C  wT                 :: vertical advective flux
527a84022c Alis*0050       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0051
                0052 C !LOCAL VARIABLES: ====================================================
8f93275e63 Jean*0053 C  i, j               :: loop indices
1b5fb69d21 Ed H*0054 C  kp1                :: =min( k+1 , Nr )
                0055 C  km1                :: =max( k-1 , 1 )
                0056 C  km2                :: =max( k-2 , 1 )
8f93275e63 Jean*0057 C  bi, bj             :: tile indices or (1,1) depending on use
1b5fb69d21 Ed H*0058 C  Cr                 :: slope ratio
8f93275e63 Jean*0059 C  Rjm, Rj, Rjp       :: differences at k-1,k,k+1
0af3073e4e Jean*0060 C  wLoc               :: velocity, vertical component
8f93275e63 Jean*0061       INTEGER i, j, kp1, km1, km2
                0062       _RL Cr, Rjm, Rj, Rjp
c6bc48bdb5 Jean*0063       _RL wLoc, wCFL
8f93275e63 Jean*0064       _RL CrMax
                0065       PARAMETER( CrMax = 1.D+6 )
                0066
527a84022c Alis*0067 C Statement function provides Limiter(Cr)
31566b6684 Alis*0068 #include "GAD_FLUX_LIMITER.h"
527a84022c Alis*0069 CEOP
31566b6684 Alis*0070
8f93275e63 Jean*0071       km2 = MAX(1,k-2)
                0072       km1 = MAX(1,k-1)
                0073       kp1 = MIN(Nr,k+1)
31566b6684 Alis*0074
                0075       IF ( k.GT.Nr) THEN
8f93275e63 Jean*0076        DO j=1-OLy,sNy+OLy
                0077         DO i=1-OLx,sNx+OLx
31566b6684 Alis*0078          wT(i,j) = 0.
                0079         ENDDO
                0080        ENDDO
                0081       ELSE
8f93275e63 Jean*0082        DO j=1-OLy,sNy+OLy
                0083         DO i=1-OLx,sNx+OLx
2b22d4bdc3 Jean*0084
99c9058df1 Jean*0085          wLoc = wFld(i,j)
c6bc48bdb5 Jean*0086          wCFL = ABS( wLoc*dTarg*recip_drC(k) )
8f93275e63 Jean*0087          Rjp = (tracer(i,j,kp1)-tracer(i,j,k))
                0088      &          *maskC(i,j,kp1,bi,bj)
                0089          Rj  = (tracer(i,j,k)  -tracer(i,j,km1))
                0090          Rjm = (tracer(i,j,km1)-tracer(i,j,km2))
                0091      &          *maskC(i,j,km2,bi,bj)
2b22d4bdc3 Jean*0092
8f93275e63 Jean*0093          IF ( rTrans(i,j).LT.zeroRL ) THEN
                0094            Cr = Rjm
                0095          ELSE
                0096            Cr = Rjp
                0097          ENDIF
                0098          IF ( ABS(Rj)*CrMax .LE. ABS(Cr) ) THEN
                0099            Cr = SIGN( CrMax, Cr )*SIGN( oneRL, Rj )
31566b6684 Alis*0100          ELSE
8f93275e63 Jean*0101            Cr = Cr/Rj
31566b6684 Alis*0102          ENDIF
8f93275e63 Jean*0103
                0104 C        calculate Limiter Function:
                0105          Cr = Limiter(Cr)
                0106
                0107          wT(i,j) = maskC(i,j,km1,bi,bj)*(
                0108      &      rTrans(i,j)*
                0109      &        (tracer(i,j,k)+tracer(i,j,km1))*0.5 _d 0
                0110      &    + ABS(rTrans(i,j))*( (oneRL-Cr) + wCFL*Cr )
                0111      &                      *Rj*0.5 _d 0 )
31566b6684 Alis*0112         ENDDO
                0113        ENDDO
                0114       ENDIF
                0115
                0116       RETURN
                0117       END