/[MITgcm]/MITgcm/model/src/calc_phi_hyd.F
ViewVC logotype

Diff of /MITgcm/model/src/calc_phi_hyd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.23 by heimbach, Fri Nov 15 03:01:21 2002 UTC revision 1.24 by jmc, Tue Dec 10 02:55:47 2002 UTC
# Line 34  C     |                 - k+1:Nr layers Line 34  C     |                 - k+1:Nr layers
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
# Line 45  C     == Global variables == Line 45  C     == Global variables ==
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"
# Line 78  CEOP Line 78  CEOP
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
# Line 121  C       P(z=eta) = P(atmospheric_loading Line 121  C       P(z=eta) = P(atmospheric_loading
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
# Line 216  CADJ GENERAL Line 212  CADJ GENERAL
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
# Line 292  C       the specific volume, analogous t Line 285  C       the specific volume, analogous t
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.
# Line 301  C             condition is simply Phi-pr Line 294  C             condition is simply Phi-pr
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)
# Line 330  Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I, Line 323  Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I,
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
# Line 341  C     is close to the Energy Cons. form Line 334  C     is close to the Energy Cons. form
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)
# Line 368  C--------- Line 361  C---------
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,
# Line 377  C    linear between 2 Tracer levels ; co Line 370  C    linear between 2 Tracer levels ; co
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)
# Line 411  C--------- Line 404  C---------
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,
# Line 422  C--------- Line 415  C---------
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)
# Line 438  C--------- Line 431  C---------
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)
# Line 459  C end: Finite Difference Form, with hFac Line 452  C end: Finite Difference Form, with hFac
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 */
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22