[MITgcm-devel] Fwd: TAF issue for "special" parameter lists
Patrick Heimbach
heimbach at MIT.EDU
Mon Feb 27 14:54:23 EST 2006
Hi Jean-Michel,
that solution works (because it avoids variable
aliasing in the parameter list).
So that's all set.
Thanks
-Patrick
On Sat, 2006-02-25 at 17:23, Jean-Michel Campin wrote:
> Hi Patrick,
>
> I propose to keep tracer & tracAB in the argument list of
> GAD_CALC_RHS, and inside GAD_CALC_RHS to use only 1 of the 2
> when GAD_CALC_RHS is called with the same argument but twice:
> CALL GAD_CALC_RHS( ... , theta, theta , ...).
>
> few questions:
> - any chance that this will work ?
> - if we keep the call like it is now (e.g., in CALC_GT),
> which argument (the 1rst one or the 2nd one ?) do I need
> to use inside GAD_CALC_RHS to get the correct AD version ?
> (may be this is stupid ?)
> - or do we need to go further and call GAD_CALC_RHS with
> ( ... , theta , dummyT, ... ) where dummyT is just unused ?
> or use gtNm1 instead of dummyT ? (I kind of like this one,
> since it's also the call we need when using AB_T ).
>
> I will check-in a modified version of gad_calc_rhs + calc_gt & calc_gs
> (and leave the ptracers for later), for the moment with plenty
> of #ifdef ALLOW_ADAMSBASHFORTH_3, and when you have time, you can try it,
> and we will see.
>
> If none of those works, we will need to have 2 versions of
> GAD_CALC_RHS (a pain !) or one version with 1 more argument when
> ALLOW_ADAMSBASHFORTH_3 is defined than when it's not (dirty !).
>
> Jean-Michel
>
> On Sat, Feb 25, 2006 at 01:44:53AM -0500, Patrick Heimbach wrote:
> >
> > Well, I agree.
> > -p.
> >
> >
> >
> > Quoting Chris Hill <cnh at MIT.EDU>:
> >
> > >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
> > >
> >
> >
> >
> > --------------------------------------------------------
> > Patrick Heimbach Massachusetts Institute of Technology
> > FON: +1/617/253-5259 EAPS, Room 54-1518
> > FAX: +1/617/253-4464 77 Massachusetts Avenue
> > mailto:heimbach at mit.edu Cambridge MA 02139
> > http://www.mit.edu/~heimbach/ USA
> >
> > _______________________________________________
> > 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
--
--------------------------------------------------------
Patrick Heimbach Massachusetts Institute of Technology
FON: +1/617/253-5259 EAPS, Room 54-1518
FAX: +1/617/253-4464 77 Massachusetts Avenue
mailto:heimbach at mit.edu Cambridge MA 02139
http://www.mit.edu/~heimbach/ USA
More information about the MITgcm-devel
mailing list