--- MITgcm/model/src/thermodynamics.F 2001/08/13 18:05:26 1.2 +++ MITgcm/model/src/thermodynamics.F 2003/10/24 05:29:35 1.54 @@ -1,31 +1,101 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/thermodynamics.F,v 1.2 2001/08/13 18:05:26 heimbach Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/thermodynamics.F,v 1.54 2003/10/24 05:29:35 edhill Exp $ C $Name: $ +#include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef ALLOW_GMREDI +# include "GMREDI_OPTIONS.h" +# endif +# ifdef ALLOW_KPP +# include "KPP_OPTIONS.h" +# endif +#ifdef ALLOW_PTRACERS +# include "PTRACERS_OPTIONS.h" +#endif +#endif /* ALLOW_AUTODIFF_TAMC */ + +CBOP +C !ROUTINE: THERMODYNAMICS +C !INTERFACE: SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid) -C /==========================================================\ -C | SUBROUTINE THERMODYNAMICS | -C | o Controlling routine for the prognostic part of the | -C | thermo-dynamics. | -C |==========================================================| -C \==========================================================/ - IMPLICIT NONE +C !DESCRIPTION: \bv +C *==========================================================* +C | SUBROUTINE THERMODYNAMICS +C | o Controlling routine for the prognostic part of the +C | thermo-dynamics. +C *=========================================================== +C | The algorithm... +C | +C | "Correction Step" +C | ================= +C | Here we update the horizontal velocities with the surface +C | pressure such that the resulting flow is either consistent +C | with the free-surface evolution or the rigid-lid: +C | U[n] = U* + dt x d/dx P +C | V[n] = V* + dt x d/dy P +C | +C | "Calculation of Gs" +C | =================== +C | This is where all the accelerations and tendencies (ie. +C | physics, parameterizations etc...) are calculated +C | rho = rho ( theta[n], salt[n] ) +C | b = b(rho, theta) +C | K31 = K31 ( rho ) +C | Gu[n] = Gu( u[n], v[n], wVel, b, ... ) +C | Gv[n] = Gv( u[n], v[n], wVel, b, ... ) +C | Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... ) +C | Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... ) +C | +C | "Time-stepping" or "Prediction" +C | ================================ +C | The models variables are stepped forward with the appropriate +C | time-stepping scheme (currently we use Adams-Bashforth II) +C | - For momentum, the result is always *only* a "prediction" +C | in that the flow may be divergent and will be "corrected" +C | later with a surface pressure gradient. +C | - Normally for tracers the result is the new field at time +C | level [n+1} *BUT* in the case of implicit diffusion the result +C | is also *only* a prediction. +C | - We denote "predictors" with an asterisk (*). +C | U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] ) +C | V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] ) +C | theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) +C | salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) +C | With implicit diffusion: +C | theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) +C | salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) +C | (1 + dt * K * d_zz) theta[n] = theta* +C | (1 + dt * K * d_zz) salt[n] = salt* +C | +C *==========================================================* +C \ev +C !USES: + IMPLICIT NONE C == Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" +#include "GAD.h" #ifdef ALLOW_PASSIVE_TRACER #include "TR1.h" #endif +#ifdef ALLOW_PTRACERS +#include "PTRACERS.h" +#endif +#ifdef ALLOW_TIMEAVE +#include "TIMEAVE_STATV.h" +#endif #ifdef ALLOW_AUTODIFF_TAMC # include "tamc.h" # include "tamc_keys.h" # include "FFIELDS.h" +# include "EOS.h" # ifdef ALLOW_KPP # include "KPP.h" # endif @@ -34,10 +104,7 @@ # endif #endif /* ALLOW_AUTODIFF_TAMC */ -#ifdef ALLOW_TIMEAVE -#include "TIMEAVE_STATV.h" -#endif - +C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C myTime - Current time in simulation C myIter - Current iteration number in simulation @@ -46,6 +113,7 @@ INTEGER myIter INTEGER myThid +C !LOCAL VARIABLES: C == Local variables C xA, yA - Per block temporaries holding face areas C uTrans, vTrans, rTrans - Per block temporaries holding flow @@ -59,22 +127,17 @@ C so we need an fVer for each C variable. C rhoK, rhoKM1 - Density at current level, and level above -C phiHyd - Hydrostatic part of the potential phiHydi. -C In z coords phiHydiHyd is the hydrostatic -C Potential (=pressure/rho0) anomaly -C In p coords phiHydiHyd is the geopotential -C surface height anomaly. C phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean) C phiSurfY or geopotentiel (atmos) in X and Y direction C KappaRT, - Total diffusion in vertical for T and S. C KappaRS (background + spatially varying, isopycnal term). +C useVariableK = T when vertical diffusion is not constant C iMin, iMax - Ranges and sub-block indices on which calculations C jMin, jMax are applied. C bi, bj C k, kup, - Index for layer above and below. kup and kDown C kDown, km1 are switched with layer to be the appropriate C index into fVerTerm. -C tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf. _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) @@ -84,7 +147,6 @@ _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) - _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -94,97 +156,29 @@ _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL tauAB - -C This is currently used by IVDC and Diagnostics +C This is currently used by IVDC and Diagnostics _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - + LOGICAL useVariableK INTEGER iMin, iMax INTEGER jMin, jMax INTEGER bi, bj INTEGER i, j INTEGER k, km1, kup, kDown + INTEGER iTracer -Cjmc : add for phiHyd output <- but not working if multi tile per CPU -c CHARACTER*(MAX_LEN_MBUF) suff -c LOGICAL DIFFERENT_MULTIPLE -c EXTERNAL DIFFERENT_MULTIPLE -Cjmc(end) - -C--- The algorithm... -C -C "Correction Step" -C ================= -C Here we update the horizontal velocities with the surface -C pressure such that the resulting flow is either consistent -C with the free-surface evolution or the rigid-lid: -C U[n] = U* + dt x d/dx P -C V[n] = V* + dt x d/dy P -C -C "Calculation of Gs" -C =================== -C This is where all the accelerations and tendencies (ie. -C physics, parameterizations etc...) are calculated -C rho = rho ( theta[n], salt[n] ) -C b = b(rho, theta) -C K31 = K31 ( rho ) -C Gu[n] = Gu( u[n], v[n], wVel, b, ... ) -C Gv[n] = Gv( u[n], v[n], wVel, b, ... ) -C Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... ) -C Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... ) -C -C "Time-stepping" or "Prediction" -C ================================ -C The models variables are stepped forward with the appropriate -C time-stepping scheme (currently we use Adams-Bashforth II) -C - For momentum, the result is always *only* a "prediction" -C in that the flow may be divergent and will be "corrected" -C later with a surface pressure gradient. -C - Normally for tracers the result is the new field at time -C level [n+1} *BUT* in the case of implicit diffusion the result -C is also *only* a prediction. -C - We denote "predictors" with an asterisk (*). -C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] ) -C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] ) -C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C With implicit diffusion: -C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] ) -C (1 + dt * K * d_zz) theta[n] = theta* -C (1 + dt * K * d_zz) salt[n] = salt* -C--- +CEOP +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_ENTER('FORWARD_STEP',myThid) +#endif + #ifdef ALLOW_AUTODIFF_TAMC C-- dummy statement to end declaration part ikey = 1 + itdkey = 1 #endif /* ALLOW_AUTODIFF_TAMC */ -C-- Set up work arrays with valid (i.e. not NaN) values -C These inital values do not alter the numerical results. They -C just ensure that all memory references are to valid floating -C point numbers. This prevents spurious hardware signals due to -C uninitialised but inert locations. - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - xA(i,j) = 0. _d 0 - yA(i,j) = 0. _d 0 - uTrans(i,j) = 0. _d 0 - vTrans(i,j) = 0. _d 0 - DO k=1,Nr - phiHyd(i,j,k) = 0. _d 0 - sigmaX(i,j,k) = 0. _d 0 - sigmaY(i,j,k) = 0. _d 0 - sigmaR(i,j,k) = 0. _d 0 - ENDDO - rhoKM1 (i,j) = 0. _d 0 - rhok (i,j) = 0. _d 0 - phiSurfX(i,j) = 0. _d 0 - phiSurfY(i,j) = 0. _d 0 - ENDDO - ENDDO - - #ifdef ALLOW_AUTODIFF_TAMC C-- HPF directive to help TAMC CHPF$ INDEPENDENT @@ -195,7 +189,7 @@ #ifdef ALLOW_AUTODIFF_TAMC C-- HPF directive to help TAMC CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS -CHPF$& ,phiHyd,utrans,vtrans,xA,yA +CHPF$& ,utrans,vtrans,xA,yA CHPF$& ,KappaRT,KappaRS CHPF$& ) #endif /* ALLOW_AUTODIFF_TAMC */ @@ -205,23 +199,32 @@ #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 - - ikey = (act1 + 1) + act2*max1 + itdkey = (act1 + 1) + act2*max1 & + act3*max1*max2 & + act4*max1*max2*max3 #endif /* ALLOW_AUTODIFF_TAMC */ -C-- Set up work arrays that need valid initial values +C-- Set up work arrays with valid (i.e. not NaN) values +C These inital values do not alter the numerical results. They +C just ensure that all memory references are to valid floating +C point numbers. This prevents spurious hardware signals due to +C uninitialised but inert locations. + DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx + xA(i,j) = 0. _d 0 + yA(i,j) = 0. _d 0 + uTrans(i,j) = 0. _d 0 + vTrans(i,j) = 0. _d 0 + rhok (i,j) = 0. _d 0 + rhoKM1 (i,j) = 0. _d 0 + phiSurfX(i,j) = 0. _d 0 + phiSurfY(i,j) = 0. _d 0 rTrans (i,j) = 0. _d 0 fVerT (i,j,1) = 0. _d 0 fVerT (i,j,2) = 0. _d 0 @@ -236,24 +239,74 @@ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C This is currently also used by IVDC and Diagnostics + sigmaX(i,j,k) = 0. _d 0 + sigmaY(i,j,k) = 0. _d 0 + sigmaR(i,j,k) = 0. _d 0 ConvectCount(i,j,k) = 0. - KappaRT(i,j,k) = 0. _d 0 - KappaRS(i,j,k) = 0. _d 0 + KappaRT(i,j,k) = 0. _d 0 + KappaRS(i,j,k) = 0. _d 0 +C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs): + gT(i,j,k,bi,bj) = 0. _d 0 + gS(i,j,k,bi,bj) = 0. _d 0 +# ifdef ALLOW_PASSIVE_TRACER +ceh3 needs an IF ( use PASSIVE_TRACER) THEN + gTr1(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef ALLOW_PTRACERS +ceh3 this should have an IF ( usePTRACERS ) THEN + DO iTracer=1,PTRACERS_numInUse + gPTr(i,j,k,bi,bj,itracer) = 0. _d 0 + ENDDO +# endif +#ifdef ALLOW_AUTODIFF_TAMC +cph all the following init. are necessary for TAF +cph although some of these are re-initialised later. +# ifdef ALLOW_GMREDI + Kwx(i,j,k,bi,bj) = 0. _d 0 + Kwy(i,j,k,bi,bj) = 0. _d 0 + Kwz(i,j,k,bi,bj) = 0. _d 0 +# ifdef GM_NON_UNITY_DIAGONAL + Kux(i,j,k,bi,bj) = 0. _d 0 + Kvy(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef GM_EXTRA_DIAGONAL + Kuz(i,j,k,bi,bj) = 0. _d 0 + Kvz(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef GM_BOLUS_ADVEC + GM_PsiX(i,j,k,bi,bj) = 0. _d 0 + GM_PsiY(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef GM_VISBECK_VARIABLE_K + VisbeckK(i,j,bi,bj) = 0. _d 0 +# endif +# endif /* ALLOW_GMREDI */ +#endif /* ALLOW_AUTODIFF_TAMC */ ENDDO ENDDO ENDDO - iMin = 1-OLx+1 + iMin = 1-OLx iMax = sNx+OLx - jMin = 1-OLy+1 + jMin = 1-OLy jMax = sNy+OLy - #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte -CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte +CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE totphihyd +CADJ & = comlev1_bibj, key=itdkey, byte=isbyte +#ifdef ALLOW_KPP +CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +#endif #endif /* ALLOW_AUTODIFF_TAMC */ +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid) +#endif + C-- Start of diagnostic loop DO k=Nr,1,-1 @@ -262,50 +315,61 @@ C? Do we still need this? cph kkey formula corrected. cph Needed for rhok, rhokm1, in the case useGMREDI. - kkey = (ikey-1)*Nr + k -CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte -CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte + kkey = (itdkey-1)*Nr + k #endif /* ALLOW_AUTODIFF_TAMC */ C-- Integrate continuity vertically for vertical velocity - CALL INTEGRATE_FOR_W( - I bi, bj, k, uVel, vVel, - O wVel, - I myThid ) +c CALL INTEGRATE_FOR_W( +c I bi, bj, k, uVel, vVel, +c O wVel, +c I myThid ) #ifdef ALLOW_OBCS #ifdef ALLOW_NONHYDROSTATIC C-- Apply OBC to W if in N-H mode - IF (useOBCS.AND.nonHydrostatic) THEN - CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid ) - ENDIF +c IF (useOBCS.AND.nonHydrostatic) THEN +c CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid ) +c ENDIF #endif /* ALLOW_NONHYDROSTATIC */ #endif /* ALLOW_OBCS */ +C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h +C-- MOST of THERMODYNAMICS will be disabled +#ifndef SINGLE_LAYER_MODE + C-- Calculate gradients of potential density for isoneutral C slope terms (e.g. GM/Redi tensor or IVDC diffusivity) c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('FIND_RHO',myThid) +#endif #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL FIND_RHO( - I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType, + I bi, bj, iMin, iMax, jMin, jMax, k, k, I theta, salt, O rhoK, I myThid ) + IF (k.GT.1) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL FIND_RHO( - I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType, + I bi, bj, iMin, iMax, jMin, jMax, k-1, k, I theta, salt, O rhoKm1, I myThid ) ENDIF +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('GRAD_SIGMA',myThid) +#endif CALL GRAD_SIGMA( I bi, bj, iMin, iMax, jMin, jMax, k, I rhoK, rhoKm1, rhoK, @@ -313,9 +377,17 @@ I myThid ) ENDIF +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte +CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ C-- Implicit Vertical Diffusion for Convection c ==> should use sigmaR !!! IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('CALC_IVDC',myThid) +#endif CALL CALC_IVDC( I bi, bj, iMin, iMax, jMin, jMax, k, I rhoKm1, rhoK, @@ -323,70 +395,109 @@ I myTime, myIter, myThid) ENDIF +#endif /* SINGLE_LAYER_MODE */ + C-- end of diagnostic k loop (Nr:1) ENDDO #ifdef ALLOW_AUTODIFF_TAMC cph avoids recomputation of integrate_for_w -CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte +CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_OBCS C-- Calculate future values on open boundaries IF (useOBCS) THEN - CALL OBCS_CALC( bi, bj, myTime+deltaT, +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('OBCS_CALC',myThid) +#endif + CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1, I uVel, vVel, wVel, theta, salt, I myThid ) ENDIF #endif /* ALLOW_OBCS */ + +#ifdef ALLOW_THERM_SEAICE + IF (useThermSeaIce) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('ICE_FORCING',myThid) +#endif C-- Determines forcing terms based on external fields -C relaxation terms, etc. - CALL EXTERNAL_FORCING_SURF( +C including effects from ice + CALL ICE_FORCING( I bi, bj, iMin, iMax, jMin, jMax, I myThid ) + ELSE +#endif /* ALLOW_THERM_SEAICE */ + +C-- Determines forcing terms based on external fields +C relaxation terms, etc. +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid) +#endif + CALL EXTERNAL_FORCING_SURF( + I bi, bj, iMin, iMax, jMin, jMax, + I myTime, myIter, myThid ) + +#ifdef ALLOW_THERM_SEAICE +C-- end of if/else block useThermSeaIce -- + ENDIF +#endif /* ALLOW_THERM_SEAICE */ + #ifdef ALLOW_AUTODIFF_TAMC cph needed for KPP CADJ STORE surfacetendencyU(:,:,bi,bj) -CADJ & = comlev1_bibj, key=ikey, byte=isbyte +CADJ & = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE surfacetendencyV(:,:,bi,bj) -CADJ & = comlev1_bibj, key=ikey, byte=isbyte +CADJ & = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE surfacetendencyS(:,:,bi,bj) -CADJ & = comlev1_bibj, key=ikey, byte=isbyte +CADJ & = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE surfacetendencyT(:,:,bi,bj) -CADJ & = comlev1_bibj, key=ikey, byte=isbyte +CADJ & = comlev1_bibj, key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ +C-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h +C-- MOST of THERMODYNAMICS will be disabled +#ifndef SINGLE_LAYER_MODE + #ifdef ALLOW_GMREDI #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte -CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte -CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte +cph storing here is needed only for one GMREDI_OPTIONS: +cph define GM_BOLUS_ADVEC +cph but I've avoided the #ifdef for now, in case more things change +CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ + C-- Calculate iso-neutral slopes for the GM/Redi parameterisation IF (useGMRedi) THEN - DO k=1,Nr - CALL GMREDI_CALC_TENSOR( - I bi, bj, iMin, iMax, jMin, jMax, k, +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid) +#endif + CALL GMREDI_CALC_TENSOR( + I bi, bj, iMin, iMax, jMin, jMax, I sigmaX, sigmaY, sigmaR, I myThid ) - ENDDO #ifdef ALLOW_AUTODIFF_TAMC ELSE - DO k=1, Nr - CALL GMREDI_CALC_TENSOR_DUMMY( - I bi, bj, iMin, iMax, jMin, jMax, k, + CALL GMREDI_CALC_TENSOR_DUMMY( + I bi, bj, iMin, iMax, jMin, jMax, I sigmaX, sigmaY, sigmaR, I myThid ) - ENDDO #endif /* ALLOW_AUTODIFF_TAMC */ ENDIF #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte -CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte -CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte +CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ #endif /* ALLOW_GMREDI */ @@ -394,6 +505,10 @@ #ifdef ALLOW_KPP C-- Compute KPP mixing coefficients IF (useKPP) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('KPP_CALC',myThid) +#endif CALL KPP_CALC( I bi, bj, myTime, myThid ) #ifdef ALLOW_AUTODIFF_TAMC @@ -405,37 +520,89 @@ #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE KPPghat (:,:,:,bi,bj) -CADJ & , KPPviscAz (:,:,:,bi,bj) CADJ & , KPPdiffKzT(:,:,:,bi,bj) CADJ & , KPPdiffKzS(:,:,:,bi,bj) CADJ & , KPPfrac (:,: ,bi,bj) -CADJ & = comlev1_bibj, key=ikey, byte=isbyte +CADJ & = comlev1_bibj, key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ #endif /* ALLOW_KPP */ #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key = ikey, byte = isbyte -CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key = ikey, byte = isbyte -CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte -CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte -CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte -CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte +CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte #ifdef ALLOW_PASSIVE_TRACER -CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte +CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte +#endif +#ifdef ALLOW_PTRACERS +cph-- moved to forward_step to avoid key computation +cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj, +cphCADJ & key=itdkey, byte=isbyte #endif #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_AIM C AIM - atmospheric intermediate model, physics package code. -C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics IF ( useAIM ) THEN - CALL TIMER_START('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid) - CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid ) - CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS [DYNAMICS]', myThid) +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid) +#endif + CALL TIMER_START('AIM_DO_PHYSICS [THERMODYNAMICS]', myThid) + CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid ) + CALL TIMER_STOP( 'AIM_DO_PHYSICS [THERMODYNAMICS]', myThid) ENDIF #endif /* ALLOW_AIM */ +#ifndef DISABLE_MULTIDIM_ADVECTION +C-- Some advection schemes are better calculated using a multi-dimensional +C method in the absence of any other terms and, if used, is done here. +C +C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h +C The default is to use multi-dimensinal advection for non-linear advection +C schemes. However, for the sake of efficiency of the adjoint it is necessary +C to be able to exclude this scheme to avoid excessive storage and +C recomputation. It *is* differentiable, if you need it. +C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to +C disable this section of code. + IF (tempMultiDimAdvec) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('GAD_ADVECTION',myThid) +#endif + CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE, + U theta,gT, + I myTime,myIter,myThid) + ENDIF + IF (saltMultiDimAdvec) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('GAD_ADVECTION',myThid) +#endif + CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY, + U salt,gS, + I myTime,myIter,myThid) + ENDIF +C Since passive tracers are configurable separately from T,S we +C call the multi-dimensional method for PTRACERS regardless +C of whether multiDimAdvection is set or not. +#ifdef ALLOW_PTRACERS + IF ( usePTRACERS ) THEN +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid) +#endif + CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid ) + ENDIF +#endif /* ALLOW_PTRACERS */ +#endif /* DISABLE_MULTIDIM_ADVECTION */ + +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid) +#endif C-- Start of thermodynamics loop DO k=Nr,1,-1 @@ -443,7 +610,7 @@ C? Patrick Is this formula correct? cph Yes, but I rewrote it. cph Also, the KappaR? need the index and subscript k! - kkey = (ikey-1)*Nr + k + kkey = (itdkey-1)*Nr + k #endif /* ALLOW_AUTODIFF_TAMC */ C-- km1 Points to level above k (=k-1) @@ -465,11 +632,26 @@ O xA,yA,uTrans,vTrans,rTrans,maskUp, I myThid) +#ifdef ALLOW_GMREDI + +C-- Residual transp = Bolus transp + Eulerian transp + IF (useGMRedi) THEN + CALL GMREDI_CALC_UVFLOW( + & uTrans, vTrans, bi, bj, k, myThid) + IF (K.GE.2) CALL GMREDI_CALC_WFLOW( + & rTrans, bi, bj, k, myThid) + ENDIF + #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte +#ifdef GM_BOLUS_ADVEC +CADJ STORE uTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE vTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte +#endif #endif /* ALLOW_AUTODIFF_TAMC */ +#endif /* ALLOW_GMREDI */ + #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL C-- Calculate the total vertical diffusivity CALL CALC_DIFFUSIVITY( @@ -477,6 +659,10 @@ I maskUp, O KappaRT,KappaRS, I myThid) +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte +# endif /* ALLOW_AUTODIFF_TAMC */ #endif iMin = 1-OLx+2 @@ -492,120 +678,139 @@ I xA,yA,uTrans,vTrans,rTrans,maskUp, I KappaRT, U fVerT, - I myTime, myThid) - tauAB = 0.5d0 + abEps + I myTime,myIter,myThid) CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,tauAB, + I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme, I theta, gT, - U gTnm1, I myIter, myThid) ENDIF + +#ifdef ALLOW_THERM_SEAICE + IF (useThermSeaIce .AND. k.EQ.1) THEN + CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid ) + ENDIF +#endif + IF ( saltStepping ) THEN CALL CALC_GS( I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown, I xA,yA,uTrans,vTrans,rTrans,maskUp, I KappaRS, U fVerS, - I myTime, myThid) - tauAB = 0.5d0 + abEps + I myTime,myIter,myThid) CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,tauAB, + I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme, I salt, gS, - U gSnm1, I myIter, myThid) ENDIF #ifdef ALLOW_PASSIVE_TRACER +ceh3 needs an IF ( usePASSIVE_TRACER ) THEN IF ( tr1Stepping ) THEN CALL CALC_GTR1( I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown, I xA,yA,uTrans,vTrans,rTrans,maskUp, I KappaRT, U fVerTr1, - I myTime, myThid) - tauAB = 0.5d0 + abEps + I myTime,myIter,myThid) CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,tauAB, + I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme, I Tr1, gTr1, - U gTr1NM1, - I myIter, myThid) + I myIter,myThid) ENDIF #endif +#ifdef ALLOW_PTRACERS + IF ( usePTRACERS ) THEN + CALL PTRACERS_INTEGRATE( + I bi,bj,k, + I xA,yA,uTrans,vTrans,rTrans,maskUp, + X KappaRS, + I myIter,myTime,myThid) + ENDIF +#endif /* ALLOW_PTRACERS */ #ifdef ALLOW_OBCS C-- Apply open boundary conditions IF (useOBCS) THEN - CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid ) + CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid ) END IF #endif /* ALLOW_OBCS */ C-- Freeze water - IF (allowFreezing) THEN + IF ( allowFreezing .AND. .NOT. useSEAICE + & .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k +CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k CADJ & , key = kkey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid ) - END IF + ENDIF C-- end of thermodynamic k loop (Nr:1) ENDDO +cswdice -- add --- +#ifdef ALLOW_THERM_SEAICE +c timeaveraging for ice model values +ceh3 This should be wrapped in an IF ( useThermSeaIce ) THEN + CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid ) +#endif +cswdice --- end add --- + + -#ifdef ALLOW_AUTODIFF_TAMC -C? Patrick? What about this one? -cph Keys iikey and idkey don't seem to be needed -cph since storing occurs on different tape for each -cph impldiff call anyways. -cph Thus, common block comlev1_impl isn't needed either. -cph Storing below needed in the case useGMREDI. - iikey = (ikey-1)*maximpl -#endif /* ALLOW_AUTODIFF_TAMC */ C-- Implicit diffusion IF (implicitDiffusion) THEN IF (tempStepping) THEN #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 1 -CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte +CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, I deltaTtracer, KappaRT, recip_HFacC, - U gTNm1, + U gT, I myThid ) ENDIF IF (saltStepping) THEN #ifdef ALLOW_AUTODIFF_TAMC - idkey = iikey + 2 -CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte +CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, I deltaTtracer, KappaRS, recip_HFacC, - U gSNm1, + U gS, I myThid ) ENDIF #ifdef ALLOW_PASSIVE_TRACER IF (tr1Stepping) THEN #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte +CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL IMPLDIFF( I bi, bj, iMin, iMax, jMin, jMax, I deltaTtracer, KappaRT, recip_HFacC, - U gTr1Nm1, + U gTr1, I myThid ) ENDIF #endif +#ifdef ALLOW_PTRACERS +C Vertical diffusion (implicit) for passive tracers + IF ( usePTRACERS ) THEN + CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid ) + ENDIF +#endif /* ALLOW_PTRACERS */ + #ifdef ALLOW_OBCS C-- Apply open boundary conditions IF (useOBCS) THEN DO K=1,Nr - CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid ) + CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid ) ENDDO END IF #endif /* ALLOW_OBCS */ @@ -613,22 +818,53 @@ C-- End If implicitDiffusion ENDIF -Ccs- +#ifdef ALLOW_TIMEAVE +ceh3 needs an IF ( useTIMEAVE ) THEN + IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN + CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount, + I Nr, deltaTclock, bi, bj, myThid) + ENDIF + useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0. + IF (taveFreq.GT.0. .AND. useVariableK ) THEN + IF (implicitDiffusion) THEN + CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT, + I Nr, 3, deltaTclock, bi, bj, myThid) + ELSE + CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT, + I Nr, 3, deltaTclock, bi, bj, myThid) + ENDIF + ENDIF +#endif /* ALLOW_TIMEAVE */ + +#endif /* SINGLE_LAYER_MODE */ + +C-- end bi,bj loops. ENDDO ENDDO -#ifdef ALLOW_AIM - IF ( useAIM ) THEN - CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid ) - ENDIF - _EXCH_XYZ_R8(gTnm1,myThid) - _EXCH_XYZ_R8(gSnm1,myThid) -#else - IF (staggerTimeStep.AND.useCubedSphereExchange) THEN - _EXCH_XYZ_R8(gTnm1,myThid) - _EXCH_XYZ_R8(gSnm1,myThid) +#ifndef DISABLE_DEBUGMODE + If (debugMode) THEN + CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid) + CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid) +#ifdef ALLOW_PTRACERS + IF ( usePTRACERS ) THEN + CALL PTRACERS_DEBUG(myThid) + ENDIF +#endif /* ALLOW_PTRACERS */ ENDIF -#endif /* ALLOW_AIM */ +#endif + +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_LEAVE('FORWARD_STEP',myThid) +#endif RETURN END