--- MITgcm/model/src/calc_phi_hyd.F	2000/11/13 16:32:57	1.8
+++ MITgcm/model/src/calc_phi_hyd.F	2001/01/12 21:02:46	1.8.2.1
@@ -1,64 +1,119 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/calc_phi_hyd.F,v 1.8 2000/11/13 16:32:57 heimbach Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/calc_phi_hyd.F,v 1.8.2.1 2001/01/12 21:02:46 adcroft Exp $
 
 #include "CPP_OPTIONS.h"
 
-      SUBROUTINE CALC_PHI_HYD( bi, bj, iMin, iMax, jMin, jMax, K,
-     I    buoyKM1, buoyKP1, phiHyd, myThid)
+      SUBROUTINE CALC_PHI_HYD(
+     I                         bi, bj, iMin, iMax, jMin, jMax, K,
+     I                         theta, salt,
+     U                         phiHyd, phiHydInterface,
+     I                         myThid)
 C     /==========================================================\
 C     | SUBROUTINE CALC_PHI_HYD                                  |
 C     | o Integrate the hydrostatic relation to find phiHyd.     |
 C     |                                                          |
+C     | On entry:                                                |
+C     |   theta,salt    are the current thermodynamics quantities|
+C     |                 (unchanged on exit)                      |
+C     |   phiHyd(i,j,1:k-1) is the hydrostatic pressure/geopot.  |
+C     |                 at cell centers (tracer points)          |
+C     |                 - 1:k-1 layers are valid                 |
+C     |                 - k:Nr layers are invalid                |
+C     |   phiHydInterface(i,j) is the hydrostatic pressure/geop. |
+C     |                 at cell the interface k (w point above)  |
+C     | On exit:                                                 |
+C     |   phiHyd(i,j,1:k) is the hydrostatic pressure/geopot.    |
+C     |                 at cell centers (tracer points)          |
+C     |                 - 1:k layers are valid                   |
+C     |                 - k+1:Nr layers are invalid              |
+C     |   phiHydInterface(i,j) is the hydrostatic pressure/geop. |
+C     |                 at cell the interface k+1 (w point below)|
+C     |                                                          |
 C     \==========================================================/
       IMPLICIT NONE
 C     == Global variables ==
 #include "SIZE.h"
-#include "DYNVARS.h"
 #include "GRID.h"
 #include "EEPARAMS.h"
 #include "PARAMS.h"
 C     == Routine arguments ==
       INTEGER bi,bj,iMin,iMax,jMin,jMax,K
-      _RL buoyKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
-      _RL buoyKP1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
+      _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
-      integer myThid
-C     == Local variables ==
-      INTEGER i,j,Km1
-      _RL halfLayer
-      _RL gamma
+      _RL phiHydInterface(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      INTEGER myThid
 
 #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
 
-      if (K.eq.1) then
-       Km1=1
-       halfLayer=0.5 _d 0
-      else
-       Km1=K-1
-       halfLayer=1.0 _d 0
-      endif
-
-C--   Scale factor for hydrostatic relation except for ocean in 
-C--   pressure coords.
-      gamma = 1. _d 0
-C--   Scale factor for hydrostatic relation for ocean in pressure 
-C--   coords.
-      IF ( buoyancyRelation .EQ. 'OCEANIC' .AND. usingPCoords ) THEN
-       gamma = recip_Gravity*recip_rhoConst
-      ENDIF
-
-C--   Contribution to phiHyd(:,:,K) from buoy(:,:,K-1) + buoy(:,:,K)
-C     (This is now the actual hydrostatic pressure|height at the T/S 
-C     points)
-      DO j=jMin,jMax
-       DO i=iMin,iMax
+C     == Local variables ==
+      INTEGER i,j
+      _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL dRloc,dRlocKp1
+
+      IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
+
+        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             *NOTE* The loading should go here but has not been implemented yet
+              phiHydInterface(i,j)=0.
+              phiHyd(i,j,k)=0.
+            ENDDO
+          ENDDO
+        ENDIF
+
+C       Calculate density
+        CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
+     &                 theta, salt,
+     &                 alphaRho, myThid)
+
+C       Hydrostatic pressure at cell centers
+        DO j=jMin,jMax
+          DO i=iMin,iMax
 #ifdef ALLOW_AUTODIFF_TAMC
+          Is this directive correct or even necessary in this new code?
 CADJ GENERAL
 #endif
-        phiHyd(i,j,K)=phiHyd(i,j,Km1)-rhoConst*halfLayer
-     &  *0.5 _d 0*( drF(Km1)+drF(K) )*gamma
-     &  *0.5 _d 0*( buoyKM1(i,j)+buoyKP1(i,j) )
-       ENDDO
-      ENDDO
+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           phiHyd(i,j,k)=phiHydInterface(i,j)+
+C    &          0.5*drF(K)*gravity*alphaRho(i,j)
+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)
+            phiHyd(i,j,k+1)=phiHyd(i,j,k)+
+     &          0.5*dRlocKp1*gravity*alphaRho(i,j)
+          ENDDO
+        ENDDO
+        
+C       Hydrostatic pressure at interface below
+        DO j=jMin,jMax
+          DO i=iMin,iMax
+            phiHydInterface(i,j)=phiHydInterface(i,j)+
+     &              drF(K)*gravity*alphaRho(i,j)
+          ENDDO
+        ENDDO
+
+      ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
+
+      ELSE
+        STOP 'CALC_PHI_HYD: We should never reach this point!'
+      ENDIF
 
 #endif