Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:41:46 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 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.)
da9f56e003 Jean*0054 #ifdef SHORTWAVE_HEATING
                0055       _RL swfracb(2)
                0056 #endif
cf336ab6c5 Ryan*0057 
                0058 C --  These factors convert the units of TFLUX and SFLUX diagnostics
                0059 C --  back to surfaceForcingT and surfaceForcingS units
                0060       fluxfac(1) = 1.0/(HeatCapacity_Cp*rUnit2mass)
                0061       fluxfac(2) = 1.0/rUnit2mass
                0062 
                0063         DO bj = myByLo(myThid), myByHi(myThid)
8d1543706e Jean*0064          DO bi = myBxLo(myThid), myBxHi(myThid)
da9f56e003 Jean*0065 
8d1543706e Jean*0066           DO iTracer = 1,2
da9f56e003 Jean*0067            k = 1
8d1543706e Jean*0068 C --       Loop for surface fluxes
cf336ab6c5 Ryan*0069            DO j=1-OLy,sNy+OLy
                0070             DO i=1-OLx,sNx+OLx
50d8304171 Ryan*0071 
                0072 #ifdef SHORTWAVE_HEATING
                0073 C --      Have to remove the shortwave from the surface flux because it is added later
                0074              IF (iTracer.EQ.1) THEN
                0075                layers_surfflux(i,j,k,iTracer,bi,bj) =
                0076      &           layers_surfflux(i,j,k,iTracer,bi,bj)
                0077 C --      Sign convention for Qsw means we have to add it to subtract it
                0078      &           +Qsw(i,j,bi,bj)
                0079              ENDIF
da9f56e003 Jean*0080 #endif /* SHORTWAVE_HEATING */
50d8304171 Ryan*0081 
cf336ab6c5 Ryan*0082              layers_surfflux(i,j,k,iTracer,bi,bj) =
                0083      &       layers_surfflux(i,j,k,iTracer,bi,bj)
8d1543706e Jean*0084      &       *recip_drF(1)*_recip_hFacC(i,j,1,bi,bj)
cf336ab6c5 Ryan*0085      &       *fluxfac(iTracer)
                0086             ENDDO
8d1543706e Jean*0087            ENDDO
da9f56e003 Jean*0088 
cf336ab6c5 Ryan*0089 C --       Loop for diffusive fluxes
                0090 C --       If done correctly, we can overwrite the flux array in place
                0091 C --       with its own divergence
                0092            DO k=1,Nr
                0093             kdown= MIN(k+1,Nr)
                0094             IF (k.EQ.Nr) THEN
                0095               downfac = 0. _d 0
                0096             ELSE
                0097               downfac = 1. _d 0
8d1543706e Jean*0098             ENDIF
50d8304171 Ryan*0099             DO j=1-OLy,sNy+OLy-1
                0100              DO i=1-OLx,sNx+OLx-1
                0101 C -- Diffusion
cf336ab6c5 Ryan*0102               tmpfac = -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0103      &          *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
                0104               layers_dfx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
8d1543706e Jean*0105      &         tmpfac * ( layers_dfx(i+1,j,k,iTracer,bi,bj) -
cf336ab6c5 Ryan*0106      &          layers_dfx(i,j,k,iTracer,bi,bj) )
                0107               layers_dfy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
8d1543706e Jean*0108      &         tmpfac * ( layers_dfy(i,j+1,k,iTracer,bi,bj) -
                0109      &          layers_dfy(i,j,k,iTracer,bi,bj) )
cf336ab6c5 Ryan*0110               layers_dfr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
                0111      &        ( layers_dfr(i,j,kdown,iTracer,bi,bj)*downfac -
                0112      &          layers_dfr(i,j,k,iTracer,bi,bj) )
50d8304171 Ryan*0113 C -- Advection
                0114               layers_afx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
                0115      &         tmpfac * ( layers_afx(i+1,j,k,iTracer,bi,bj) -
                0116      &          layers_afx(i,j,k,iTracer,bi,bj) )
                0117               layers_afy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
                0118      &         tmpfac * ( layers_afy(i,j+1,k,iTracer,bi,bj) -
                0119      &          layers_afy(i,j,k,iTracer,bi,bj) )
                0120               layers_afr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
                0121      &        ( layers_afr(i,j,kdown,iTracer,bi,bj)*downfac -
                0122      &          layers_afr(i,j,k,iTracer,bi,bj) )
                0123 
                0124 #ifdef SHORTWAVE_HEATING
da9f56e003 Jean*0125               IF (iTracer.EQ.1) THEN
50d8304171 Ryan*0126                 swfracb(1)=abs(rF(k))
                0127                 swfracb(2)=abs(rF(k+1))
                0128                 CALL SWFRAC(
                0129      I             2, minusOne,
                0130      U             swfracb,
                0131      I             1.0, 1, myThid )
                0132 C  ----- debuggin
                0133 C                IF ((i.EQ.0).AND.(j.EQ.0)) THEN
                0134 C                  WRITE(msgBuf,'(2A,I3,A,F6.2,A,F6.2)')
                0135 C     &             'S/R LAYERS_THERMODYNAMICS:',
                0136 C     &             ' k=', k,
                0137 C     &             ' swfracb(1)=', swfracb(1),
                0138 C     &             ' swfracb(2)=', swfracb(2)
                0139 C                  CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0140 C     &                         SQUEEZE_RIGHT, myThid )
                0141 C                ENDIF
                0142 C            --- kdown == kp1
                0143 C                kp1 = k+1
                0144                 IF (k.EQ.Nr) THEN
                0145 C                  kp1 = k
                0146                   swfracb(2)=0. _d 0
                0147                 ENDIF
                0148                 layers_sw(i,j,k,iTracer,bi,bj) =
                0149      &            layers_sw(i,j,k,iTracer,bi,bj)
                0150      &            -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
                0151      &                      -swfracb(2)*maskC(i,j,kdown,bi,bj))
                0152      &            *fluxfac(1)
                0153      &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0154               ENDIF
da9f56e003 Jean*0155 #endif /* SHORTWAVE_HEATING */
                0156 
cf336ab6c5 Ryan*0157              ENDDO
8d1543706e Jean*0158             ENDDO
cf336ab6c5 Ryan*0159            ENDDO
                0160           ENDDO
                0161          ENDDO
                0162         ENDDO
                0163 
                0164 C-    TFLUX (=total heat flux, match heat-content variations, [W/m2])
                0165 C       IF ( fluidIsWater .AND.
                0166 C      &     DIAGNOSTICS_IS_ON('TFLUX   ',myThid) ) THEN
                0167 C        DO bj = myByLo(myThid), myByHi(myThid)
                0168 C         DO bi = myBxLo(myThid), myBxHi(myThid)
                0169 C          DO j = 1,sNy
                0170 C           DO i = 1,sNx
                0171 C            tmp1k(i,j,bi,bj) =
                0172 C #ifdef SHORTWAVE_HEATING
                0173 C      &      -Qsw(i,j,bi,bj)+
                0174 C #endif
                0175 C      &      (surfaceForcingT(i,j,bi,bj)+surfaceForcingTice(i,j,bi,bj))
                0176 C      &      *HeatCapacity_Cp*rUnit2mass
                0177 C           ENDDO
                0178 C          ENDDO
                0179 C #ifdef NONLIN_FRSURF
                0180 C          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0181 C      &        .AND. useRealFreshWaterFlux ) THEN
                0182 C           DO j=1,sNy
                0183 C            DO i=1,sNx
                0184 C             tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
                0185 C      &       + PmEpR(i,j,bi,bj)*theta(i,j,ks,bi,bj)*HeatCapacity_Cp
                0186 C            ENDDO
                0187 C           ENDDO
                0188 C          ENDIF
                0189 C #endif /* NONLIN_FRSURF */
                0190 C         ENDDO
                0191 C        ENDDO
                0192 C        CALL DIAGNOSTICS_FILL( tmp1k,'TFLUX   ',0,1,0,1,1,myThid )
                0193 C       ENDIF
8d1543706e Jean*0194 C
cf336ab6c5 Ryan*0195 C C-    SFLUX (=total salt flux, match salt-content variations [g/m2/s])
                0196 C       IF ( fluidIsWater .AND.
                0197 C      &     DIAGNOSTICS_IS_ON('SFLUX   ',myThid) ) THEN
                0198 C        DO bj = myByLo(myThid), myByHi(myThid)
                0199 C         DO bi = myBxLo(myThid), myBxHi(myThid)
                0200 C          DO j = 1,sNy
                0201 C           DO i = 1,sNx
                0202 C            tmp1k(i,j,bi,bj) =
                0203 C      &      surfaceForcingS(i,j,bi,bj)*rUnit2mass
                0204 C           ENDDO
                0205 C          ENDDO
8d1543706e Jean*0206 C
cf336ab6c5 Ryan*0207 C #ifdef NONLIN_FRSURF
                0208 C          IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0209 C      &        .AND. useRealFreshWaterFlux ) THEN
                0210 C           DO j=1,sNy
                0211 C            DO i=1,sNx
                0212 C             tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
                0213 C      &       + PmEpR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
                0214 C            ENDDO
                0215 C           ENDDO
                0216 C          ENDIF
                0217 C #endif /* NONLIN_FRSURF */
8d1543706e Jean*0218 C
cf336ab6c5 Ryan*0219 C         ENDDO
                0220 C        ENDDO
                0221 C        CALL DIAGNOSTICS_FILL( tmp1k,'SFLUX   ',0,1,0,1,1,myThid )
                0222 C       ENDIF
                0223 
                0224 C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
                0225 C      IF ( kLev .EQ. kSurface ) THEN
                0226 C       DO j=1,sNy
                0227 C        DO i=1,sNx
                0228 C          gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
                0229 C     &      +surfaceForcingT(i,j,bi,bj)
                0230 C     &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
                0231 C        ENDDO
                0232 C       ENDDO
                0233 C      ELSEIF ( kSurface.EQ.-1 ) THEN
                0234 C       DO j=1,sNy
                0235 C        DO i=1,sNx
                0236 C         IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
                0237 C          gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
                0238 C     &      +surfaceForcingT(i,j,bi,bj)
                0239 C     &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
                0240 C         ENDIF
                0241 C        ENDDO
                0242 C       ENDDO
                0243 C      ENDIF
                0244 
                0245 C--   Divergence of fluxes
                0246 C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
                0247 C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
                0248 C     DO j=1-OLy,sNy+OLy-1
                0249 C      DO i=1-OLx,sNx+OLx-1
                0250 C       gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
                0251 C    &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0252 C    &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
                0253 C    &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
                0254 C    &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
                0255 C    &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
                0256 C    &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
                0257 C    &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
                0258 C    &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
                0259 C    &                  )*maskInC(i,j,bi,bj)
                0260 C    &    )
                0261 C      ENDDO
8d1543706e Jean*0262 C     ENDDO
cf336ab6c5 Ryan*0263 
                0264 #endif /* LAYERS_THERMODYNAMICS */
ee16a2cae4 Ryan*0265 #endif /* USE_LAYERS */
cf336ab6c5 Ryan*0266       RETURN
                0267       END
                0268 C -- end of S/R LAYERS_CALC_RHS
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!