C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst2u1_adv_r.F,v 1.4 2006/06/19 14:40:43 jmc Exp $ C $Name: $ #include "GAD_OPTIONS.h" CBOP C !ROUTINE: GAD_DST2U1_ADV_R C !INTERFACE: ========================================================== SUBROUTINE GAD_DST2U1_ADV_R( I bi,bj,k, advectionScheme, I deltaTloc, rTrans, wFld, I tracer, O wT, I myThid ) C !DESCRIPTION: C Calculates the area integrated vertical flux due to advection C of a tracer using second-order Direct Space and Time (DST-2) C interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme. C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "GRID.h" #include "GAD.h" C !INPUT PARAMETERS: =================================================== C bi,bj :: tile indices C k :: vertical level C advectionScheme :: advection scheme to use: either 2nd Order DST C or 1rst Order Upwind C deltaTloc :: local time-step (s) C rTrans :: vertical volume transport C wFld :: vertical flow C tracer :: tracer field C myThid :: thread number INTEGER bi,bj,k INTEGER advectionScheme _RL deltaTloc _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) INTEGER myThid C !OUTPUT PARAMETERS: ================================================== C wT :: vertical advective flux _RL wT (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== C i,j :: loop indices C km1 :: =max( k-1 , 1 ) C rLimit :: centered (vs upwind) fraction C wLoc :: velocity, vertical component C wCFL :: Courant-Friedrich-Levy number INTEGER i,j,km1 _RL wLoc, wCFL, rLimit, wAbs CEOP rLimit = 0. _d 0 IF ( advectionScheme.EQ.ENUM_DST2 ) rLimit = 1. _d 0 km1=MAX(1,k-1) IF ( k.LE.1 .OR. k.GT.Nr) THEN DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx wT(i,j) = 0. ENDDO ENDDO ELSE DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx wLoc = wFld(i,j) c wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj) wCFL = ABS(wLoc*deltaTloc*recip_drC(k)) c wT(i,j) = maskC(i,j,km1,bi,bj)*( c & rTrans(i,j)*(tracer(i,j,km1)+tracer(i,j,k))*0.5 _d 0 c & + ( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )*ABS(rTrans(i,j)) c & *(tracer(i,j,km1)-tracer(i,j,k))*0.5 _d 0*rkSign c & ) C-- above formulation produces large truncation error when: C 1rst.O upWind and w*rkSign > 0 & |tracer(k-1)| << |tracer(k)| C or w*rkSign < 0 & |tracer(k-1)| >> |tracer(k)| C-- change to a more robust expression: wAbs = ABS(rTrans(i,j))*rkSign & *( 1. _d 0 - rLimit*(1. _d 0 - wCFL) ) wT(i,j) = maskC(i,j,km1,bi,bj)*( & ( rTrans(i,j)+wAbs )* 0.5 _d 0 * tracer(i,j,km1) & + ( rTrans(i,j)-wAbs )* 0.5 _d 0 * tracer(i,j,k) & ) ENDDO ENDDO ENDIF RETURN END