File indexing completed on 2024-04-23 05:10:04 UTC
view on githubraw file Latest commit 8fabf868 on 2024-04-22 22:38:11 UTC
9bf3441755 Jean*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
517dbdc414 Jean*0003 #ifdef ALLOW_AUTODIFF
0004 # include "AUTODIFF_OPTIONS.h"
0005 #endif
924557e60a Chri*0006
9366854e02 Chri*0007
0008
0009
fb481a83c2 Alis*0010 SUBROUTINE CALC_PHI_HYD(
f1f873e18c Jean*0011 I bi, bj, iMin, iMax, jMin, jMax, k,
0012 U phiHydF,
0013 O phiHydC, dPhiHydX, dPhiHydY,
842349fcb7 Jean*0014 I myTime, myIter, myThid )
9366854e02 Chri*0015
0016
11734d097d Chri*0017
0f5ba23f97 Jean*0018
9366854e02 Chri*0019
f1f873e18c Jean*0020
0021
0022
0f5ba23f97 Jean*0023
f1f873e18c Jean*0024
0025
0026
0027
0f5ba23f97 Jean*0028
f1f873e18c Jean*0029
0030
0031
0032
9366854e02 Chri*0033
0034
0035
924557e60a Chri*0036 IMPLICIT NONE
0037
0038 #include "SIZE.h"
0039 #include "GRID.h"
81bc00c2f0 Chri*0040 #include "EEPARAMS.h"
924557e60a Chri*0041 #include "PARAMS.h"
7c50f07931 Mart*0042 #ifdef ALLOW_AUTODIFF_TAMC
612f1eab29 Patr*0043 #include "tamc.h"
7c50f07931 Mart*0044 #endif /* ALLOW_AUTODIFF_TAMC */
9e5f2093d3 Alis*0045 #include "SURFACE.h"
60c223928f Mart*0046 #include "DYNVARS.h"
612f1eab29 Patr*0047
9366854e02 Chri*0048
924557e60a Chri*0049
0f5ba23f97 Jean*0050
f1f873e18c Jean*0051
0052
0053
0054
0055
0056
0057
0058
0059 INTEGER bi,bj,iMin,iMax,jMin,jMax,k
0060 _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0061 _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
2b8326b7dd Jean*0062 _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0063 _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7e1bd93d4f Jean*0064 _RL myTime
0065 INTEGER myIter, myThid
0f5ba23f97 Jean*0066
1dbaea09ee Chri*0067 #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
0068
9366854e02 Chri*0069
fb481a83c2 Alis*0070
3d9080e0a5 Jean*0071
0072
0073
0074
f1f873e18c Jean*0075 INTEGER i,j
fb481a83c2 Alis*0076 _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
3d9080e0a5 Jean*0077 #ifndef DISABLE_SIGMA_CODE
0078 _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0079 _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0080 _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0081 _RL pKappaC, locDepth, fullDepth
0082 #endif /* DISABLE_SIGMA_CODE */
8fabf86828 Jean*0083 _RL thetaRef, locBuoy
cfbe180681 Jean*0084 _RL dRlocM,dRlocP, ddRloc
f1f873e18c Jean*0085 _RL ddPIm, ddPIp, rec_dRm, rec_dRp
0086 _RL surfPhiFac
0087 LOGICAL useDiagPhiRlow, addSurfPhiAnom
3d9080e0a5 Jean*0088 LOGICAL useFVgradPhi
d9efc637ba Mart*0089 #ifdef ALLOW_AUTODIFF_TAMC
0090
0091
0092 INTEGER tkey, kkey
0093 #endif
9366854e02 Chri*0094
1adfdf6281 Jean*0095 useDiagPhiRlow = .TRUE.
2b8326b7dd Jean*0096 addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
3d9080e0a5 Jean*0097 useFVgradPhi = selectSigmaCoord.NE.0
0098
f1f873e18c Jean*0099 surfPhiFac = 0.
0100 IF (addSurfPhiAnom) surfPhiFac = 1.
45adaae56c Jean*0101
0102
0f5ba23f97 Jean*0103
463053c692 Jean*0104
f1f873e18c Jean*0105
0106
0f5ba23f97 Jean*0107
f1f873e18c Jean*0108
0f5ba23f97 Jean*0109
f1f873e18c Jean*0110
45adaae56c Jean*0111
fb481a83c2 Alis*0112
612f1eab29 Patr*0113 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0114 tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
7448700841 Mart*0115 kkey = k + (tkey-1)*Nr
612f1eab29 Patr*0116 #endif /* ALLOW_AUTODIFF_TAMC */
0117
0f5ba23f97 Jean*0118
f1f873e18c Jean*0119
0120
0121 IF (k.EQ.1) THEN
2b8326b7dd Jean*0122 DO j=1-OLy,sNy+OLy
0123 DO i=1-OLx,sNx+OLx
f1f873e18c Jean*0124 phiHydF(i,j) = 0.
0125 ENDDO
0126 ENDDO
0127 ENDIF
7e1bd93d4f Jean*0128
0129
f1f873e18c Jean*0130 IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
fb481a83c2 Alis*0131
0132
0133
842349fcb7 Jean*0134 #ifdef ALLOW_AUTODIFF_TAMC
7e1bd93d4f Jean*0135
842349fcb7 Jean*0136 #endif /* ALLOW_AUTODIFF_TAMC */
fb481a83c2 Alis*0137
842349fcb7 Jean*0138 IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
f1f873e18c Jean*0139
612f1eab29 Patr*0140 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0141
0142
612f1eab29 Patr*0143 #endif /* ALLOW_AUTODIFF_TAMC */
842349fcb7 Jean*0144 CALL FIND_RHO_2D(
0145 I iMin, iMax, jMin, jMax, k,
b98355c8c5 jm-c 0146 I theta(1-OLx,1-OLy,k,bi,bj),
0147 I salt(1-OLx,1-OLy,k,bi,bj),
842349fcb7 Jean*0148 O alphaRho,
0149 I k, bi, bj, myThid )
0150 ELSE
0151 DO j=jMin,jMax
0152 DO i=iMin,iMax
0153 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
0154 ENDDO
0155 ENDDO
0156 ENDIF
0f5ba23f97 Jean*0157
a6cbc7a360 Mart*0158 #ifdef ALLOW_SHELFICE
0f5ba23f97 Jean*0159
a6cbc7a360 Mart*0160
807d4a3bf5 Jean*0161 IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
0162
0163
0164
0165 DO j=jMin,jMax
0166 DO i=iMin,iMax
0167 IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
0168 ENDDO
0169 ENDDO
0170 ELSEIF ( useShelfIce ) THEN
a6cbc7a360 Mart*0171 DO j=jMin,jMax
0172 DO i=iMin,iMax
0173 alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
0174 ENDDO
0175 ENDDO
0176 ENDIF
0177 #endif /* ALLOW_SHELFICE */
fb481a83c2 Alis*0178
af53f2701c Jean*0179 #ifdef ALLOW_MOM_COMMON
0d373cd03b Jean*0180
ee19c4a296 Jean*0181 IF ( quasiHydrostatic ) THEN
0182 CALL MOM_QUASIHYDROSTATIC( bi, bj, k, uVel, vVel,
0183 U alphaRho,
0184 I myTime, myIter, myThid )
775ad00e06 Alis*0185 ENDIF
af53f2701c Jean*0186 #endif /* ALLOW_MOM_COMMON */
7448700841 Mart*0187 #ifdef ALLOW_AUTODIFF_TAMC
0188
0189 #endif
775ad00e06 Alis*0190
f1f873e18c Jean*0191 #ifdef NONLIN_FRSURF
2b8326b7dd Jean*0192 IF ( addSurfPhiAnom .AND.
0193 & uniformFreeSurfLev .AND. k.EQ.1 ) THEN
f1f873e18c Jean*0194 DO j=jMin,jMax
0195 DO i=iMin,iMax
0196 phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
0197 & *gravity*alphaRho(i,j)*recip_rhoConst
0198 ENDDO
0199 ENDDO
0200 ENDIF
0201 #endif /* NONLIN_FRSURF */
1adfdf6281 Jean*0202
f1f873e18c Jean*0203
7e1bd93d4f Jean*0204
0205 IF (integr_GeoPot.EQ.1) THEN
0206
0207
0208
0209
0210
2b8326b7dd Jean*0211
0212 IF ( uniformFreeSurfLev ) THEN
0213 DO j=jMin,jMax
0214 DO i=iMin,iMax
0215 phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0216 & + halfRL*drF(k)*gravFacC(k)*gravity
0217 & *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0218 phiHydF(i,j) = phiHydF(i,j)
bd053599dc Jean*0219 & + drF(k)*gravFacC(k)*gravity
0220 & *alphaRho(i,j)*recip_rhoConst
7e1bd93d4f Jean*0221 ENDDO
0222 ENDDO
2b8326b7dd Jean*0223 ELSE
0224 DO j=jMin,jMax
0225 DO i=iMin,iMax
0226 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
0227 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
0228 #ifdef NONLIN_FRSURF
0229 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
0230 #endif
bd053599dc Jean*0231 phiHydC(i,j) = ddRloc*gravFacC(k)*gravity
0232 & *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0233 ELSE
0234 phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0235 & + halfRL*drF(k)*gravFacC(k)*gravity
0236 & *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0237 ENDIF
bd053599dc Jean*0238 phiHydF(i,j) = phiHydC(i,j)
0239 & + halfRL*drF(k)*gravFacC(k)*gravity
0240 & *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0241 ENDDO
0242 ENDDO
0243 ENDIF
7e1bd93d4f Jean*0244
0245 ELSE
0246
0247
2b8326b7dd Jean*0248
0249
f1f873e18c Jean*0250
bd053599dc Jean*0251 dRlocM = halfRL*drC(k)*gravFacF(k)
0252 IF (k.EQ.1) dRlocM = (rF(k)-rC(k))*gravFacF(k)
2b8326b7dd Jean*0253 IF (k.EQ.Nr) THEN
bd053599dc Jean*0254 dRlocP = (rC(k)-rF(k+1))*gravFacF(k+1)
2b8326b7dd Jean*0255 ELSE
bd053599dc Jean*0256 dRlocP = halfRL*drC(k+1)*gravFacF(k+1)
2b8326b7dd Jean*0257 ENDIF
0258 IF ( uniformFreeSurfLev ) THEN
7e1bd93d4f Jean*0259 DO j=jMin,jMax
0260 DO i=iMin,iMax
2b8326b7dd Jean*0261 phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0262 & + dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0263 phiHydF(i,j) = phiHydC(i,j)
bd053599dc Jean*0264 & + dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
fb481a83c2 Alis*0265 ENDDO
7e1bd93d4f Jean*0266 ENDDO
2b8326b7dd Jean*0267 ELSE
0d373cd03b Jean*0268 rec_dRm = oneRL/(rF(k)-rC(k))
0269 rec_dRp = oneRL/(rC(k)-rF(k+1))
2b8326b7dd Jean*0270 DO j=jMin,jMax
0271 DO i=iMin,iMax
0272 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
0273 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
0274 #ifdef NONLIN_FRSURF
0275 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
0276 #endif
0d373cd03b Jean*0277 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
0278 & +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
0279 & )*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0280 ELSE
0281 phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0282 & + dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0283 ENDIF
bd053599dc Jean*0284 phiHydF(i,j) = phiHydC(i,j)
0285 & + dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0286 ENDDO
0287 ENDDO
0288 ENDIF
7e1bd93d4f Jean*0289
0290
0291 ENDIF
0f5ba23f97 Jean*0292
7e1bd93d4f Jean*0293
f1f873e18c Jean*0294 ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
9e5f2093d3 Alis*0295
ff02675122 Jean*0296
0297
e305438401 Mart*0298 #ifdef ALLOW_AUTODIFF_TAMC
0299
0300 #endif /* ALLOW_AUTODIFF_TAMC */
9e5f2093d3 Alis*0301
842349fcb7 Jean*0302 IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
1adfdf6281 Jean*0303
9e5f2093d3 Alis*0304 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0305
0306
9e5f2093d3 Alis*0307 #endif /* ALLOW_AUTODIFF_TAMC */
842349fcb7 Jean*0308 CALL FIND_RHO_2D(
0309 I iMin, iMax, jMin, jMax, k,
b98355c8c5 jm-c 0310 I theta(1-OLx,1-OLy,k,bi,bj),
0311 I salt(1-OLx,1-OLy,k,bi,bj),
842349fcb7 Jean*0312 O alphaRho,
0313 I k, bi, bj, myThid )
0314 ELSE
0315 DO j=jMin,jMax
0316 DO i=iMin,iMax
0317 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
0318 ENDDO
0319 ENDDO
0320 ENDIF
c5c831ccb0 Jean*0321
7448700841 Mart*0322 #ifdef ALLOW_AUTODIFF_TAMC
0323
0324 #endif
8fabf86828 Jean*0325
0326
0327
1adfdf6281 Jean*0328 DO j=jMin,jMax
0329 DO i=iMin,iMax
8fabf86828 Jean*0330 locBuoy = alphaRho(i,j)*recip_rhoConst
0331 alphaRho(i,j) = -maskC(i,j,k,bi,bj)*recip_rhoConst
0332 & * locBuoy/( oneRL + locBuoy )
1adfdf6281 Jean*0333 ENDDO
0334 ENDDO
0335
0d373cd03b Jean*0336 #ifdef ALLOW_MOM_COMMON
0337
ee19c4a296 Jean*0338 IF ( quasiHydrostatic ) THEN
0339 CALL MOM_QUASIHYDROSTATIC( bi, bj, k, uVel, vVel,
0340 U alphaRho,
0341 I myTime, myIter, myThid )
0d373cd03b Jean*0342 ENDIF
0343 #endif /* ALLOW_MOM_COMMON */
7448700841 Mart*0344 #ifdef ALLOW_AUTODIFF_TAMC
0345
0346 #endif
0d373cd03b Jean*0347
7e1bd93d4f Jean*0348
0349
0350 IF (integr_GeoPot.EQ.1) THEN
0351
0352
0353 DO j=jMin,jMax
9e5f2093d3 Alis*0354 DO i=iMin,iMax
fb481a83c2 Alis*0355
7e1bd93d4f Jean*0356
0357
0358
2b8326b7dd Jean*0359
807d4a3bf5 Jean*0360 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0361 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
0362 #ifdef NONLIN_FRSURF
0363 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
0364 #endif
0365 phiHydC(i,j) = ddRloc*alphaRho(i,j)
0f5ba23f97 Jean*0366
f1f873e18c Jean*0367
0d373cd03b Jean*0368
f1f873e18c Jean*0369
0370
0371 ELSE
0d373cd03b Jean*0372 phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
0373
f1f873e18c Jean*0374 ENDIF
0375
0d373cd03b Jean*0376 phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
f1f873e18c Jean*0377
7e1bd93d4f Jean*0378 ENDDO
0379 ENDDO
0380
0381 ELSE
f1f873e18c Jean*0382
0383
0d373cd03b Jean*0384 dRlocM = halfRL*drC(k)
f1f873e18c Jean*0385 IF (k.EQ.1) dRlocM=rF(k)-rC(k)
0386 IF (k.EQ.Nr) THEN
0387 dRlocP=rC(k)-rF(k+1)
0388 ELSE
0d373cd03b Jean*0389 dRlocP=halfRL*drC(k+1)
f1f873e18c Jean*0390 ENDIF
0d373cd03b Jean*0391 rec_dRm = oneRL/(rF(k)-rC(k))
0392 rec_dRp = oneRL/(rC(k)-rF(k+1))
7e1bd93d4f Jean*0393
0394 DO j=jMin,jMax
0395 DO i=iMin,iMax
0396
0397
0398
807d4a3bf5 Jean*0399 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0400 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
0401 #ifdef NONLIN_FRSURF
0402 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
0403 #endif
0d373cd03b Jean*0404 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
0405 & +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
f1f873e18c Jean*0406 & )*alphaRho(i,j)
0407 ELSE
0408 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
0409 ENDIF
0410 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
9e5f2093d3 Alis*0411 ENDDO
7e1bd93d4f Jean*0412 ENDDO
0413
0414
0415 ENDIF
fb481a83c2 Alis*0416
f1f873e18c Jean*0417 ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
45adaae56c Jean*0418
fb481a83c2 Alis*0419
0420
0421
0422
5edeab65c1 Jean*0423 IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
6cdaaf6acc Jean*0424
5edeab65c1 Jean*0425 IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
cfbe180681 Jean*0426
5edeab65c1 Jean*0427 thetaRef = thetaConst
0428 ELSE
cfbe180681 Jean*0429
5edeab65c1 Jean*0430 thetaRef = tRef(k)
0431 ENDIF
0432 DO j=jMin,jMax
0433 DO i=iMin,iMax
b98355c8c5 jm-c 0434 alphaRho(i,j) = ( theta(i,j,k,bi,bj)
0435 & *( salt(i,j,k,bi,bj)*atm_Rq + oneRL )
5edeab65c1 Jean*0436 & - thetaRef )*maskC(i,j,k,bi,bj)
0437 ENDDO
0438 ENDDO
0439 ELSE
0440 DO j=jMin,jMax
0441 DO i=iMin,iMax
0442 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
0443 ENDDO
0444 ENDDO
cfbe180681 Jean*0445 ENDIF
6cdaaf6acc Jean*0446
0d373cd03b Jean*0447 #ifdef ALLOW_MOM_COMMON
0448
ee19c4a296 Jean*0449 IF ( quasiHydrostatic ) THEN
0450 CALL MOM_QUASIHYDROSTATIC( bi, bj, k, uVel, vVel,
0451 U alphaRho,
0452 I myTime, myIter, myThid )
0d373cd03b Jean*0453 ENDIF
0454 #endif /* ALLOW_MOM_COMMON */
7448700841 Mart*0455 #ifdef ALLOW_AUTODIFF_TAMC
0456
0457 #endif
0d373cd03b Jean*0458
f1f873e18c Jean*0459
fb481a83c2 Alis*0460
3d9080e0a5 Jean*0461 #ifndef DISABLE_SIGMA_CODE
0462
0463 IF ( useFVgradPhi ) THEN
0464
0465 fullDepth = rF(1)-rF(Nr+1)
0466 DO j=jMin,jMax
0467 DO i=iMin,iMax
0468 locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
0469 #ifdef NONLIN_FRSURF
0470 locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
0471 #endif
0472 pKappaF(i,j) = (
0473 & ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
0474 & + bHybSigmF( k )*locDepth
0475 & )/atm_Po )**atm_kappa
0476 pKappaC = (
0477 & ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
0478 & + bHybSigmC( k )*locDepth
0479 & )/atm_Po )**atm_kappa
0480 pKappaU(i,j) = (
0481 & ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
0482 & + bHybSigmF(k+1)*locDepth
0483 & )/atm_Po )**atm_kappa
0484 phiHydC(i,j) = phiHydF(i,j)
0485 & + atm_Cp*( pKappaF(i,j) - pKappaC )*alphaRho(i,j)
0486 phiHydU(i,j) = phiHydF(i,j)
0487 & + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
0488 ENDDO
0489 ENDDO
0490
0491
0492 ELSEIF (integr_GeoPot.EQ.0) THEN
0493 #else /* DISABLE_SIGMA_CODE */
f1f873e18c Jean*0494 IF (integr_GeoPot.EQ.0) THEN
3d9080e0a5 Jean*0495 #endif /* DISABLE_SIGMA_CODE */
f1f873e18c Jean*0496
45adaae56c Jean*0497
0498
0f5ba23f97 Jean*0499
e9b27c9813 Alis*0500
45adaae56c Jean*0501
0502
f1f873e18c Jean*0503 IF (k.EQ.1) THEN
0504 ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
0505 & -((rC( k )/atm_Po)**atm_kappa) )
0506 ELSE
0507 ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
0d373cd03b Jean*0508 & -((rC( k )/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0509 ENDIF
0510 IF (k.EQ.Nr) THEN
0511 ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0512 & -((rF(k+1)/atm_Po)**atm_kappa) )
0513 ELSE
0514 ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0d373cd03b Jean*0515 & -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0516 ENDIF
45adaae56c Jean*0517
f1f873e18c Jean*0518 DO j=jMin,jMax
0519 DO i=iMin,iMax
6cdaaf6acc Jean*0520 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
0521 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
45adaae56c Jean*0522 ENDDO
f1f873e18c Jean*0523 ENDDO
45adaae56c Jean*0524
fb481a83c2 Alis*0525
0526
f1f873e18c Jean*0527 ELSEIF (integr_GeoPot.EQ.1) THEN
0528
45adaae56c Jean*0529
0530
0531
0532
0533
0f5ba23f97 Jean*0534
45adaae56c Jean*0535
0536
f1f873e18c Jean*0537 ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
0538 & -((rC( k )/atm_Po)**atm_kappa) )
0539 ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0540 & -((rF(k+1)/atm_Po)**atm_kappa) )
0541 DO j=jMin,jMax
0542 DO i=iMin,iMax
807d4a3bf5 Jean*0543 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0544 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
0545 #ifdef NONLIN_FRSURF
0546 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
0547 #endif
0548 phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
6cdaaf6acc Jean*0549 & *ddPIm*alphaRho(i,j)
f1f873e18c Jean*0550 ELSE
6cdaaf6acc Jean*0551 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
f1f873e18c Jean*0552 ENDIF
6cdaaf6acc Jean*0553 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
45adaae56c Jean*0554 ENDDO
f1f873e18c Jean*0555 ENDDO
0556
fb481a83c2 Alis*0557
0558
f1f873e18c Jean*0559 ELSEIF ( integr_GeoPot.EQ.2
0560 & .OR. integr_GeoPot.EQ.3 ) THEN
0f5ba23f97 Jean*0561
f1f873e18c Jean*0562
0563
45adaae56c Jean*0564
0565
0566
0f5ba23f97 Jean*0567
45adaae56c Jean*0568
f1f873e18c Jean*0569 IF (k.EQ.1) THEN
0570 ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
0571 & -((rC( k )/atm_Po)**atm_kappa) )
0572 ELSE
0573 ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
0d373cd03b Jean*0574 & -((rC( k )/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0575 ENDIF
0576 IF (k.EQ.Nr) THEN
0577 ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0578 & -((rF(k+1)/atm_Po)**atm_kappa) )
0579 ELSE
0580 ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0d373cd03b Jean*0581 & -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0582 ENDIF
0d373cd03b Jean*0583 rec_dRm = oneRL/(rF(k)-rC(k))
0584 rec_dRp = oneRL/(rC(k)-rF(k+1))
f1f873e18c Jean*0585 DO j=jMin,jMax
0586 DO i=iMin,iMax
807d4a3bf5 Jean*0587 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0588 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
0589 #ifdef NONLIN_FRSURF
0590 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
0591 #endif
0d373cd03b Jean*0592 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
0593 & +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
0594 & )*alphaRho(i,j)
f1f873e18c Jean*0595 ELSE
6cdaaf6acc Jean*0596 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
f1f873e18c Jean*0597 ENDIF
6cdaaf6acc Jean*0598 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
fb481a83c2 Alis*0599 ENDDO
f1f873e18c Jean*0600 ENDDO
0601
45adaae56c Jean*0602
fb481a83c2 Alis*0603
f1f873e18c Jean*0604 ELSE
0605 STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
0606 ENDIF
fb481a83c2 Alis*0607
45adaae56c Jean*0608
fb481a83c2 Alis*0609 ELSE
463053c692 Jean*0610 STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
fb481a83c2 Alis*0611 ENDIF
924557e60a Chri*0612
7448700841 Mart*0613 #ifdef ALLOW_AUTODIFF_TAMC
0614 # ifdef NONLIN_FRSURF
0615
0616
0617
0618 # endif
0619 #endif
0620
3d9080e0a5 Jean*0621 IF ( .NOT. useFVgradPhi ) THEN
0622
0623
0624 IF ( momPressureForcing ) THEN
0625 CALL CALC_GRAD_PHI_HYD(
0626 I k, bi, bj, iMin,iMax, jMin,jMax,
b98355c8c5 jm-c 0627 I phiHydC, alphaRho,
3d9080e0a5 Jean*0628 O dPhiHydX, dPhiHydY,
0629 I myTime, myIter, myThid)
0630 ENDIF
0631
0632 #ifndef DISABLE_SIGMA_CODE
0633 ELSE
0634
0635
0636 IF ( fluidIsWater ) THEN
0637 STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
0638 ENDIF
0639 IF ( momPressureForcing ) THEN
0640 CALL CALC_GRAD_PHI_FV(
0641 I k, bi, bj, iMin,iMax, jMin,jMax,
0642 I phiHydF, phiHydU, pKappaF, pKappaU,
0643 O dPhiHydX, dPhiHydY,
0644 I myTime, myIter, myThid)
0645 ENDIF
0646 DO j=jMin,jMax
0647 DO i=iMin,iMax
0648 phiHydF(i,j) = phiHydU(i,j)
0649 ENDDO
0650 ENDDO
0651
0652 #endif /* DISABLE_SIGMA_CODE */
0653
0654 ENDIF
0655
f1f873e18c Jean*0656
0657
0658
0659
0660 IF (useDiagPhiRlow) THEN
0661 CALL DIAGS_PHI_RLOW(
0662 I k, bi, bj, iMin,iMax, jMin,jMax,
b98355c8c5 jm-c 0663 I phiHydF, phiHydC, alphaRho,
0f5ba23f97 Jean*0664 I myTime, myIter, myThid)
f1f873e18c Jean*0665 ENDIF
0666
0667
0668 CALL DIAGS_PHI_HYD(
0669 I k, bi, bj, iMin,iMax, jMin,jMax,
0670 I phiHydC,
0671 I myTime, myIter, myThid)
0672
45adaae56c Jean*0673 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
1dbaea09ee Chri*0674
09f06f8597 Jean*0675 RETURN
0676 END