8 |
C !INTERFACE: |
C !INTERFACE: |
9 |
SUBROUTINE CALC_PHI_HYD( |
SUBROUTINE CALC_PHI_HYD( |
10 |
I bi, bj, iMin, iMax, jMin, jMax, K, |
I bi, bj, iMin, iMax, jMin, jMax, K, |
11 |
I theta, salt, |
I tFld, sFld, |
12 |
U phiHyd, |
U phiHyd, |
13 |
I myThid) |
I myThid) |
14 |
C !DESCRIPTION: \bv |
C !DESCRIPTION: \bv |
18 |
C *==========================================================* |
C *==========================================================* |
19 |
C | Potential (ocean: Pressure/rho ; atmos = geopotential)| |
C | Potential (ocean: Pressure/rho ; atmos = geopotential)| |
20 |
C | On entry: | |
C | On entry: | |
21 |
C | theta,salt are the current thermodynamics quantities| |
C | tFld,sFld are the current thermodynamics quantities| |
22 |
C | (unchanged on exit) | |
C | (unchanged on exit) | |
23 |
C | phiHyd(i,j,1:k-1) is the hydrostatic Potential | |
C | phiHyd(i,j,1:k-1) is the hydrostatic Potential | |
24 |
C | at cell centers (tracer points) | |
C | at cell centers (tracer points) | |
51 |
#include "tamc_keys.h" |
#include "tamc_keys.h" |
52 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
53 |
#include "SURFACE.h" |
#include "SURFACE.h" |
54 |
|
#include "DYNVARS.h" |
55 |
|
|
56 |
C !INPUT/OUTPUT PARAMETERS: |
C !INPUT/OUTPUT PARAMETERS: |
57 |
C == Routine arguments == |
C == Routine arguments == |
58 |
INTEGER bi,bj,iMin,iMax,jMin,jMax,K |
INTEGER bi,bj,iMin,iMax,jMin,jMax,K |
59 |
_RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
_RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
60 |
_RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
_RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) |
61 |
_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
62 |
INTEGER myThid |
INTEGER myThid |
63 |
|
|
133 |
C Calculate density |
C Calculate density |
134 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
135 |
kkey = (ikey-1)*Nr + k |
kkey = (ikey-1)*Nr + k |
136 |
CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
137 |
CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
138 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
139 |
CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType, |
CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, |
140 |
& theta, salt, |
& tFld, sFld, |
141 |
& alphaRho, myThid) |
& alphaRho, myThid) |
142 |
|
|
143 |
C Hydrostatic pressure at cell centers |
C Hydrostatic pressure at cell centers |
164 |
C---------- This discretization is the "energy conserving" form |
C---------- This discretization is the "energy conserving" form |
165 |
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 |
166 |
C |
C |
167 |
|
|
168 |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
169 |
& 0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst |
& 0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst |
170 |
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 |
& 0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst |
& 0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst |
172 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
173 |
|
|
174 |
|
C---------- Compute bottom pressure deviation from gravity*rho0*H |
175 |
|
C This has to be done starting from phiHyd at the current |
176 |
|
C tracer point and .5 of the cell's thickness has to be |
177 |
|
C substracted from hFacC |
178 |
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
179 |
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
180 |
|
& + (hFacC(i,j,k,bi,bj)-0.5)*drF(K) |
181 |
|
& *gravity*alphaRho(i,j)*recip_rhoConst |
182 |
|
& + gravity*etaN(i,j,bi,bj) |
183 |
|
ENDIF |
184 |
|
C----------------------------------------------------------------------- |
185 |
|
|
186 |
ENDDO |
ENDDO |
187 |
ENDDO |
ENDDO |
188 |
|
|
214 |
C Calculate density |
C Calculate density |
215 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
216 |
kkey = (ikey-1)*Nr + k |
kkey = (ikey-1)*Nr + k |
217 |
CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
218 |
CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
219 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
220 |
CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType, |
CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, |
221 |
& theta, salt, |
& tFld, sFld, |
222 |
& alphaRho, myThid) |
& alphaRho, myThid) |
223 |
|
|
224 |
C Hydrostatic pressure at cell centers |
C Hydrostatic pressure at cell centers |
245 |
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)+ |
246 |
& 0.5*dRlocKp1*locAlpha |
& 0.5*dRlocKp1*locAlpha |
247 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
248 |
|
|
249 |
|
C---------- Compute gravity*(sea surface elevation) |
250 |
|
C This has to be done starting from phiHyd at the current |
251 |
|
C tracer point and .5 of the cell's thickness has to be |
252 |
|
C substracted from hFacC |
253 |
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
254 |
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
255 |
|
& + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha |
256 |
|
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
257 |
|
ENDIF |
258 |
|
C----------------------------------------------------------------------- |
259 |
|
|
260 |
ENDDO |
ENDDO |
261 |
ENDDO |
ENDDO |
262 |
|
|
283 |
DO i=iMin,iMax |
DO i=iMin,iMax |
284 |
phiHyd(i,j,K)= |
phiHyd(i,j,K)= |
285 |
& ddPIp*maskC(i,j,K,bi,bj) |
& ddPIp*maskC(i,j,K,bi,bj) |
286 |
& *(theta(I,J,K,bi,bj)-tRef(K)) |
& *(tFld(I,J,K,bi,bj)-tRef(K)) |
287 |
ENDDO |
ENDDO |
288 |
ENDDO |
ENDDO |
289 |
ELSE |
ELSE |
294 |
DO i=iMin,iMax |
DO i=iMin,iMax |
295 |
phiHyd(i,j,K)=phiHyd(i,j,K-1) |
phiHyd(i,j,K)=phiHyd(i,j,K-1) |
296 |
& +ddPI*maskC(i,j,K-1,bi,bj) |
& +ddPI*maskC(i,j,K-1,bi,bj) |
297 |
& *(theta(I,J,K-1,bi,bj)-tRef(K-1)) |
& *(tFld(I,J,K-1,bi,bj)-tRef(K-1)) |
298 |
& +ddPI*maskC(i,j, K ,bi,bj) |
& +ddPI*maskC(i,j, K ,bi,bj) |
299 |
& *(theta(I,J, K ,bi,bj)-tRef( K )) |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
300 |
C Old code (atmos-exact) looked like this |
C Old code (atmos-exact) looked like this |
301 |
Cold phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI* |
Cold phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI* |
302 |
Cold & (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K)) |
Cold & (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K)) |
303 |
ENDDO |
ENDDO |
304 |
ENDDO |
ENDDO |
305 |
ENDIF |
ENDIF |
323 |
DO i=iMin,iMax |
DO i=iMin,iMax |
324 |
phiHyd(i,j,K) = |
phiHyd(i,j,K) = |
325 |
& ddPIp*_hFacC(I,J, K ,bi,bj) |
& ddPIp*_hFacC(I,J, K ,bi,bj) |
326 |
& *(theta(I,J, K ,bi,bj)-tRef( K )) |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
327 |
ENDDO |
ENDDO |
328 |
ENDDO |
ENDDO |
329 |
ELSE |
ELSE |
335 |
DO i=iMin,iMax |
DO i=iMin,iMax |
336 |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
337 |
& +ddPIm*_hFacC(I,J,K-1,bi,bj) |
& +ddPIm*_hFacC(I,J,K-1,bi,bj) |
338 |
& *(theta(I,J,K-1,bi,bj)-tRef(K-1)) |
& *(tFld(I,J,K-1,bi,bj)-tRef(K-1)) |
339 |
& +ddPIp*_hFacC(I,J, K ,bi,bj) |
& +ddPIp*_hFacC(I,J, K ,bi,bj) |
340 |
& *(theta(I,J, K ,bi,bj)-tRef( K )) |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
341 |
ENDDO |
ENDDO |
342 |
ENDDO |
ENDDO |
343 |
ENDIF |
ENDIF |
362 |
phiHyd(i,j,K) = |
phiHyd(i,j,K) = |
363 |
& ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half) |
& ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half) |
364 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) ) |
365 |
& *(theta(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
366 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
367 |
ENDDO |
ENDDO |
368 |
ENDDO |
ENDDO |
375 |
DO i=iMin,iMax |
DO i=iMin,iMax |
376 |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
377 |
& + ddPIm*0.5 |
& + ddPIm*0.5 |
378 |
& *(theta(i,j,K-1,bi,bj)-tRef(K-1)) |
& *(tFld(i,j,K-1,bi,bj)-tRef(K-1)) |
379 |
& * maskC(i,j,K-1,bi,bj) |
& * maskC(i,j,K-1,bi,bj) |
380 |
& +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half) |
& +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half) |
381 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) ) |
382 |
& *(theta(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
383 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
384 |
ENDDO |
ENDDO |
385 |
ENDDO |
ENDDO |
407 |
phiHyd(i,j,K) = |
phiHyd(i,j,K) = |
408 |
& ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
& ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
409 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
410 |
& *(theta(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
411 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
412 |
ENDDO |
ENDDO |
413 |
ENDDO |
ENDDO |
422 |
DO i=iMin,iMax |
DO i=iMin,iMax |
423 |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
424 |
& + ddPIm*0.5 |
& + ddPIm*0.5 |
425 |
& *(theta(i,j,K-1,bi,bj)-tRef(K-1)) |
& *(tFld(i,j,K-1,bi,bj)-tRef(K-1)) |
426 |
& * maskC(i,j,K-1,bi,bj) |
& * maskC(i,j,K-1,bi,bj) |
427 |
& +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
& +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
428 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
429 |
& *(theta(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
430 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
431 |
ENDDO |
ENDDO |
432 |
ENDDO |
ENDDO |