Back to home page

MITgcm

 
 

    


File indexing completed on 2025-07-08 05:11:10 UTC

view on githubraw file Latest commit 00c7090d on 2025-07-07 16:10:22 UTC
cf336ab6c5 Ryan*0001 #include "LAYERS_OPTIONS.h"
                0002 C--  File layers_thermodynamics.F:
                0003 C--   Contents
                0004 C--   o LAYERS_CALC_RHS
                0005 
                0006 CBOP 0
                0007 C     !ROUTINE: LAYERS_CALC_RHS
                0008 C     !INTERFACE:
                0009       SUBROUTINE LAYERS_CALC_RHS(
                0010      I                  myThid )
                0011 
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | SUBROUTINE LAYERS_CALC_RHS
                0015 C     | Recalculate the divergence of the RHS terms in T and S eqns.
                0016 C     | Replaces the values of layers_surfflux, layers_df? IN PLACE
                0017 C     | with the corresponding tendencies (same units as GT and GS)
                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 "GRID.h"
50d8304171 Ryan*0028 #include "FFIELDS.h"
cf336ab6c5 Ryan*0029 #include "LAYERS_SIZE.h"
                0030 #include "LAYERS.h"
                0031 
                0032 C !INPUT PARAMETERS:
                0033 C     myThid    :: my Thread Id number
                0034       INTEGER myThid
                0035 CEOP
                0036 
                0037 #ifdef ALLOW_LAYERS
ee16a2cae4 Ryan*0038 #ifdef LAYERS_THERMODYNAMICS
cf336ab6c5 Ryan*0039 C !LOCAL VARIABLES:
                0040 C     bi, bj   :: tile indices
                0041 C     i,j      :: horizontal indices
8d1543706e Jean*0042 C     k        :: vertical index for model grid
cf336ab6c5 Ryan*0043 C     kdown    :: temporary placeholder
                0044 C     fluxfac  :: scaling factor for converting surface flux to tendency
8d1543706e Jean*0045 C     fluxfac  :: scaling factor for converting diffusive flux to tendency
cf336ab6c5 Ryan*0046 C     downfac  :: mask for lower point
                0047 
                0048       INTEGER bi, bj
da9f56e003 Jean*0049       INTEGER i,j,k,kdown,iTracer
                0050       _RL fluxfac(2), downfac, tmpfac
                0051 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
50d8304171 Ryan*0052       _RL minusone
                0053       PARAMETER (minusOne=-1.)
cf336ab6c5 Ryan*0054 
                0055 C --  These factors convert the units of TFLUX and SFLUX diagnostics
                0056 C --  back to surfaceForcingT and surfaceForcingS units
                0057       fluxfac(1) = 1.0/(HeatCapacity_Cp*rUnit2mass)
                0058       fluxfac(2) = 1.0/rUnit2mass
                0059 
                0060         DO bj = myByLo(myThid), myByHi(myThid)
8d1543706e Jean*0061          DO bi = myBxLo(myThid), myBxHi(myThid)
da9f56e003 Jean*0062 
8d1543706e Jean*0063           DO iTracer = 1,2
da9f56e003 Jean*0064            k = 1
8d1543706e Jean*0065 C --       Loop for surface fluxes
cf336ab6c5 Ryan*0066            DO j=1-OLy,sNy+OLy
                0067             DO i=1-OLx,sNx+OLx
50d8304171 Ryan*0068 
                0069 #ifdef SHORTWAVE_HEATING
                0070 C --      Have to remove the shortwave from the surface flux because it is added later
                0071              IF (iTracer.EQ.1) THEN
                0072                layers_surfflux(i,j,k,iTracer,bi,bj) =
                0073      &           layers_surfflux(i,j,k,iTracer,bi,bj)
                0074 C --      Sign convention for Qsw means we have to add it to subtract it
                0075      &           +Qsw(i,j,bi,bj)
                0076              ENDIF
da9f56e003 Jean*0077 #endif /* SHORTWAVE_HEATING */
50d8304171 Ryan*0078 
cf336ab6c5 Ryan*0079              layers_surfflux(i,j,k,iTracer,bi,bj) =
                0080      &       layers_surfflux(i,j,k,iTracer,bi,bj)
8d1543706e Jean*0081      &       *recip_drF(1)*_recip_hFacC(i,j,1,bi,bj)
cf336ab6c5 Ryan*0082      &       *fluxfac(iTracer)
                0083             ENDDO
8d1543706e Jean*0084            ENDDO
da9f56e003 Jean*0085 
cf336ab6c5 Ryan*0086 C --       Loop for diffusive fluxes
                0087 C --       If done correctly, we can overwrite the flux array in place
                0088 C --       with its own divergence
                0089            DO k=1,Nr
00c7090dc0 Mart*0090             downFac = 1. _d 0
                0091 C     note: here kdown is for the mask below current level k
                0092             IF ( usingZCoords ) THEN
                0093              kdown = MIN(k+1,Nr)
                0094              IF ( k.EQ.Nr ) downFac = 0. _d 0
cf336ab6c5 Ryan*0095             ELSE
00c7090dc0 Mart*0096 C     this is the oceanic pressure coordinate case
                0097              kdown = MAX(k-1,1)
                0098              IF ( k.EQ.1  ) downFac = 0. _d 0
8d1543706e Jean*0099             ENDIF
50d8304171 Ryan*0100             DO j=1-OLy,sNy+OLy-1
                0101              DO i=1-OLx,sNx+OLx-1
                0102 C -- Diffusion
cf336ab6c5 Ryan*0103               tmpfac = -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0104      &          *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
                0105               layers_dfx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
8d1543706e Jean*0106      &         tmpfac * ( layers_dfx(i+1,j,k,iTracer,bi,bj) -
cf336ab6c5 Ryan*0107      &          layers_dfx(i,j,k,iTracer,bi,bj) )
                0108               layers_dfy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
8d1543706e Jean*0109      &         tmpfac * ( layers_dfy(i,j+1,k,iTracer,bi,bj) -
                0110      &          layers_dfy(i,j,k,iTracer,bi,bj) )
cf336ab6c5 Ryan*0111               layers_dfr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
                0112      &        ( layers_dfr(i,j,kdown,iTracer,bi,bj)*downfac -
                0113      &          layers_dfr(i,j,k,iTracer,bi,bj) )
50d8304171 Ryan*0114 C -- Advection
                0115               layers_afx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
                0116      &         tmpfac * ( layers_afx(i+1,j,k,iTracer,bi,bj) -
                0117      &          layers_afx(i,j,k,iTracer,bi,bj) )
                0118               layers_afy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
                0119      &         tmpfac * ( layers_afy(i,j+1,k,iTracer,bi,bj) -
                0120      &          layers_afy(i,j,k,iTracer,bi,bj) )
                0121               layers_afr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
                0122      &        ( layers_afr(i,j,kdown,iTracer,bi,bj)*downfac -
                0123      &          layers_afr(i,j,k,iTracer,bi,bj) )
                0124 
                0125 #ifdef SHORTWAVE_HEATING
da9f56e003 Jean*0126               IF (iTracer.EQ.1) THEN
50d8304171 Ryan*0127                 layers_sw(i,j,k,iTracer,bi,bj) =
                0128      &            layers_sw(i,j,k,iTracer,bi,bj)
00c7090dc0 Mart*0129      &            + Qsw(i,j,bi,bj)*gravitySign
                0130      &            *( SWFrac3D(i,j,k,bi,bj) - SWFrac3D(i,j,k+1,bi,bj) )
50d8304171 Ryan*0131      &            *fluxfac(1)
                0132      &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0133               ENDIF
da9f56e003 Jean*0134 #endif /* SHORTWAVE_HEATING */
                0135 
cf336ab6c5 Ryan*0136              ENDDO
8d1543706e Jean*0137             ENDDO
cf336ab6c5 Ryan*0138            ENDDO
                0139           ENDDO
                0140          ENDDO
                0141         ENDDO
                0142 
                0143 C-    TFLUX (=total heat flux, match heat-content variations, [W/m2])
                0144 C       IF ( fluidIsWater .AND.
                0145 C      &     DIAGNOSTICS_IS_ON('TFLUX   ',myThid) ) THEN
                0146 C        DO bj = myByLo(myThid), myByHi(myThid)
                0147 C         DO bi = myBxLo(myThid), myBxHi(myThid)
                0148 C          DO j = 1,sNy
                0149 C           DO i = 1,sNx
                0150 C            tmp1k(i,j,bi,bj) =
                0151 C #ifdef SHORTWAVE_HEATING
                0152 C      &      -Qsw(i,j,bi,bj)+
                0153 C #endif
                0154 C      &      (surfaceForcingT(i,j,bi,bj)+surfaceForcingTice(i,j,bi,bj))
                0155 C      &      *HeatCapacity_Cp*rUnit2mass
                0156 C           ENDDO
                0157 C          ENDDO
                0158 C #ifdef NONLIN_FRSURF
                0159 C          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0160 C      &        .AND. useRealFreshWaterFlux ) THEN
                0161 C           DO j=1,sNy
                0162 C            DO i=1,sNx
                0163 C             tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
                0164 C      &       + PmEpR(i,j,bi,bj)*theta(i,j,ks,bi,bj)*HeatCapacity_Cp
                0165 C            ENDDO
                0166 C           ENDDO
                0167 C          ENDIF
                0168 C #endif /* NONLIN_FRSURF */
                0169 C         ENDDO
                0170 C        ENDDO
                0171 C        CALL DIAGNOSTICS_FILL( tmp1k,'TFLUX   ',0,1,0,1,1,myThid )
                0172 C       ENDIF
8d1543706e Jean*0173 C
cf336ab6c5 Ryan*0174 C C-    SFLUX (=total salt flux, match salt-content variations [g/m2/s])
                0175 C       IF ( fluidIsWater .AND.
                0176 C      &     DIAGNOSTICS_IS_ON('SFLUX   ',myThid) ) THEN
                0177 C        DO bj = myByLo(myThid), myByHi(myThid)
                0178 C         DO bi = myBxLo(myThid), myBxHi(myThid)
                0179 C          DO j = 1,sNy
                0180 C           DO i = 1,sNx
                0181 C            tmp1k(i,j,bi,bj) =
                0182 C      &      surfaceForcingS(i,j,bi,bj)*rUnit2mass
                0183 C           ENDDO
                0184 C          ENDDO
8d1543706e Jean*0185 C
cf336ab6c5 Ryan*0186 C #ifdef NONLIN_FRSURF
                0187 C          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0188 C      &        .AND. useRealFreshWaterFlux ) THEN
                0189 C           DO j=1,sNy
                0190 C            DO i=1,sNx
                0191 C             tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
                0192 C      &       + PmEpR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
                0193 C            ENDDO
                0194 C           ENDDO
                0195 C          ENDIF
                0196 C #endif /* NONLIN_FRSURF */
8d1543706e Jean*0197 C
cf336ab6c5 Ryan*0198 C         ENDDO
                0199 C        ENDDO
                0200 C        CALL DIAGNOSTICS_FILL( tmp1k,'SFLUX   ',0,1,0,1,1,myThid )
                0201 C       ENDIF
                0202 
                0203 C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
                0204 C      IF ( kLev .EQ. kSurface ) THEN
                0205 C       DO j=1,sNy
                0206 C        DO i=1,sNx
                0207 C          gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
                0208 C     &      +surfaceForcingT(i,j,bi,bj)
                0209 C     &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
                0210 C        ENDDO
                0211 C       ENDDO
                0212 C      ELSEIF ( kSurface.EQ.-1 ) THEN
                0213 C       DO j=1,sNy
                0214 C        DO i=1,sNx
                0215 C         IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
                0216 C          gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
                0217 C     &      +surfaceForcingT(i,j,bi,bj)
                0218 C     &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
                0219 C         ENDIF
                0220 C        ENDDO
                0221 C       ENDDO
                0222 C      ENDIF
                0223 
                0224 C--   Divergence of fluxes
                0225 C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
                0226 C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
                0227 C     DO j=1-OLy,sNy+OLy-1
                0228 C      DO i=1-OLx,sNx+OLx-1
                0229 C       gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
                0230 C    &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0231 C    &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
                0232 C    &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
                0233 C    &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
                0234 C    &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
                0235 C    &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
                0236 C    &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
                0237 C    &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
                0238 C    &                  )*maskInC(i,j,bi,bj)
                0239 C    &    )
                0240 C      ENDDO
8d1543706e Jean*0241 C     ENDDO
cf336ab6c5 Ryan*0242 
                0243 #endif /* LAYERS_THERMODYNAMICS */
ee16a2cae4 Ryan*0244 #endif /* USE_LAYERS */
cf336ab6c5 Ryan*0245       RETURN
                0246       END
                0247 C -- end of S/R LAYERS_CALC_RHS