pkg/generic_advdiff/gad_dst2u1_adv_y.F

76b8386f6e Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: GAD_DST2U1_ADV_Y
                0005
                0006 C !INTERFACE: ==========================================================
0ef5a3ea3b Jean*0007       SUBROUTINE GAD_DST2U1_ADV_Y(
692dd30681 Jean*0008      I           bi,bj,k, advectionScheme, calcCFL,
                0009      I           deltaTloc, vTrans, vFld,
76b8386f6e Jean*0010      I           tracer,
                0011      O           vT,
                0012      I           myThid )
                0013
                0014 C !DESCRIPTION:
                0015 C  Calculates the area integrated meridional 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
692dd30681 Jean*0030 C  calcCFL           :: =T: calculate CFL number ; =F: take vFld as CFL
                0031 C  deltaTloc         :: local time-step (s)
76b8386f6e Jean*0032 C  vTrans            :: meridional volume transport
692dd30681 Jean*0033 C  vFld              :: meridional flow / CFL number
76b8386f6e Jean*0034 C  tracer            :: tracer field
                0035 C  myThid            :: thread number
                0036       INTEGER bi,bj, k, advectionScheme
692dd30681 Jean*0037       LOGICAL calcCFL
76b8386f6e Jean*0038       _RL deltaTloc
                0039       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0040       _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76b8386f6e Jean*0041       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       INTEGER myThid
                0043
                0044 C !OUTPUT PARAMETERS: ==================================================
                0045 C  vT                :: meridional advective flux
                0046       _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047
                0048 C !LOCAL VARIABLES: ====================================================
                0049 C  i,j               :: loop indices
                0050 C  yLimit            :: centered (vs upwind) fraction
                0051 C  vCFL              :: Courant-Friedrich-Levy number
                0052       INTEGER i,j
692dd30681 Jean*0053       _RL vCFL, yLimit, vAbs
76b8386f6e Jean*0054 CEOP
                0055
                0056       yLimit = 0. _d 0
                0057       IF ( advectionScheme.EQ.ENUM_DST2 ) yLimit = 1. _d 0
                0058
                0059       DO i=1-Olx,sNx+Olx
                0060        vT(i,1-Oly)=0.
                0061       ENDDO
                0062       DO j=1-Oly+1,sNy+Oly
                0063        DO i=1-Olx,sNx+Olx
                0064
692dd30681 Jean*0065         vCFL = vFld(i,j)
                0066         IF ( calcCFL ) vCFL = ABS( vFld(i,j)*deltaTloc
62743db5c5 Jean*0067      &                  *recip_dyC(i,j,bi,bj)*recip_deepFacC(k) )
76b8386f6e Jean*0068
0ef5a3ea3b Jean*0069 c       vT(i,j) =
                0070 c    &     vTrans(i,j)*(tracer(i,j-1)+tracer(i,j))*0.5 _d 0
                0071 c    &   + ( 1. _d 0 - yLimit*(1. _d 0 - vCFL) )*ABS(vTrans(i,j))
                0072 c    &                *(tracer(i,j-1)-tracer(i,j))*0.5 _d 0
                0073 C-- above formulation produces large truncation error when:
                0074 C    1rst.O upWind and   v > 0 & |tracer(i,j-1)| << |tracer(i,j)|
                0075 C                   or   v < 0 & |tracer(i,j-1)| >> |tracer(i,j)|
                0076 C-- change to a more robust expression:
                0077         vAbs = ABS(vTrans(i,j))
                0078      &       *( 1. _d 0 - yLimit*(1. _d 0 - vCFL) )
                0079         vT(i,j) = ( vTrans(i,j)+vAbs )* 0.5 _d 0 * tracer(i,j-1)
                0080      &          + ( vTrans(i,j)-vAbs )* 0.5 _d 0 * tracer(i,j)
                0081
76b8386f6e Jean*0082        ENDDO
                0083       ENDDO
                0084
                0085       RETURN
                0086       END