C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.77 2016/03/13 01:44:02 jmc Exp $ C $Name: $ #include "GAD_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: GAD_ADVECTION C !INTERFACE: ========================================================== SUBROUTINE GAD_ADVECTION( I implicitAdvection, advectionScheme, vertAdvecScheme, I trIdentity, deltaTLev, I uFld, vFld, wFld, tracer, O gTracer, I bi,bj, myTime,myIter,myThid) C !DESCRIPTION: C Calculates the tendency of a tracer due to advection. C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection} C and can only be used for the non-linear advection schemes such as the C direct-space-time method and flux-limiters. C C The algorithm is as follows: C \begin{itemize} C \item{$\theta^{(n+1/3)} = \theta^{(n)} C - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$} C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)} C - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$} C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)} C - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$} C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$} C \end{itemize} C C The tendency (output) is over-written by this routine. C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "GAD.h" #ifdef ALLOW_AUTODIFF # include "tamc.h" # include "tamc_keys.h" # ifdef ALLOW_PTRACERS # include "PTRACERS_SIZE.h" # endif #endif #ifdef ALLOW_EXCH2 #include "W2_EXCH2_SIZE.h" #include "W2_EXCH2_TOPOLOGY.h" #endif /* ALLOW_EXCH2 */ C !INPUT PARAMETERS: =================================================== C implicitAdvection :: implicit vertical advection (later on) C advectionScheme :: advection scheme to use (Horizontal plane) C vertAdvecScheme :: advection scheme to use (vertical direction) C trIdentity :: tracer identifier C uFld :: Advection velocity field, zonal component C vFld :: Advection velocity field, meridional component C wFld :: Advection velocity field, vertical component C tracer :: tracer field C bi,bj :: tile indices C myTime :: current time C myIter :: iteration number C myThid :: thread number LOGICAL implicitAdvection INTEGER advectionScheme, vertAdvecScheme INTEGER trIdentity _RL deltaTLev(Nr) _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) INTEGER bi,bj _RL myTime INTEGER myIter INTEGER myThid C !OUTPUT PARAMETERS: ================================================== C gTracer :: tendency array _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) C !FUNCTIONS: ========================================================== #ifdef ALLOW_DIAGNOSTICS CHARACTER*4 GAD_DIAG_SUFX EXTERNAL GAD_DIAG_SUFX LOGICAL DIAGNOSTICS_IS_ON EXTERNAL DIAGNOSTICS_IS_ON #endif C !LOCAL VARIABLES: ==================================================== C maskUp :: 2-D array for mask at W points C maskLocW :: 2-D array for mask at West points C maskLocS :: 2-D array for mask at South points C [iMin,iMax]Upd :: loop range to update tracer field C [jMin,jMax]Upd :: loop range to update tracer field C i,j,k :: loop indices C kUp :: index into 2 1/2D array, toggles between 1 and 2 C kDown :: index into 2 1/2D array, toggles between 2 and 1 C xA,yA :: areas of X and Y face of tracer cells C uTrans,vTrans :: 2-D arrays of volume transports at U,V points C rTrans :: 2-D arrays of volume transports at W points C rTransKp :: vertical volume transport at interface k+1 C rTran3d :: 3-D array of volume transport at W points C afr :: 3-D array of vertical advective flux C af :: 2-D array for horizontal advective flux C afx :: 2-D array for horizontal advective flux, x direction C afy :: 2-D array for horizontal advective flux, y direction C fVerT :: 2 1/2D arrays for vertical advective flux C localTij :: 2-D array, temporary local copy of tracer field C localT3d :: 3-D array, temporary local copy of tracer field C kp1Msk :: flag (0,1) for over-riding mask for W levels C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir C interiorOnly :: only update the interior of myTile, but not the edges C overlapOnly :: only update the edges of myTile, but not the interior C npass :: number of passes in multi-dimensional method C ipass :: number of the current pass being made C myTile :: variables used to determine which cube face C nCFace :: owns a tile for cube grid runs using C :: multi-dim advection. C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube c _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd INTEGER i,j,k,kUp,kDown _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rTran3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL afr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL localT3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #ifdef GAD_MULTIDIM_COMPRESSIBLE _RL tmpTrac _RL localVol(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL locVol3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) #endif _RL kp1Msk LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns LOGICAL interiorOnly, overlapOnly INTEGER npass, ipass INTEGER nCFace LOGICAL N_edge, S_edge, E_edge, W_edge #ifdef ALLOW_AUTODIFF_TAMC C msgBuf :: Informational/error message buffer CHARACTER*(MAX_LEN_MBUF) msgBuf #endif #ifdef ALLOW_EXCH2 INTEGER myTile #endif #ifdef ALLOW_DIAGNOSTICS CHARACTER*8 diagName CHARACTER*4 diagSufx LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR #endif /* ALLOW_DIAGNOSTICS */ CEOP #ifdef ALLOW_AUTODIFF_TAMC act0 = trIdentity max0 = maxpass act1 = bi - myBxLo(myThid) max1 = myBxHi(myThid) - myBxLo(myThid) + 1 act2 = bj - myByLo(myThid) max2 = myByHi(myThid) - myByLo(myThid) + 1 act3 = myThid - 1 max3 = nTx*nTy act4 = ikey_dynamics - 1 igadkey = act0 & + act1*max0 & + act2*max0*max1 & + act3*max0*max1*max2 & + act4*max0*max1*max2*max3 IF (trIdentity.GT.maxpass) THEN WRITE(msgBuf,'(A,2I3)') & 'GAD_ADVECTION: maxpass < trIdentity ', & maxpass, trIdentity CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R GAD_ADVECTION' ENDIF #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_DIAGNOSTICS C-- Set diagnostics flags and suffix for the current tracer doDiagAdvX = .FALSE. doDiagAdvY = .FALSE. doDiagAdvR = .FALSE. IF ( useDiagnostics ) THEN diagSufx = GAD_DIAG_SUFX( trIdentity, myThid ) diagName = 'ADVx'//diagSufx doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid ) diagName = 'ADVy'//diagSufx doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid ) diagName = 'ADVr'//diagSufx doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid ) ENDIF #endif /* ALLOW_DIAGNOSTICS */ C-- Set up work arrays with valid (i.e. not NaN) values C These inital values do not alter the numerical results. They C just ensure that all memory references are to valid floating C point numbers. This prevents spurious hardware signals due to C uninitialised but inert locations. DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C- xA,yA,vFld,uTrans,vTrans are set over the full domain C => no need for extra initialisation c xA(i,j) = 0. _d 0 c yA(i,j) = 0. _d 0 c uTrans(i,j) = 0. _d 0 c vTrans(i,j) = 0. _d 0 C- rTransKp is set over the full domain: no need for extra initialisation c rTransKp(i,j)= 0. _d 0 C- rTrans and fVerT need to be initialised to zero: rTrans(i,j) = 0. _d 0 fVerT(i,j,1) = 0. _d 0 fVerT(i,j,2) = 0. _d 0 #ifdef ALLOW_AUTODIFF # ifdef GAD_MULTIDIM_COMPRESSIBLE localVol(i,j) = 0. _d 0 # endif localTij(i,j) = 0. _d 0 #endif /* ALLOW_AUTODIFF */ ENDDO ENDDO C-- Set tile-specific parameters for horizontal fluxes IF (useCubedSphereExchange) THEN npass = 3 #ifdef ALLOW_AUTODIFF_TAMC IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3' #endif #ifdef ALLOW_EXCH2 myTile = W2_myTileList(bi,bj) nCFace = exch2_myFace(myTile) N_edge = exch2_isNedge(myTile).EQ.1 S_edge = exch2_isSedge(myTile).EQ.1 E_edge = exch2_isEedge(myTile).EQ.1 W_edge = exch2_isWedge(myTile).EQ.1 #else nCFace = bi N_edge = .TRUE. S_edge = .TRUE. E_edge = .TRUE. W_edge = .TRUE. #endif ELSE npass = 2 nCFace = 0 N_edge = .FALSE. S_edge = .FALSE. E_edge = .FALSE. W_edge = .FALSE. ENDIF C-- Start of k loop for horizontal fluxes DO k=1,Nr #ifdef ALLOW_AUTODIFF_TAMC kkey = (igadkey-1)*Nr + k CADJ STORE tracer(:,:,k,bi,bj) = CADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C-- Get temporary terms used by tendency routines DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k) & *drF(k)*_hFacW(i,j,k,bi,bj) yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k) & *drF(k)*_hFacS(i,j,k,bi,bj) ENDDO ENDDO C-- Calculate "volume transports" through tracer cell faces. C anelastic: scaled by rhoFacC (~ mass transport) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uTrans(i,j) = uFld(i,j,k)*xA(i,j)*rhoFacC(k) vTrans(i,j) = vFld(i,j,k)*yA(i,j)*rhoFacC(k) ENDDO ENDDO C-- Make local copy of tracer array and mask West & South DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx localTij(i,j) = tracer(i,j,k,bi,bj) #ifdef GAD_MULTIDIM_COMPRESSIBLE localVol(i,j) = rA(i,j,bi,bj)*deepFac2C(k) & *rhoFacC(k)*drF(k)*hFacC(i,j,k,bi,bj) & + ( oneRS - maskC(i,j,k,bi,bj) ) #endif #ifdef ALLOW_OBCS maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj) maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj) #else /* ALLOW_OBCS */ maskLocW(i,j) = _maskW(i,j,k,bi,bj) maskLocS(i,j) = _maskS(i,j,k,bi,bj) #endif /* ALLOW_OBCS */ ENDDO ENDDO IF (useCubedSphereExchange) THEN withSigns = .FALSE. CALL FILL_CS_CORNER_UV_RS( & withSigns, maskLocW,maskLocS, bi,bj, myThid ) ENDIF C-- Multiple passes for different directions on different tiles C-- For cube need one pass for each of red, green and blue axes. DO ipass=1,npass #ifdef ALLOW_AUTODIFF_TAMC passkey = ipass & + (k-1) *maxpass & + (igadkey-1)*maxpass*Nr IF (npass .GT. maxpass) THEN STOP 'GAD_ADVECTION: npass > maxcube. check tamc.h' ENDIF #endif /* ALLOW_AUTODIFF_TAMC */ interiorOnly = .FALSE. overlapOnly = .FALSE. IF (useCubedSphereExchange) THEN C- CubedSphere : pass 3 times, with partial update of local tracer field IF (ipass.EQ.1) THEN overlapOnly = MOD(nCFace,3).EQ.0 interiorOnly = MOD(nCFace,3).NE.0 calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2 calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5 ELSEIF (ipass.EQ.2) THEN overlapOnly = MOD(nCFace,3).EQ.2 interiorOnly = MOD(nCFace,3).EQ.1 calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4 calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1 ELSE interiorOnly = .TRUE. calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6 calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3 ENDIF ELSE C- not CubedSphere calc_fluxes_X = MOD(ipass,2).EQ.1 calc_fluxes_Y = .NOT.calc_fluxes_X ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- X direction #ifdef ALLOW_AUTODIFF C- Always reset advective flux in X DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx af(i,j) = 0. ENDDO ENDDO # ifndef DISABLE_MULTIDIM_ADVECTION CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte CADJ STORE af(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte # endif #endif /* ALLOW_AUTODIFF */ IF (calc_fluxes_X) THEN C- Do not compute fluxes if C a) needed in overlap only C and b) the overlap of myTile are not cube-face Edges IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN C- Internal exchange for calculations in X IF ( overlapOnly ) THEN CALL FILL_CS_CORNER_TR_RL( 1, .FALSE., & localTij, bi,bj, myThid ) ENDIF C- Advective flux in X #ifndef ALLOW_AUTODIFF DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx af(i,j) = 0. ENDDO ENDDO #else /* ALLOW_AUTODIFF */ # ifndef DISABLE_MULTIDIM_ADVECTION CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte # endif #endif /* ALLOW_AUTODIFF */ IF ( advectionScheme.EQ.ENUM_UPWIND_1RST & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE., I deltaTLev(k),uTrans,uFld(1-OLx,1-OLy,k), localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k), I uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k), I uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k), I uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij, O af, myThid ) #ifndef ALLOW_AUTODIFF ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k), I uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_PPM_NULL_LIMIT .OR. & advectionScheme.EQ.ENUM_PPM_MONO_LIMIT .OR. & advectionScheme.EQ.ENUM_PPM_WENO_LIMIT) THEN CALL GAD_PPM_ADV_X( advectionScheme, bi, bj, k , .TRUE., I deltaTLev(k), uFld(1-OLx,1-OLy,k), uTrans, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_PQM_NULL_LIMIT .OR. & advectionScheme.EQ.ENUM_PQM_MONO_LIMIT .OR. & advectionScheme.EQ.ENUM_PQM_WENO_LIMIT) THEN CALL GAD_PQM_ADV_X( advectionScheme, bi, bj, k , .TRUE., I deltaTLev(k), uFld(1-OLx,1-OLy,k), uTrans, localTij, O af, myThid ) #endif /* ndef ALLOW_AUTODIFF */ ELSE STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim' ENDIF #ifdef ALLOW_OBCS IF ( useOBCS ) THEN C- replace advective flux with 1st order upwind scheme estimate CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k, I maskLocW, uTrans, localTij, U af, myThid ) ENDIF #endif /* ALLOW_OBCS */ C- Internal exchange for next calculations in Y IF ( overlapOnly .AND. ipass.EQ.1 ) THEN CALL FILL_CS_CORNER_TR_RL( 2, .FALSE., & localTij, bi,bj, myThid ) ENDIF C- Advective flux in X : done ENDIF C- Update the local tracer field where needed: C use "maksInC" to prevent updating tracer field in OB regions #ifdef ALLOW_AUTODIFF_TAMC # ifdef GAD_MULTIDIM_COMPRESSIBLE CADJ STORE localVol(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte # endif #endif /* ALLOW_AUTODIFF_TAMC */ C update in overlap-Only IF ( overlapOnly ) THEN iMinUpd = 1-OLx+1 iMaxUpd = sNx+OLx-1 C- notes: these 2 lines below have no real effect (because recip_hFac=0 C in corner region) but safer to keep them. IF ( W_edge ) iMinUpd = 1 IF ( E_edge ) iMaxUpd = sNx IF ( S_edge ) THEN DO j=1-OLy,0 DO i=iMinUpd,iMaxUpd #ifdef GAD_MULTIDIM_COMPRESSIBLE tmpTrac = localTij(i,j)*localVol(i,j) & -deltaTLev(k)*( af(i+1,j) - af(i,j) ) & *maskInC(i,j,bi,bj) localVol(i,j) = localVol(i,j) & -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) ) & *maskInC(i,j,bi,bj) localTij(i,j) = tmpTrac/localVol(i,j) #else /* GAD_MULTIDIM_COMPRESSIBLE */ localTij(i,j) = localTij(i,j) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( af(i+1,j)-af(i,j) & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) & )*maskInC(i,j,bi,bj) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ ENDDO ENDDO ENDIF IF ( N_edge ) THEN DO j=sNy+1,sNy+OLy DO i=iMinUpd,iMaxUpd #ifdef GAD_MULTIDIM_COMPRESSIBLE tmpTrac = localTij(i,j)*localVol(i,j) & -deltaTLev(k)*( af(i+1,j) - af(i,j) ) & *maskInC(i,j,bi,bj) localVol(i,j) = localVol(i,j) & -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) ) & *maskInC(i,j,bi,bj) localTij(i,j) = tmpTrac/localVol(i,j) #else /* GAD_MULTIDIM_COMPRESSIBLE */ localTij(i,j) = localTij(i,j) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( af(i+1,j)-af(i,j) & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) & )*maskInC(i,j,bi,bj) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ ENDDO ENDDO ENDIF ELSE C do not only update the overlap jMinUpd = 1-OLy jMaxUpd = sNy+OLy IF ( interiorOnly .AND. S_edge ) jMinUpd = 1 IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy DO j=jMinUpd,jMaxUpd DO i=1-OLx+1,sNx+OLx-1 #ifdef GAD_MULTIDIM_COMPRESSIBLE tmpTrac = localTij(i,j)*localVol(i,j) & -deltaTLev(k)*( af(i+1,j) - af(i,j) ) & *maskInC(i,j,bi,bj) localVol(i,j) = localVol(i,j) & -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) ) & *maskInC(i,j,bi,bj) localTij(i,j) = tmpTrac/localVol(i,j) #else /* GAD_MULTIDIM_COMPRESSIBLE */ localTij(i,j) = localTij(i,j) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( af(i+1,j)-af(i,j) & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) & )*maskInC(i,j,bi,bj) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ ENDDO ENDDO C- keep advective flux (for diagnostics) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx afx(i,j) = af(i,j) ENDDO ENDDO C- end if/else update overlap-Only ENDIF C-- End of X direction ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Y direction #ifdef ALLOW_AUTODIFF C- Always reset advective flux in Y DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx af(i,j) = 0. ENDDO ENDDO # ifndef DISABLE_MULTIDIM_ADVECTION CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte CADJ STORE af(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte # endif #endif /* ALLOW_AUTODIFF */ IF (calc_fluxes_Y) THEN C- Do not compute fluxes if C a) needed in overlap only C and b) the overlap of myTile are not cube-face edges IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN C- Internal exchange for calculations in Y IF ( overlapOnly ) THEN CALL FILL_CS_CORNER_TR_RL( 2, .FALSE., & localTij, bi,bj, myThid ) ENDIF C- Advective flux in Y #ifndef ALLOW_AUTODIFF DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx af(i,j) = 0. ENDDO ENDDO #else /* ALLOW_AUTODIFF */ #ifndef DISABLE_MULTIDIM_ADVECTION CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte #endif #endif /* ALLOW_AUTODIFF */ IF ( advectionScheme.EQ.ENUM_UPWIND_1RST & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE., I deltaTLev(k),vTrans,vFld(1-OLx,1-OLy,k), localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k), I vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k), I vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k), I vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij, O af, myThid ) #ifndef ALLOW_AUTODIFF ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k), I vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_PPM_NULL_LIMIT .OR. & advectionScheme.EQ.ENUM_PPM_MONO_LIMIT .OR. & advectionScheme.EQ.ENUM_PPM_WENO_LIMIT) THEN CALL GAD_PPM_ADV_Y(advectionScheme, bi, bj, k , .TRUE., I deltaTLev(k), vFld(1-OLX,1-OLy,k), vTrans, localTij, O af, myThid ) ELSEIF (advectionScheme.EQ.ENUM_PQM_NULL_LIMIT .OR. & advectionScheme.EQ.ENUM_PQM_MONO_LIMIT .OR. & advectionScheme.EQ.ENUM_PQM_WENO_LIMIT) THEN CALL GAD_PQM_ADV_Y(advectionScheme, bi, bj, k , .TRUE., I deltaTLev(k), vFld(1-OLX,1-OLy,k), vTrans, localTij, O af, myThid ) #endif /* ndef ALLOW_AUTODIFF */ ELSE STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' ENDIF #ifdef ALLOW_OBCS IF ( useOBCS ) THEN C- replace advective flux with 1st order upwind scheme estimate CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k, I maskLocS, vTrans, localTij, U af, myThid ) ENDIF #endif /* ALLOW_OBCS */ C- Internal exchange for next calculations in X IF ( overlapOnly .AND. ipass.EQ.1 ) THEN CALL FILL_CS_CORNER_TR_RL( 1, .FALSE., & localTij, bi,bj, myThid ) ENDIF C- Advective flux in Y : done ENDIF C- Update the local tracer field where needed: C use "maksInC" to prevent updating tracer field in OB regions #ifdef ALLOW_AUTODIFF_TAMC # ifdef GAD_MULTIDIM_COMPRESSIBLE CADJ STORE localVol(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte # endif #endif /* ALLOW_AUTODIFF_TAMC */ C update in overlap-Only IF ( overlapOnly ) THEN jMinUpd = 1-OLy+1 jMaxUpd = sNy+OLy-1 C- notes: these 2 lines below have no real effect (because recip_hFac=0 C in corner region) but safer to keep them. IF ( S_edge ) jMinUpd = 1 IF ( N_edge ) jMaxUpd = sNy IF ( W_edge ) THEN DO j=jMinUpd,jMaxUpd DO i=1-OLx,0 #ifdef GAD_MULTIDIM_COMPRESSIBLE tmpTrac = localTij(i,j)*localVol(i,j) & -deltaTLev(k)*( af(i,j+1) - af(i,j) ) & *maskInC(i,j,bi,bj) localVol(i,j) = localVol(i,j) & -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) ) & *maskInC(i,j,bi,bj) localTij(i,j) = tmpTrac/localVol(i,j) #else /* GAD_MULTIDIM_COMPRESSIBLE */ localTij(i,j) = localTij(i,j) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( af(i,j+1)-af(i,j) & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) & )*maskInC(i,j,bi,bj) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ ENDDO ENDDO ENDIF IF ( E_edge ) THEN DO j=jMinUpd,jMaxUpd DO i=sNx+1,sNx+OLx #ifdef GAD_MULTIDIM_COMPRESSIBLE tmpTrac = localTij(i,j)*localVol(i,j) & -deltaTLev(k)*( af(i,j+1) - af(i,j) ) & *maskInC(i,j,bi,bj) localVol(i,j) = localVol(i,j) & -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) ) & *maskInC(i,j,bi,bj) localTij(i,j) = tmpTrac/localVol(i,j) #else /* GAD_MULTIDIM_COMPRESSIBLE */ localTij(i,j) = localTij(i,j) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( af(i,j+1)-af(i,j) & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) & )*maskInC(i,j,bi,bj) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ ENDDO ENDDO ENDIF ELSE C do not only update the overlap iMinUpd = 1-OLx iMaxUpd = sNx+OLx IF ( interiorOnly .AND. W_edge ) iMinUpd = 1 IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx DO j=1-OLy+1,sNy+OLy-1 DO i=iMinUpd,iMaxUpd #ifdef GAD_MULTIDIM_COMPRESSIBLE tmpTrac = localTij(i,j)*localVol(i,j) & -deltaTLev(k)*( af(i,j+1) - af(i,j) ) & *maskInC(i,j,bi,bj) localVol(i,j) = localVol(i,j) & -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) ) & *maskInC(i,j,bi,bj) localTij(i,j) = tmpTrac/localVol(i,j) #else /* GAD_MULTIDIM_COMPRESSIBLE */ localTij(i,j) = localTij(i,j) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( af(i,j+1)-af(i,j) & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) & )*maskInC(i,j,bi,bj) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ ENDDO ENDDO C- keep advective flux (for diagnostics) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx afy(i,j) = af(i,j) ENDDO ENDDO C end if/else update overlap-Only ENDIF C-- End of Y direction ENDIF C-- End of ipass loop ENDDO IF ( implicitAdvection ) THEN C- explicit advection is done ; store tendency in gTracer: #ifdef GAD_MULTIDIM_COMPRESSIBLE STOP 'GAD_ADVECTION: missing code for implicitAdvection' #endif /* GAD_MULTIDIM_COMPRESSIBLE */ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx gTracer(i,j,k) = & ( localTij(i,j) - tracer(i,j,k,bi,bj) )/deltaTLev(k) ENDDO ENDDO ELSE C- horizontal advection done; store intermediate result in 3D array: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx #ifdef GAD_MULTIDIM_COMPRESSIBLE locVol3d(i,j,k) = localVol(i,j) #endif /* GAD_MULTIDIM_COMPRESSIBLE */ localT3d(i,j,k) = localTij(i,j) ENDDO ENDDO ENDIF #ifdef ALLOW_DIAGNOSTICS IF ( doDiagAdvX ) THEN diagName = 'ADVx'//diagSufx CALL DIAGNOSTICS_FILL( afx, diagName, k,1, 2,bi,bj, myThid ) ENDIF IF ( doDiagAdvY ) THEN diagName = 'ADVy'//diagSufx CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid ) ENDIF #ifdef ALLOW_LAYERS IF ( useLayers ) THEN CALL LAYERS_FILL(afx,trIdentity,'AFX',k,1,2,bi,bj,myThid) CALL LAYERS_FILL(afy,trIdentity,'AFY',k,1,2,bi,bj,myThid) ENDIF #endif /* ALLOW_LAYERS */ #endif /* ALLOW_DIAGNOSTICS */ #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevC & .AND. trIdentity.EQ.GAD_TEMPERATURE & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0 & .AND. nPx.EQ.1 .AND. nPy.EQ.1 & .AND. useCubedSphereExchange ) THEN CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION', & afx,afy, k, standardMessageUnit,bi,bj,myThid ) ENDIF #endif /* ALLOW_DEBUG */ C-- End of K loop for horizontal fluxes ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( .NOT.implicitAdvection ) THEN #ifndef ALLOW_AUTODIFF IF (vertAdvecScheme.EQ.ENUM_PPM_NULL_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PPM_MONO_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PPM_WENO_LIMIT) THEN C-- PPM-style vertical advection DO k=1,Nr IF (k.EQ.1) THEN C-- vertical transport: surface DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTran3d(i,j,k) = 0. _d 0 ENDDO ENDDO ELSE C-- vertical transport: interior DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTran3d(i,j,k) = wFld(i,j,k)*rA(i,j,bi,bj) & *deepFac2F(k)*rhoFacF(k) & *maskC(i,j,k-1,bi,bj) ENDDO ENDDO ENDIF ENDDO C-- calc. PPM vertical flux data CALL GAD_PPM_ADV_R( vertAdvecScheme, bi, bj, I deltaTLev, wFld, rTran3d, localT3d, O afr, myThid ) ENDIF IF (vertAdvecScheme.EQ.ENUM_PQM_NULL_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PQM_MONO_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PQM_WENO_LIMIT) THEN C-- PQM-style vertical advection DO k=1,Nr IF (k.EQ.1) THEN C-- vertical transport: surface DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTran3d(i,j,k) = 0. _d 0 ENDDO ENDDO ELSE C-- vertical transport: interior DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTran3d(i,j,k) = wFld(i,j,k)*rA(i,j,bi,bj) & *deepFac2F(k)*rhoFacF(k) & *maskC(i,j,k-1,bi,bj) ENDDO ENDDO ENDIF ENDDO C-- calc. PQM vertical flux data CALL GAD_PQM_ADV_R( vertAdvecScheme, bi, bj, I deltaTLev, wFld, rTran3d, localT3d, O afr, myThid ) ENDIF #endif /* ndef ALLOW_AUTODIFF */ C-- Start of k loop for vertical flux DO k=Nr,1,-1 #ifdef ALLOW_AUTODIFF_TAMC kkey = (igadkey-1)*Nr + (Nr-k+1) #endif /* ALLOW_AUTODIFF_TAMC */ C-- kUp Cycles through 1,2 to point to w-layer above C-- kDown Cycles through 2,1 to point to w-layer below kUp = 1+MOD(k+1,2) kDown= 1+MOD(k,2) kp1Msk=1. IF (k.EQ.Nr) kp1Msk=0. #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE rtrans(:,:) = CADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte #endif C-- Compute Vertical transport #ifdef ALLOW_AIM C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr IF ( k.EQ.1 .OR. & (useAIM .AND. trIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr) & ) THEN #else IF ( k.EQ.1 ) THEN #endif C- Surface interface : DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTransKp(i,j) = kp1Msk*rTrans(i,j) rTrans(i,j) = 0. fVerT(i,j,kUp) = 0. ENDDO ENDDO ELSE C- Interior interface : DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTransKp(i,j) = kp1Msk*rTrans(i,j) rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj) & *deepFac2F(k)*rhoFacF(k) & *maskC(i,j,k-1,bi,bj) fVerT(i,j,kUp) = 0. ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC cphmultiCADJ STORE localT3d(:,:,k) cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte cphmultiCADJ STORE rTrans(:,:) cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C- Compute vertical advective flux in the interior: IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST & .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme, I deltaTLev(k),rTrans,wFld(1-OLx,1-OLy,k),localT3d, O fVerT(1-OLx,1-OLy,kUp), myThid ) ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k), I rTrans, wFld(1-OLx,1-OLy,k), localT3d, O fVerT(1-OLx,1-OLy,kUp), myThid ) ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k), I rTrans, wFld(1-OLx,1-OLy,k), localT3d, O fVerT(1-OLx,1-OLy,kUp), myThid ) ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k), I rTrans, wFld(1-OLx,1-OLy,k), localT3d, O fVerT(1-OLx,1-OLy,kUp), myThid ) #ifndef ALLOW_AUTODIFF ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k), I rTrans, wFld(1-OLx,1-OLy,k), localT3d, O fVerT(1-OLx,1-OLy,kUp), myThid ) ELSEIF (vertAdvecScheme.EQ.ENUM_PPM_NULL_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PPM_MONO_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PPM_WENO_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PQM_NULL_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PQM_MONO_LIMIT .OR. & vertAdvecScheme.EQ.ENUM_PQM_WENO_LIMIT) THEN C- copy level from 3d flux data DO j = 1-OLy,sNy+OLy DO i = 1-OLx,sNx+OLx fVerT(i,j,kUp) = afr(i,j,k) ENDDO ENDDO #endif /* ndef ALLOW_AUTODIFF */ ELSE STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' ENDIF C- end Surface/Interior if bloc ENDIF #ifdef ALLOW_AUTODIFF_TAMC cphmultiCADJ STORE rTrans(:,:) cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte cphmultiCADJ STORE rTranskp(:,:) cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte cph --- following storing of fVerT is critical for correct cph --- gradient with multiDimAdvection cph --- Without it, kDown component is not properly recomputed cph --- This is a TAF bug (and no warning available) CADJ STORE fVerT(:,:,:) CADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C-- Divergence of vertical fluxes #ifdef GAD_MULTIDIM_COMPRESSIBLE DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx tmpTrac = localT3d(i,j,k)*locVol3d(i,j,k) & -deltaTLev(k)*( fVerT(i,j,kDown)-fVerT(i,j,kUp) ) & *rkSign*maskInC(i,j,bi,bj) localVol(i,j) = locVol3d(i,j,k) & -deltaTLev(k)*( rTransKp(i,j) - rTrans(i,j) ) & *rkSign*maskInC(i,j,bi,bj) C- localTij only needed for Variance Bugget: can be move there localTij(i,j) = tmpTrac/localVol(i,j) C-- without rescaling of tendencies: c gTracer(i,j,k) = c & ( localTij(i,j) - tracer(i,j,k,bi,bj) )/deltaTLev(k) C-- Non-Lin Free-Surf: consistent with rescaling of tendencies C (in FREESURF_RESCALE_G) and RealFreshFlux/addMass. C Also valid for linear Free-Surf (r & r* coords) w/wout RealFreshFlux C and consistent with linFSConserveTr and "surfExpan_" monitor. gTracer(i,j,k) = & ( tmpTrac - tracer(i,j,k,bi,bj)*localVol(i,j) ) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj) & *recip_rhoFacC(k) & /deltaTLev(k) ENDDO ENDDO #else /* GAD_MULTIDIM_COMPRESSIBLE */ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx localTij(i,j) = localT3d(i,j,k) & -deltaTLev(k)*recip_rhoFacC(k) & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *( fVerT(i,j,kDown)-fVerT(i,j,kUp) & -tracer(i,j,k,bi,bj)*(rTransKp(i,j)-rTrans(i,j)) & )*rkSign*maskInC(i,j,bi,bj) gTracer(i,j,k) = & ( localTij(i,j) - tracer(i,j,k,bi,bj) )/deltaTLev(k) ENDDO ENDDO #endif /* GAD_MULTIDIM_COMPRESSIBLE */ #ifdef ALLOW_DIAGNOSTICS IF ( doDiagAdvR ) THEN diagName = 'ADVr'//diagSufx CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp), & diagName, k,1, 2,bi,bj, myThid ) ENDIF #ifdef ALLOW_LAYERS IF ( useLayers ) THEN CALL LAYERS_FILL( fVerT(1-OLx,1-OLy,kUp), trIdentity, & 'AFR', k, 1, 2,bi,bj, myThid) ENDIF #endif /* ALLOW_LAYERS */ #endif /* ALLOW_DIAGNOSTICS */ C-- End of K loop for vertical flux ENDDO C-- end of if not.implicitAdvection block ENDIF RETURN END