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
0005
0006
0007 SUBROUTINE FORCING_SURF_RELAX(
0008 I iMin, iMax, jMin, jMax,
0009 I myTime, myIter, myThid )
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 IMPLICIT NONE
0020
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
0035
0036
0037
0038
0039
0040 INTEGER iMin, iMax
0041 INTEGER jMin, jMax
0042 _RL myTime
0043 INTEGER myIter
0044 INTEGER myThid
0045
0046
0047
0048
0049
0050
0051 INTEGER bi,bj
0052 INTEGER i,j
0053 INTEGER ks
0054
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
0073
0074 #ifdef ALLOW_SEAICE
0075 IF ( useSEAICE .AND. (.NOT. SEAICErestoreUnderIce) ) THEN
0076
0077 DO j = jMin, jMax
0078 DO i = iMin, iMax
0079
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
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
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
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
0112 #ifdef NONLIN_FRSURF
0113
0114
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
0160 ENDDO
0161 ENDDO
0162
0163
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
0245 tmpFac = HeatCapacity_Cp*rUnit2mass
0246 CALL DIAGNOSTICS_SCALE_FILL( surfaceForcingT, tmpFac, 1,
0247 & 'TRELAX ', 0, 1, 0,1,1, myThid )
0248
0249
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