10 |
I bi, bj, iMin, iMax, jMin, jMax, K, |
I bi, bj, iMin, iMax, jMin, jMax, K, |
11 |
I tFld, sFld, |
I tFld, sFld, |
12 |
U phiHyd, |
U phiHyd, |
13 |
I myThid) |
O dPhiHydX, dPhiHydY, |
14 |
|
I myTime, myIter, myThid) |
15 |
C !DESCRIPTION: \bv |
C !DESCRIPTION: \bv |
16 |
C *==========================================================* |
C *==========================================================* |
17 |
C | SUBROUTINE CALC_PHI_HYD | |
C | SUBROUTINE CALC_PHI_HYD | |
60 |
_RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
_RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
61 |
_RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
_RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
62 |
_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
63 |
INTEGER myThid |
_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
64 |
|
_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
65 |
|
_RL myTime |
66 |
|
INTEGER myIter, myThid |
67 |
|
|
68 |
#ifdef INCLUDE_PHIHYD_CALCULATION_CODE |
#ifdef INCLUDE_PHIHYD_CALCULATION_CODE |
69 |
|
|
74 |
_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
75 |
_RL dRloc,dRlocKp1,locAlpha |
_RL dRloc,dRlocKp1,locAlpha |
76 |
_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm |
_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm |
77 |
|
INTEGER iMnLoc,jMnLoc |
78 |
|
PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 ) |
79 |
CEOP |
CEOP |
80 |
|
|
|
zero = 0. _d 0 |
|
|
one = 1. _d 0 |
|
|
half = .5 _d 0 |
|
|
|
|
81 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
82 |
C Atmosphere: |
C Atmosphere: |
83 |
C integr_GeoPot => select one option for the integration of the Geopotential: |
C integr_GeoPot => select one option for the integration of the Geopotential: |
105 |
& + act4*max1*max2*max3 |
& + act4*max1*max2*max3 |
106 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
107 |
|
|
108 |
|
|
109 |
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
110 |
IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN |
IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN |
111 |
C This is the hydrostatic pressure calculation for the Ocean |
C This is the hydrostatic pressure calculation for the Ocean |
112 |
C which uses the FIND_RHO() routine to calculate density |
C which uses the FIND_RHO() routine to calculate density |
113 |
C before integrating g*rho over the current layer/interface |
C before integrating g*rho over the current layer/interface |
114 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
115 |
|
CADJ GENERAL |
116 |
|
#endif /* ALLOW_AUTODIFF_TAMC */ |
117 |
|
|
118 |
dRloc=drC(k) |
dRloc=drC(k) |
119 |
IF (k.EQ.1) dRloc=drF(1) |
IF (k.EQ.1) dRloc=drF(1) |
129 |
DO j=jMin,jMax |
DO j=jMin,jMax |
130 |
DO i=iMin,iMax |
DO i=iMin,iMax |
131 |
phiHyd(i,j,k) = phi0surf(i,j,bi,bj) |
phiHyd(i,j,k) = phi0surf(i,j,bi,bj) |
132 |
|
c phiHyd(i,j,k) = 0. |
133 |
ENDDO |
ENDDO |
134 |
ENDDO |
ENDDO |
135 |
ENDIF |
ENDIF |
149 |
CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid) |
CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid) |
150 |
ENDIF |
ENDIF |
151 |
|
|
152 |
C Hydrostatic pressure at cell centers |
C--- Hydrostatic pressure at cell centers |
153 |
DO j=jMin,jMax |
|
154 |
|
IF (integr_GeoPot.EQ.1) THEN |
155 |
|
C -- Finite Volume Form |
156 |
|
|
157 |
|
DO j=jMin,jMax |
158 |
DO i=iMin,iMax |
DO i=iMin,iMax |
|
#ifdef ALLOW_AUTODIFF_TAMC |
|
|
c Patrick, is this directive correct or even necessary in |
|
|
c this new code? |
|
|
c Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+... |
|
|
c within the k-loop. |
|
|
CADJ GENERAL |
|
|
#endif /* ALLOW_AUTODIFF_TAMC */ |
|
159 |
|
|
160 |
CmlC---------- This discretization is the "finite volume" form |
C---------- This discretization is the "finite volume" form |
161 |
CmlC which has not been used to date since it does not |
C which has not been used to date since it does not |
162 |
CmlC conserve KE+PE exactly even though it is more natural |
C conserve KE+PE exactly even though it is more natural |
163 |
CmlC |
C |
164 |
Cml IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
165 |
Cml phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
166 |
Cml & + hFacC(i,j,k,bi,bj) |
& + hFacC(i,j,k,bi,bj) |
167 |
Cml & *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
& *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
168 |
Cml & + gravity*etaN(i,j,bi,bj) |
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
169 |
Cml ENDIF |
ENDIF |
170 |
Cml IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
171 |
Cml & drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
& + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
172 |
Cml phiHyd(i,j,k)=phiHyd(i,j,k)+ |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
173 |
Cml & 0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
& + half*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
174 |
CmlC----------------------------------------------------------------------- |
|
175 |
|
ENDDO |
176 |
|
ENDDO |
177 |
|
|
178 |
|
ELSE |
179 |
|
C -- Finite Difference Form |
180 |
|
|
181 |
|
DO j=jMin,jMax |
182 |
|
DO i=iMin,iMax |
183 |
|
|
184 |
C---------- This discretization is the "energy conserving" form |
C---------- This discretization is the "energy conserving" form |
185 |
C which has been used since at least Adcroft et al., MWR 1997 |
C which has been used since at least Adcroft et al., MWR 1997 |
186 |
C |
C |
187 |
|
|
188 |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
phiHyd(i,j,k)=phiHyd(i,j,k) |
189 |
& 0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst |
& +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst |
190 |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
191 |
& 0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst |
& +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst |
|
C----------------------------------------------------------------------- |
|
192 |
|
|
193 |
C---------- Compute bottom pressure deviation from gravity*rho0*H |
C---------- Compute bottom pressure deviation from gravity*rho0*H |
194 |
C This has to be done starting from phiHyd at the current |
C This has to be done starting from phiHyd at the current |
196 |
C substracted from hFacC |
C substracted from hFacC |
197 |
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
198 |
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
199 |
& + (hFacC(i,j,k,bi,bj)-.5)*drF(K) |
& + (hFacC(i,j,k,bi,bj)-half)*drF(K) |
200 |
& *gravity*alphaRho(i,j)*recip_rhoConst |
& *gravity*alphaRho(i,j)*recip_rhoConst |
201 |
& + gravity*etaN(i,j,bi,bj) |
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
202 |
ENDIF |
ENDIF |
|
C----------------------------------------------------------------------- |
|
203 |
|
|
204 |
ENDDO |
ENDDO |
205 |
ENDDO |
ENDDO |
206 |
|
|
207 |
|
C -- end if integr_GeoPot = ... |
208 |
|
ENDIF |
209 |
|
|
210 |
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
211 |
ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN |
ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN |
212 |
C This is the hydrostatic pressure calculation for the Ocean |
C This is the hydrostatic pressure calculation for the Ocean |
213 |
C which uses the FIND_RHO() routine to calculate density |
C which uses the FIND_RHO() routine to calculate density |
214 |
C before integrating g*rho over the current layer/interface |
C before integrating (1/rho)'*dp over the current layer/interface |
215 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
216 |
CADJ GENERAL |
CADJ GENERAL |
217 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
228 |
DO j=jMin,jMax |
DO j=jMin,jMax |
229 |
DO i=iMin,iMax |
DO i=iMin,iMax |
230 |
phiHyd(i,j,k) = phi0surf(i,j,bi,bj) |
phiHyd(i,j,k) = phi0surf(i,j,bi,bj) |
231 |
|
c phiHyd(i,j,k) = 0. |
232 |
ENDDO |
ENDDO |
233 |
ENDDO |
ENDDO |
234 |
ENDIF |
ENDIF |
247 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
248 |
|
|
249 |
|
|
250 |
C Hydrostatic pressure at cell centers |
C---- Hydrostatic pressure at cell centers |
251 |
DO j=jMin,jMax |
|
252 |
|
IF (integr_GeoPot.EQ.1) THEN |
253 |
|
C -- Finite Volume Form |
254 |
|
|
255 |
|
DO j=jMin,jMax |
256 |
DO i=iMin,iMax |
DO i=iMin,iMax |
257 |
locAlpha=alphaRho(i,j)+rhoConst |
locAlpha=alphaRho(i,j)+rhoConst |
258 |
IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha |
locAlpha=maskC(i,j,k,bi,bj)* |
259 |
|
& (one/locAlpha - recip_rhoConst) |
260 |
|
c IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha |
261 |
|
|
262 |
|
C---------- This discretization is the "finite volume" form |
263 |
|
C which has not been used to date since it does not |
264 |
|
C conserve KE+PE exactly even though it is more natural |
265 |
|
C |
266 |
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
267 |
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
268 |
|
& + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha |
269 |
|
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
270 |
|
ENDIF |
271 |
|
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
272 |
|
& + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha |
273 |
|
phiHyd(i,j,k)=phiHyd(i,j,k) |
274 |
|
& +(hFacC(i,j,k,bi,bj)-half)*drF(K)*locAlpha |
275 |
|
|
276 |
CmlC---------- This discretization is the "finite volume" form |
ENDDO |
277 |
CmlC which has not been used to date since it does not |
ENDDO |
278 |
CmlC conserve KE+PE exactly even though it is more natural |
|
279 |
CmlC |
ELSE |
280 |
Cml IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
C -- Finite Difference Form |
281 |
Cml phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
|
282 |
Cml & + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha |
DO j=jMin,jMax |
283 |
Cml & + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
DO i=iMin,iMax |
284 |
Cml ENDIF |
locAlpha=alphaRho(i,j)+rhoConst |
285 |
Cml IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
locAlpha=maskC(i,j,k,bi,bj)* |
286 |
Cml & drF(K)*locAlpha |
& (one/locAlpha - recip_rhoConst) |
287 |
Cml phiHyd(i,j,k)=phiHyd(i,j,k)+ |
c IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha |
|
Cml & 0.5*drF(K)*locAlpha |
|
|
CmlC----------------------------------------------------------------------- |
|
288 |
|
|
289 |
C---------- This discretization is the "energy conserving" form |
C---------- This discretization is the "energy conserving" form |
|
C which has been used since at least Adcroft et al., MWR 1997 |
|
|
C |
|
290 |
|
|
291 |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
phiHyd(i,j,k)=phiHyd(i,j,k) |
292 |
& 0.5*dRloc*locAlpha |
& + half*dRloc*locAlpha |
293 |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
294 |
& 0.5*dRlocKp1*locAlpha |
& + half*dRlocKp1*locAlpha |
295 |
|
|
|
C----------------------------------------------------------------------- |
|
296 |
|
|
297 |
C---------- Compute gravity*(sea surface elevation) first |
C---------- Compute gravity*(sea surface elevation) first |
298 |
C This has to be done starting from phiHyd at the current |
C This has to be done starting from phiHyd at the current |
300 |
C substracted from hFacC |
C substracted from hFacC |
301 |
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
302 |
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
303 |
& + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha |
& + (hFacC(i,j,k,bi,bj)-half)*drF(k)*locAlpha |
304 |
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
305 |
ENDIF |
ENDIF |
|
C----------------------------------------------------------------------- |
|
306 |
|
|
307 |
ENDDO |
ENDDO |
308 |
ENDDO |
ENDDO |
309 |
|
|
310 |
|
C -- end if integr_GeoPot = ... |
311 |
|
ENDIF |
312 |
|
|
313 |
ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN |
ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN |
314 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
331 |
& -((rC(K)/atm_Po)**atm_kappa) ) |
& -((rC(K)/atm_Po)**atm_kappa) ) |
332 |
DO j=jMin,jMax |
DO j=jMin,jMax |
333 |
DO i=iMin,iMax |
DO i=iMin,iMax |
334 |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
c phiHyd(i,j,K)= |
335 |
& +ddPIp*maskC(i,j,K,bi,bj) |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+ |
336 |
|
& ddPIp*maskC(i,j,K,bi,bj) |
337 |
& *(tFld(I,J,K,bi,bj)-tRef(K)) |
& *(tFld(I,J,K,bi,bj)-tRef(K)) |
338 |
ENDDO |
ENDDO |
339 |
ENDDO |
ENDDO |
372 |
& -((rC(K)/atm_Po)**atm_kappa) ) |
& -((rC(K)/atm_Po)**atm_kappa) ) |
373 |
DO j=jMin,jMax |
DO j=jMin,jMax |
374 |
DO i=iMin,iMax |
DO i=iMin,iMax |
375 |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
c phiHyd(i,j,K)= |
376 |
& +ddPIp*_hFacC(I,J, K ,bi,bj) |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+ |
377 |
|
& ddPIp*_hFacC(I,J, K ,bi,bj) |
378 |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
379 |
ENDDO |
ENDDO |
380 |
ENDDO |
ENDDO |
411 |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
412 |
DO j=jMin,jMax |
DO j=jMin,jMax |
413 |
DO i=iMin,iMax |
DO i=iMin,iMax |
414 |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
c phiHyd(i,j,K)= |
415 |
& +( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half) |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+ |
416 |
|
& ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half) |
417 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) ) |
418 |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
419 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
457 |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
458 |
DO j=jMin,jMax |
DO j=jMin,jMax |
459 |
DO i=iMin,iMax |
DO i=iMin,iMax |
460 |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
c phiHyd(i,j,K)= |
461 |
& +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+ |
462 |
|
& ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
463 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
464 |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
465 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
497 |
STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !' |
STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !' |
498 |
ENDIF |
ENDIF |
499 |
|
|
500 |
|
IF (momPressureForcing) THEN |
501 |
|
iMnLoc = MAX(1-Olx+1,iMin) |
502 |
|
jMnLoc = MAX(1-Oly+1,jMin) |
503 |
|
CALL CALC_GRAD_PHI_HYD( |
504 |
|
I k, bi, bj, iMnLoc,iMax, jMnLoc,jMax, |
505 |
|
I phiHyd, alphaRho, tFld, sFld, |
506 |
|
O dPhiHydX, dPhiHydY, |
507 |
|
I myTime, myIter, myThid) |
508 |
|
ENDIF |
509 |
|
|
510 |
#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */ |
#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */ |
511 |
|
|
512 |
RETURN |
RETURN |