[MITgcm-devel] Fwd: TAF issue for "special" parameter lists
Chris Hill
cnh at mit.edu
Fri Feb 24 21:03:33 EST 2006
P. (JM),
It would be better not do it this way even for forward code.
The reason is that certain compiler optimizations will assume "no
aliasing" which means assuming that localT and locaABT refer to
different pieces of memory and therefore certain transformations that
safely permute the ordering of computations are allowed. Aliased
arguments can lead to extremely hard to track down bugs suddenly
appearing with some compiler all of a sudden. Some compilers provide
assume_alising flags but that tends to harm performance since it assumes
just about everything might be aliased with everything and the compiler
gets paranoid about lots of useful optimizations.
So far we have avoided any aliasing in the core code (not sure about
sea-ice, fizhi etc....) for just this reason. It would be much better to
continue the practice of not having aliased dummy arguments even if it
means adding some more logic in gad_calc_rhs.
Chris
Patrick Heimbach wrote:
>
> ----- Forwarded message from heimbach at MIT.EDU -----
> Date: Fri, 24 Feb 2006 17:57:15 -0500
> From: Patrick Heimbach <heimbach at MIT.EDU>
> Reply-To: heimbach at MIT.EDU
> Subject: TAF issue for "special" parameter lists
> To: Ralf Giering <Ralf.Giering at FastOpt.com>
>
>
> Ralf,
>
> some of my colleagues have recently introduced a construct
> into the model which breaks the adjoint.
> It's not very clean programming, but works well for them
> in pure forward to be able to switch between different
> variations of Adams-Bashforth.
>
> The general structure is as following (the underlying code
> is calls of S/R gad_calc_rhs by S/R calc_gt, S/R calc_gs
> for the case NOT AB_3)
>
> call foo( ..., var, var, ... )
>
> subroutine foo( ..., var1, var2, ... )
>
> u = f1(var1)
> v = f2(var2)
>
> var is pure input to S/R foo, not modified.
> The reason for doing this is to be able to use the same
> argument list for a variation of AB (the AB_3) for which
> call gad_calc_rhs( ..., theta, theta, ... ) ! AB_2
> is modified into
> call gad_calc_rhs( ..., gTnm1, theta, ... ) ! AB_3
>
> Sensitivities for this construct are not
> properly accumulated, here's why:
>
> call adfoo( ..., advar, advar, ... )
>
> subroutine adfoo( ..., advar1, advar2, ...)
>
> advar2 = advar2 + adf2(advar2)
> advar1 = advar1 + adf1(advar1)
>
> The link of accumulating advar1 and advar2 onto advar
> (because they're the same) is missing; strictly:
> advar = advar1 + advar2
>
> I have actually hand-coded this variation, i.e.
> introduce an auxiliary variable advarh in the call
> and then accumulate advarh onto advar, something like
> call adfoo( ..., advar, advarh, ... )
> advar = advar + advarh
> advarh = 0.
> It works.
> This seems to show a strategy for TAF to properly
> handle the construct.
>
> The effect is subtle, not visible at first (and in
> gradient checks), but after 13yr integration, the "ECCO"
> gradient is messed up.
>
> I am appending the routines where this problem occurs:
>
> calc_gt.F (has the call to "foo" = gad_calc_rhs with
> the problematic parameter list)
> gad_calc_rhs.F ( = "foo" where theta=Tracer, and theta=TracAB)
> adgad_calc_rhs (just to show how adTracer and adTracAB
> are treated separately)
> adcalc_gt_broken (the code generated by TAF,
> it has adtheta twice in param. list)
> adcalc_gt_fixed (I introduced on auxiliary variable adthetah,
> which after the call gets added onto adtheta).
>
> Hope this is enough info.
> I think you'll quickly recognize the general problem.
> It's not urgent, we'll re-write with a workaround,
> but I think, you can modify TAF in above way to
> handle this, provided TAF is able to recognize that
> theta is twice in the argument list
> (which might not be easy/general though).
>
> Cheers
> -Patrick
>
>
>
> ------------------------------------------------------------------------
>
> C $Header: /u/gcmpack/MITgcm/model/src/calc_gt.F,v 1.46 2005/12/08 15:44:33 heimbach Exp $
> C $Name: $
>
> #include "PACKAGES_CONFIG.h"
> #include "CPP_OPTIONS.h"
>
> CBOP
> C !ROUTINE: CALC_GT
> C !INTERFACE:
> SUBROUTINE CALC_GT(
> I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
> I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
> I KappaRT,
> U fVerT,
> I myTime,myIter,myThid )
> C !DESCRIPTION: \bv
> C *==========================================================*
> C | SUBROUTINE CALC_GT
> C | o Calculate the temperature tendency terms.
> C *==========================================================*
> C | A procedure called EXTERNAL_FORCING_T is called from
> C | here. These procedures can be used to add per problem
> C | heat flux source terms.
> C | Note: Although it is slightly counter-intuitive the
> C | EXTERNAL_FORCING routine is not the place to put
> C | file I/O. Instead files that are required to
> C | calculate the external source terms are generally
> C | read during the model main loop. This makes the
> C | logisitics of multi-processing simpler and also
> C | makes the adjoint generation simpler. It also
> C | allows for I/O to overlap computation where that
> C | is supported by hardware.
> C | Aside from the problem specific term the code here
> C | forms the tendency terms due to advection and mixing
> C | The baseline implementation here uses a centered
> C | difference form for the advection term and a tensorial
> C | divergence of a flux form for the diffusive term. The
> C | diffusive term is formulated so that isopycnal mixing and
> C | GM-style subgrid-scale terms can be incorporated b simply
> C | setting the diffusion tensor terms appropriately.
> C *==========================================================*
> C \ev
>
> C !USES:
> IMPLICIT NONE
> C == GLobal variables ==
> #include "SIZE.h"
> #include "DYNVARS.h"
> #include "EEPARAMS.h"
> #include "PARAMS.h"
> #ifdef ALLOW_GENERIC_ADVDIFF
> #include "GAD.h"
> #endif
> #ifdef ALLOW_AUTODIFF_TAMC
> # include "tamc.h"
> # include "tamc_keys.h"
> #endif
>
> C !INPUT/OUTPUT PARAMETERS:
> C == Routine arguments ==
> C fVerT :: Flux of temperature (T) in the vertical
> C direction at the upper(U) and lower(D) faces of a cell.
> C maskUp :: Land mask used to denote base of the domain.
> C xA :: Tracer cell face area normal to X
> C yA :: Tracer cell face area normal to X
> C uTrans :: Zonal volume transport through cell face
> C vTrans :: Meridional volume transport through cell face
> C rTrans :: Vertical volume transport at interface k
> C rTransKp1 :: Vertical volume transport at inteface k+1
> C bi, bj, iMin, iMax, jMin, jMax :: Range of points for which calculation
> C results will be set.
> C myThid :: Instance number for this innvocation of CALC_GT
> _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
> _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 rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> INTEGER k,kUp,kDown,kM1
> INTEGER bi,bj,iMin,iMax,jMin,jMax
> _RL myTime
> INTEGER myIter
> INTEGER myThid
> CEOP
>
> #ifdef ALLOW_GENERIC_ADVDIFF
> C === Local variables ===
> LOGICAL calcAdvection
> INTEGER iterNb
> #ifdef ALLOW_ADAMSBASHFORTH_3
> INTEGER m1, m2
> #endif
>
> #ifdef ALLOW_AUTODIFF_TAMC
> 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
> itdkey = (act1 + 1) + act2*max1
> & + act3*max1*max2
> & + act4*max1*max2*max3
> kkey = (itdkey-1)*Nr + k
> #endif /* ALLOW_AUTODIFF_TAMC */
>
> #ifdef ALLOW_AUTODIFF_TAMC
> C-- only the kUp part of fverT is set in this subroutine
> C-- the kDown is still required
> fVerT(1,1,kDown) = fVerT(1,1,kDown)
> # ifdef NONLIN_FRSURF
> CADJ STORE fVerT(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
> # endif
> #endif
>
> C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
>
> calcAdvection = tempAdvection .AND. .NOT.tempMultiDimAdvec
> iterNb = myIter
> IF (staggerTimeStep) iterNb = myIter -1
>
> #ifdef ALLOW_ADAMSBASHFORTH_3
> IF ( AdamsBashforth_T ) THEN
> m1 = 1 + MOD(iterNb+1,2)
> m2 = 1 + MOD( iterNb ,2)
> CALL GAD_CALC_RHS(
> I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
> I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
> I uVel, vVel, wVel,
> I diffKhT, diffK4T, KappaRT,
> I gtNm(1-Olx,1-Oly,1,1,1,m2), theta,
> I GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
> I calcAdvection, tempImplVertAdv,
> U fVerT, gT,
> I myTime, myIter, myThid )
> ELSE
> #endif /* ALLOW_ADAMSBASHFORTH_3 */
> CALL GAD_CALC_RHS(
> I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
> I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
> I uVel, vVel, wVel,
> I diffKhT, diffK4T, KappaRT, theta, theta,
> I GAD_TEMPERATURE, tempAdvScheme, tempVertAdvScheme,
> I calcAdvection, tempImplVertAdv,
> U fVerT, gT,
> I myTime, myIter, myThid )
> #ifdef ALLOW_ADAMSBASHFORTH_3
> ENDIF
> #endif
>
> C-- External thermal forcing term(s) inside Adams-Bashforth:
> IF ( tempForcing .AND. forcing_In_AB )
> & CALL EXTERNAL_FORCING_T(
> I iMin,iMax,jMin,jMax,bi,bj,k,
> I myTime,myThid)
>
> IF ( AdamsBashforthGt ) THEN
> #ifdef ALLOW_ADAMSBASHFORTH_3
> CALL ADAMS_BASHFORTH3(
> I bi, bj, k,
> U gT, gtNm,
> I tempStartAB, iterNb, myThid )
> #else
> CALL ADAMS_BASHFORTH2(
> I bi, bj, k,
> U gT, gtNm1,
> I iterNb, myThid )
> #endif
> ENDIF
>
> C-- External thermal forcing term(s) outside Adams-Bashforth:
> IF ( tempForcing .AND. .NOT.forcing_In_AB )
> & CALL EXTERNAL_FORCING_T(
> I iMin,iMax,jMin,jMax,bi,bj,k,
> I myTime,myThid)
>
> #ifdef NONLIN_FRSURF
> IF (nonlinFreeSurf.GT.0) THEN
> CALL FREESURF_RESCALE_G(
> I bi, bj, k,
> U gT,
> I myThid )
> IF ( AdamsBashforthGt ) THEN
> #ifdef ALLOW_ADAMSBASHFORTH_3
> CALL FREESURF_RESCALE_G(
> I bi, bj, k,
> U gtNm(1-Olx,1-Oly,1,1,1,1),
> I myThid )
> CALL FREESURF_RESCALE_G(
> I bi, bj, k,
> U gtNm(1-Olx,1-Oly,1,1,1,2),
> I myThid )
> #else
> CALL FREESURF_RESCALE_G(
> I bi, bj, k,
> U gtNm1,
> I myThid )
> #endif
> ENDIF
> ENDIF
> #endif /* NONLIN_FRSURF */
>
> #endif /* ALLOW_GENERIC_ADVDIFF */
>
> RETURN
> END
>
>
> ------------------------------------------------------------------------
>
> C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.38 2005/11/06 22:14:02 jmc Exp $
> C $Name: $
>
> #include "GAD_OPTIONS.h"
>
> CBOP
> C !ROUTINE: GAD_CALC_RHS
>
> C !INTERFACE: ==========================================================
> SUBROUTINE GAD_CALC_RHS(
> I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
> I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
> I uVel, vVel, wVel,
> I diffKh, diffK4, KappaR, Tracer, TracAB,
> I tracerIdentity, advectionScheme, vertAdvecScheme,
> I calcAdvection, implicitAdvection,
> U fVerT, gTracer,
> I myTime, myIter, myThid )
>
> C !DESCRIPTION:
> C Calculates the tendency of a tracer due to advection and diffusion.
> C It calculates the fluxes in each direction indepentently and then
> C sets the tendency to the divergence of these fluxes. The advective
> C fluxes are only calculated here when using the linear advection schemes
> C otherwise only the diffusive and parameterized fluxes are calculated.
> C
> C Contributions to the flux are calculated and added:
> C \begin{equation*}
> C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
> C \end{equation*}
> C
> C The tendency is the divergence of the fluxes:
> C \begin{equation*}
> C G_\theta = G_\theta + \nabla \cdot {\bf F}
> C \end{equation*}
> C
> C The tendency is assumed to contain data on entry.
>
> C !USES: ===============================================================
> IMPLICIT NONE
> #include "SIZE.h"
> #include "EEPARAMS.h"
> #include "PARAMS.h"
> #include "GRID.h"
> #include "SURFACE.h"
> #include "GAD.h"
>
> #ifdef ALLOW_AUTODIFF_TAMC
> #include "tamc.h"
> #include "tamc_keys.h"
> #endif /* ALLOW_AUTODIFF_TAMC */
>
> C !INPUT PARAMETERS: ===================================================
> C bi,bj :: tile indices
> C iMin,iMax :: loop range for called routines
> C jMin,jMax :: loop range for called routines
> C kup :: index into 2 1/2D array, toggles between 1|2
> C kdown :: index into 2 1/2D array, toggles between 2|1
> C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
> 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 rTransKp1 :: 2-D array of volume trans at W pts, interf k+1
> C maskUp :: 2-D array for mask at W points
> C uVel,vVel,wVel :: 3 components of the velcity field (3-D array)
> C diffKh :: horizontal diffusion coefficient
> C diffK4 :: bi-harmonic diffusion coefficient
> C KappaR :: 2-D array for vertical diffusion coefficient, interf k
> C Tracer :: tracer field
> C TracAB :: tracer field but extrapolated fwd in time (to nIter+1/2)
> C if applying AB on tracer field (instead of on tendency)
> C tracerIdentity :: tracer identifier (required for KPP,GM)
> C advectionScheme :: advection scheme to use (Horizontal plane)
> C vertAdvecScheme :: advection scheme to use (Vertical direction)
> C calcAdvection :: =False if Advec computed with multiDim scheme
> C implicitAdvection:: =True if vertical Advec computed implicitly
> C myTime :: current time
> C myIter :: iteration number
> C myThid :: thread number
> INTEGER bi,bj,iMin,iMax,jMin,jMax
> INTEGER k,kUp,kDown,kM1
> _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 rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> _RL diffKh, diffK4
> _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> _RL TracAB(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> INTEGER tracerIdentity
> INTEGER advectionScheme, vertAdvecScheme
> LOGICAL calcAdvection
> LOGICAL implicitAdvection
> _RL myTime
> INTEGER myIter, myThid
>
> C !OUTPUT PARAMETERS: ==================================================
> C gTracer :: tendency array
> C fVerT :: 2 1/2D arrays for vertical advective flux
> _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
>
> C !LOCAL VARIABLES: ====================================================
> C i,j :: loop indices
> C df4 :: used for storing del^2 T for bi-harmonic term
> C fZon :: zonal flux
> C fMer :: meridional flux
> C af :: advective flux
> C df :: diffusive flux
> C localT :: local copy of tracer field
> C locABT :: local copy of (AB-extrapolated) tracer field
> #ifdef ALLOW_DIAGNOSTICS
> CHARACTER*8 diagName
> CHARACTER*4 GAD_DIAG_SUFX, diagSufx
> EXTERNAL GAD_DIAG_SUFX
> #endif
> INTEGER i,j
> _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RL advFac, rAdvFac
> CEOP
>
> #ifdef ALLOW_AUTODIFF_TAMC
> C-- only the kUp part of fverT is set in this subroutine
> C-- the kDown is still required
> fVerT(1,1,kDown) = fVerT(1,1,kDown)
> #endif
>
> #ifdef ALLOW_DIAGNOSTICS
> C-- Set diagnostic suffix for the current tracer
> IF ( useDiagnostics ) THEN
> diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
> ENDIF
> #endif
>
> advFac = 0. _d 0
> IF (calcAdvection) advFac = 1. _d 0
> rAdvFac = rkSign*advFac
> IF (implicitAdvection) rAdvFac = 0. _d 0
>
> DO j=1-OLy,sNy+OLy
> DO i=1-OLx,sNx+OLx
> fZon(i,j) = 0. _d 0
> fMer(i,j) = 0. _d 0
> fVerT(i,j,kUp) = 0. _d 0
> df(i,j) = 0. _d 0
> df4(i,j) = 0. _d 0
> ENDDO
> ENDDO
>
> C-- Make local copy of tracer array
> DO j=1-OLy,sNy+OLy
> DO i=1-OLx,sNx+OLx
> localT(i,j)=Tracer(i,j,k,bi,bj)
> locABT(i,j)=TracAB(i,j,k,bi,bj)
> ENDDO
> ENDDO
>
> C-- Unless we have already calculated the advection terms we initialize
> C the tendency to zero.
> C <== now done earlier at the beginning of thermodynamics.
> c IF (calcAdvection) THEN
> c DO j=1-Oly,sNy+Oly
> c DO i=1-Olx,sNx+Olx
> c gTracer(i,j,k,bi,bj)=0. _d 0
> c ENDDO
> c ENDDO
> c ENDIF
>
> C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
> IF (diffK4 .NE. 0.) THEN
> CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
> CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
> CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
> ENDIF
>
> C-- Initialize net flux in X direction
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fZon(i,j) = 0. _d 0
> ENDDO
> ENDDO
>
> C- Advective flux in X
> IF (calcAdvection) THEN
> IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
> CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
> ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
> & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
> CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
> I dTtracerLev(k), uTrans, uVel, locABT,
> O af, myThid )
> ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
> CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
> I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
> CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
> ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
> CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
> ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
> CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
> I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
> IF ( inAdMode ) THEN
> cph This block is to trick the adjoint:
> cph IF inAdExact=.FALSE., we want to use DST3
> cph with limiters in forward, but without limiters in reverse.
> CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
> I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ELSE
> CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
> I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ENDIF
> ELSE
> STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
> ENDIF
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fZon(i,j) = fZon(i,j) + af(i,j)
> ENDDO
> ENDDO
> #ifdef ALLOW_DIAGNOSTICS
> IF ( useDiagnostics ) THEN
> diagName = 'ADVx'//diagSufx
> CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
> ENDIF
> #endif
> ENDIF
>
> C- Diffusive flux in X
> IF (diffKh.NE.0.) THEN
> CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
> ELSE
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> df(i,j) = 0. _d 0
> ENDDO
> ENDDO
> ENDIF
>
> C- Add bi-harmonic diffusive flux in X
> IF (diffK4 .NE. 0.) THEN
> CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
> ENDIF
>
> #ifdef ALLOW_GMREDI
> C- GM/Redi flux in X
> IF (useGMRedi) THEN
> C *note* should update GMREDI_XTRANSPORT to set df *aja*
> CALL GMREDI_XTRANSPORT(
> I iMin,iMax,jMin,jMax,bi,bj,K,
> I xA,Tracer,tracerIdentity,
> U df,
> I myThid)
> ENDIF
> #endif
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fZon(i,j) = fZon(i,j) + df(i,j)
> ENDDO
> ENDDO
>
> #ifdef ALLOW_DIAGNOSTICS
> C- Diagnostics of Tracer flux in X dir (mainly Diffusive term),
> C excluding advective terms:
> IF ( useDiagnostics .AND.
> & (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
> diagName = 'DIFx'//diagSufx
> CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
> ENDIF
> #endif
>
> C-- Initialize net flux in Y direction
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fMer(i,j) = 0. _d 0
> ENDDO
> ENDDO
>
> C- Advective flux in Y
> IF (calcAdvection) THEN
> IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
> CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
> ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
> & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
> CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
> I dTtracerLev(k), vTrans, vVel, locABT,
> O af, myThid )
> ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
> CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
> I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
> CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
> ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
> CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
> ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
> CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
> I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
> IF ( inAdMode ) THEN
> cph This block is to trick the adjoint:
> cph IF inAdExact=.FALSE., we want to use DST3
> cph with limiters in forward, but without limiters in reverse.
> CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
> I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ELSE
> CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
> I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
> O af, myThid )
> ENDIF
> ELSE
> STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
> ENDIF
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fMer(i,j) = fMer(i,j) + af(i,j)
> ENDDO
> ENDDO
> #ifdef ALLOW_DIAGNOSTICS
> IF ( useDiagnostics ) THEN
> diagName = 'ADVy'//diagSufx
> CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
> ENDIF
> #endif
> ENDIF
>
> C- Diffusive flux in Y
> IF (diffKh.NE.0.) THEN
> CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
> ELSE
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> df(i,j) = 0. _d 0
> ENDDO
> ENDDO
> ENDIF
>
> C- Add bi-harmonic flux in Y
> IF (diffK4 .NE. 0.) THEN
> CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
> ENDIF
>
> #ifdef ALLOW_GMREDI
> C- GM/Redi flux in Y
> IF (useGMRedi) THEN
> C *note* should update GMREDI_YTRANSPORT to set df *aja*
> CALL GMREDI_YTRANSPORT(
> I iMin,iMax,jMin,jMax,bi,bj,K,
> I yA,Tracer,tracerIdentity,
> U df,
> I myThid)
> ENDIF
> #endif
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fMer(i,j) = fMer(i,j) + df(i,j)
> ENDDO
> ENDDO
>
> #ifdef ALLOW_DIAGNOSTICS
> C- Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
> C excluding advective terms:
> IF ( useDiagnostics .AND.
> & (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
> diagName = 'DIFy'//diagSufx
> CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
> ENDIF
> #endif
>
> C-- Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
> C- Advective flux in R
> #ifdef ALLOW_AIM
> C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
> IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
> & (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
> & ) THEN
> #else
> IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
> #endif
> C- Compute vertical advective flux in the interior:
> IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
> CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
> ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
> & .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
> CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
> I dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
> O af, myThid )
> ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
> CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
> I dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
> O af, myThid )
> ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
> CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
> ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
> CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
> ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
> CALL GAD_DST3_ADV_R( bi,bj,k,
> I dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
> O af, myThid )
> ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
> cph This block is to trick the adjoint:
> cph IF inAdExact=.FALSE., we want to use DST3
> cph with limiters in forward, but without limiters in reverse.
> IF ( inAdMode ) THEN
> CALL GAD_DST3_ADV_R( bi,bj,k,
> I dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
> O af, myThid )
> ELSE
> CALL GAD_DST3FL_ADV_R( bi,bj,k,
> I dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
> O af, myThid )
> ENDIF
> ELSE
> STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
> ENDIF
> C- add the advective flux to fVerT
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
> ENDDO
> ENDDO
> #ifdef ALLOW_DIAGNOSTICS
> IF ( useDiagnostics ) THEN
> diagName = 'ADVr'//diagSufx
> CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
> C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
> C does it only if k=1 (never the case here)
> IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
> ENDIF
> #endif
> ENDIF
>
> C- Diffusive flux in R
> C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
> C boundary condition.
> IF (implicitDiffusion) THEN
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> df(i,j) = 0. _d 0
> ENDDO
> ENDDO
> ELSE
> CALL GAD_DIFF_R(bi,bj,k,KappaR,Tracer,df,myThid)
> ENDIF
>
> #ifdef ALLOW_GMREDI
> C- GM/Redi flux in R
> IF (useGMRedi) THEN
> C *note* should update GMREDI_RTRANSPORT to set df *aja*
> CALL GMREDI_RTRANSPORT(
> I iMin,iMax,jMin,jMax,bi,bj,K,
> I Tracer,tracerIdentity,
> U df,
> I myThid)
> ENDIF
> #endif
>
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
> ENDDO
> ENDDO
>
> #ifdef ALLOW_DIAGNOSTICS
> C- Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
> C Explicit terms only & excluding advective terms:
> IF ( useDiagnostics .AND.
> & (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
> diagName = 'DFrE'//diagSufx
> CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
> ENDIF
> #endif
>
> #ifdef ALLOW_KPP
> C- Set non local KPP transport term (ghat):
> IF ( useKPP .AND. k.GE.2 ) THEN
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> df(i,j) = 0. _d 0
> ENDDO
> ENDDO
> IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
> CALL KPP_TRANSPORT_T(
> I iMin,iMax,jMin,jMax,bi,bj,k,km1,
> O df )
> ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
> CALL KPP_TRANSPORT_S(
> I iMin,iMax,jMin,jMax,bi,bj,k,km1,
> O df )
> #ifdef ALLOW_PTRACERS
> ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
> CALL KPP_TRANSPORT_PTR(
> I iMin,iMax,jMin,jMax,bi,bj,k,km1,
> I tracerIdentity-GAD_TR1+1,
> O df )
> #endif
> ELSE
> PRINT*,'invalid tracer indentity: ', tracerIdentity
> STOP 'GAD_CALC_RHS: Ooops'
> ENDIF
> DO j=1-Oly,sNy+Oly
> DO i=1-Olx,sNx+Olx
> fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
> ENDDO
> ENDDO
> ENDIF
> #endif
>
> C-- Divergence of fluxes
> DO j=1-Oly,sNy+Oly-1
> DO i=1-Olx,sNx+Olx-1
> gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
> & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
> & *( (fZon(i+1,j)-fZon(i,j))
> & +(fMer(i,j+1)-fMer(i,j))
> & +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
> & -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
> & +(vTrans(i,j+1)-vTrans(i,j))
> & +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
> & )*advFac
> & )
> ENDDO
> ENDDO
>
> #ifdef ALLOW_DEBUG
> IF ( debugLevel .GE. debLevB
> & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
> & .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
> & .AND. nPx.EQ.1 .AND. nPy.EQ.1
> & .AND. useCubedSphereExchange ) THEN
> CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
> & fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
> ENDIF
> #endif /* ALLOW_DEBUG */
>
> RETURN
> END
>
>
> ------------------------------------------------------------------------
>
>
> subroutine adgad_calc_rhs( bi, bj, imin, imax, jmin, jmax, k, km1,
> $ kup, kdown, xa, ya, utrans, vtrans, rtrans, rtranskp1, maskup,
> $diffkh, diffk4, kappar, tracer, tracab, traceridentity,
> $advectionscheme, vertadvecscheme, calcadvection,
> $implicitadvection, mythid, adutrans, advtrans, adrtrans,
> $adrtranskp1, adkappar, adtracer, adtracab, adfvert, adgtracer )
>
> ...
>
> C----------------------------------------------
> C RESET LOCAL ADJOINT VARIABLES
> C----------------------------------------------
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> adaf(ip1,ip2) = 0.d0
> end do
> end do
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> addf(ip1,ip2) = 0.d0
> end do
> end do
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> addf4(ip1,ip2) = 0.d0
> end do
> end do
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> adfmer(ip1,ip2) = 0.d0
> end do
> end do
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> adfzon(ip1,ip2) = 0.d0
> end do
> end do
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> adlocabt(ip1,ip2) = 0.d0
> end do
> end do
> do ip2 = 1-oly, sny+oly
> do ip1 = 1-olx, snx+olx
> adlocalt(ip1,ip2) = 0.d0
> end do
> end do
>
> C----------------------------------------------
> C ROUTINE BODY
> C----------------------------------------------
> advfac = 0.d0
> if (calcadvection) then
> advfac = 1.d0
> endif
> radvfac = rksign*advfac
> if (implicitadvection) then
> radvfac = 0.d0
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> localt(i,j) = tracer(i,j,k,bi,bj)
> locabt(i,j) = tracab(i,j,k,bi,bj)
> end do
> end do
> do j = 1-oly, sny+oly-1
> do i = 1-olx, snx+olx-1
> adfmer(i,j+1) = adfmer(i,j+1)-adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)
> adfmer(i,j) = adfmer(i,j)+adgtracer(i,j,k,bi,bj)*recip_hfacc(
> $i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)
> adfvert(i,j,kdown) = adfvert(i,j,kdown)-adgtracer(i,j,k,bi,bj)
> $*recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*rksign
> adfvert(i,j,kup) = adfvert(i,j,kup)+adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*rksign
> adfzon(i+1,j) = adfzon(i+1,j)-adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)
> adfzon(i,j) = adfzon(i,j)+adgtracer(i,j,k,bi,bj)*recip_hfacc(
> $i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)
> adlocalt(i,j) = adlocalt(i,j)+adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*(utrans(
> $i+1,j)-utrans(i,j)+vtrans(i,j+1)-vtrans(i,j)+(rtranskp1(i,j)-
> $rtrans(i,j))*radvfac)*advfac
> adrtrans(i,j) = adrtrans(i,j)-adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*localt(
> $i,j)*radvfac*advfac
> adrtranskp1(i,j) = adrtranskp1(i,j)+adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*localt(
> $i,j)*radvfac*advfac
> adutrans(i+1,j) = adutrans(i+1,j)+adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*localt(
> $i,j)*advfac
> adutrans(i,j) = adutrans(i,j)-adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*localt(
> $i,j)*advfac
> advtrans(i,j+1) = advtrans(i,j+1)+adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*localt(
> $i,j)*advfac
> advtrans(i,j) = advtrans(i,j)-adgtracer(i,j,k,bi,bj)*
> $recip_hfacc(i,j,k,bi,bj)*recip_drf(k)*recip_ra(i,j,bi,bj)*localt(
> $i,j)*advfac
> end do
> end do
> if (usekpp .and. k .ge. 2) then
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = addf(i,j)+adfvert(i,j,kup)*maskup(i,j)
> end do
> end do
> if (traceridentity .eq. gad_temperature) then
> call adkpp_transport_t( imin,imax,jmin,jmax,bi,bj,k,km1,addf )
> else if (traceridentity .eq. gad_salinity) then
> call adkpp_transport_s( imin,imax,jmin,jmax,bi,bj,k,km1,addf )
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = 0.d0
> end do
> end do
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = addf(i,j)+adfvert(i,j,kup)*maskup(i,j)
> end do
> end do
> if (usegmredi) then
> call adgmredi_rtransport( imin,imax,jmin,jmax,bi,bj,k,tracer,
> $adtracer,addf )
> endif
> if (implicitdiffusion) then
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = 0.d0
> end do
> end do
> else
> call adgad_diff_r( bi,bj,k,kappar,tracer,adkappar,adtracer,addf
> $)
> endif
> if (calcadvection .and. ( .not. implicitadvection) .and. k .ge. 2)
> $ then
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adaf(i,j) = adaf(i,j)+adfvert(i,j,kup)
> end do
> end do
> if (vertadvecscheme .eq. enum_centered_2nd) then
> call adgad_c2_adv_r( bi,bj,k,rtrans,tracab,adrtrans,adtracab,
> $adaf )
> else if (vertadvecscheme .eq. enum_upwind_1rst .or.
> $vertadvecscheme .eq. enum_dst2) then
> call adgad_dst2u1_adv_r( bi,bj,k,vertadvecscheme,dttracerlev(
> $k),rtrans,tracab(1-olx,1-oly,1,bi,bj),adrtrans,adtracab(1-olx,1-
> $oly,1,bi,bj),adaf )
> else if (vertadvecscheme .eq. enum_flux_limit) then
> call adgad_fluxlimit_adv_r( bi,bj,k,dttracerlev(k),rtrans,
> $tracab(1-olx,1-oly,1,bi,bj),adrtrans,adtracab(1-olx,1-oly,1,bi,bj)
> $,adaf )
> else if (vertadvecscheme .eq. enum_upwind_3rd) then
> call adgad_u3_adv_r( bi,bj,k,rtrans,tracab,adrtrans,adtracab,
> $adaf )
> else if (vertadvecscheme .eq. enum_centered_4th) then
> call adgad_c4_adv_r( bi,bj,k,rtrans,tracab,adrtrans,adtracab,
> $adaf )
> else if (vertadvecscheme .eq. enum_dst3) then
> call adgad_dst3_adv_r( bi,bj,k,dttracerlev(k),rtrans,tracab(1-
> $olx,1-oly,1,bi,bj),adrtrans,adtracab(1-olx,1-oly,1,bi,bj),adaf )
> else if (vertadvecscheme .eq. enum_dst3_flux_limit) then
> if (inadmode) then
> call adgad_dst3_adv_r( bi,bj,k,dttracerlev(k),rtrans,tracab(
> $1-olx,1-oly,1,bi,bj),adrtrans,adtracab(1-olx,1-oly,1,bi,bj),adaf )
> else
> call adgad_dst3fl_adv_r( bi,bj,k,dttracerlev(k),rtrans,
> $tracab(1-olx,1-oly,1,bi,bj),adrtrans,adtracab(1-olx,1-oly,1,bi,bj)
> $,adaf )
> endif
> endif
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = addf(i,j)+adfmer(i,j)
> end do
> end do
> if (usegmredi) then
> call adgmredi_ytransport( imin,imax,jmin,jmax,bi,bj,k,ya,tracer,
> $traceridentity,mythid,adtracer,addf )
> endif
> if (diffk4 .ne. 0.) then
> call adgad_biharm_y( bi,bj,ya,diffk4,addf4,addf )
> endif
> if (diffkh .ne. 0.) then
> call adgad_diff_y( bi,bj,ya,diffkh,adlocalt,addf )
> else
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = 0.d0
> end do
> end do
> endif
> if (calcadvection) then
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adaf(i,j) = adaf(i,j)+adfmer(i,j)
> end do
> end do
> if (advectionscheme .eq. enum_centered_2nd) then
> call adgad_c2_adv_y( vtrans,locabt,advtrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_upwind_1rst .or.
> $advectionscheme .eq. enum_dst2) then
> call adgad_dst2u1_adv_y( bi,bj,k,advectionscheme,dttracerlev(
> $k),vtrans,locabt,advtrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_flux_limit) then
> call adgad_fluxlimit_adv_y( bi,bj,k,dttracerlev(k),vtrans,
> $masks(1-olx,1-oly,k,bi,bj),locabt,advtrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_upwind_3rd) then
> call adgad_u3_adv_y( bi,bj,k,vtrans,locabt,advtrans,adlocabt,
> $adaf )
> else if (advectionscheme .eq. enum_centered_4th) then
> call adgad_c4_adv_y( bi,bj,k,vtrans,locabt,advtrans,adlocabt,
> $adaf )
> else if (advectionscheme .eq. enum_dst3) then
> call adgad_dst3_adv_y( bi,bj,k,dttracerlev(k),vtrans,masks(1-
> $olx,1-oly,k,bi,bj),locabt,advtrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_dst3_flux_limit) then
> if (inadmode) then
> call adgad_dst3_adv_y( bi,bj,k,dttracerlev(k),vtrans,masks(
> $1-olx,1-oly,k,bi,bj),locabt,advtrans,adlocabt,adaf )
> else
> call adgad_dst3fl_adv_y( bi,bj,k,dttracerlev(k),vtrans,
> $masks(1-olx,1-oly,k,bi,bj),locabt,advtrans,adlocabt,adaf )
> endif
> endif
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adfmer(i,j) = 0.d0
> end do
> end do
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> addf(i,j) = addf(i,j)+adfzon(i,j)
> end do
> end do
> if (usegmredi) then
> call adgmredi_xtransport( imin,imax,jmin,jmax,bi,bj,k,xa,tracer,
> $traceridentity,mythid,adtracer,addf )
> endif
> if (diffk4 .ne. 0.) then
> call adgad_biharm_x( bi,bj,xa,diffk4,addf4,addf )
> endif
> if (diffkh .ne. 0.) then
> call adgad_diff_x( bi,bj,xa,diffkh,adlocalt,addf )
> endif
> if (calcadvection) then
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adaf(i,j) = adaf(i,j)+adfzon(i,j)
> end do
> end do
> if (advectionscheme .eq. enum_centered_2nd) then
> call adgad_c2_adv_x( utrans,locabt,adutrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_upwind_1rst .or.
> $advectionscheme .eq. enum_dst2) then
> call adgad_dst2u1_adv_x( bi,bj,k,advectionscheme,dttracerlev(
> $k),utrans,locabt,adutrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_flux_limit) then
> call adgad_fluxlimit_adv_x( bi,bj,k,dttracerlev(k),utrans,
> $maskw(1-olx,1-oly,k,bi,bj),locabt,adutrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_upwind_3rd) then
> call adgad_u3_adv_x( bi,bj,k,utrans,locabt,adutrans,adlocabt,
> $adaf )
> else if (advectionscheme .eq. enum_centered_4th) then
> call adgad_c4_adv_x( bi,bj,k,utrans,locabt,adutrans,adlocabt,
> $adaf )
> else if (advectionscheme .eq. enum_dst3) then
> call adgad_dst3_adv_x( bi,bj,k,dttracerlev(k),utrans,maskw(1-
> $olx,1-oly,k,bi,bj),locabt,adutrans,adlocabt,adaf )
> else if (advectionscheme .eq. enum_dst3_flux_limit) then
> if (inadmode) then
> call adgad_dst3_adv_x( bi,bj,k,dttracerlev(k),utrans,maskw(
> $1-olx,1-oly,k,bi,bj),locabt,adutrans,adlocabt,adaf )
> else
> call adgad_dst3fl_adv_x( bi,bj,k,dttracerlev(k),utrans,
> $maskw(1-olx,1-oly,k,bi,bj),locabt,adutrans,adlocabt,adaf )
> endif
> endif
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adfzon(i,j) = 0.d0
> end do
> end do
> if (diffk4 .ne. 0.) then
> call adgad_del2( bi,bj,k,adfzon,adfmer,addf4 )
> call adgad_grad_y( bi,bj,ya,adlocalt,adfmer )
> call adgad_grad_x( bi,bj,xa,adlocalt,adfzon )
> endif
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adtracab(i,j,k,bi,bj) = adtracab(i,j,k,bi,bj)+adlocabt(i,j)
> adlocabt(i,j) = 0.d0
> adtracer(i,j,k,bi,bj) = adtracer(i,j,k,bi,bj)+adlocalt(i,j)
> adlocalt(i,j) = 0.d0
> end do
> end do
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adfvert(i,j,kup) = 0.d0
> end do
> end do
>
> end
>
>
>
> ------------------------------------------------------------------------
>
>
> subroutine adcalc_gt( bi, bj, imin, imax, jmin, jmax, k, km1, kup,
> $ kdown, xa, ya, utrans, vtrans, rtrans, rtranskp1, maskup,
> $kappart, myiter, mythid, adutrans, advtrans, adrtrans,
> $adrtranskp1, adkappart, adfvert )
>
> ...
>
>
> C----------------------------------------------
> C ROUTINE BODY
> C----------------------------------------------
> calcadvection = tempadvection .and. ( .not. tempmultidimadvec)
> iternb = myiter
> if (staggertimestep) then
> iternb = myiter-1
> endif
> if (tempforcing .and. ( .not. forcing_in_ab)) then
> call adexternal_forcing_t( bi,bj,k )
> endif
> if (adamsbashforthgt) then
> call adadams_bashforth2( bi,bj,k,iternb,adgt,adgtnm1 )
> endif
> if (tempforcing .and. forcing_in_ab) then
> call adexternal_forcing_t( bi,bj,k )
> endif
> call adgad_calc_rhs( bi,bj,imin,imax,jmin,jmax,k,km1,kup,kdown,xa,
> $ya,utrans,vtrans,rtrans,rtranskp1,maskup,diffkht,diffk4t,kappart,
> $theta,theta,gad_temperature,tempadvscheme,tempvertadvscheme,
> $calcadvection,tempimplvertadv,mythid,adutrans,advtrans,adrtrans,
> $adrtranskp1,adkappart,adtheta,adtheta,adfvert,adgt )
>
> end
>
>
> ------------------------------------------------------------------------
>
>
> subroutine adcalc_gt( bi, bj, imin, imax, jmin, jmax, k, km1, kup,
> $ kdown, xa, ya, utrans, vtrans, rtrans, rtranskp1, maskup,
> $kappart, myiter, mythid, adutrans, advtrans, adrtrans,
> $adrtranskp1, adkappart, adfvert )
>
> ...
>
> C----------------------------------------------
> C ROUTINE BODY
> C----------------------------------------------
> cph(
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adthetah(i,j,k,bi,bj) = 0.D0
> adthetah(i,j,km1,bi,bj) = 0.D0
> end do
> end do
> cph)
> calcadvection = tempadvection .and. ( .not. tempmultidimadvec)
> iternb = myiter
> if (staggertimestep) then
> iternb = myiter-1
> endif
> if (tempforcing .and. ( .not. forcing_in_ab)) then
> call adexternal_forcing_t( bi,bj,k )
> endif
> if (adamsbashforthgt) then
> call adadams_bashforth2( bi,bj,k,iternb,adgt,adgtnm1 )
> endif
> if (tempforcing .and. forcing_in_ab) then
> call adexternal_forcing_t( bi,bj,k )
> endif
> call adgad_calc_rhs( bi,bj,imin,imax,jmin,jmax,k,km1,kup,kdown,xa,
> $ya,utrans,vtrans,rtrans,rtranskp1,maskup,diffkht,diffk4t,kappart,
> $theta,theta,gad_temperature,tempadvscheme,tempvertadvscheme,
> $calcadvection,tempimplvertadv,mythid,adutrans,advtrans,adrtrans,
> $adrtranskp1,adkappart,adtheta,adthetah,adfvert,adgt )
> cph(
> do j = 1-oly, sny+oly
> do i = 1-olx, snx+olx
> adtheta(i,j,k,bi,bj) = adtheta(i,j,k,bi,bj) +
> & adthetah(i,j,k,bi,bj)
> adthetah(i,j,k,bi,bj) = 0.D0
> adtheta(i,j,km1,bi,bj) = adtheta(i,j,km1,bi,bj) +
> & adthetah(i,j,km1,bi,bj)
> adthetah(i,j,km1,bi,bj) = 0.D0
> end do
> end do
> cph)
>
> end
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list