pkg/layers/layers_calc_divergence.F

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