model/src/calc_phi_hyd.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.24 2002/12/10 02:55:47 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_PHI_HYD
C     !INTERFACE:
      SUBROUTINE CALC_PHI_HYD(
     I                         bi, bj, iMin, iMax, jMin, jMax, K,
     I                         tFld, sFld,
     U                         phiHyd,
     O                         dPhiHydX, dPhiHydY,
     I                         myTime, myIter, myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE CALC_PHI_HYD                                  |
C     | o Integrate the hydrostatic relation to find the Hydros. | 
C     *==========================================================*
C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|
C     | On entry:                                                |
C     |   tFld,sFld     are the current thermodynamics quantities|
C     |                 (unchanged on exit)                      |
C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |
C     |                 at cell centers (tracer points)          |
C     |                 - 1:k-1 layers are valid                 |
C     |                 - k:Nr layers are invalid                |
C     |   phiHyd(i,j,k) is the hydrostatic Potential             |
C     |  (ocean only_^) at cell the interface k (w point above)  |
C     | On exit:                                                 |
C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |
C     |                 at cell centers (tracer points)          |
C     |                 - 1:k layers are valid                   |
C     |                 - k+1:Nr layers are invalid              |
C     |   phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho)   |
C     |  (ocean only-^) at cell the interface k+1 (w point below)|
C     | Atmosphere:                                              |
C     |   integr_GeoPot allows to select one integration method  |
C     |    (see the list below)                                  |
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
c #include "FFIELDS.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */
#include "SURFACE.h"
#include "DYNVARS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
      INTEGER bi,bj,iMin,iMax,jMin,jMax,K
      _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL myTime
      INTEGER myIter, myThid
      
#ifdef INCLUDE_PHIHYD_CALCULATION_CODE

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j, Kp1
      _RL zero, one, half
      _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dRloc,dRlocKp1,locAlpha
      _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
      INTEGER iMnLoc,jMnLoc
      PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
CEOP

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C  Atmosphere:  
C integr_GeoPot => select one option for the integration of the Geopotential:
C   = 0 : Energy Conserving Form, No hFac ;
C   = 1 : Finite Volume Form, with hFac, linear in P by Half level;
C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
C     2 : case Tracer level at the middle of InterFace_W; 
C     3 : case InterFace_W  at the middle of Tracer levels;
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1

          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1

          act3 = myThid - 1
          max3 = nTx*nTy

          act4 = ikey_dynamics - 1

          ikey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
C       This is the hydrostatic pressure calculation for the Ocean
C       which uses the FIND_RHO() routine to calculate density
C       before integrating g*rho over the current layer/interface
#ifdef      ALLOW_AUTODIFF_TAMC
CADJ GENERAL
#endif      /* ALLOW_AUTODIFF_TAMC */

        dRloc=drC(k)
        IF (k.EQ.1) dRloc=drF(1)
        IF (k.EQ.Nr) THEN
          dRlocKp1=0.
        ELSE
          dRlocKp1=drC(k+1)
        ENDIF

C--     If this is the top layer we impose the boundary condition
C       P(z=eta) = P(atmospheric_loading)
        IF (k.EQ.1) THEN
          DO j=jMin,jMax
            DO i=iMin,iMax
              phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
c             phiHyd(i,j,k) = 0.
            ENDDO
          ENDDO
        ENDIF

C       Calculate density
#ifdef ALLOW_AUTODIFF_TAMC
        kkey = (ikey-1)*Nr + k
CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
        CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
     &                 tFld, sFld,
     &                 alphaRho, myThid)

C Quasi-hydrostatic terms are added in as if they modify the buoyancy
        IF (quasiHydrostatic) THEN
         CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
        ENDIF

C---   Hydrostatic pressure at cell centers

       IF (integr_GeoPot.EQ.1) THEN
C  --  Finite Volume Form

         DO j=jMin,jMax
          DO i=iMin,iMax

C---------- This discretization is the "finite volume" form
C           which has not been used to date since it does not
C           conserve KE+PE exactly even though it is more natural
C
          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
     &          + hFacC(i,j,k,bi,bj)
     &             *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
          ENDIF
           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
     &            + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
           phiHyd(i,j,k)=phiHyd(i,j,k)+
     &       + half*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst

          ENDDO
         ENDDO

       ELSE
C  --  Finite Difference Form

         DO j=jMin,jMax
          DO i=iMin,iMax

C---------- This discretization is the "energy conserving" form
C           which has been used since at least Adcroft et al., MWR 1997
C
             
            phiHyd(i,j,k)=phiHyd(i,j,k)
     &        +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
     &        +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst

C---------- Compute bottom pressure deviation from gravity*rho0*H
C           This has to be done starting from phiHyd at the current
C           tracer point and .5 of the cell's thickness has to be
C           substracted from hFacC
            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
     &              + (hFacC(i,j,k,bi,bj)-half)*drF(K)
     &                   *gravity*alphaRho(i,j)*recip_rhoConst
     &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
            ENDIF

          ENDDO
         ENDDO

C  --  end if integr_GeoPot = ...
       ENDIF
        
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
C       This is the hydrostatic pressure calculation for the Ocean
C       which uses the FIND_RHO() routine to calculate density
C       before integrating (1/rho)'*dp over the current layer/interface
#ifdef      ALLOW_AUTODIFF_TAMC
CADJ GENERAL
#endif      /* ALLOW_AUTODIFF_TAMC */

        dRloc=drC(k)
        IF (k.EQ.1) dRloc=drF(1)
        IF (k.EQ.Nr) THEN
          dRlocKp1=0.
        ELSE
          dRlocKp1=drC(k+1)
        ENDIF

        IF (k.EQ.1) THEN
          DO j=jMin,jMax
            DO i=iMin,iMax
              phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
c             phiHyd(i,j,k) = 0.
            ENDDO
          ENDDO
        ENDIF

C       Calculate density
#ifdef ALLOW_AUTODIFF_TAMC
            kkey = (ikey-1)*Nr + k
CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
        CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
     &                 tFld, sFld,
     &                 alphaRho, myThid)
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */


C----  Hydrostatic pressure at cell centers

       IF (integr_GeoPot.EQ.1) THEN
C  --  Finite Volume Form

         DO j=jMin,jMax
          DO i=iMin,iMax
            locAlpha=alphaRho(i,j)+rhoConst
            locAlpha=maskC(i,j,k,bi,bj)*
     &              (one/locAlpha - recip_rhoConst)
c           IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha

C---------- This discretization is the "finite volume" form
C           which has not been used to date since it does not
C           conserve KE+PE exactly even though it is more natural
C
            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
            ENDIF
            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
            phiHyd(i,j,k)=phiHyd(i,j,k)
     &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*locAlpha

          ENDDO
         ENDDO

       ELSE
C  --  Finite Difference Form

         DO j=jMin,jMax
          DO i=iMin,iMax
            locAlpha=alphaRho(i,j)+rhoConst
            locAlpha=maskC(i,j,k,bi,bj)*
     &              (one/locAlpha - recip_rhoConst)
c           IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha

C---------- This discretization is the "energy conserving" form

            phiHyd(i,j,k)=phiHyd(i,j,k)
     &          + half*dRloc*locAlpha
            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
     &          + half*dRlocKp1*locAlpha


C---------- Compute gravity*(sea surface elevation) first
C           This has to be done starting from phiHyd at the current
C           tracer point and .5 of the cell's thickness has to be
C           substracted from hFacC
            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
     &              + (hFacC(i,j,k,bi,bj)-half)*drF(k)*locAlpha
     &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
            ENDIF

          ENDDO
         ENDDO

C  --  end if integr_GeoPot = ...
       ENDIF

      ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C       This is the hydrostatic geopotential calculation for the Atmosphere
C       The ideal gas law is used implicitly here rather than calculating
C       the specific volume, analogous to the oceanic case.

C       Integrate d Phi / d pi

      IF (integr_GeoPot.EQ.0) THEN
C  --  Energy Conserving Form, No hFac  --
C------------ The integration for the first level phi(k=1) is the same
C             for both the "finite volume" and energy conserving methods.
Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary 
C             condition is simply Phi-prime(Ro_surf)=0.
C           o convention ddPI > 0 (same as drF & drC)
C-----------------------------------------------------------------------
        IF (K.EQ.1) THEN
          ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
     &                  -((rC(K)/atm_Po)**atm_kappa) )
          DO j=jMin,jMax
           DO i=iMin,iMax
c            phiHyd(i,j,K)= 
             phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
     &          ddPIp*maskC(i,j,K,bi,bj)
     &               *(tFld(I,J,K,bi,bj)-tRef(K))
           ENDDO
          ENDDO
        ELSE
C-------- This discretization is the energy conserving form
          ddPI=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
     &                 -((rC( K )/atm_Po)**atm_kappa) )*0.5
          DO j=jMin,jMax
           DO i=iMin,iMax
              phiHyd(i,j,K)=phiHyd(i,j,K-1)
     &           +ddPI*maskC(i,j,K-1,bi,bj)
     &                *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
     &           +ddPI*maskC(i,j, K ,bi,bj)
     &                *(tFld(I,J, K ,bi,bj)-tRef( K ))
C             Old code (atmos-exact) looked like this
Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K))
           ENDDO
          ENDDO
        ENDIF
C end: Energy Conserving Form, No hFac  --
C-----------------------------------------------------------------------

      ELSEIF (integr_GeoPot.EQ.1) THEN
C  --  Finite Volume Form, with hFac, linear in P by Half level  --
C---------
C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
C   Note: a true Finite Volume form should be linear between 2 Interf_W :
C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
C   also: if Interface_W at the middle between tracer levels, this form
C     is close to the Energy Cons. form in the Interior, except for the 
C     non-linearity in PI(p)
C---------
        IF (K.EQ.1) THEN
          ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
     &                  -((rC(K)/atm_Po)**atm_kappa) )
          DO j=jMin,jMax
           DO i=iMin,iMax
c            phiHyd(i,j,K)=
             phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
     &          ddPIp*_hFacC(I,J, K ,bi,bj)
     &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
           ENDDO
          ENDDO
        ELSE
          ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
     &                  -((rF( K )/atm_Po)**atm_kappa) )
          ddPIp=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
     &                  -((rC( K )/atm_Po)**atm_kappa) )
          DO j=jMin,jMax
           DO i=iMin,iMax
             phiHyd(i,j,K) = phiHyd(i,j,K-1)
     &         +ddPIm*_hFacC(I,J,K-1,bi,bj)
     &               *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
     &         +ddPIp*_hFacC(I,J, K ,bi,bj)
     &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
           ENDDO
          ENDDO
        ENDIF
C end: Finite Volume Form, with hFac, linear in P by Half level  --
C-----------------------------------------------------------------------

      ELSEIF (integr_GeoPot.EQ.2) THEN
C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C  Finite Difference formulation consistent with Partial Cell,
C    case Tracer level at the middle of InterFace_W 
C    linear between 2 Tracer levels ; conserve energy in the Interior 
C---------
        Kp1 = min(Nr,K+1)
        IF (K.EQ.1) THEN
          ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
     &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0
          ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
     &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
          DO j=jMin,jMax
           DO i=iMin,iMax
c            phiHyd(i,j,K)=
             phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
     &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
     &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
     &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
     &               * maskC(i,j, K ,bi,bj)
           ENDDO
          ENDDO
        ELSE
          ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
     &                  -((rC( K )/atm_Po)**atm_kappa) )
          ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
     &                  -((rC(Kp1)/atm_Po)**atm_kappa) ) 
          DO j=jMin,jMax
           DO i=iMin,iMax
             phiHyd(i,j,K) = phiHyd(i,j,K-1)
     &        + ddPIm*0.5
     &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
     &               * maskC(i,j,K-1,bi,bj)
     &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
     &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
     &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
     &               * maskC(i,j, K ,bi,bj)
           ENDDO
          ENDDO
        ENDIF
C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
C-----------------------------------------------------------------------

      ELSEIF (integr_GeoPot.EQ.3) THEN
C  --  Finite Difference Form, with hFac, Interface_W = middle  --
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C  Finite Difference formulation consistent with Partial Cell,
C   Valid & accurate if Interface_W at middle between tracer levels
C   linear in p between 2 Tracer levels ; conserve energy in the Interior 
C---------
        Kp1 = min(Nr,K+1)
        IF (K.EQ.1) THEN
          ratioRm=0.5*drF(K)/(rF(k)-rC(K))
          ratioRp=drF(K)*recip_drC(Kp1)
          ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
     &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0
          ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
     &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
          DO j=jMin,jMax
           DO i=iMin,iMax
c            phiHyd(i,j,K)=
             phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
     &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
     &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
     &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
     &               * maskC(i,j, K ,bi,bj)
           ENDDO
          ENDDO
        ELSE
          ratioRm=drF(K)*recip_drC(K)
          ratioRp=drF(K)*recip_drC(Kp1)
          ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
     &                  -((rC( K )/atm_Po)**atm_kappa) )
          ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
     &                  -((rC(Kp1)/atm_Po)**atm_kappa) ) 
          DO j=jMin,jMax
           DO i=iMin,iMax
             phiHyd(i,j,K) = phiHyd(i,j,K-1)
     &        + ddPIm*0.5
     &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
     &               * maskC(i,j,K-1,bi,bj)
     &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
     &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
     &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
     &               * maskC(i,j, K ,bi,bj)
           ENDDO
          ENDDO
        ENDIF
C end: Finite Difference Form, with hFac, Interface_W = middle  --
C-----------------------------------------------------------------------

      ELSE
        STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      ELSE
        STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
      ENDIF

      IF (momPressureForcing) THEN 
        iMnLoc = MAX(1-Olx+1,iMin)
        jMnLoc = MAX(1-Oly+1,jMin)
        CALL CALC_GRAD_PHI_HYD(
     I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
     I                         phiHyd, alphaRho, tFld, sFld,
     O                         dPhiHydX, dPhiHydY,
     I                         myTime, myIter, myThid)  
      ENDIF

#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */

      RETURN
      END
1	jmc	1.25	C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.24 2002/12/10 02:55:47 jmc Exp $
2	cnh	1.16	C $Name: $
3	cnh	1.1
4	cnh	1.6	#include "CPP_OPTIONS.h"
5	cnh	1.1
6	cnh	1.16	CBOP
7			C !ROUTINE: CALC_PHI_HYD
8			C !INTERFACE:
9	adcroft	1.9	SUBROUTINE CALC_PHI_HYD(
10			I bi, bj, iMin, iMax, jMin, jMax, K,
11	mlosch	1.20	I tFld, sFld,
12	adcroft	1.9	U phiHyd,
13	jmc	1.25	O dPhiHydX, dPhiHydY,
14			I myTime, myIter, myThid)
15	cnh	1.16	C !DESCRIPTION: \bv
16			C ==========================================================
17	cnh	1.1	C \| SUBROUTINE CALC_PHI_HYD \|
18	jmc	1.11	C \| o Integrate the hydrostatic relation to find the Hydros. \|
19	cnh	1.16	C ==========================================================
20	jmc	1.11	C \| Potential (ocean: Pressure/rho ; atmos = geopotential)\|
21	adcroft	1.9	C \| On entry: \|
22	mlosch	1.20	C \| tFld,sFld are the current thermodynamics quantities\|
23	adcroft	1.9	C \| (unchanged on exit) \|
24	jmc	1.11	C \| phiHyd(i,j,1:k-1) is the hydrostatic Potential \|
25	adcroft	1.9	C \| at cell centers (tracer points) \|
26			C \| - 1:k-1 layers are valid \|
27			C \| - k:Nr layers are invalid \|
28	jmc	1.11	C \| phiHyd(i,j,k) is the hydrostatic Potential \|
29	jmc	1.14	C \| (ocean only_^) at cell the interface k (w point above) \|
30	adcroft	1.9	C \| On exit: \|
31	jmc	1.11	C \| phiHyd(i,j,1:k) is the hydrostatic Potential \|
32	adcroft	1.9	C \| at cell centers (tracer points) \|
33			C \| - 1:k layers are valid \|
34			C \| - k+1:Nr layers are invalid \|
35	jmc	1.11	C \| phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho) \|
36	jmc	1.14	C \| (ocean only-^) at cell the interface k+1 (w point below)\|
37			C \| Atmosphere: \|
38	jmc	1.24	C \| integr_GeoPot allows to select one integration method \|
39	jmc	1.14	C \| (see the list below) \|
40	cnh	1.16	C ==========================================================
41			C \ev
42			C !USES:
43	cnh	1.1	IMPLICIT NONE
44			C == Global variables ==
45			#include "SIZE.h"
46			#include "GRID.h"
47			#include "EEPARAMS.h"
48			#include "PARAMS.h"
49	jmc	1.24	c #include "FFIELDS.h"
50	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
51			#include "tamc.h"
52			#include "tamc_keys.h"
53			#endif /* ALLOW_AUTODIFF_TAMC */
54	adcroft	1.19	#include "SURFACE.h"
55	mlosch	1.20	#include "DYNVARS.h"
56	heimbach	1.13
57	cnh	1.16	C !INPUT/OUTPUT PARAMETERS:
58	cnh	1.1	C == Routine arguments ==
59			INTEGER bi,bj,iMin,iMax,jMin,jMax,K
60	mlosch	1.20	_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)
62	cnh	1.2	_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63	jmc	1.25	_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	jmc	1.14
68	adcroft	1.9	#ifdef INCLUDE_PHIHYD_CALCULATION_CODE
69
70	cnh	1.16	C !LOCAL VARIABLES:
71	cnh	1.1	C == Local variables ==
72	jmc	1.14	INTEGER i,j, Kp1
73			_RL zero, one, half
74	adcroft	1.9	_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75	adcroft	1.19	_RL dRloc,dRlocKp1,locAlpha
76	jmc	1.14	_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
77	jmc	1.25	INTEGER iMnLoc,jMnLoc
78			PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
79	cnh	1.16	CEOP
80	jmc	1.14
81			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
82			C Atmosphere:
83	jmc	1.24	C integr_GeoPot => select one option for the integration of the Geopotential:
84	jmc	1.14	C = 0 : Energy Conserving Form, No hFac ;
85			C = 1 : Finite Volume Form, with hFac, linear in P by Half level;
86			C =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
87			C 2 : case Tracer level at the middle of InterFace_W;
88			C 3 : case InterFace_W at the middle of Tracer levels;
89			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
90	adcroft	1.9
91	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
92			act1 = bi - myBxLo(myThid)
93			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
94
95			act2 = bj - myByLo(myThid)
96			max2 = myByHi(myThid) - myByLo(myThid) + 1
97
98			act3 = myThid - 1
99			max3 = nTx*nTy
100
101			act4 = ikey_dynamics - 1
102
103			ikey = (act1 + 1) + act2*max1
104			& + act3max1max2
105			& + act4max1max2*max3
106			#endif /* ALLOW_AUTODIFF_TAMC */
107
108	jmc	1.25
109			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
110	adcroft	1.9	IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
111			C This is the hydrostatic pressure calculation for the Ocean
112			C which uses the FIND_RHO() routine to calculate density
113			C before integrating g*rho over the current layer/interface
114	jmc	1.25	#ifdef ALLOW_AUTODIFF_TAMC
115			CADJ GENERAL
116			#endif /* ALLOW_AUTODIFF_TAMC */
117	adcroft	1.9
118			dRloc=drC(k)
119			IF (k.EQ.1) dRloc=drF(1)
120			IF (k.EQ.Nr) THEN
121			dRlocKp1=0.
122			ELSE
123			dRlocKp1=drC(k+1)
124			ENDIF
125
126			C-- If this is the top layer we impose the boundary condition
127			C P(z=eta) = P(atmospheric_loading)
128			IF (k.EQ.1) THEN
129			DO j=jMin,jMax
130			DO i=iMin,iMax
131	jmc	1.24	phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
132	jmc	1.25	c phiHyd(i,j,k) = 0.
133	adcroft	1.9	ENDDO
134			ENDDO
135			ENDIF
136
137			C Calculate density
138	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
139	heimbach	1.23	kkey = (ikey-1)*Nr + k
140			CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
141	mlosch	1.20	CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
142	heimbach	1.13	#endif /* ALLOW_AUTODIFF_TAMC */
143	mlosch	1.20	CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
144			& tFld, sFld,
145	adcroft	1.9	& alphaRho, myThid)
146	adcroft	1.22
147			C Quasi-hydrostatic terms are added in as if they modify the buoyancy
148			IF (quasiHydrostatic) THEN
149			CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
150			ENDIF
151	adcroft	1.9
152	jmc	1.25	C--- Hydrostatic pressure at cell centers
153
154			IF (integr_GeoPot.EQ.1) THEN
155			C -- Finite Volume Form
156
157			DO j=jMin,jMax
158	adcroft	1.9	DO i=iMin,iMax
159
160	jmc	1.25	C---------- This discretization is the "finite volume" form
161			C which has not been used to date since it does not
162			C conserve KE+PE exactly even though it is more natural
163			C
164			IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
165			phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
166			& + hFacC(i,j,k,bi,bj)
167			& drF(K)gravityalphaRho(i,j)recip_rhoConst
168			& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
169			ENDIF
170			IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
171			& + drF(K)gravityalphaRho(i,j)*recip_rhoConst
172			phiHyd(i,j,k)=phiHyd(i,j,k)+
173			& + halfdrF(K)gravityalphaRho(i,j)recip_rhoConst
174
175			ENDDO
176			ENDDO
177
178			ELSE
179			C -- Finite Difference Form
180
181			DO j=jMin,jMax
182			DO i=iMin,iMax
183	adcroft	1.9
184			C---------- This discretization is the "energy conserving" form
185			C which has been used since at least Adcroft et al., MWR 1997
186			C
187	mlosch	1.20
188	jmc	1.25	phiHyd(i,j,k)=phiHyd(i,j,k)
189			& +halfdRlocgravityalphaRho(i,j)recip_rhoConst
190			IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
191			& +halfdRlocKp1gravityalphaRho(i,j)recip_rhoConst
192	mlosch	1.20
193			C---------- Compute bottom pressure deviation from gravityrho0H
194			C This has to be done starting from phiHyd at the current
195			C tracer point and .5 of the cell's thickness has to be
196			C substracted from hFacC
197			IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
198			phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
199	jmc	1.25	& + (hFacC(i,j,k,bi,bj)-half)*drF(K)
200	mlosch	1.20	& gravityalphaRho(i,j)*recip_rhoConst
201	jmc	1.25	& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
202	mlosch	1.20	ENDIF
203
204	adcroft	1.9	ENDDO
205	jmc	1.25	ENDDO
206
207			C -- end if integr_GeoPot = ...
208			ENDIF
209	adcroft	1.9
210	jmc	1.25	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
211	adcroft	1.19	ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
212			C This is the hydrostatic pressure calculation for the Ocean
213			C which uses the FIND_RHO() routine to calculate density
214	jmc	1.25	C before integrating (1/rho)'*dp over the current layer/interface
215	mlosch	1.21	#ifdef ALLOW_AUTODIFF_TAMC
216			CADJ GENERAL
217			#endif /* ALLOW_AUTODIFF_TAMC */
218	adcroft	1.19
219			dRloc=drC(k)
220			IF (k.EQ.1) dRloc=drF(1)
221			IF (k.EQ.Nr) THEN
222			dRlocKp1=0.
223			ELSE
224			dRlocKp1=drC(k+1)
225			ENDIF
226
227			IF (k.EQ.1) THEN
228			DO j=jMin,jMax
229			DO i=iMin,iMax
230	jmc	1.24	phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
231	jmc	1.25	c phiHyd(i,j,k) = 0.
232	adcroft	1.19	ENDDO
233			ENDDO
234			ENDIF
235
236			C Calculate density
237			#ifdef ALLOW_AUTODIFF_TAMC
238			kkey = (ikey-1)*Nr + k
239	heimbach	1.23	CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
240	mlosch	1.20	CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
241	adcroft	1.19	#endif /* ALLOW_AUTODIFF_TAMC */
242	mlosch	1.20	CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
243			& tFld, sFld,
244	adcroft	1.19	& alphaRho, myThid)
245	heimbach	1.23	#ifdef ALLOW_AUTODIFF_TAMC
246			CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
247			#endif /* ALLOW_AUTODIFF_TAMC */
248
249	adcroft	1.19
250	jmc	1.25	C---- Hydrostatic pressure at cell centers
251
252			IF (integr_GeoPot.EQ.1) THEN
253			C -- Finite Volume Form
254
255			DO j=jMin,jMax
256	adcroft	1.19	DO i=iMin,iMax
257	mlosch	1.21	locAlpha=alphaRho(i,j)+rhoConst
258	jmc	1.25	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	adcroft	1.19
276	jmc	1.25	ENDDO
277			ENDDO
278
279			ELSE
280			C -- Finite Difference Form
281
282			DO j=jMin,jMax
283			DO i=iMin,iMax
284			locAlpha=alphaRho(i,j)+rhoConst
285			locAlpha=maskC(i,j,k,bi,bj)*
286			& (one/locAlpha - recip_rhoConst)
287			c IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
288	adcroft	1.9
289	adcroft	1.19	C---------- This discretization is the "energy conserving" form
290	mlosch	1.21
291	jmc	1.25	phiHyd(i,j,k)=phiHyd(i,j,k)
292			& + halfdRloclocAlpha
293			IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
294			& + halfdRlocKp1locAlpha
295	mlosch	1.21
296	mlosch	1.20
297	mlosch	1.21	C---------- Compute gravity*(sea surface elevation) first
298	mlosch	1.20	C This has to be done starting from phiHyd at the current
299			C tracer point and .5 of the cell's thickness has to be
300			C substracted from hFacC
301			IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
302			phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
303	jmc	1.25	& + (hFacC(i,j,k,bi,bj)-half)drF(k)locAlpha
304	mlosch	1.20	& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
305			ENDIF
306
307	adcroft	1.19	ENDDO
308	jmc	1.25	ENDDO
309
310			C -- end if integr_GeoPot = ...
311			ENDIF
312	adcroft	1.9
313			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
314	jmc	1.14	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
315	adcroft	1.9	C This is the hydrostatic geopotential calculation for the Atmosphere
316			C The ideal gas law is used implicitly here rather than calculating
317			C the specific volume, analogous to the oceanic case.
318
319			C Integrate d Phi / d pi
320
321	jmc	1.24	IF (integr_GeoPot.EQ.0) THEN
322	jmc	1.14	C -- Energy Conserving Form, No hFac --
323			C------------ The integration for the first level phi(k=1) is the same
324			C for both the "finite volume" and energy conserving methods.
325	adcroft	1.17	Ci NOTE o Working with geopotential Anomaly, the geopotential boundary
326			C condition is simply Phi-prime(Ro_surf)=0.
327	jmc	1.14	C o convention ddPI > 0 (same as drF & drC)
328			C-----------------------------------------------------------------------
329	adcroft	1.9	IF (K.EQ.1) THEN
330	jmc	1.24	ddPIp=atm_Cp( ((rF(K)/atm_Po)*atm_kappa)
331			& -((rC(K)/atm_Po)**atm_kappa) )
332	adcroft	1.9	DO j=jMin,jMax
333	jmc	1.14	DO i=iMin,iMax
334	jmc	1.25	c phiHyd(i,j,K)=
335			phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
336			& ddPIp*maskC(i,j,K,bi,bj)
337	mlosch	1.20	& *(tFld(I,J,K,bi,bj)-tRef(K))
338	jmc	1.14	ENDDO
339			ENDDO
340			ELSE
341			C-------- This discretization is the energy conserving form
342	jmc	1.24	ddPI=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
343			& -((rC( K )/atm_Po)*atm_kappa) )0.5
344	jmc	1.14	DO j=jMin,jMax
345			DO i=iMin,iMax
346			phiHyd(i,j,K)=phiHyd(i,j,K-1)
347			& +ddPI*maskC(i,j,K-1,bi,bj)
348	mlosch	1.20	& *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
349	jmc	1.14	& +ddPI*maskC(i,j, K ,bi,bj)
350	mlosch	1.20	& *(tFld(I,J, K ,bi,bj)-tRef( K ))
351	jmc	1.14	C Old code (atmos-exact) looked like this
352			Cold phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
353	mlosch	1.20	Cold & (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K))
354	jmc	1.14	ENDDO
355			ENDDO
356			ENDIF
357			C end: Energy Conserving Form, No hFac --
358	adcroft	1.9	C-----------------------------------------------------------------------
359	jmc	1.14
360	jmc	1.24	ELSEIF (integr_GeoPot.EQ.1) THEN
361	jmc	1.14	C -- Finite Volume Form, with hFac, linear in P by Half level --
362			C---------
363			C Finite Volume formulation consistent with Partial Cell, linear in p by piece
364			C Note: a true Finite Volume form should be linear between 2 Interf_W :
365			C phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
366			C also: if Interface_W at the middle between tracer levels, this form
367			C is close to the Energy Cons. form in the Interior, except for the
368			C non-linearity in PI(p)
369			C---------
370			IF (K.EQ.1) THEN
371	jmc	1.24	ddPIp=atm_Cp( ((rF(K)/atm_Po)*atm_kappa)
372			& -((rC(K)/atm_Po)**atm_kappa) )
373	jmc	1.14	DO j=jMin,jMax
374			DO i=iMin,iMax
375	jmc	1.25	c phiHyd(i,j,K)=
376			phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
377			& ddPIp*_hFacC(I,J, K ,bi,bj)
378	mlosch	1.20	& *(tFld(I,J, K ,bi,bj)-tRef( K ))
379	adcroft	1.9	ENDDO
380			ENDDO
381			ELSE
382	jmc	1.24	ddPIm=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
383			& -((rF( K )/atm_Po)**atm_kappa) )
384			ddPIp=atm_Cp( ((rF( K )/atm_Po)*atm_kappa)
385			& -((rC( K )/atm_Po)**atm_kappa) )
386	jmc	1.14	DO j=jMin,jMax
387			DO i=iMin,iMax
388			phiHyd(i,j,K) = phiHyd(i,j,K-1)
389	mlosch	1.18	& +ddPIm*_hFacC(I,J,K-1,bi,bj)
390	mlosch	1.20	& *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
391	mlosch	1.18	& +ddPIp*_hFacC(I,J, K ,bi,bj)
392	mlosch	1.20	& *(tFld(I,J, K ,bi,bj)-tRef( K ))
393	jmc	1.14	ENDDO
394			ENDDO
395			ENDIF
396			C end: Finite Volume Form, with hFac, linear in P by Half level --
397	adcroft	1.9	C-----------------------------------------------------------------------
398
399	jmc	1.24	ELSEIF (integr_GeoPot.EQ.2) THEN
400	jmc	1.14	C -- Finite Difference Form, with hFac, Tracer Lev. = middle --
401			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
402			C Finite Difference formulation consistent with Partial Cell,
403			C case Tracer level at the middle of InterFace_W
404			C linear between 2 Tracer levels ; conserve energy in the Interior
405			C---------
406			Kp1 = min(Nr,K+1)
407			IF (K.EQ.1) THEN
408	jmc	1.24	ddPIm=atm_Cp( ((rF( K )/atm_Po)*atm_kappa)
409			& -((rC( K )/atm_Po)*atm_kappa) ) 2. _d 0
410			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
411			& -((rC(Kp1)/atm_Po)**atm_kappa) )
412	jmc	1.14	DO j=jMin,jMax
413			DO i=iMin,iMax
414	jmc	1.25	c phiHyd(i,j,K)=
415			phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
416			& ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
417	mlosch	1.18	& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
418	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
419	jmc	1.14	& * maskC(i,j, K ,bi,bj)
420			ENDDO
421			ENDDO
422			ELSE
423	jmc	1.24	ddPIm=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
424			& -((rC( K )/atm_Po)**atm_kappa) )
425			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
426			& -((rC(Kp1)/atm_Po)**atm_kappa) )
427	jmc	1.14	DO j=jMin,jMax
428			DO i=iMin,iMax
429			phiHyd(i,j,K) = phiHyd(i,j,K-1)
430			& + ddPIm*0.5
431	mlosch	1.20	& *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
432	jmc	1.14	& * maskC(i,j,K-1,bi,bj)
433	mlosch	1.18	& +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
434			& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
435	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
436	jmc	1.14	& * maskC(i,j, K ,bi,bj)
437			ENDDO
438			ENDDO
439			ENDIF
440			C end: Finite Difference Form, with hFac, Tracer Lev. = middle --
441	adcroft	1.9	C-----------------------------------------------------------------------
442
443	jmc	1.24	ELSEIF (integr_GeoPot.EQ.3) THEN
444	jmc	1.14	C -- Finite Difference Form, with hFac, Interface_W = middle --
445			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
446			C Finite Difference formulation consistent with Partial Cell,
447			C Valid & accurate if Interface_W at middle between tracer levels
448			C linear in p between 2 Tracer levels ; conserve energy in the Interior
449			C---------
450			Kp1 = min(Nr,K+1)
451			IF (K.EQ.1) THEN
452			ratioRm=0.5*drF(K)/(rF(k)-rC(K))
453			ratioRp=drF(K)*recip_drC(Kp1)
454	jmc	1.24	ddPIm=atm_Cp( ((rF( K )/atm_Po)*atm_kappa)
455			& -((rC( K )/atm_Po)*atm_kappa) ) 2. _d 0
456			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
457			& -((rC(Kp1)/atm_Po)**atm_kappa) )
458	jmc	1.14	DO j=jMin,jMax
459			DO i=iMin,iMax
460	jmc	1.25	c phiHyd(i,j,K)=
461			phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
462			& ( ddPImmax(zero,(_hFacC(i,j,K,bi,bj)-one)ratioRm+half)
463	mlosch	1.18	& +ddPIpmin(zero, _hFacC(i,j,K,bi,bj)ratioRp -half) )
464	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
465	jmc	1.14	& * maskC(i,j, K ,bi,bj)
466			ENDDO
467			ENDDO
468			ELSE
469			ratioRm=drF(K)*recip_drC(K)
470			ratioRp=drF(K)*recip_drC(Kp1)
471	jmc	1.24	ddPIm=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
472			& -((rC( K )/atm_Po)**atm_kappa) )
473			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
474			& -((rC(Kp1)/atm_Po)**atm_kappa) )
475	adcroft	1.9	DO j=jMin,jMax
476	jmc	1.14	DO i=iMin,iMax
477			phiHyd(i,j,K) = phiHyd(i,j,K-1)
478			& + ddPIm*0.5
479	mlosch	1.20	& *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
480	jmc	1.14	& * maskC(i,j,K-1,bi,bj)
481	mlosch	1.18	& +(ddPImmax(zero,(_hFacC(i,j,K,bi,bj)-one)ratioRm+half)
482			& +ddPIpmin(zero, _hFacC(i,j,K,bi,bj)ratioRp -half) )
483	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
484	jmc	1.14	& * maskC(i,j, K ,bi,bj)
485			ENDDO
486	adcroft	1.9	ENDDO
487			ENDIF
488	jmc	1.14	C end: Finite Difference Form, with hFac, Interface_W = middle --
489			C-----------------------------------------------------------------------
490	cnh	1.1
491	jmc	1.14	ELSE
492	jmc	1.24	STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
493	jmc	1.14	ENDIF
494	cnh	1.6
495	jmc	1.14	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
496	adcroft	1.9	ELSE
497	jmc	1.24	STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
498	jmc	1.25	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	cnh	1.5	ENDIF
509	cnh	1.1
510	jmc	1.14	#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
511	cnh	1.6
512	jmc	1.11	RETURN
513			END