--- MITgcm_contrib/jscott/code_rafmod/thermodynamics.F 2007/08/21 16:34:17 1.1 +++ MITgcm_contrib/jscott/code_rafmod/thermodynamics.F 2009/09/03 20:40:01 1.2 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/jscott/code_rafmod/thermodynamics.F,v 1.1 2007/08/21 16:34:17 jscott Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/jscott/code_rafmod/thermodynamics.F,v 1.2 2009/09/03 20:40:01 jscott Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" @@ -6,6 +6,9 @@ #ifdef ALLOW_GENERIC_ADVDIFF # include "GAD_OPTIONS.h" #endif +#ifdef ALLOW_LONGSTEP +#include "LONGSTEP_OPTIONS.h" +#endif #ifdef ALLOW_AUTODIFF_TAMC # ifdef ALLOW_GMREDI @@ -78,6 +81,7 @@ #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" +#include "RESTART.h" #include "DYNVARS.h" #include "GRID.h" #ifdef ALLOW_GENERIC_ADVDIFF @@ -85,8 +89,11 @@ # include "GAD_SOM_VARS.h" #endif #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP # include "PTRACERS_SIZE.h" -# include "PTRACERS.h" +# include "PTRACERS_PARAMS.h" +# include "PTRACERS_FIELDS.h" +#endif #endif #ifdef ALLOW_TIMEAVE # include "TIMEAVE_STATV.h" @@ -96,6 +103,7 @@ # include "tamc.h" # include "tamc_keys.h" # include "FFIELDS.h" +# include "SURFACE.h" # include "EOS.h" # ifdef ALLOW_KPP # include "KPP.h" @@ -134,7 +142,7 @@ C variable. C kappaRT, - Total diffusion in vertical at level k, for T and S C kappaRS (background + spatially varying, isopycnal term). -C kappaRTr - Total diffusion in vertical at level k, +C kappaRTr - Total diffusion in vertical at level k, C for each passive Tracer C kappaRk - Total diffusion in vertical, all levels, 1 tracer C useVariableK = T when vertical diffusion is not constant @@ -142,7 +150,7 @@ 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 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) @@ -159,13 +167,15 @@ _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly) #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP _RL fVerP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num) _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num) #endif +#endif _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) _RL diffKh3d_x (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) _RL diffKh3d_y (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) - INTEGER iMin, iMax + INTEGER iMin, iMax INTEGER jMin, jMax INTEGER bi, bj INTEGER i, j @@ -177,8 +187,10 @@ LOGICAL useVariableK #endif #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP INTEGER iTracer, ip #endif +#endif CEOP @@ -199,9 +211,17 @@ #endif /* ALLOW_AUTODIFF_TAMC */ C-- Compute correction at the surface for Lin Free Surf. - IF (linFSConserveTr) - & CALL CALC_WSURF_TR(theta,salt,wVel, - & myTime,myIter,myThid) +#ifdef ALLOW_AUTODIFF_TAMC + TsurfCor = 0. _d 0 + SsurfCor = 0. _d 0 +#endif + IF (linFSConserveTr) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics, byte=isbyte +#endif + CALL CALC_WSURF_TR(theta,salt,wVel, + & myTime,myIter,myThid) + ENDIF DO bj=myByLo(myThid),myByHi(myThid) @@ -212,8 +232,10 @@ CHPF$& ,kappaRT,kappaRS CHPF$& ) # ifdef ALLOW_PTRACERS +# ifndef ALLOW_LONGSTEP CHPF$ INDEPENDENT, NEW (fVerP,kappaRTr) # endif +# endif #endif /* ALLOW_AUTODIFF_TAMC */ DO bi=myBxLo(myThid),myBxHi(myThid) @@ -267,6 +289,7 @@ ENDDO #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP IF ( usePTRACERS ) THEN DO ip=1,PTRACERS_num DO j=1-OLy,sNy+OLy @@ -289,6 +312,7 @@ ENDDO ENDIF #endif +#endif #ifdef ALLOW_ADAMSBASHFORTH_3 C- Apply AB on T,S : @@ -318,7 +342,7 @@ 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-- Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h C-- MOST of THERMODYNAMICS will be disabled #ifndef SINGLE_LAYER_MODE @@ -344,7 +368,7 @@ 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. -#ifdef GAD_ALLOW_SOM_ADVECT +#ifdef GAD_ALLOW_TS_SOM_ADV IF ( tempSOM_Advection ) THEN #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) @@ -352,27 +376,27 @@ #endif CALL GAD_SOM_ADVECT( I tempImplVertAdv, tempAdvScheme, tempVertAdvScheme, - I GAD_TEMPERATURE, + I GAD_TEMPERATURE, dTtracerLev, I uVel, vVel, wVel, theta, U som_T, O gT, I bi,bj,myTime,myIter,myThid) ELSEIF (tempMultiDimAdvec) THEN -#else /* GAD_ALLOW_SOM_ADVECT */ +#else /* GAD_ALLOW_TS_SOM_ADV */ IF (tempMultiDimAdvec) THEN -#endif /* GAD_ALLOW_SOM_ADVECT */ +#endif /* GAD_ALLOW_TS_SOM_ADV */ #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 GAD_TEMPERATURE, dTtracerLev, I uVel, vVel, wVel, theta, O gT, I bi,bj,myTime,myIter,myThid) ENDIF -#ifdef GAD_ALLOW_SOM_ADVECT +#ifdef GAD_ALLOW_TS_SOM_ADV IF ( saltSOM_Advection ) THEN #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) @@ -380,22 +404,22 @@ #endif CALL GAD_SOM_ADVECT( I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme, - I GAD_SALINITY, + I GAD_SALINITY, dTtracerLev, I uVel, vVel, wVel, salt, U som_S, O gS, I bi,bj,myTime,myIter,myThid) ELSEIF (saltMultiDimAdvec) THEN -#else /* GAD_ALLOW_SOM_ADVECT */ +#else /* GAD_ALLOW_TS_SOM_ADV */ IF (saltMultiDimAdvec) THEN -#endif /* GAD_ALLOW_SOM_ADVECT */ +#endif /* GAD_ALLOW_TS_SOM_ADV */ #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 GAD_SALINITY, dTtracerLev, I uVel, vVel, wVel, salt, O gS, I bi,bj,myTime,myIter,myThid) @@ -405,6 +429,7 @@ C call the multi-dimensional method for PTRACERS regardless C of whether multiDimAdvection is set or not. #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP IF ( usePTRACERS ) THEN #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) @@ -412,6 +437,7 @@ #endif CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid ) ENDIF +#endif /* ALLOW_LONGSTEP */ #endif /* ALLOW_PTRACERS */ #endif /* DISABLE_MULTIDIM_ADVECTION */ @@ -421,13 +447,14 @@ #endif Cjrs - IF (diffKhT.ne.0) THEN + IF (diffKhT.NE.0) THEN CALL CALC_MLD(bi,bj,diffkh3d_x,diffkh3d_y,myThid) ENDIF + C-- Start of thermodynamics loop DO k=Nr,1,-1 -#ifdef ALLOW_AUTODIFF_TAMC +#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! @@ -508,8 +535,8 @@ ENDIF ENDIF # ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte -CADJ STORE wfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte +CADJ STORE wfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte # ifdef GM_BOLUS_ADVEC CADJ STORE ufld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE vfld(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte @@ -528,7 +555,7 @@ O kappaRT,kappaRS, I myThid) ENDIF -# ifdef ALLOW_AUTODIFF_TAMC +# ifdef ALLOW_AUTODIFF_TAMC CADJ STORE kappaRT(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE kappaRS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte # endif /* ALLOW_AUTODIFF_TAMC */ @@ -542,10 +569,12 @@ C-- Calculate active tracer tendencies (gT,gS,...) C and step forward storing result in gT, gS, etc. C-- -# ifdef ALLOW_AUTODIFF_TAMC +# ifdef ALLOW_AUTODIFF_TAMC # if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI) +# ifdef GM_NON_UNITY_DIAGONAL CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte +# endif # ifdef GM_EXTRA_DIAGONAL CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte @@ -553,7 +582,7 @@ # endif # endif /* ALLOW_AUTODIFF_TAMC */ C -#ifdef ALLOW_AUTODIFF_TAMC +#ifdef ALLOW_AUTODIFF_TAMC # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL) cph-test CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:) @@ -566,7 +595,7 @@ #endif /* ALLOW_AUTODIFF_TAMC */ C IF ( tempStepping ) THEN -#ifdef ALLOW_AUTODIFF_TAMC +#ifdef ALLOW_AUTODIFF_TAMC CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL) CADJ STORE gt(:,:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte @@ -583,14 +612,14 @@ #ifdef ALLOW_ADAMSBASHFORTH_3 IF ( AdamsBashforth_T ) THEN CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme, + I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,dTtracerLev(k), I gtNm(1-Olx,1-Oly,1,1,1,m2), U gT, I myIter, myThid) ELSE #endif CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme, + I bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,dTtracerLev(k), I theta, U gT, I myIter, myThid) @@ -617,14 +646,14 @@ #ifdef ALLOW_ADAMSBASHFORTH_3 IF ( AdamsBashforth_S ) THEN CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme, + I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,dTtracerLev(k), I gsNm(1-Olx,1-Oly,1,1,1,m2), U gS, I myIter, myThid) ELSE #endif CALL TIMESTEP_TRACER( - I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme, + I bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,dTtracerLev(k), I salt, U gS, I myIter, myThid) @@ -634,6 +663,7 @@ ENDIF #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP IF ( usePTRACERS ) THEN IF ( .NOT.implicitDiffusion ) THEN CALL PTRACERS_CALC_DIFF( @@ -642,17 +672,18 @@ O kappaRTr, I myThid) ENDIF -# ifdef ALLOW_AUTODIFF_TAMC +# ifdef ALLOW_AUTODIFF_TAMC CADJ STORE kappaRTr(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte # endif /* ALLOW_AUTODIFF_TAMC */ CALL PTRACERS_INTEGRATE( I bi,bj,k, I xA, yA, maskUp, uFld, vFld, wFld, I uTrans, vTrans, rTrans, rTransKp1, - I kappaRTr, + I kappaRTr,diffKh3d_x,diffKh3d_y, U fVerP, I myTime,myIter,myThid) ENDIF +#endif /*ALLOW_LONGSTEP */ #endif /* ALLOW_PTRACERS */ #ifdef ALLOW_OBCS @@ -664,7 +695,7 @@ 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 +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 @@ -677,22 +708,71 @@ C-- end of thermodynamic k loop (Nr:1) ENDDO -C All explicit advection/diffusion/sources should now be -C done. The updated tracer field is in gPtr. Accumalate -C explicit tendency and also reset gPtr to initial tracer +#ifdef ALLOW_DOWN_SLOPE + IF ( tempStepping .AND. useDOWN_SLOPE ) THEN + IF ( usingPCoords ) THEN + CALL DWNSLP_APPLY( + I GAD_TEMPERATURE, bi, bj, kSurfC, + I recip_drF, recip_hFacC, recip_rA, + I dTtracerLev, + I theta, + U gT, + I myTime, myIter, myThid ) + ELSE + CALL DWNSLP_APPLY( + I GAD_TEMPERATURE, bi, bj, kLowC, + I recip_drF, recip_hFacC, recip_rA, + I dTtracerLev, + I theta, + U gT, + I myTime, myIter, myThid ) + ENDIF + ENDIF + IF ( saltStepping .AND. useDOWN_SLOPE ) THEN + IF ( usingPCoords ) THEN + CALL DWNSLP_APPLY( + I GAD_SALINITY, bi, bj, kSurfC, + I recip_drF, recip_hFacC, recip_rA, + I dTtracerLev, + I salt, + U gS, + I myTime, myIter, myThid ) + ELSE + CALL DWNSLP_APPLY( + I GAD_SALINITY, bi, bj, kLowC, + I recip_drF, recip_hFacC, recip_rA, + I dTtracerLev, + I salt, + U gS, + I myTime, myIter, myThid ) + ENDIF + ENDIF +#ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP + IF ( usePTRACERS .AND. useDOWN_SLOPE ) THEN + CALL PTRACERS_DWNSLP_APPLY( + I bi, bj, myTime, myIter, myThid ) + ENDIF +#endif /*ALLOW_LONGSTEP */ +#endif /* ALLOW_PTRACERS */ +#endif /* ALLOW_DOWN_SLOPE */ + +C All explicit advection/diffusion/sources should now be +C done. The updated tracer field is in gPtr. Accumalate +C explicit tendency and also reset gPtr to initial tracer C field for implicit matrix calculation #ifdef ALLOW_MATRIX - IF (useMATRIX) + IF (useMATRIX) & CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid) -#endif +#endif iMin = 1 iMax = sNx jMin = 1 jMax = sNy -C-- Implicit vertical advection & diffusion +C-- Implicit vertical advection & diffusion IF ( tempStepping .AND. implicitDiffusion ) THEN CALL CALC_3D_DIFFUSIVITY( I bi,bj,iMin,iMax,jMin,jMax, @@ -704,6 +784,7 @@ IF ( tempImplVertAdv ) THEN CALL GAD_IMPLICIT_R( I tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE, + I dTtracerLev, I kappaRk, wVel, theta, U gT, I bi, bj, myTime, myIter, myThid ) @@ -734,7 +815,7 @@ c I Nr, 3, deltaTclock, bi, bj, myThid) ENDIF ENDIF -#endif /* ALLOW_TIMEAVE */ +#endif /* ALLOW_TIMEAVE */ IF ( saltStepping .AND. implicitDiffusion ) THEN CALL CALC_3D_DIFFUSIVITY( @@ -748,6 +829,7 @@ IF ( saltImplVertAdv ) THEN CALL GAD_IMPLICIT_R( I saltImplVertAdv, saltAdvScheme, GAD_SALINITY, + I dTtracerLev, I kappaRk, wVel, salt, U gS, I bi, bj, myTime, myIter, myThid ) @@ -767,12 +849,14 @@ ENDIF #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP IF ( usePTRACERS ) THEN C-- Vertical advection/diffusion (implicit) for passive tracers CALL PTRACERS_IMPLICIT( U kappaRk, I bi, bj, myTime, myIter, myThid ) ENDIF +#endif /* ALLOW_LONGSTEP */ #endif /* ALLOW_PTRACERS */ #ifdef ALLOW_OBCS @@ -807,9 +891,11 @@ CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid) #endif #ifdef ALLOW_PTRACERS +#ifndef ALLOW_LONGSTEP IF ( usePTRACERS ) THEN CALL PTRACERS_DEBUG(myThid) ENDIF +#endif /* ALLOW_LONGSTEP */ #endif /* ALLOW_PTRACERS */ ENDIF #endif /* ALLOW_DEBUG */