model/src/calc_phi_hyd.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.27 2003/02/09 02:58:39 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 )
      LOGICAL useDiagPhiRlow
CEOP
      useDiagPhiRlow = .TRUE.

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
c             phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
              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---   Diagnose Hydrostatic pressure at the bottom:
       IF (useDiagPhiRlow) THEN
          CALL DIAGS_PHI_RLOW(
     I                        k, bi, bj, iMin,iMax, jMin,jMax,
     I                        phiHyd, alphaRho, tFld, sFld,
     I                        myTime, myIter, 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.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

          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
c             phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
              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--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
        DO j=jMin,jMax
          DO i=iMin,iMax
            locAlpha=alphaRho(i,j)+rhoConst
            alphaRho(i,j)=maskC(i,j,k,bi,bj)*
     &              (one/locAlpha - recip_rhoConst)
          ENDDO
        ENDDO

C---   Diagnose Sea-surface height (Hydrostatic geopotential at r=Rlow):
       IF (useDiagPhiRlow) THEN
          CALL DIAGS_PHI_RLOW(
     I                        k, bi, bj, iMin,iMax, jMin,jMax,
     I                        phiHyd, alphaRho, tFld, sFld,
     I                        myTime, myIter, 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.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
     &          + hFacC(i,j,k,bi,bj)*drF(K)*alphaRho(i,j)
            phiHyd(i,j,k)=phiHyd(i,j,k)
     &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*alphaRho(i,j)

          ENDDO
         ENDDO

       ELSE
C  --  Finite Difference Form

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

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

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

          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)= phi0surf(i,j,bi,bj)+
             phiHyd(i,j,K)= 
     &          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)= phi0surf(i,j,bi,bj)+
             phiHyd(i,j,K)=
     &          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)= phi0surf(i,j,bi,bj)+
             phiHyd(i,j,K)=
     &        ( 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)= phi0surf(i,j,bi,bj)+
             phiHyd(i,j,K)=
     &        ( 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.28	C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.27 2003/02/09 02:58:39 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	jmc	1.27	LOGICAL useDiagPhiRlow
80	cnh	1.16	CEOP
81	jmc	1.27	useDiagPhiRlow = .TRUE.
82	jmc	1.14
83			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
84			C Atmosphere:
85	jmc	1.24	C integr_GeoPot => select one option for the integration of the Geopotential:
86	jmc	1.14	C = 0 : Energy Conserving Form, No hFac ;
87			C = 1 : Finite Volume Form, with hFac, linear in P by Half level;
88			C =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
89			C 2 : case Tracer level at the middle of InterFace_W;
90			C 3 : case InterFace_W at the middle of Tracer levels;
91			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
92	adcroft	1.9
93	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
94			act1 = bi - myBxLo(myThid)
95			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
96
97			act2 = bj - myByLo(myThid)
98			max2 = myByHi(myThid) - myByLo(myThid) + 1
99
100			act3 = myThid - 1
101			max3 = nTx*nTy
102
103			act4 = ikey_dynamics - 1
104
105			ikey = (act1 + 1) + act2*max1
106			& + act3max1max2
107			& + act4max1max2*max3
108			#endif /* ALLOW_AUTODIFF_TAMC */
109
110	jmc	1.25
111			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
112	adcroft	1.9	IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
113			C This is the hydrostatic pressure calculation for the Ocean
114			C which uses the FIND_RHO() routine to calculate density
115			C before integrating g*rho over the current layer/interface
116	jmc	1.25	#ifdef ALLOW_AUTODIFF_TAMC
117			CADJ GENERAL
118			#endif /* ALLOW_AUTODIFF_TAMC */
119	adcroft	1.9
120			dRloc=drC(k)
121			IF (k.EQ.1) dRloc=drF(1)
122			IF (k.EQ.Nr) THEN
123			dRlocKp1=0.
124			ELSE
125			dRlocKp1=drC(k+1)
126			ENDIF
127
128			C-- If this is the top layer we impose the boundary condition
129			C P(z=eta) = P(atmospheric_loading)
130			IF (k.EQ.1) THEN
131			DO j=jMin,jMax
132			DO i=iMin,iMax
133	jmc	1.26	c phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
134			phiHyd(i,j,k) = 0.
135	adcroft	1.9	ENDDO
136			ENDDO
137			ENDIF
138
139			C Calculate density
140	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
141	heimbach	1.23	kkey = (ikey-1)*Nr + k
142			CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
143	mlosch	1.20	CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
144	heimbach	1.13	#endif /* ALLOW_AUTODIFF_TAMC */
145	mlosch	1.20	CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
146			& tFld, sFld,
147	adcroft	1.9	& alphaRho, myThid)
148	adcroft	1.22
149			C Quasi-hydrostatic terms are added in as if they modify the buoyancy
150			IF (quasiHydrostatic) THEN
151			CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
152			ENDIF
153	adcroft	1.9
154	jmc	1.27	C--- Diagnose Hydrostatic pressure at the bottom:
155			IF (useDiagPhiRlow) THEN
156			CALL DIAGS_PHI_RLOW(
157			I k, bi, bj, iMin,iMax, jMin,jMax,
158			I phiHyd, alphaRho, tFld, sFld,
159			I myTime, myIter, myThid)
160			ENDIF
161
162	jmc	1.25	C--- Hydrostatic pressure at cell centers
163
164			IF (integr_GeoPot.EQ.1) THEN
165			C -- Finite Volume Form
166
167			DO j=jMin,jMax
168	adcroft	1.9	DO i=iMin,iMax
169
170	jmc	1.25	C---------- This discretization is the "finite volume" form
171			C which has not been used to date since it does not
172			C conserve KE+PE exactly even though it is more natural
173			C
174			IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
175			& + drF(K)gravityalphaRho(i,j)*recip_rhoConst
176	jmc	1.28	phiHyd(i,j,k)=phiHyd(i,j,k)
177	jmc	1.25	& + halfdrF(K)gravityalphaRho(i,j)recip_rhoConst
178
179			ENDDO
180			ENDDO
181
182			ELSE
183			C -- Finite Difference Form
184
185			DO j=jMin,jMax
186			DO i=iMin,iMax
187	adcroft	1.9
188			C---------- This discretization is the "energy conserving" form
189			C which has been used since at least Adcroft et al., MWR 1997
190			C
191	jmc	1.25	phiHyd(i,j,k)=phiHyd(i,j,k)
192			& +halfdRlocgravityalphaRho(i,j)recip_rhoConst
193			IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
194			& +halfdRlocKp1gravityalphaRho(i,j)recip_rhoConst
195	mlosch	1.20
196	adcroft	1.9	ENDDO
197	jmc	1.25	ENDDO
198
199			C -- end if integr_GeoPot = ...
200			ENDIF
201	adcroft	1.9
202	jmc	1.25	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
203	adcroft	1.19	ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
204			C This is the hydrostatic pressure calculation for the Ocean
205			C which uses the FIND_RHO() routine to calculate density
206	jmc	1.25	C before integrating (1/rho)'*dp over the current layer/interface
207	mlosch	1.21	#ifdef ALLOW_AUTODIFF_TAMC
208			CADJ GENERAL
209			#endif /* ALLOW_AUTODIFF_TAMC */
210	adcroft	1.19
211			dRloc=drC(k)
212			IF (k.EQ.1) dRloc=drF(1)
213			IF (k.EQ.Nr) THEN
214			dRlocKp1=0.
215			ELSE
216			dRlocKp1=drC(k+1)
217			ENDIF
218
219			IF (k.EQ.1) THEN
220			DO j=jMin,jMax
221			DO i=iMin,iMax
222	jmc	1.26	c phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
223			phiHyd(i,j,k) = 0.
224	adcroft	1.19	ENDDO
225			ENDDO
226			ENDIF
227
228	jmc	1.27	C-- Calculate density
229	adcroft	1.19	#ifdef ALLOW_AUTODIFF_TAMC
230			kkey = (ikey-1)*Nr + k
231	heimbach	1.23	CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
232	mlosch	1.20	CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
233	adcroft	1.19	#endif /* ALLOW_AUTODIFF_TAMC */
234	mlosch	1.20	CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
235			& tFld, sFld,
236	adcroft	1.19	& alphaRho, myThid)
237	heimbach	1.23	#ifdef ALLOW_AUTODIFF_TAMC
238			CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
239			#endif /* ALLOW_AUTODIFF_TAMC */
240
241	jmc	1.27	C-- Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
242			DO j=jMin,jMax
243			DO i=iMin,iMax
244			locAlpha=alphaRho(i,j)+rhoConst
245			alphaRho(i,j)=maskC(i,j,k,bi,bj)*
246			& (one/locAlpha - recip_rhoConst)
247			ENDDO
248			ENDDO
249
250			C--- Diagnose Sea-surface height (Hydrostatic geopotential at r=Rlow):
251			IF (useDiagPhiRlow) THEN
252			CALL DIAGS_PHI_RLOW(
253			I k, bi, bj, iMin,iMax, jMin,jMax,
254			I phiHyd, alphaRho, tFld, sFld,
255			I myTime, myIter, myThid)
256			ENDIF
257	adcroft	1.19
258	jmc	1.25	C---- Hydrostatic pressure at cell centers
259
260			IF (integr_GeoPot.EQ.1) THEN
261			C -- Finite Volume Form
262
263			DO j=jMin,jMax
264	adcroft	1.19	DO i=iMin,iMax
265	jmc	1.25
266			C---------- This discretization is the "finite volume" form
267			C which has not been used to date since it does not
268			C conserve KE+PE exactly even though it is more natural
269			C
270			IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
271	jmc	1.27	& + hFacC(i,j,k,bi,bj)drF(K)alphaRho(i,j)
272	jmc	1.25	phiHyd(i,j,k)=phiHyd(i,j,k)
273	jmc	1.27	& +(hFacC(i,j,k,bi,bj)-half)drF(K)alphaRho(i,j)
274	adcroft	1.19
275	jmc	1.25	ENDDO
276			ENDDO
277
278			ELSE
279			C -- Finite Difference Form
280
281			DO j=jMin,jMax
282			DO i=iMin,iMax
283	adcroft	1.9
284	adcroft	1.19	C---------- This discretization is the "energy conserving" form
285	mlosch	1.21
286	jmc	1.25	phiHyd(i,j,k)=phiHyd(i,j,k)
287	jmc	1.27	& + halfdRlocalphaRho(i,j)
288	jmc	1.25	IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
289	jmc	1.27	& + halfdRlocKp1alphaRho(i,j)
290	mlosch	1.20
291	adcroft	1.19	ENDDO
292	jmc	1.25	ENDDO
293
294			C -- end if integr_GeoPot = ...
295			ENDIF
296	adcroft	1.9
297			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
298	jmc	1.14	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
299	adcroft	1.9	C This is the hydrostatic geopotential calculation for the Atmosphere
300			C The ideal gas law is used implicitly here rather than calculating
301			C the specific volume, analogous to the oceanic case.
302
303			C Integrate d Phi / d pi
304
305	jmc	1.24	IF (integr_GeoPot.EQ.0) THEN
306	jmc	1.14	C -- Energy Conserving Form, No hFac --
307			C------------ The integration for the first level phi(k=1) is the same
308			C for both the "finite volume" and energy conserving methods.
309	adcroft	1.17	Ci NOTE o Working with geopotential Anomaly, the geopotential boundary
310			C condition is simply Phi-prime(Ro_surf)=0.
311	jmc	1.14	C o convention ddPI > 0 (same as drF & drC)
312			C-----------------------------------------------------------------------
313	adcroft	1.9	IF (K.EQ.1) THEN
314	jmc	1.24	ddPIp=atm_Cp( ((rF(K)/atm_Po)*atm_kappa)
315			& -((rC(K)/atm_Po)**atm_kappa) )
316	adcroft	1.9	DO j=jMin,jMax
317	jmc	1.14	DO i=iMin,iMax
318	jmc	1.26	c phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
319			phiHyd(i,j,K)=
320	jmc	1.25	& ddPIp*maskC(i,j,K,bi,bj)
321	mlosch	1.20	& *(tFld(I,J,K,bi,bj)-tRef(K))
322	jmc	1.14	ENDDO
323			ENDDO
324			ELSE
325			C-------- This discretization is the energy conserving form
326	jmc	1.24	ddPI=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
327			& -((rC( K )/atm_Po)*atm_kappa) )0.5
328	jmc	1.14	DO j=jMin,jMax
329			DO i=iMin,iMax
330			phiHyd(i,j,K)=phiHyd(i,j,K-1)
331			& +ddPI*maskC(i,j,K-1,bi,bj)
332	mlosch	1.20	& *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
333	jmc	1.14	& +ddPI*maskC(i,j, K ,bi,bj)
334	mlosch	1.20	& *(tFld(I,J, K ,bi,bj)-tRef( K ))
335	jmc	1.14	C Old code (atmos-exact) looked like this
336			Cold phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
337	mlosch	1.20	Cold & (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K))
338	jmc	1.14	ENDDO
339			ENDDO
340			ENDIF
341			C end: Energy Conserving Form, No hFac --
342	adcroft	1.9	C-----------------------------------------------------------------------
343	jmc	1.14
344	jmc	1.24	ELSEIF (integr_GeoPot.EQ.1) THEN
345	jmc	1.14	C -- Finite Volume Form, with hFac, linear in P by Half level --
346			C---------
347			C Finite Volume formulation consistent with Partial Cell, linear in p by piece
348			C Note: a true Finite Volume form should be linear between 2 Interf_W :
349			C phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
350			C also: if Interface_W at the middle between tracer levels, this form
351			C is close to the Energy Cons. form in the Interior, except for the
352			C non-linearity in PI(p)
353			C---------
354			IF (K.EQ.1) THEN
355	jmc	1.24	ddPIp=atm_Cp( ((rF(K)/atm_Po)*atm_kappa)
356			& -((rC(K)/atm_Po)**atm_kappa) )
357	jmc	1.14	DO j=jMin,jMax
358			DO i=iMin,iMax
359	jmc	1.26	c phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
360			phiHyd(i,j,K)=
361	jmc	1.25	& ddPIp*_hFacC(I,J, K ,bi,bj)
362	mlosch	1.20	& *(tFld(I,J, K ,bi,bj)-tRef( K ))
363	adcroft	1.9	ENDDO
364			ENDDO
365			ELSE
366	jmc	1.24	ddPIm=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
367			& -((rF( K )/atm_Po)**atm_kappa) )
368			ddPIp=atm_Cp( ((rF( K )/atm_Po)*atm_kappa)
369			& -((rC( K )/atm_Po)**atm_kappa) )
370	jmc	1.14	DO j=jMin,jMax
371			DO i=iMin,iMax
372			phiHyd(i,j,K) = phiHyd(i,j,K-1)
373	mlosch	1.18	& +ddPIm*_hFacC(I,J,K-1,bi,bj)
374	mlosch	1.20	& *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
375	mlosch	1.18	& +ddPIp*_hFacC(I,J, K ,bi,bj)
376	mlosch	1.20	& *(tFld(I,J, K ,bi,bj)-tRef( K ))
377	jmc	1.14	ENDDO
378			ENDDO
379			ENDIF
380			C end: Finite Volume Form, with hFac, linear in P by Half level --
381	adcroft	1.9	C-----------------------------------------------------------------------
382
383	jmc	1.24	ELSEIF (integr_GeoPot.EQ.2) THEN
384	jmc	1.14	C -- Finite Difference Form, with hFac, Tracer Lev. = middle --
385			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
386			C Finite Difference formulation consistent with Partial Cell,
387			C case Tracer level at the middle of InterFace_W
388			C linear between 2 Tracer levels ; conserve energy in the Interior
389			C---------
390			Kp1 = min(Nr,K+1)
391			IF (K.EQ.1) THEN
392	jmc	1.24	ddPIm=atm_Cp( ((rF( K )/atm_Po)*atm_kappa)
393			& -((rC( K )/atm_Po)*atm_kappa) ) 2. _d 0
394			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
395			& -((rC(Kp1)/atm_Po)**atm_kappa) )
396	jmc	1.14	DO j=jMin,jMax
397			DO i=iMin,iMax
398	jmc	1.26	c phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
399			phiHyd(i,j,K)=
400	jmc	1.25	& ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
401	mlosch	1.18	& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
402	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
403	jmc	1.14	& * maskC(i,j, K ,bi,bj)
404			ENDDO
405			ENDDO
406			ELSE
407	jmc	1.24	ddPIm=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
408			& -((rC( K )/atm_Po)**atm_kappa) )
409			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
410			& -((rC(Kp1)/atm_Po)**atm_kappa) )
411	jmc	1.14	DO j=jMin,jMax
412			DO i=iMin,iMax
413			phiHyd(i,j,K) = phiHyd(i,j,K-1)
414			& + ddPIm*0.5
415	mlosch	1.20	& *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
416	jmc	1.14	& * maskC(i,j,K-1,bi,bj)
417	mlosch	1.18	& +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
418			& +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
419	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
420	jmc	1.14	& * maskC(i,j, K ,bi,bj)
421			ENDDO
422			ENDDO
423			ENDIF
424			C end: Finite Difference Form, with hFac, Tracer Lev. = middle --
425	adcroft	1.9	C-----------------------------------------------------------------------
426
427	jmc	1.24	ELSEIF (integr_GeoPot.EQ.3) THEN
428	jmc	1.14	C -- Finite Difference Form, with hFac, Interface_W = middle --
429			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
430			C Finite Difference formulation consistent with Partial Cell,
431			C Valid & accurate if Interface_W at middle between tracer levels
432			C linear in p between 2 Tracer levels ; conserve energy in the Interior
433			C---------
434			Kp1 = min(Nr,K+1)
435			IF (K.EQ.1) THEN
436			ratioRm=0.5*drF(K)/(rF(k)-rC(K))
437			ratioRp=drF(K)*recip_drC(Kp1)
438	jmc	1.24	ddPIm=atm_Cp( ((rF( K )/atm_Po)*atm_kappa)
439			& -((rC( K )/atm_Po)*atm_kappa) ) 2. _d 0
440			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
441			& -((rC(Kp1)/atm_Po)**atm_kappa) )
442	jmc	1.14	DO j=jMin,jMax
443			DO i=iMin,iMax
444	jmc	1.26	c phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
445			phiHyd(i,j,K)=
446	jmc	1.25	& ( ddPImmax(zero,(_hFacC(i,j,K,bi,bj)-one)ratioRm+half)
447	mlosch	1.18	& +ddPIpmin(zero, _hFacC(i,j,K,bi,bj)ratioRp -half) )
448	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
449	jmc	1.14	& * maskC(i,j, K ,bi,bj)
450			ENDDO
451			ENDDO
452			ELSE
453			ratioRm=drF(K)*recip_drC(K)
454			ratioRp=drF(K)*recip_drC(Kp1)
455	jmc	1.24	ddPIm=atm_Cp( ((rC(K-1)/atm_Po)*atm_kappa)
456			& -((rC( K )/atm_Po)**atm_kappa) )
457			ddPIp=atm_Cp( ((rC( K )/atm_Po)*atm_kappa)
458			& -((rC(Kp1)/atm_Po)**atm_kappa) )
459	adcroft	1.9	DO j=jMin,jMax
460	jmc	1.14	DO i=iMin,iMax
461			phiHyd(i,j,K) = phiHyd(i,j,K-1)
462			& + ddPIm*0.5
463	mlosch	1.20	& *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
464	jmc	1.14	& * maskC(i,j,K-1,bi,bj)
465	mlosch	1.18	& +(ddPImmax(zero,(_hFacC(i,j,K,bi,bj)-one)ratioRm+half)
466			& +ddPIpmin(zero, _hFacC(i,j,K,bi,bj)ratioRp -half) )
467	mlosch	1.20	& *(tFld(i,j, K ,bi,bj)-tRef( K ))
468	jmc	1.14	& * maskC(i,j, K ,bi,bj)
469			ENDDO
470	adcroft	1.9	ENDDO
471			ENDIF
472	jmc	1.14	C end: Finite Difference Form, with hFac, Interface_W = middle --
473			C-----------------------------------------------------------------------
474	cnh	1.1
475	jmc	1.14	ELSE
476	jmc	1.24	STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
477	jmc	1.14	ENDIF
478	cnh	1.6
479	jmc	1.14	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
480	adcroft	1.9	ELSE
481	jmc	1.24	STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
482	jmc	1.25	ENDIF
483
484			IF (momPressureForcing) THEN
485			iMnLoc = MAX(1-Olx+1,iMin)
486			jMnLoc = MAX(1-Oly+1,jMin)
487			CALL CALC_GRAD_PHI_HYD(
488			I k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
489			I phiHyd, alphaRho, tFld, sFld,
490			O dPhiHydX, dPhiHydY,
491			I myTime, myIter, myThid)
492	cnh	1.5	ENDIF
493	cnh	1.1
494	jmc	1.14	#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
495	cnh	1.6
496	jmc	1.11	RETURN
497			END