C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/thermodynamics.F,v 1.72 2004/07/06 00:58:40 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" #ifdef ALLOW_PTRACERS # include "PTRACERS_OPTIONS.h" #endif #ifdef ALLOW_AUTODIFF_TAMC # ifdef ALLOW_GMREDI # include "GMREDI_OPTIONS.h" # endif # ifdef ALLOW_KPP # include "KPP_OPTIONS.h" # endif #endif /* ALLOW_AUTODIFF_TAMC */ CBOP C !ROUTINE: THERMODYNAMICS C !INTERFACE: SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid) 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 # ifdef ALLOW_GMREDI # include "GMREDI.h" # endif # ifdef ALLOW_EBM # include "EBM.h" # endif #endif /* ALLOW_AUTODIFF_TAMC */ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C myTime - Current time in simulation C myIter - Current iteration number in simulation C myThid - Thread number for this instance of the routine. _RL myTime 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 C transport C o uTrans: Zonal transport C o vTrans: Meridional transport C o rTrans: Vertical transport C rTransKp1 o vertical volume transp. at interface k+1 C maskUp o maskUp: land/water mask for W points C fVer[STUV] o fVer: Vertical flux term - note fVer C is "pipelined" in the vertical C so we need an fVer for each C variable. 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. _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 fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) #ifdef ALLOW_PASSIVE_TRACER _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) #endif #ifdef ALLOW_PTRACERS _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num) #endif _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) _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 kp1Msk LOGICAL useVariableK INTEGER iMin, iMax INTEGER jMin, jMax INTEGER bi, bj INTEGER i, j INTEGER k, km1, kup, kDown INTEGER iTracer, ip CEOP #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) & CALL DEBUG_ENTER('THERMODYNAMICS',myThid) #endif #ifdef ALLOW_AUTODIFF_TAMC C-- dummy statement to end declaration part ikey = 1 itdkey = 1 #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_AUTODIFF_TAMC C-- HPF directive to help TAMC CHPF$ INDEPENDENT #endif /* ALLOW_AUTODIFF_TAMC */ DO bj=myByLo(myThid),myByHi(myThid) #ifdef ALLOW_AUTODIFF_TAMC C-- HPF directive to help TAMC CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS CHPF$& ,utrans,vtrans,xA,yA CHPF$& ,KappaRT,KappaRS CHPF$& ) #endif /* ALLOW_AUTODIFF_TAMC */ DO bi=myBxLo(myThid),myBxHi(myThid) #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 #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 rTrans (i,j) = 0. _d 0 rTransKp1(i,j) = 0. _d 0 fVerT (i,j,1) = 0. _d 0 fVerT (i,j,2) = 0. _d 0 fVerS (i,j,1) = 0. _d 0 fVerS (i,j,2) = 0. _d 0 #ifdef ALLOW_PASSIVE_TRACER fVerTr1(i,j,1) = 0. _d 0 fVerTr1(i,j,2) = 0. _d 0 #endif #ifdef ALLOW_PTRACERS DO ip=1,PTRACERS_num fVerP (i,j,1,ip) = 0. _d 0 fVerP (i,j,2,ip) = 0. _d 0 ENDDO #endif ENDDO ENDDO DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C This is currently also used by IVDC and Diagnostics 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 ENDDO ENDDO ENDDO c iMin = 1-OLx c iMax = sNx+OLx c jMin = 1-OLy c jMax = sNy+OLy #ifdef ALLOW_AUTODIFF_TAMC 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 */ #ifdef ALLOW_AUTODIFF_TAMC cph avoids recomputation of integrate_for_w CADJ STORE wvel (:,:,:,bi,bj) = 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_AUTODIFF_TAMC 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=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 */ #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 #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) & CALL DEBUG_CALL('GAD_ADVECTION',myThid) #endif CALL GAD_ADVECTION( I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme, I GAD_TEMPERATURE, I uVel, vVel, wVel, theta, O gT, I bi,bj,myTime,myIter,myThid) ENDIF IF (saltMultiDimAdvec) THEN #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) & CALL DEBUG_CALL('GAD_ADVECTION',myThid) #endif CALL GAD_ADVECTION( I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme, I GAD_SALINITY, I uVel, vVel, wVel, salt, O gS, I bi,bj,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 #ifdef ALLOW_DEBUG 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 */ #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) & CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid) #endif C-- Start of thermodynamics loop DO k=Nr,1,-1 #ifdef ALLOW_AUTODIFF_TAMC C? Patrick Is this formula correct? cph Yes, but I rewrote it. cph Also, the KappaR? need the index and subscript k! kkey = (itdkey-1)*Nr + k #endif /* ALLOW_AUTODIFF_TAMC */ C-- km1 Points to level above k (=k-1) C-- kup Cycles through 1,2 to point to layer above C-- kDown Cycles through 2,1 to point to current layer km1 = MAX(1,k-1) kup = 1+MOD(k+1,2) kDown= 1+MOD(k,2) iMin = 1-OLx iMax = sNx+OLx jMin = 1-OLy jMax = sNy+OLy kp1Msk=1. IF (k.EQ.Nr) kp1Msk=0. DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx rTransKp1(i,j) = kp1Msk*rTrans(i,j) ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte #endif C-- Get temporary terms used by tendency routines CALL CALC_COMMON_FACTORS ( I bi,bj,iMin,iMax,jMin,jMax,k, O xA,yA,uTrans,vTrans,rTrans,maskUp, I myThid) IF (k.EQ.1) THEN C- Surface interface : DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx rTrans(i,j) = 0. ENDDO ENDDO ELSE C- Interior interface : DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj) ENDDO ENDDO ENDIF #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 rTrans(:,:) = 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 #endif #endif /* ALLOW_AUTODIFF_TAMC */ #endif /* ALLOW_GMREDI */ #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL C-- Calculate the total vertical diffusivity CALL CALC_DIFFUSIVITY( I bi,bj,iMin,iMax,jMin,jMax,k, 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 iMax = sNx+OLx-1 jMin = 1-OLy+2 jMax = sNy+OLy-1 C-- Calculate active tracer tendencies (gT,gS,...) C and step forward storing result in gTnm1, gSnm1, etc. IF ( tempStepping ) THEN CALL 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) CALL TIMESTEP_TRACER( I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme, I theta, gT, I myIter, myThid) ENDIF IF ( saltStepping ) THEN CALL CALC_GS( I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown, I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp, I KappaRS, U fVerS, I myTime,myIter,myThid) CALL TIMESTEP_TRACER( I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme, I salt, gS, 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,rTransKp1,maskUp, I KappaRT, U fVerTr1, I myTime,myIter,myThid) CALL TIMESTEP_TRACER( I bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme, I Tr1, gTr1, I myIter,myThid) ENDIF #endif #ifdef ALLOW_PTRACERS IF ( usePTRACERS ) THEN CALL PTRACERS_INTEGRATE( I bi,bj,k, I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp, X fVerP, 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, gT, gS, myThid ) END IF #endif /* ALLOW_OBCS */ C-- Freeze water C this bit of code is left here for backward compatibility. C freezing at surface level has been moved to FORWARD_STEP IF ( useOldFreezing .AND. .NOT. useSEAICE & .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN #ifdef ALLOW_AUTODIFF_TAMC 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 ) ENDIF C-- end of thermodynamic k loop (Nr:1) ENDDO C-- Implicit vertical advection & diffusion #ifdef INCLUDE_IMPLVERTADV_CODE IF ( tempImplVertAdv ) THEN CALL GAD_IMPLICIT_R( I tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE, I kappaRT, wVel, theta, U gT, I bi, bj, myTime, myIter, myThid ) ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN #else /* INCLUDE_IMPLVERTADV_CODE */ IF ( tempStepping .AND. implicitDiffusion ) THEN #endif /* INCLUDE_IMPLVERTADV_CODE */ #ifdef ALLOW_AUTODIFF_TAMC 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 gT, I myThid ) ENDIF #ifdef INCLUDE_IMPLVERTADV_CODE IF ( saltImplVertAdv ) THEN CALL GAD_IMPLICIT_R( I saltImplVertAdv, saltAdvScheme, GAD_SALINITY, I kappaRS, wVel, salt, U gS, I bi, bj, myTime, myIter, myThid ) ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN #else /* INCLUDE_IMPLVERTADV_CODE */ IF ( saltStepping .AND. implicitDiffusion ) THEN #endif /* INCLUDE_IMPLVERTADV_CODE */ #ifdef ALLOW_AUTODIFF_TAMC 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 gS, I myThid ) ENDIF #ifdef ALLOW_PASSIVE_TRACER IF ( tr1Stepping .AND. implicitDiffusion ) THEN #ifdef ALLOW_AUTODIFF_TAMC 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 gTr1, I myThid ) ENDIF #endif #ifdef ALLOW_PTRACERS c #ifdef INCLUDE_IMPLVERTADV_CODE c IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN c ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN c #else IF ( usePTRACERS .AND. implicitDiffusion ) THEN C-- Vertical diffusion (implicit) for passive tracers CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid ) ENDIF #endif /* ALLOW_PTRACERS */ #ifdef ALLOW_OBCS C-- Apply open boundary conditions IF ( ( implicitDiffusion & .OR. tempImplVertAdv & .OR. saltImplVertAdv & ) .AND. useOBCS ) THEN DO K=1,Nr CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid ) ENDDO ENDIF #endif /* ALLOW_OBCS */ #ifdef ALLOW_TIMEAVE #ifndef HRCUBE IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount, 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 /* ndef HRCUBE */ #endif /* ALLOW_TIMEAVE */ #endif /* SINGLE_LAYER_MODE */ C-- end bi,bj loops. ENDDO ENDDO #ifdef ALLOW_DEBUG 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 #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) & CALL DEBUG_LEAVE('THERMODYNAMICS',myThid) #endif RETURN END