Back to home page

MITgcm

 
 

    


File indexing completed on 2021-06-27 05:11:26 UTC

view on githubraw file Latest commit 4e4ad91a on 2021-06-26 16:30:07 UTC
4e66ab0b67 Oliv*0001 #include "LONGSTEP_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: LONGSTEP_THERMODYNAMICS
                0005 C     !INTERFACE:
7908448559 Jean*0006       SUBROUTINE LONGSTEP_THERMODYNAMICS( myTime, myIter, myThid )
4e66ab0b67 Oliv*0007 C     !DESCRIPTION: \bv
                0008 C     *==========================================================*
                0009 C     | SUBROUTINE LONGSTEP_THERMODYNAMICS
                0010 C     | o Controlling routine for the prognostics of passive tracers
                0011 C     |   with longer time step.
                0012 C     *===========================================================
                0013 C     | This is a copy of THERMODYNAMICS, but only with the
                0014 C     | parts relevant to ptracers, and dynamics fields replaced
                0015 C     | by their longstep averages.
                0016 C     | When THERMODYNAMICS is changed, this routine probably has
                0017 C     | to be changed too :(
                0018 C     *==========================================================*
                0019 C     \ev
                0020 
                0021 C     !USES:
                0022       IMPLICIT NONE
                0023 C     == Global variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "RESTART.h"
                0028 #include "DYNVARS.h"
                0029 #include "GRID.h"
36487395b9 Jean*0030 #include "SURFACE.h"
4e66ab0b67 Oliv*0031 #ifdef ALLOW_GENERIC_ADVDIFF
                0032 # include "GAD.h"
                0033 #endif
                0034 #include "LONGSTEP_PARAMS.h"
                0035 #include "LONGSTEP.h"
                0036 #ifdef ALLOW_PTRACERS
                0037 # include "PTRACERS_SIZE.h"
                0038 # include "PTRACERS_PARAMS.h"
                0039 # include "PTRACERS_FIELDS.h"
                0040 #endif
                0041 
                0042 C     !INPUT/OUTPUT PARAMETERS:
                0043 C     == Routine arguments ==
7908448559 Jean*0044 C     myTime :: Current time in simulation
                0045 C     myIter :: Current iteration number in simulation
                0046 C     myThid :: Thread number for this instance of the routine.
4e66ab0b67 Oliv*0047       _RL myTime
                0048       INTEGER myIter
                0049       INTEGER myThid
                0050 
                0051 #ifdef ALLOW_LONGSTEP
                0052 C     !LOCAL VARIABLES:
                0053 C     == Local variables
7908448559 Jean*0054 C     uFld,vFld,wFld :: Local copy of velocity field (3 components)
                0055 C     kappaRk        :: Total diffusion in vertical, all levels, 1 tracer
                0056 C     useVariableK   :: T when vertical diffusion is not constant
                0057 C     iMin, iMax     :: Ranges and sub-block indices on which calculations
                0058 C     jMin, jMax        are applied.
                0059 C     bi, bj         :: Tile indices
                0060 C     i, j, k        :: loop indices
                0061       _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0062       _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0063       _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
36487395b9 Jean*0064       _RL kappaRk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0065       _RS recip_hFacNew(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
4e66ab0b67 Oliv*0066       INTEGER iMin, iMax
                0067       INTEGER jMin, jMax
                0068       INTEGER bi, bj
7908448559 Jean*0069       INTEGER i, j, k
4e66ab0b67 Oliv*0070 CEOP
                0071 
                0072 #ifdef ALLOW_DEBUG
1378a26402 Jean*0073       IF (debugMode) CALL DEBUG_ENTER('LONGSTEP_THERMODYNAMICS',myThid)
4e66ab0b67 Oliv*0074 #endif
                0075 
                0076 C     time for a ptracer time step?
                0077       IF ( LS_doTimeStep ) THEN
                0078 
                0079       DO bj=myByLo(myThid),myByHi(myThid)
                0080        DO bi=myBxLo(myThid),myBxHi(myThid)
                0081 
                0082 C--   Set up work arrays with valid (i.e. not NaN) values
                0083 C     These inital values do not alter the numerical results. They
                0084 C     just ensure that all memory references are to valid floating
                0085 C     point numbers. This prevents spurious hardware signals due to
                0086 C     uninitialised but inert locations.
                0087 
                0088         DO k=1,Nr
                0089          DO j=1-OLy,sNy+OLy
                0090           DO i=1-OLx,sNx+OLx
                0091 C This is currently also used by IVDC and Diagnostics
                0092            kappaRk(i,j,k)    = 0. _d 0
                0093           ENDDO
                0094          ENDDO
                0095         ENDDO
                0096 
36487395b9 Jean*0097 C--     Compute new reciprocal hFac for implicit calculation
                0098 #ifdef NONLIN_FRSURF
                0099         IF ( nonlinFreeSurf.GT.0 ) THEN
                0100          IF ( select_rStar.GT.0 ) THEN
                0101 # ifndef DISABLE_RSTAR_CODE
                0102           DO k=1,Nr
                0103            DO j=1-OLy,sNy+OLy
                0104             DO i=1-OLx,sNx+OLx
                0105              recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj)
                0106      &                            / rStarExpC(i,j,bi,bj)
                0107             ENDDO
                0108            ENDDO
                0109           ENDDO
                0110 # endif /* DISABLE_RSTAR_CODE */
                0111          ELSEIF ( selectSigmaCoord.NE.0 ) THEN
                0112 # ifndef DISABLE_SIGMA_CODE
                0113           DO k=1,Nr
                0114            DO j=1-OLy,sNy+OLy
                0115             DO i=1-OLx,sNx+OLx
                0116              recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj)
7908448559 Jean*0117      &        /( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
36487395b9 Jean*0118      &                    *dBHybSigF(k)*recip_drF(k)
                0119      &                    *recip_hFacC(i,j,k,bi,bj)
                0120      &         )
                0121             ENDDO
                0122            ENDDO
                0123           ENDDO
                0124 # endif /* DISABLE_RSTAR_CODE */
                0125          ELSE
                0126           DO k=1,Nr
                0127            DO j=1-OLy,sNy+OLy
                0128             DO i=1-OLx,sNx+OLx
                0129              IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
                0130               recip_hFacNew(i,j,k) = 1. _d 0 / hFac_surfC(i,j,bi,bj)
                0131              ELSE
                0132               recip_hFacNew(i,j,k) = recip_hFacC(i,j,k,bi,bj)
                0133              ENDIF
                0134             ENDDO
                0135            ENDDO
                0136           ENDDO
                0137          ENDIF
                0138         ELSE
                0139 #endif /* NONLIN_FRSURF */
                0140           DO k=1,Nr
                0141            DO j=1-OLy,sNy+OLy
                0142             DO i=1-OLx,sNx+OLx
                0143              recip_hFacNew(i,j,k) = _recip_hFacC(i,j,k,bi,bj)
                0144             ENDDO
                0145            ENDDO
                0146           ENDDO
                0147 #ifdef NONLIN_FRSURF
                0148         ENDIF
                0149 #endif /* NONLIN_FRSURF */
                0150 
f941e53e0a Jean*0151         iMin = 1-OLx
                0152         iMax = sNx+OLx
                0153         jMin = 1-OLy
                0154         jMax = sNy+OLy
4e66ab0b67 Oliv*0155 
f941e53e0a Jean*0156 C     need to recompute surfaceForcingPtr using LS_fwFlux
                0157         CALL LONGSTEP_FORCING_SURF(
                0158      I        bi, bj, iMin, iMax, jMin, jMax,
                0159      I        myTime,myIter,myThid )
4e66ab0b67 Oliv*0160 
f941e53e0a Jean*0161 C--   Set up 3-D velocity field that we use to advect tracers:
                0162 C-    just do a local copy:
                0163         DO k=1,Nr
                0164          DO j=1-OLy,sNy+OLy
                0165           DO i=1-OLx,sNx+OLx
                0166            uFld(i,j,k) = LS_uVel(i,j,k,bi,bj)
                0167            vFld(i,j,k) = LS_vVel(i,j,k,bi,bj)
                0168            wFld(i,j,k) = LS_wVel(i,j,k,bi,bj)
                0169           ENDDO
                0170          ENDDO
                0171         ENDDO
                0172 #ifdef ALLOW_GMREDI
                0173 C-    add Bolus velocity to Eulerian-mean velocity:
                0174         IF (useGMRedi) THEN
                0175           CALL  GMREDI_RESIDUAL_FLOW(
                0176      U                  uFld, vFld, wFld,
                0177      I                  bi, bj, myIter, myThid )
                0178         ENDIF
                0179 #endif /* ALLOW_GMREDI */
4e66ab0b67 Oliv*0180 
                0181 #ifdef ALLOW_PTRACERS
f941e53e0a Jean*0182 C--     Calculate passive tracer tendencies
                0183 C       and step forward, storing result in gPtr
                0184 C       Also apply open boundary conditions for each passive tracer
                0185         IF ( usePTRACERS ) THEN
                0186 #ifdef ALLOW_DEBUG
                0187            IF (debugMode) CALL DEBUG_CALL('PTRACERS_INTEGRATE',myThid)
                0188 #endif
                0189            CALL PTRACERS_INTEGRATE(
                0190      I          bi, bj, recip_hFacNew,
                0191      I          uFld, vFld, wFld,
                0192      U          kappaRk,
                0193      I          myTime, myIter, myThid )
4e66ab0b67 Oliv*0194         ENDIF
                0195 #endif /* ALLOW_PTRACERS */
                0196 
                0197 C--   end bi,bj loops.
                0198        ENDDO
                0199       ENDDO
                0200 
                0201 #ifdef ALLOW_DEBUG
8830b8f970 Jean*0202       IF ( debugLevel.GE.debLevD ) THEN
4e66ab0b67 Oliv*0203        CALL DEBUG_STATS_RL(Nr,LS_uVel,'LS_Uvel (THERMODYNAMICS)',myThid)
                0204        CALL DEBUG_STATS_RL(Nr,LS_vVel,'LS_Vvel (THERMODYNAMICS)',myThid)
                0205        CALL DEBUG_STATS_RL(Nr,LS_wVel,'LS_Wvel (THERMODYNAMICS)',myThid)
                0206        CALL DEBUG_STATS_RL(Nr,LS_theta,'LS_Theta (THERMODYNAMICS)',
                0207      &                     myThid)
                0208        CALL DEBUG_STATS_RL(Nr,LS_salt,'LS_Salt (THERMODYNAMICS)',myThid)
                0209 #ifdef ALLOW_PTRACERS
                0210        IF ( usePTRACERS ) THEN
                0211          CALL PTRACERS_DEBUG(myThid)
                0212        ENDIF
                0213 #endif /* ALLOW_PTRACERS */
                0214       ENDIF
                0215 #endif /* ALLOW_DEBUG */
                0216 
                0217 C     LS_doTimeStep
                0218       ENDIF
                0219 
                0220 #ifdef ALLOW_DEBUG
1378a26402 Jean*0221       IF (debugMode) CALL DEBUG_LEAVE('LONGSTEP_THERMODYNAMICS',myThid)
4e66ab0b67 Oliv*0222 #endif
                0223 
                0224 #endif /* ALLOW_LONGSTEP */
                0225 
                0226       RETURN
                0227       END