pkg/generic_advdiff/gad_dst2u1_adv_r.F

76b8386f6e Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: GAD_DST2U1_ADV_R
                0005
                0006 C !INTERFACE: ==========================================================
0ef5a3ea3b Jean*0007       SUBROUTINE GAD_DST2U1_ADV_R(
76b8386f6e Jean*0008      I           bi,bj,k, advectionScheme,
0af3073e4e Jean*0009      I           deltaTloc, rTrans, wFld,
76b8386f6e Jean*0010      I           tracer,
                0011      O           wT,
                0012      I           myThid )
                0013
                0014 C !DESCRIPTION:
                0015 C  Calculates the area integrated vertical flux due to advection
0ef5a3ea3b Jean*0016 C  of a tracer using second-order Direct Space and Time (DST-2)
76b8386f6e Jean*0017 C  interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme.
                0018
                0019 C !USES: ===============================================================
                0020       IMPLICIT NONE
                0021 #include "SIZE.h"
                0022 #include "GRID.h"
                0023 #include "GAD.h"
                0024
                0025 C !INPUT PARAMETERS: ===================================================
                0026 C  bi,bj             :: tile indices
                0027 C  k                 :: vertical level
                0028 C  advectionScheme   :: advection scheme to use: either 2nd Order DST
                0029 C                                                or 1rst Order Upwind
                0030 C  deltaTloc         :: local time-step (s)
                0031 C  rTrans            :: vertical volume transport
0af3073e4e Jean*0032 C  wFld              :: vertical flow
76b8386f6e Jean*0033 C  tracer            :: tracer field
                0034 C  myThid            :: thread number
                0035       INTEGER bi,bj,k
                0036       INTEGER advectionScheme
                0037       _RL deltaTloc
                0038       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0039       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76b8386f6e Jean*0040       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0041       INTEGER myThid
                0042
                0043 C !OUTPUT PARAMETERS: ==================================================
                0044 C  wT                :: vertical advective flux
                0045       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046
                0047 C !LOCAL VARIABLES: ====================================================
                0048 C  i,j               :: loop indices
                0049 C  km1               :: =max( k-1 , 1 )
                0050 C  rLimit            :: centered (vs upwind) fraction
0af3073e4e Jean*0051 C  wLoc              :: velocity, vertical component
76b8386f6e Jean*0052 C  wCFL              :: Courant-Friedrich-Levy number
                0053       INTEGER i,j,km1
0af3073e4e Jean*0054       _RL wLoc, wCFL, rLimit, wAbs
76b8386f6e Jean*0055 CEOP
                0056
                0057       rLimit = 0. _d 0
                0058       IF ( advectionScheme.EQ.ENUM_DST2 ) rLimit = 1. _d 0
                0059
                0060       km1=MAX(1,k-1)
                0061
                0062       IF ( k.LE.1 .OR. k.GT.Nr) THEN
                0063        DO j=1-Oly,sNy+Oly
                0064         DO i=1-Olx,sNx+Olx
                0065          wT(i,j) = 0.
                0066         ENDDO
                0067        ENDDO
                0068       ELSE
                0069        DO j=1-Oly,sNy+Oly
                0070         DO i=1-Olx,sNx+Olx
                0071
99c9058df1 Jean*0072          wLoc = wFld(i,j)
                0073 c        wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
0af3073e4e Jean*0074          wCFL = ABS(wLoc*deltaTloc*recip_drC(k))
76b8386f6e Jean*0075
0ef5a3ea3b Jean*0076 c        wT(i,j) = maskC(i,j,km1,bi,bj)*(
                0077 c    &     rTrans(i,j)*(tracer(i,j,km1)+tracer(i,j,k))*0.5 _d 0
                0078 c    &   + ( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )*ABS(rTrans(i,j))
                0079 c    &                *(tracer(i,j,km1)-tracer(i,j,k))*0.5 _d 0*rkSign
                0080 c    &                                  )
                0081 C-- above formulation produces large truncation error when:
                0082 C    1rst.O upWind and   w*rkSign > 0 & |tracer(k-1)| << |tracer(k)|
                0083 C                   or   w*rkSign < 0 & |tracer(k-1)| >> |tracer(k)|
                0084 C-- change to a more robust expression:
                0085          wAbs = ABS(rTrans(i,j))*rkSign
                0086      &       *( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )
                0087          wT(i,j) = maskC(i,j,km1,bi,bj)*(
                0088      &             ( rTrans(i,j)+wAbs )* 0.5 _d 0 * tracer(i,j,km1)
                0089      &           + ( rTrans(i,j)-wAbs )* 0.5 _d 0 * tracer(i,j,k)
76b8386f6e Jean*0090      &                                  )
                0091         ENDDO
                0092        ENDDO
                0093       ENDIF
                0094
                0095       RETURN
                0096       END