Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F122652881
thermodynamics.F
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Jul 21, 09:32
Size
12 KB
Mime Type
text/x-c
Expires
Wed, Jul 23, 09:32 (2 d)
Engine
blob
Format
Raw Data
Handle
27535350
Attached To
R2795 mitgcm_lac_leman_abirani
thermodynamics.F
View Options
C $Header: /u/gcmpack/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
Event Timeline
Log In to Comment