C $Header: /home/ubuntu/mnt/e9_copy/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