34 |
C | phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho) | |
C | phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho) | |
35 |
C | (ocean only-^) at cell the interface k+1 (w point below)| |
C | (ocean only-^) at cell the interface k+1 (w point below)| |
36 |
C | Atmosphere: | |
C | Atmosphere: | |
37 |
C | Integr_GeoPot allows to select one integration method | |
C | integr_GeoPot allows to select one integration method | |
38 |
C | (see the list below) | |
C | (see the list below) | |
39 |
C *==========================================================* |
C *==========================================================* |
40 |
C \ev |
C \ev |
45 |
#include "GRID.h" |
#include "GRID.h" |
46 |
#include "EEPARAMS.h" |
#include "EEPARAMS.h" |
47 |
#include "PARAMS.h" |
#include "PARAMS.h" |
48 |
#include "FFIELDS.h" |
c #include "FFIELDS.h" |
49 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
50 |
#include "tamc.h" |
#include "tamc.h" |
51 |
#include "tamc_keys.h" |
#include "tamc_keys.h" |
78 |
|
|
79 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
80 |
C Atmosphere: |
C Atmosphere: |
81 |
C Integr_GeoPot => select one option for the integration of the Geopotential: |
C integr_GeoPot => select one option for the integration of the Geopotential: |
82 |
C = 0 : Energy Conserving Form, No hFac ; |
C = 0 : Energy Conserving Form, No hFac ; |
83 |
C = 1 : Finite Volume Form, with hFac, linear in P by Half level; |
C = 1 : Finite Volume Form, with hFac, linear in P by Half level; |
84 |
C =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels |
C =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels |
121 |
IF (k.EQ.1) THEN |
IF (k.EQ.1) THEN |
122 |
DO j=jMin,jMax |
DO j=jMin,jMax |
123 |
DO i=iMin,iMax |
DO i=iMin,iMax |
124 |
#ifdef ATMOSPHERIC_LOADING |
phiHyd(i,j,k) = phi0surf(i,j,bi,bj) |
|
phiHyd(i,j,k)=pload(i,j,bi,bj)*recip_rhoConst |
|
|
#else |
|
|
phiHyd(i,j,k)=0. _d 0 |
|
|
#endif |
|
125 |
ENDDO |
ENDDO |
126 |
ENDDO |
ENDDO |
127 |
ENDIF |
ENDIF |
212 |
IF (k.EQ.1) THEN |
IF (k.EQ.1) THEN |
213 |
DO j=jMin,jMax |
DO j=jMin,jMax |
214 |
DO i=iMin,iMax |
DO i=iMin,iMax |
215 |
phiHyd(i,j,k)=0. |
phiHyd(i,j,k) = phi0surf(i,j,bi,bj) |
|
#ifdef ATMOSPHERIC_LOADING |
|
|
phiHyd(i,j,k)=pload(i,j,bi,bj) |
|
|
#endif |
|
216 |
ENDDO |
ENDDO |
217 |
ENDDO |
ENDDO |
218 |
ENDIF |
ENDIF |
285 |
|
|
286 |
C Integrate d Phi / d pi |
C Integrate d Phi / d pi |
287 |
|
|
288 |
IF (Integr_GeoPot.EQ.0) THEN |
IF (integr_GeoPot.EQ.0) THEN |
289 |
C -- Energy Conserving Form, No hFac -- |
C -- Energy Conserving Form, No hFac -- |
290 |
C------------ The integration for the first level phi(k=1) is the same |
C------------ The integration for the first level phi(k=1) is the same |
291 |
C for both the "finite volume" and energy conserving methods. |
C for both the "finite volume" and energy conserving methods. |
294 |
C o convention ddPI > 0 (same as drF & drC) |
C o convention ddPI > 0 (same as drF & drC) |
295 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
296 |
IF (K.EQ.1) THEN |
IF (K.EQ.1) THEN |
297 |
ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa) |
298 |
& -((rC(K)/atm_po)**atm_kappa) ) |
& -((rC(K)/atm_Po)**atm_kappa) ) |
299 |
DO j=jMin,jMax |
DO j=jMin,jMax |
300 |
DO i=iMin,iMax |
DO i=iMin,iMax |
301 |
phiHyd(i,j,K)= |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
302 |
& ddPIp*maskC(i,j,K,bi,bj) |
& +ddPIp*maskC(i,j,K,bi,bj) |
303 |
& *(tFld(I,J,K,bi,bj)-tRef(K)) |
& *(tFld(I,J,K,bi,bj)-tRef(K)) |
304 |
ENDDO |
ENDDO |
305 |
ENDDO |
ENDDO |
306 |
ELSE |
ELSE |
307 |
C-------- This discretization is the energy conserving form |
C-------- This discretization is the energy conserving form |
308 |
ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa) |
ddPI=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa) |
309 |
& -((rC( K )/atm_po)**atm_kappa) )*0.5 |
& -((rC( K )/atm_Po)**atm_kappa) )*0.5 |
310 |
DO j=jMin,jMax |
DO j=jMin,jMax |
311 |
DO i=iMin,iMax |
DO i=iMin,iMax |
312 |
phiHyd(i,j,K)=phiHyd(i,j,K-1) |
phiHyd(i,j,K)=phiHyd(i,j,K-1) |
323 |
C end: Energy Conserving Form, No hFac -- |
C end: Energy Conserving Form, No hFac -- |
324 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
325 |
|
|
326 |
ELSEIF (Integr_GeoPot.EQ.1) THEN |
ELSEIF (integr_GeoPot.EQ.1) THEN |
327 |
C -- Finite Volume Form, with hFac, linear in P by Half level -- |
C -- Finite Volume Form, with hFac, linear in P by Half level -- |
328 |
C--------- |
C--------- |
329 |
C Finite Volume formulation consistent with Partial Cell, linear in p by piece |
C Finite Volume formulation consistent with Partial Cell, linear in p by piece |
334 |
C non-linearity in PI(p) |
C non-linearity in PI(p) |
335 |
C--------- |
C--------- |
336 |
IF (K.EQ.1) THEN |
IF (K.EQ.1) THEN |
337 |
ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa) |
338 |
& -((rC(K)/atm_po)**atm_kappa) ) |
& -((rC(K)/atm_Po)**atm_kappa) ) |
339 |
DO j=jMin,jMax |
DO j=jMin,jMax |
340 |
DO i=iMin,iMax |
DO i=iMin,iMax |
341 |
phiHyd(i,j,K) = |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
342 |
& ddPIp*_hFacC(I,J, K ,bi,bj) |
& +ddPIp*_hFacC(I,J, K ,bi,bj) |
343 |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
& *(tFld(I,J, K ,bi,bj)-tRef( K )) |
344 |
ENDDO |
ENDDO |
345 |
ENDDO |
ENDDO |
346 |
ELSE |
ELSE |
347 |
ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa) |
ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa) |
348 |
& -((rF( K )/atm_po)**atm_kappa) ) |
& -((rF( K )/atm_Po)**atm_kappa) ) |
349 |
ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa) |
350 |
& -((rC( K )/atm_po)**atm_kappa) ) |
& -((rC( K )/atm_Po)**atm_kappa) ) |
351 |
DO j=jMin,jMax |
DO j=jMin,jMax |
352 |
DO i=iMin,iMax |
DO i=iMin,iMax |
353 |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
361 |
C end: Finite Volume Form, with hFac, linear in P by Half level -- |
C end: Finite Volume Form, with hFac, linear in P by Half level -- |
362 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
363 |
|
|
364 |
ELSEIF (Integr_GeoPot.EQ.2) THEN |
ELSEIF (integr_GeoPot.EQ.2) THEN |
365 |
C -- Finite Difference Form, with hFac, Tracer Lev. = middle -- |
C -- Finite Difference Form, with hFac, Tracer Lev. = middle -- |
366 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
367 |
C Finite Difference formulation consistent with Partial Cell, |
C Finite Difference formulation consistent with Partial Cell, |
370 |
C--------- |
C--------- |
371 |
Kp1 = min(Nr,K+1) |
Kp1 = min(Nr,K+1) |
372 |
IF (K.EQ.1) THEN |
IF (K.EQ.1) THEN |
373 |
ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa) |
ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa) |
374 |
& -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0 |
& -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0 |
375 |
ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa) |
376 |
& -((rC(Kp1)/atm_po)**atm_kappa) ) |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
377 |
DO j=jMin,jMax |
DO j=jMin,jMax |
378 |
DO i=iMin,iMax |
DO i=iMin,iMax |
379 |
phiHyd(i,j,K) = |
phiHyd(i,j,K)= phi0surf(i,j,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 |
& *(tFld(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 |
386 |
ELSE |
ELSE |
387 |
ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa) |
ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa) |
388 |
& -((rC( K )/atm_po)**atm_kappa) ) |
& -((rC( K )/atm_Po)**atm_kappa) ) |
389 |
ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa) |
390 |
& -((rC(Kp1)/atm_po)**atm_kappa) ) |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
391 |
DO j=jMin,jMax |
DO j=jMin,jMax |
392 |
DO i=iMin,iMax |
DO i=iMin,iMax |
393 |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
404 |
C end: Finite Difference Form, with hFac, Tracer Lev. = middle -- |
C end: Finite Difference Form, with hFac, Tracer Lev. = middle -- |
405 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
406 |
|
|
407 |
ELSEIF (Integr_GeoPot.EQ.3) THEN |
ELSEIF (integr_GeoPot.EQ.3) THEN |
408 |
C -- Finite Difference Form, with hFac, Interface_W = middle -- |
C -- Finite Difference Form, with hFac, Interface_W = middle -- |
409 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
410 |
C Finite Difference formulation consistent with Partial Cell, |
C Finite Difference formulation consistent with Partial Cell, |
415 |
IF (K.EQ.1) THEN |
IF (K.EQ.1) THEN |
416 |
ratioRm=0.5*drF(K)/(rF(k)-rC(K)) |
ratioRm=0.5*drF(K)/(rF(k)-rC(K)) |
417 |
ratioRp=drF(K)*recip_drC(Kp1) |
ratioRp=drF(K)*recip_drC(Kp1) |
418 |
ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa) |
ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa) |
419 |
& -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0 |
& -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0 |
420 |
ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa) |
421 |
& -((rC(Kp1)/atm_po)**atm_kappa) ) |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
422 |
DO j=jMin,jMax |
DO j=jMin,jMax |
423 |
DO i=iMin,iMax |
DO i=iMin,iMax |
424 |
phiHyd(i,j,K) = |
phiHyd(i,j,K)= phi0surf(i,j,bi,bj) |
425 |
& ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
& +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half) |
426 |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp -half) ) |
427 |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
& *(tFld(i,j, K ,bi,bj)-tRef( K )) |
428 |
& * maskC(i,j, K ,bi,bj) |
& * maskC(i,j, K ,bi,bj) |
431 |
ELSE |
ELSE |
432 |
ratioRm=drF(K)*recip_drC(K) |
ratioRm=drF(K)*recip_drC(K) |
433 |
ratioRp=drF(K)*recip_drC(Kp1) |
ratioRp=drF(K)*recip_drC(Kp1) |
434 |
ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa) |
ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa) |
435 |
& -((rC( K )/atm_po)**atm_kappa) ) |
& -((rC( K )/atm_Po)**atm_kappa) ) |
436 |
ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa) |
437 |
& -((rC(Kp1)/atm_po)**atm_kappa) ) |
& -((rC(Kp1)/atm_Po)**atm_kappa) ) |
438 |
DO j=jMin,jMax |
DO j=jMin,jMax |
439 |
DO i=iMin,iMax |
DO i=iMin,iMax |
440 |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
phiHyd(i,j,K) = phiHyd(i,j,K-1) |
452 |
C----------------------------------------------------------------------- |
C----------------------------------------------------------------------- |
453 |
|
|
454 |
ELSE |
ELSE |
455 |
STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !' |
STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !' |
456 |
ENDIF |
ENDIF |
457 |
|
|
458 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
459 |
ELSE |
ELSE |
460 |
STOP 'CALC_PHI_HYD: We should never reach this point!' |
STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !' |
461 |
ENDIF |
ENDIF |
462 |
|
|
463 |
#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */ |
#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */ |