C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.47 2016/03/10 20:55:56 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOP
C !ROUTINE: CALC_PHI_HYD
C !INTERFACE:
SUBROUTINE CALC_PHI_HYD(
I bi, bj, iMin, iMax, jMin, jMax, k,
I tFld, sFld,
U phiHydF,
O phiHydC, 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 | phiHydF(i,j) is the hydrostatic Potential anomaly
C | at middle between tracer points k-1,k
C | On exit:
C | phiHydC(i,j) is the hydrostatic Potential anomaly
C | at cell centers (tracer points), level k
C | phiHydF(i,j) is the hydrostatic Potential anomaly
C | at middle between tracer points k,k+1
C | dPhiHydX,Y hydrostatic Potential gradient (X&Y dir)
C | at cell centers (tracer points), level k
C | integr_GeoPot allows to select one integration method
C | 1= Finite volume form ; else= Finite difference form
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#ifdef ALLOW_AUTODIFF
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF */
#include "SURFACE.h"
#include "DYNVARS.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C bi, bj, k :: tile and level indices
C iMin,iMax,jMin,jMax :: computational domain
C tFld :: potential temperature
C sFld :: salinity
C phiHydF :: hydrostatic potential anomaly at middle between
C 2 centers (entry: Interf_k ; output: Interf_k+1)
C phiHydC :: hydrostatic potential anomaly at cell center
C dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
C myTime :: current time
C myIter :: current iteration number
C myThid :: thread number for this instance of the routine.
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 phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_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 ==
C phiHydU :: hydrostatic potential anomaly at interface k+1 (Upper k)
C pKappaF :: (p/Po)^kappa at interface k
C pKappaU :: (p/Po)^kappa at interface k+1 (Upper k)
C pKappaC :: (p/Po)^kappa at cell center k
INTEGER i,j
_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifndef DISABLE_SIGMA_CODE
_RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL pKappaC, locDepth, fullDepth
#endif /* DISABLE_SIGMA_CODE */
_RL thetaRef, locAlpha
_RL dRlocM,dRlocP, ddRloc
_RL ddPIm, ddPIp, rec_dRm, rec_dRp
_RL surfPhiFac
LOGICAL useDiagPhiRlow, addSurfPhiAnom
LOGICAL useFVgradPhi
CEOP
useDiagPhiRlow = .TRUE.
addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
useFVgradPhi = selectSigmaCoord.NE.0
surfPhiFac = 0.
IF (addSurfPhiAnom) surfPhiFac = 1.
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, accurate with Topo full cell;
C = 1 : Finite Volume Form, with Part-Cell, linear in P by Half level;
C =2,3: Finite Difference Form, with Part-Cell,
C linear in P between 2 Tracer levels.
C can handle both cases: Tracer lev at the middle of InterFace_W
C and InterFace_W at the middle of Tracer lev;
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-- Initialize phiHydF to zero :
C note: atmospheric_loading or Phi_topo anomaly are incorporated
C later in S/R calc_grad_phi_hyd
IF (k.EQ.1) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
phiHydF(i,j) = 0.
ENDDO
ENDDO
ENDIF
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 */
IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
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 & kind = isbyte
CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
CADJ & kind = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
CALL FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I tFld(1-OLx,1-OLy,k,bi,bj),
I sFld(1-OLx,1-OLy,k,bi,bj),
O alphaRho,
I k, bi, bj, myThid )
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_SHELFICE
C mask rho, so that there is no contribution of phiHyd from
C overlying shelfice (whose density we do not know)
IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
C- note: does not work for down_slope pkg which needs rho below the bottom.
C setting rho=0 above the ice-shelf base is enough (and works in both cases)
C but might be slower (--> keep original masking if not using down_slope pkg)
DO j=jMin,jMax
DO i=iMin,iMax
IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
ENDDO
ENDDO
ELSEIF ( useShelfIce ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* ALLOW_SHELFICE */
#ifdef ALLOW_MOM_COMMON
C-- Quasi-hydrostatic terms are added in as if they modify the buoyancy
IF (quasiHydrostatic) THEN
CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
ENDIF
#endif /* ALLOW_MOM_COMMON */
#ifdef NONLIN_FRSURF
IF ( addSurfPhiAnom .AND.
& uniformFreeSurfLev .AND. k.EQ.1 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
& *gravity*alphaRho(i,j)*recip_rhoConst
ENDDO
ENDDO
ENDIF
#endif /* NONLIN_FRSURF */
C---- Hydrostatic pressure at cell centers
IF (integr_GeoPot.EQ.1) THEN
C -- Finite Volume Form
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
IF ( uniformFreeSurfLev ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
phiHydC(i,j) = phiHydF(i,j)
& + halfRL*drF(k)*gravFacC(k)*gravity
& *alphaRho(i,j)*recip_rhoConst
phiHydF(i,j) = phiHydF(i,j)
& + drF(k)*gravFacC(k)*gravity
& *alphaRho(i,j)*recip_rhoConst
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
#ifdef NONLIN_FRSURF
ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
#endif
phiHydC(i,j) = ddRloc*gravFacC(k)*gravity
& *alphaRho(i,j)*recip_rhoConst
ELSE
phiHydC(i,j) = phiHydF(i,j)
& + halfRL*drF(k)*gravFacC(k)*gravity
& *alphaRho(i,j)*recip_rhoConst
ENDIF
phiHydF(i,j) = phiHydC(i,j)
& + halfRL*drF(k)*gravFacC(k)*gravity
& *alphaRho(i,j)*recip_rhoConst
ENDDO
ENDDO
ENDIF
ELSE
C -- Finite Difference Form
C---------- This discretization is the "energy conserving" form
C which has been used since at least Adcroft et al., MWR 1997
dRlocM = halfRL*drC(k)*gravFacF(k)
IF (k.EQ.1) dRlocM = (rF(k)-rC(k))*gravFacF(k)
IF (k.EQ.Nr) THEN
dRlocP = (rC(k)-rF(k+1))*gravFacF(k+1)
ELSE
dRlocP = halfRL*drC(k+1)*gravFacF(k+1)
ENDIF
IF ( uniformFreeSurfLev ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
phiHydC(i,j) = phiHydF(i,j)
& + dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
phiHydF(i,j) = phiHydC(i,j)
& + dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
ENDDO
ENDDO
ELSE
rec_dRm = oneRL/(rF(k)-rC(k))
rec_dRp = oneRL/(rC(k)-rF(k+1))
DO j=jMin,jMax
DO i=iMin,iMax
IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
#ifdef NONLIN_FRSURF
ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
#endif
phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
& +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
& )*gravity*alphaRho(i,j)*recip_rhoConst
ELSE
phiHydC(i,j) = phiHydF(i,j)
& + dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
ENDIF
phiHydF(i,j) = phiHydC(i,j)
& + dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
ENDDO
ENDDO
ENDIF
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 before
C integrating (1/rho)_prime*dp over the current layer/interface
#ifdef ALLOW_AUTODIFF_TAMC
CADJ GENERAL
#endif /* ALLOW_AUTODIFF_TAMC */
IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
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 & kind = isbyte
CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
CADJ & kind = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
CALL FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I tFld(1-OLx,1-OLy,k,bi,bj),
I sFld(1-OLx,1-OLy,k,bi,bj),
O alphaRho,
I k, bi, bj, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte,
CADJ & kind = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
ENDDO
ENDDO
ENDIF
C-- Calculate specific volume anomaly : alpha_prime = 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)*
& (oneRL/locAlpha - recip_rhoConst)
ENDDO
ENDDO
#ifdef ALLOW_MOM_COMMON
C-- Quasi-hydrostatic terms are added as if they modify the specific-volume
IF (quasiHydrostatic) THEN
CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
ENDIF
#endif /* ALLOW_MOM_COMMON */
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
IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
#ifdef NONLIN_FRSURF
ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
#endif
phiHydC(i,j) = ddRloc*alphaRho(i,j)
c--to reproduce results of c48d_post: uncomment those 4+1 lines
c phiHydC(i,j)=phiHydF(i,j)
c & +(hFacC(i,j,k,bi,bj)-halfRL)*drF(k)*alphaRho(i,j)
c phiHydF(i,j)=phiHydF(i,j)
c & + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
ELSE
phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
c phiHydF(i,j) = phiHydF(i,j) + drF(k)*alphaRho(i,j)
ENDIF
c-- and comment this last one:
phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
c-----
ENDDO
ENDDO
ELSE
C -- Finite Difference Form, with Part-Cell Bathy
dRlocM = halfRL*drC(k)
IF (k.EQ.1) dRlocM=rF(k)-rC(k)
IF (k.EQ.Nr) THEN
dRlocP=rC(k)-rF(k+1)
ELSE
dRlocP=halfRL*drC(k+1)
ENDIF
rec_dRm = oneRL/(rF(k)-rC(k))
rec_dRp = oneRL/(rC(k)-rF(k+1))
DO j=jMin,jMax
DO i=iMin,iMax
C---------- This discretization is the "energy conserving" form
IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
#ifdef NONLIN_FRSURF
ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
#endif
phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
& +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
& )*alphaRho(i,j)
ELSE
phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
ENDIF
phiHydF(i,j) = phiHydC(i,j) + dRlocP*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.
IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
C-- virtual potential temperature anomaly (including water vapour effect)
IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
C- isothermal (theta=const) reference state
thetaRef = thetaConst
ELSE
C- horizontally uniform (tRef) reference state
thetaRef = tRef(k)
ENDIF
DO j=jMin,jMax
DO i=iMin,iMax
alphaRho(i,j) = ( tFld(i,j,k,bi,bj)
& *( sFld(i,j,k,bi,bj)*atm_Rq + oneRL )
& - thetaRef )*maskC(i,j,k,bi,bj)
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_MOM_COMMON
C-- Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
IF (quasiHydrostatic) THEN
CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
ENDIF
#endif /* ALLOW_MOM_COMMON */
C--- Integrate d Phi / d pi
#ifndef DISABLE_SIGMA_CODE
C -- Finite Volume Form, integrated to r-level (cell center & upper interface)
IF ( useFVgradPhi ) THEN
fullDepth = rF(1)-rF(Nr+1)
DO j=jMin,jMax
DO i=iMin,iMax
locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
#ifdef NONLIN_FRSURF
locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
#endif
pKappaF(i,j) = (
& ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
& + bHybSigmF( k )*locDepth
& )/atm_Po )**atm_kappa
pKappaC = (
& ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
& + bHybSigmC( k )*locDepth
& )/atm_Po )**atm_kappa
pKappaU(i,j) = (
& ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
& + bHybSigmF(k+1)*locDepth
& )/atm_Po )**atm_kappa
phiHydC(i,j) = phiHydF(i,j)
& + atm_Cp*( pKappaF(i,j) - pKappaC )*alphaRho(i,j)
phiHydU(i,j) = phiHydF(i,j)
& + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
ENDDO
ENDDO
C end: Finite Volume Form, integrated to r-level.
ELSEIF (integr_GeoPot.EQ.0) THEN
#else /* DISABLE_SIGMA_CODE */
IF (integr_GeoPot.EQ.0) THEN
#endif /* DISABLE_SIGMA_CODE */
C -- Energy Conserving Form, accurate with Full cell topo --
C------------ The integration for the first level phi(k=1) is the same
C for both the "finite volume" and energy conserving methods.
C *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
ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
& -((rC( k )/atm_Po)**atm_kappa) )
ELSE
ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
& -((rC( k )/atm_Po)**atm_kappa) )*halfRL
ENDIF
IF (k.EQ.Nr) THEN
ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rF(k+1)/atm_Po)**atm_kappa) )
ELSE
ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
ENDIF
C-------- This discretization is the energy conserving form
DO j=jMin,jMax
DO i=iMin,iMax
phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
ENDDO
ENDDO
C end: Energy Conserving Form, No hFac --
C-----------------------------------------------------------------------
ELSEIF (integr_GeoPot.EQ.1) THEN
C -- Finite Volume Form, with Part-Cell Topo, 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---------
ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
& -((rC( k )/atm_Po)**atm_kappa) )
ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rF(k+1)/atm_Po)**atm_kappa) )
DO j=jMin,jMax
DO i=iMin,iMax
IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
#ifdef NONLIN_FRSURF
ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
#endif
phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
& *ddPIm*alphaRho(i,j)
ELSE
phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
ENDIF
phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
ENDDO
ENDDO
C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
C-----------------------------------------------------------------------
ELSEIF ( integr_GeoPot.EQ.2
& .OR. integr_GeoPot.EQ.3 ) THEN
C -- Finite Difference Form, with Part-Cell Topo,
C works with Interface_W at the middle between 2.Tracer_Level
C and with Tracer_Level at the middle between 2.Interface_W.
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---------
IF (k.EQ.1) THEN
ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
& -((rC( k )/atm_Po)**atm_kappa) )
ELSE
ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
& -((rC( k )/atm_Po)**atm_kappa) )*halfRL
ENDIF
IF (k.EQ.Nr) THEN
ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rF(k+1)/atm_Po)**atm_kappa) )
ELSE
ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
ENDIF
rec_dRm = oneRL/(rF(k)-rC(k))
rec_dRp = oneRL/(rC(k)-rF(k+1))
DO j=jMin,jMax
DO i=iMin,iMax
IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
#ifdef NONLIN_FRSURF
ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
#endif
phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
& +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
& )*alphaRho(i,j)
ELSE
phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
ENDIF
phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
ENDDO
ENDDO
C end: Finite Difference Form, with Part-Cell Topo
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 ( .NOT. useFVgradPhi ) THEN
C-- r-coordinate and r*-coordinate cases:
IF ( momPressureForcing ) THEN
CALL CALC_GRAD_PHI_HYD(
I k, bi, bj, iMin,iMax, jMin,jMax,
I phiHydC, alphaRho, tFld, sFld,
O dPhiHydX, dPhiHydY,
I myTime, myIter, myThid)
ENDIF
#ifndef DISABLE_SIGMA_CODE
ELSE
C-- else (SigmaCoords part)
IF ( fluidIsWater ) THEN
STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
ENDIF
IF ( momPressureForcing ) THEN
CALL CALC_GRAD_PHI_FV(
I k, bi, bj, iMin,iMax, jMin,jMax,
I phiHydF, phiHydU, pKappaF, pKappaU,
O dPhiHydX, dPhiHydY,
I myTime, myIter, myThid)
ENDIF
DO j=jMin,jMax
DO i=iMin,iMax
phiHydF(i,j) = phiHydU(i,j)
ENDDO
ENDDO
#endif /* DISABLE_SIGMA_CODE */
C-- end if-not/else useFVgradPhi
ENDIF
C--- Diagnose Phi at boundary r=R_low :
C = Ocean bottom pressure (Ocean, Z-coord.)
C = Sea-surface height (Ocean, P-coord.)
C = Top atmosphere height (Atmos, P-coord.)
IF (useDiagPhiRlow) THEN
CALL DIAGS_PHI_RLOW(
I k, bi, bj, iMin,iMax, jMin,jMax,
I phiHydF, phiHydC, alphaRho, tFld, sFld,
I myTime, myIter, myThid)
ENDIF
C--- Diagnose Full Hydrostatic Potential at cell center level
CALL DIAGS_PHI_HYD(
I k, bi, bj, iMin,iMax, jMin,jMax,
I phiHydC,
I myTime, myIter, myThid)
#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
RETURN
END