Back to home page

MITgcm

 
 

    


File indexing completed on 2021-08-12 05:10:57 UTC

view on githubraw file Latest commit 0320e252 on 2021-08-11 16:08:52 UTC
afdb812123 Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003 
                0004 CBOP
                0005 C     !ROUTINE: FORCING_SURF_RELAX
                0006 C     !INTERFACE:
                0007       SUBROUTINE FORCING_SURF_RELAX(
                0008      I                   iMin, iMax, jMin, jMax,
                0009      I                   myTime, myIter, myThid )
                0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE FORCING_SURF_RELAX
                0013 C     | o Calculate relaxation surface forcing terms
                0014 C     |   for temperature and salinity
                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 "FFIELDS.h"
                0025 #include "DYNVARS.h"
                0026 #include "GRID.h"
                0027 #include "SURFACE.h"
                0028 #ifdef ALLOW_SEAICE
                0029 # include "SEAICE_SIZE.h"
                0030 # include "SEAICE_PARAMS.h"
                0031 # include "SEAICE.h"
                0032 #endif /* ALLOW_SEAICE */
                0033 
                0034 C     !INPUT/OUTPUT PARAMETERS:
                0035 C     === Routine arguments ===
                0036 C     iMin,iMax, jMin,jMax :: Range of points for calculation
                0037 C     myTime  :: Current time in simulation
                0038 C     myIter  :: Current iteration number in simulation
                0039 C     myThid  :: Thread no. that called this routine.
                0040       INTEGER iMin, iMax
                0041       INTEGER jMin, jMax
                0042       _RL myTime
                0043       INTEGER myIter
                0044       INTEGER myThid
                0045 
                0046 C     !LOCAL VARIABLES:
                0047 C     === Local variables ===
                0048 C     bi,bj  :: tile indices
                0049 C     i,j    :: loop indices
                0050 C     ks     :: index of surface interface layer
                0051       INTEGER bi,bj
                0052       INTEGER i,j
                0053       INTEGER ks
                0054 CEOP
                0055 #ifdef ALLOW_DIAGNOSTICS
                0056       _RL tmpFac
                0057 #endif /* ALLOW_DIAGNOSTICS */
                0058 #ifdef ALLOW_BALANCE_RELAX
                0059       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0060       _RL sumTile(nSx,nSy), sumGlob, globAver
                0061 #endif /* ALLOW_BALANCE_RELAX */
                0062 
                0063       IF ( usingPCoords ) THEN
                0064        ks        = Nr
                0065       ELSE
                0066        ks        = 1
                0067       ENDIF
                0068 
                0069       DO bj=myByLo(myThid),myByHi(myThid)
                0070        DO bi=myBxLo(myThid),myBxHi(myThid)
                0071 
                0072 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0073 
                0074 #ifdef ALLOW_SEAICE
                0075        IF ( useSEAICE .AND. (.NOT. SEAICErestoreUnderIce) ) THEN
                0076 C     Do not restore under sea-ice
                0077         DO j = jMin, jMax
                0078          DO i = iMin, iMax
                0079 C     Heat Flux (restoring term) :
                0080           surfaceForcingT(i,j,bi,bj) =
                0081      &      -lambdaThetaClimRelax(i,j,bi,bj)*(1.-AREA(i,j,bi,bj))
                0082      &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
                0083      &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
                0084 C     Salt Flux (restoring term) :
                0085           surfaceForcingS(i,j,bi,bj) =
                0086      &      -lambdaSaltClimRelax(i,j,bi,bj) *(1.-AREA(i,j,bi,bj))
                0087      &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
                0088      &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
                0089          ENDDO
                0090         ENDDO
                0091        ELSE
                0092 #endif /* ALLOW_SEAICE */
                0093         DO j = jMin, jMax
                0094          DO i = iMin, iMax
                0095 C     Heat Flux (restoring term) :
                0096           surfaceForcingT(i,j,bi,bj) =
                0097      &      -lambdaThetaClimRelax(i,j,bi,bj)
                0098      &         *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
                0099      &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
                0100 C     Salt Flux (restoring term) :
                0101           surfaceForcingS(i,j,bi,bj) =
                0102      &      -lambdaSaltClimRelax(i,j,bi,bj)
                0103      &         *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
                0104      &         *drF(ks)*_hFacC(i,j,ks,bi,bj)
                0105          ENDDO
                0106         ENDDO
                0107 #ifdef ALLOW_SEAICE
                0108        ENDIF
                0109 #endif /* ALLOW_SEAICE */
                0110 
                0111 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0112 #ifdef NONLIN_FRSURF
                0113 C-    T,S surface forcing will be applied (thermodynamics) after the update
                0114 C     of surf.thickness (hFac): account for change in surf.thickness
                0115        IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
                0116         IF ( select_rStar.GT.0 ) THEN
                0117 # ifndef DISABLE_RSTAR_CODE
                0118          DO j=jMin,jMax
                0119           DO i=iMin,iMax
                0120             surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
                0121      &                                  * rStarExpC(i,j,bi,bj)
                0122             surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
                0123      &                                  * rStarExpC(i,j,bi,bj)
                0124           ENDDO
                0125          ENDDO
                0126 # endif /* DISABLE_RSTAR_CODE */
                0127         ELSEIF ( selectSigmaCoord.NE.0 ) THEN
                0128 # ifndef DISABLE_SIGMA_CODE
                0129          DO j=jMin,jMax
                0130           DO i=iMin,iMax
                0131             surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
                0132      &        *( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
                0133      &                    *dBHybSigF(ks)*recip_drF(ks)
                0134      &                    *recip_hFacC(i,j,ks,bi,bj)
                0135      &         )
                0136             surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
                0137      &        *( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
                0138      &                    *dBHybSigF(ks)*recip_drF(ks)
                0139      &                    *recip_hFacC(i,j,ks,bi,bj)
                0140      &         )
                0141           ENDDO
                0142          ENDDO
                0143 # endif /* DISABLE_SIGMA_CODE */
                0144         ELSE
                0145          DO j=jMin,jMax
                0146           DO i=iMin,iMax
                0147            IF (ks.EQ.kSurfC(i,j,bi,bj)) THEN
                0148             surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
                0149      &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
                0150             surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
                0151      &             *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
                0152            ENDIF
                0153           ENDDO
                0154          ENDDO
                0155         ENDIF
                0156        ENDIF
                0157 #endif /* NONLIN_FRSURF */
                0158 
                0159 C--   end bi,bj loops.
                0160        ENDDO
                0161       ENDDO
                0162 
                0163 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0164 
                0165 #ifdef ALLOW_BALANCE_RELAX
                0166 
                0167       IF ( balanceThetaClimRelax ) THEN
                0168         DO bj=myByLo(myThid),myByHi(myThid)
                0169          DO bi=myBxLo(myThid),myBxHi(myThid)
                0170           sumTile(bi,bj) = 0. _d 0
                0171           DO j=1,sNy
                0172            DO i=1,sNx
                0173              sumTile(bi,bj) = sumTile(bi,bj)
                0174      &                      + surfaceForcingT(i,j,bi,bj)
                0175      &                       *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0176            ENDDO
                0177           ENDDO
                0178          ENDDO
                0179         ENDDO
                0180         CALL GLOBAL_SUM_TILE_RL( sumTile, sumGlob, myThid )
                0181         globAver = sumGlob
                0182         IF ( globalArea.GT.zeroRL ) globAver = globAver / globalArea
                0183         DO bj=myByLo(myThid),myByHi(myThid)
                0184          DO bi=myBxLo(myThid),myBxHi(myThid)
                0185           DO j=1-OLy,sNy+OLy
                0186            DO i=1-OLx,sNx+OLx
                0187              surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
                0188      &                                  - globAver
                0189            ENDDO
                0190           ENDDO
                0191          ENDDO
                0192         ENDDO
                0193         IF ( balancePrintMean ) THEN
                0194          _BEGIN_MASTER( myThid )
0320e25227 Mart*0195          WRITE(msgBuf,'(2A,1PE21.14,A,I10)') 'rm Global mean of',
                0196      &                ' SSTrelax= ', globAver, '  @ it=', myIter
afdb812123 Jean*0197          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0198      &                       SQUEEZE_RIGHT, myThid )
                0199          _END_MASTER( myThid )
                0200         ENDIF
                0201       ENDIF
                0202 
                0203       IF ( balanceSaltClimRelax ) THEN
                0204         DO bj=myByLo(myThid),myByHi(myThid)
                0205          DO bi=myBxLo(myThid),myBxHi(myThid)
                0206           sumTile(bi,bj) = 0. _d 0
                0207           DO j=1,sNy
                0208            DO i=1,sNx
                0209              sumTile(bi,bj) = sumTile(bi,bj)
                0210      &                      + surfaceForcingS(i,j,bi,bj)
                0211      &                       *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0212            ENDDO
                0213           ENDDO
                0214          ENDDO
                0215         ENDDO
                0216         CALL GLOBAL_SUM_TILE_RL( sumTile, sumGlob, myThid )
                0217         globAver = sumGlob
                0218         IF ( globalArea.GT.zeroRL ) globAver = globAver / globalArea
                0219         DO bj=myByLo(myThid),myByHi(myThid)
                0220          DO bi=myBxLo(myThid),myBxHi(myThid)
                0221           DO j=1-OLy,sNy+OLy
                0222            DO i=1-OLx,sNx+OLx
                0223              surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
                0224      &                                  - globAver
                0225            ENDDO
                0226           ENDDO
                0227          ENDDO
                0228         ENDDO
                0229         IF ( balancePrintMean ) THEN
                0230          _BEGIN_MASTER( myThid )
0320e25227 Mart*0231          WRITE(msgBuf,'(2A,1PE21.14,A,I10)') 'rm Global mean of',
                0232      &                ' SSSrelax= ', globAver, '  @ it=', myIter
afdb812123 Jean*0233          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0234      &                       SQUEEZE_RIGHT, myThid )
                0235          _END_MASTER( myThid )
                0236         ENDIF
                0237       ENDIF
                0238 
                0239 #endif /* ALLOW_BALANCE_RELAX */
                0240 
                0241 #ifdef ALLOW_DIAGNOSTICS
                0242       IF ( useDiagnostics ) THEN
                0243 
                0244 C     tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
                0245         tmpFac = HeatCapacity_Cp*rUnit2mass
                0246         CALL DIAGNOSTICS_SCALE_FILL( surfaceForcingT, tmpFac, 1,
                0247      &                              'TRELAX  ', 0, 1, 0,1,1, myThid )
                0248 
                0249 C     sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
                0250         tmpFac = rUnit2mass
                0251         CALL DIAGNOSTICS_SCALE_FILL( surfaceForcingS, tmpFac, 1,
                0252      &                              'SRELAX  ', 0, 1, 0,1,1, myThid )
                0253 
                0254       ENDIF
                0255 #endif /* ALLOW_DIAGNOSTICS */
                0256 
                0257       RETURN
                0258       END