C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/thermodynamics.F,v 1.161 2016/08/24 15:00:51 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif #if (defined ALLOW_PTRACERS) && (!defined ALLOW_LONGSTEP) # define DO_PTRACERS_HERE #endif #ifdef ALLOW_AUTODIFF # ifdef ALLOW_GMREDI # include "GMREDI_OPTIONS.h" # endif # ifdef ALLOW_KPP # include "KPP_OPTIONS.h" # endif #endif /* ALLOW_AUTODIFF */ 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 \ev C !USES: IMPLICIT NONE C == Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "RESTART.h" #include "DYNVARS.h" #include "GRID.h" #include "SURFACE.h" #include "FFIELDS.h" #ifdef ALLOW_GENERIC_ADVDIFF # include "GAD.h" #endif #ifdef DO_PTRACERS_HERE # include "PTRACERS_SIZE.h" # include "PTRACERS_PARAMS.h" # include "PTRACERS_FIELDS.h" #endif #ifdef ALLOW_AUTODIFF # include "tamc.h" # include "tamc_keys.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 # ifdef ALLOW_SALT_PLUME # include "SALT_PLUME.h" # endif #endif /* ALLOW_AUTODIFF */ 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 #ifdef ALLOW_GENERIC_ADVDIFF # ifdef ALLOW_MONITOR C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE # endif /* ALLOW_MONITOR */ C !LOCAL VARIABLES: C == Local variables C uFld,vFld,wFld :: Local copy of velocity field (3 components) C kappaRk :: Total diffusion in vertical, all levels, 1 tracer C recip_hFacNew :: reciprocal of futur time-step hFacC C bi, bj :: Tile indices C i, j, k :: loop indices _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL kappaRk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RS recip_hFacNew(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) INTEGER bi, bj INTEGER i, j, k # ifdef ALLOW_MONITOR LOGICAL monOutputCFL _RL wrTime _RL trAdvCFL(3,nSx,nSy) # endif /* ALLOW_MONITOR */ CEOP #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_ENTER('THERMODYNAMICS',myThid) #endif # ifdef ALLOW_MONITOR monOutputCFL = .FALSE. IF ( monitorSelect.GE.2 ) THEN wrTime = myTime IF ( .NOT.staggerTimeStep ) wrTime = myTime + deltaTClock monOutputCFL = & DIFFERENT_MULTIPLE( monitorFreq, wrTime, deltaTClock ) ENDIF # endif /* ALLOW_MONITOR */ #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 */ C-- Compute correction at the surface for Lin Free Surf. #ifdef ALLOW_AUTODIFF 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 #ifdef ALLOW_LAYERS IF ( useLayers ) THEN CALL LAYERS_WSURF_TR( theta, salt, wVel, & myTime, myIter, myThid ) ENDIF #endif /* ALLOW_LAYERS */ #ifdef DO_PTRACERS_HERE #ifdef ALLOW_AUTODIFF DO k=1,PTRACERS_numInUse meanSurfCorPTr(k) = 0.0 _d 0 ENDDO #endif IF ( PTRACERS_calcSurfCor ) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE ptracer = comlev1, key = ikey_dynamics, byte=isbyte #endif CALL PTRACERS_CALC_WSURF_TR(wVel,myTime,myIter,myThid) ENDIF #endif /* DO_PTRACERS_HERE */ DO bj=myByLo(myThid),myByHi(myThid) 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 k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx recip_hFacNew(i,j,k) = 0. _d 0 C This is currently also used by IVDC and Diagnostics kappaRk(i,j,k) = 0. _d 0 ENDDO ENDDO ENDDO C-- Compute new reciprocal hFac for implicit calculation #ifdef NONLIN_FRSURF IF ( nonlinFreeSurf.GT.0 ) THEN IF ( select_rStar.GT.0 ) THEN # ifndef DISABLE_RSTAR_CODE DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj) & / rStarExpC(i,j,bi,bj) ENDDO ENDDO ENDDO # endif /* DISABLE_RSTAR_CODE */ ELSEIF ( selectSigmaCoord.NE.0 ) THEN # ifndef DISABLE_SIGMA_CODE DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj) & /( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf & *dBHybSigF(k)*recip_drF(k) & *recip_hFacC(i,j,k,bi,bj) & ) ENDDO ENDDO ENDDO # endif /* DISABLE_RSTAR_CODE */ ELSE DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN recip_hFacNew(i,j,k) = 1. _d 0 / hFac_surfC(i,j,bi,bj) ELSE recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDIF ELSE #endif /* NONLIN_FRSURF */ DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx recip_hFacNew(i,j,k) = _recip_hFacC(i,j,k,bi,bj) ENDDO ENDDO ENDDO #ifdef NONLIN_FRSURF ENDIF #endif /* NONLIN_FRSURF */ C-- Set up 3-D velocity field that we use to advect tracers: C- just do a local copy: DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uFld(i,j,k) = uVel(i,j,k,bi,bj) vFld(i,j,k) = vVel(i,j,k,bi,bj) wFld(i,j,k) = wVel(i,j,k,bi,bj) ENDDO ENDDO ENDDO #ifdef ALLOW_GMREDI C- add Bolus velocity to Eulerian-mean velocity: IF (useGMRedi) THEN CALL GMREDI_RESIDUAL_FLOW( U uFld, vFld, wFld, I bi, bj, myIter, myThid ) ENDIF #endif /* ALLOW_GMREDI */ #ifdef ALLOW_MONITOR IF ( monOutputCFL ) THEN CALL MON_CALC_ADVCFL_TILE( Nr, bi, bj, I uFld, vFld, wFld, dTtracerLev, O trAdvCFL(1,bi,bj), I myIter, myThid ) c WRITE(standardMessageUnit,'(A,I8,2I3,A,1P3E14.6)') c & ' trAdv_CFL: it,bi,bj=', myIter,bi,bj, c & ' , CFL =', (trAdvCFL(i,bi,bj),i=1,3) ENDIF #endif /* ALLOW_MONITOR */ C- Apply AB on T,S : moved inside TEMP/SALT_INTEGRATE #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE recip_hFacNew(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte CADJ STORE uFld (:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE vFld (:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE wFld (:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte # ifdef ALLOW_SALT_PLUME CADJ STORE saltPlumeFlux(:,:,bi,bj) = CADJ & comlev1_bibj, key=itdkey,kind = isbyte CADJ STORE saltPlumeDepth(:,:,bi,bj) = CADJ & comlev1_bibj, key=itdkey,kind = isbyte # endif # if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI) # ifdef GM_NON_UNITY_DIAGONAL CADJ STORE kux(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE kvy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte # endif # ifdef GM_EXTRA_DIAGONAL CADJ STORE kuz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte CADJ STORE kvz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte # endif # endif #endif /* ALLOW_AUTODIFF_TAMC */ C-- Calculate active tracer tendencies and step forward in time. C Active tracer arrays are updated while adjustments (filters, C conv.adjustment) are applied later in TRACERS_CORRECTION_STEP. IF ( tempStepping ) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('TEMP_INTEGRATE',myThid) #endif CALL TEMP_INTEGRATE( I bi, bj, recip_hFacNew, I uFld, vFld, wFld, U kappaRk, I myTime, myIter, myThid ) ENDIF IF ( saltStepping ) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('SALT_INTEGRATE',myThid) #endif CALL SALT_INTEGRATE( I bi, bj, recip_hFacNew, I uFld, vFld, wFld, U kappaRk, I myTime, myIter, myThid ) ENDIF #ifdef DO_PTRACERS_HERE C-- Calculate passive tracer tendencies and step forward in time. C Passive tracer arrays are updated while adjustments (filters, C conv.adjustment) are applied later in TRACERS_CORRECTION_STEP. C Also apply open boundary conditions for each passive tracer IF ( usePTRACERS ) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('PTRACERS_INTEGRATE',myThid) #endif CALL PTRACERS_INTEGRATE( I bi, bj, recip_hFacNew, I uFld, vFld, wFld, U kappaRk, I myTime, myIter, myThid ) ENDIF #endif /* DO_PTRACERS_HERE */ #ifdef ALLOW_OBCS C-- Apply open boundary conditions IF ( useOBCS ) THEN CALL OBCS_APPLY_TS( bi, bj, 0, theta, salt, myThid ) ENDIF #endif /* ALLOW_OBCS */ #ifdef ALLOW_FRICTION_HEATING #ifdef ALLOW_DIAGNOSTICS IF ( addFrictionHeating .AND. useDiagnostics ) THEN CALL DIAGNOSTICS_FILL_RS( frictionHeating, 'HeatDiss', & 0, Nr, 1, bi, bj, myThid ) ENDIF #endif /* ALLOW_DIAGNOSTICS */ C- Reset frictionHeating to zero IF ( addFrictionHeating .AND. .NOT.staggerTimeStep ) THEN DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx frictionHeating(i,j,k,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDIF #endif /* ALLOW_FRICTION_HEATING */ C-- end bi,bj loops. ENDDO ENDDO #ifdef ALLOW_MONITOR IF ( monOutputCFL ) THEN CALL MON_CALC_ADVCFL_GLOB( I trAdvCFL, myIter, myThid ) ENDIF #endif /* ALLOW_MONITOR */ #ifdef ALLOW_DEBUG IF ( debugLevel.GE.debLevD ) 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) #ifndef ALLOW_ADAMSBASHFORTH_3 CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid) CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid) #endif #ifdef DO_PTRACERS_HERE IF ( usePTRACERS ) THEN CALL PTRACERS_DEBUG(myThid) ENDIF #endif /* DO_PTRACERS_HERE */ ENDIF #endif /* ALLOW_DEBUG */ #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_LEAVE('THERMODYNAMICS',myThid) #endif #endif /* ALLOW_GENERIC_ADVDIFF */ RETURN END