model/src/calc_phi_hyd.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.19 2002/08/15 17:25:31 adcroft 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,
     I                         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"
#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)
      INTEGER 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
CEOP

      zero = 0. _d 0
      one  = 1. _d 0
      half = .5 _d 0

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 */

      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

        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
#ifdef ATMOSPHERIC_LOADING
              phiHyd(i,j,k)=pload(i,j,bi,bj)*recip_rhoConst
#else
              phiHyd(i,j,k)=0. _d 0
#endif
            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       Hydrostatic pressure at cell centers
        DO j=jMin,jMax
          DO i=iMin,iMax
#ifdef      ALLOW_AUTODIFF_TAMC
c           Patrick, is this directive correct or even necessary in
c           this new code?
c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...
c           within the k-loop.
CADJ GENERAL
#endif      /* ALLOW_AUTODIFF_TAMC */

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
c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
c    &              drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
c           phiHyd(i,j,k)=phiHyd(i,j,k)+
c    &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
C-----------------------------------------------------------------------

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)+
     &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
     &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
C-----------------------------------------------------------------------

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)-0.5)*drF(K)
     &                   *gravity*alphaRho(i,j)*recip_rhoConst
     &              + gravity*etaN(i,j,bi,bj)
            ENDIF
C-----------------------------------------------------------------------

          ENDDO
        ENDDO
        
      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 g*rho over the current layer/interface

        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)=0.
              phiHyd(i,j,k)=pload(i,j,bi,bj)
c    &  -Ro_surf(i,j,bi,bj)*recip_rhoNil
c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSurfC(i,j,bi,bj) ))/1000.
c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSurfC(i,j,bi,bj) ))*recip_rhoNil
            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       Hydrostatic pressure at cell centers
        DO j=jMin,jMax
          DO i=iMin,iMax
            locAlpha=alphaRho(i,j)+rhoNil
            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
c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
c    &              drF(K)*locAlpha
c           phiHyd(i,j,k)=phiHyd(i,j,k)+
c    &          0.5*drF(K)*locAlpha
C-----------------------------------------------------------------------

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)+
     &          0.5*dRloc*locAlpha
            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
     &          0.5*dRlocKp1*locAlpha
C-----------------------------------------------------------------------

C---------- Compute gravity*(sea surface elevation)
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)-0.5)*drF(k)*locAlpha
     &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
            ENDIF
C-----------------------------------------------------------------------

          ENDDO
        ENDDO

      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
             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
             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
             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
             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: We should never reach this point!'
      ENDIF

#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */

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