[MITgcm-devel] Fwd: TAF issue for "special" parameter lists
Jean-Michel Campin
jmc at ocean.mit.edu
Sat Feb 25 13:42:47 EST 2006
Hi Chris,
You will need to tell me more about this "no aliasing" thing,
because I don't quiet understand the unsafe part of it
for this particular case (where both arguments are input-only
and therefore not modified in the S/R).
But this definitively help in deciding how to modify the code
to fix the adjoint.
I am looking at this right now.
See you,
Jean-Michel
On Fri, Feb 24, 2006 at 04:03:33PM -1000, Chris Hill wrote:
> 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
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list