C $Header: /u/gcmpack/MITgcm/model/src/diags_phi_rlow.F,v 1.12 2016/03/10 20:55:56 jmc Exp $
C $Name: $
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: DIAGS_PHI_RLOW
C !INTERFACE:
SUBROUTINE DIAGS_PHI_RLOW(
I k, bi, bj, iMin,iMax, jMin,jMax,
I phiHydF, phiHydC, alphRho, tFld, sFld,
I myTime, myIter, myThid)
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R DIAGS_PHI_RLOW
C | o Diagnose Phi-Hydrostatic at r-lower boundary
C | = bottom pressure (ocean in z-coord) ;
C | = sea surface elevation (ocean in p-coord) ;
C | = height at the top of atmosphere (in p-coord) ;
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C bi,bj :: tile index
C iMin,iMax,jMin,jMax :: Loop counters
C phiHydF :: hydrostatic potential anomaly at middle between
C 2 centers k & k+1 (interface k+1)
C phiHydC :: hydrostatic potential anomaly at cell center
C (atmos: =Geopotential ; ocean-z: =Pressure/rho)
C alphRho :: Density (z-coord) or specific volume (p-coord)
C tFld :: Potential temp.
C sFld :: Salinity
C myTime :: Current time
C myIter :: Current iteration number
C myThid :: my Thread Id number
INTEGER k, bi,bj, iMin,iMax, jMin,jMax
_RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_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 myTime
INTEGER myIter, myThid
#ifdef INCLUDE_PHIHYD_CALCULATION_CODE
C !LOCAL VARIABLES:
C == Local variables ==
C i,j :: Loop counters
INTEGER i,j
_RL ddRloc, ratioRm, ratioRp
#ifdef NONLIN_FRSURF
_RL facP, dPhiRef
#endif /* NONLIN_FRSURF */
CEOP
IF ( usingZCoords ) THEN
C----- Compute bottom pressure deviation from gravity*rho0*H
C Start from phiHyd at the (bottom) tracer point and add Del_h*g*rho_prime
C with Del_h = distance from the bottom up to tracer point
C-- Initialise to zero (otherwise phi0surf accumulate over land)
IF ( k.EQ.1 ) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
phiHydLow(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDIF
IF (integr_GeoPot.EQ.1) THEN
C -- Finite Volume Form
DO j=jMin,jMax
DO i=iMin,iMax
IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
ddRloc = rC(k)-R_low(i,j,bi,bj)
phiHydLow(i,j,bi,bj) = phiHydC(i,j)
& + ddRloc*gravFacC(k)*gravity*alphRho(i,j)*recip_rhoConst
ENDIF
ENDDO
ENDDO
ELSE
C -- Finite Difference Form
ratioRm = oneRL
ratioRp = oneRL
IF (k.GT.1 ) ratioRm = halfRL*drC(k)/(rF(k)-rC(k))
IF (k.LT.Nr) ratioRp = halfRL*drC(k+1)/(rC(k)-rF(k+1))
ratioRm = ratioRm*gravFacF(k)
ratioRp = ratioRp*gravFacF(k+1)
DO j=jMin,jMax
DO i=iMin,iMax
IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
ddRloc = rC(k)-R_low(i,j,bi,bj)
phiHydLow(i,j,bi,bj) = phiHydC(i,j)
& +( MIN(zeroRL,ddRloc)*ratioRm
& +MAX(zeroRL,ddRloc)*ratioRp
& )*gravity*alphRho(i,j)*recip_rhoConst
ENDIF
ENDDO
ENDDO
C -- end if integr_GeoPot = ...
ENDIF
C -- end usingZCoords
ENDIF
IF ( k.EQ.Nr ) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C -- last level (bottom): rescale (r*) and add surface contribution
IF ( usingPCoords ) THEN
C -- P coordinate : Phi(R_low) is simply at the top :
DO j=jMin,jMax
DO i=iMin,iMax
phiHydLow(i,j,bi,bj) = phiHydF(i,j)
ENDDO
ENDDO
ENDIF
#ifdef NONLIN_FRSURF
c IF ( select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
IF ( select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
C- Integral of b.dr = rStarFac * Integral of b.dr* :
IF ( fluidIsAir ) THEN
C- Consistent with Phi'= Integr[ theta'.dPi ] :
DO j=jMin,jMax
DO i=iMin,iMax
facP = pStarFacK(i,j,bi,bj)
dPhiRef = phiRef(2*k+1) - gravity*topoZ(i,j,bi,bj)
& - phi0surf(i,j,bi,bj)
phiHydLow(i,j,bi,bj) =
& phiHydLow(i,j,bi,bj)*facP
& + MAX( dPhiRef, zeroRL )*( facP - 1. _d 0 )
& + phi0surf(i,j,bi,bj)
ENDDO
ENDDO
ELSEIF ( usingPCoords ) THEN
C- Note: dPhiRef*(rStarFacC -1) = 1/rho*PSo*((eta+PSo)/PSo -1) = Bo_surf*etaN
C so that this expression is the same as before (same as in "ELSE" block below)
DO j=jMin,jMax
DO i=iMin,iMax
dPhiRef = ( Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) )
& *recip_rhoConst
phiHydLow(i,j,bi,bj) =
& phiHydLow(i,j,bi,bj)*rStarFacC(i,j,bi,bj)
& + dPhiRef*( rStarFacC(i,j,bi,bj) - 1. _d 0 )
& + phi0surf(i,j,bi,bj)
ENDDO
ENDDO
ELSE
C- Note: dPhiRef*(rStarFacC -1) = g*H*((eta+H)/H -1) = Bo_surf*etaN
C so that this expression is the same as before (same as in "ELSE" block below)
DO j=jMin,jMax
DO i=iMin,iMax
dPhiRef = ( Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj) )
& *gravity
phiHydLow(i,j,bi,bj) =
& phiHydLow(i,j,bi,bj)*rStarFacC(i,j,bi,bj)
& + dPhiRef*( rStarFacC(i,j,bi,bj) - 1. _d 0 )
& + phi0surf(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
ELSE
#else /* NONLIN_FRSURF */
IF ( .TRUE. ) THEN
#endif /* NONLIN_FRSURF */
DO j=jMin,jMax
DO i=iMin,iMax
phiHydLow(i,j,bi,bj) = phiHydLow(i,j,bi,bj)
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
& + phi0surf(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C -- end if k=Nr.
ENDIF
#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
RETURN
END