model/src/diags_phi_rlow.F

C $Header: /u/gcmpack/MITgcm/model/src/diags_phi_rlow.F,v 1.11 2014/05/02 16:57:43 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
1	jmc	1.12	C $Header: /u/gcmpack/MITgcm/model/src/diags_phi_rlow.F,v 1.11 2014/05/02 16:57:43 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "CPP_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: DIAGS_PHI_RLOW
8			C !INTERFACE:
9	jmc	1.5	SUBROUTINE DIAGS_PHI_RLOW(
10	jmc	1.1	I k, bi, bj, iMin,iMax, jMin,jMax,
11	jmc	1.2	I phiHydF, phiHydC, alphRho, tFld, sFld,
12	jmc	1.1	I myTime, myIter, myThid)
13			C !DESCRIPTION: \bv
14			C ==========================================================
15	jmc	1.5	C \| S/R DIAGS_PHI_RLOW
16	jmc	1.1	C \| o Diagnose Phi-Hydrostatic at r-lower boundary
17	jmc	1.5	C \| = bottom pressure (ocean in z-coord) ;
18	jmc	1.1	C \| = sea surface elevation (ocean in p-coord) ;
19			C \| = height at the top of atmosphere (in p-coord) ;
20			C ==========================================================
21			C \ev
22
23			C !USES:
24			IMPLICIT NONE
25			C == Global variables ==
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30			#include "SURFACE.h"
31			#include "DYNVARS.h"
32
33			C !INPUT/OUTPUT PARAMETERS:
34			C == Routine Arguments ==
35	jmc	1.5	C bi,bj :: tile index
36	jmc	1.1	C iMin,iMax,jMin,jMax :: Loop counters
37	jmc	1.5	C phiHydF :: hydrostatic potential anomaly at middle between
38	jmc	1.2	C 2 centers k & k+1 (interface k+1)
39			C phiHydC :: hydrostatic potential anomaly at cell center
40	jmc	1.1	C (atmos: =Geopotential ; ocean-z: =Pressure/rho)
41			C alphRho :: Density (z-coord) or specific volume (p-coord)
42			C tFld :: Potential temp.
43	jmc	1.5	C sFld :: Salinity
44			C myTime :: Current time
45			C myIter :: Current iteration number
46			C myThid :: my Thread Id number
47	jmc	1.1	INTEGER k, bi,bj, iMin,iMax, jMin,jMax
48	jmc	1.2	_RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49			_RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50	jmc	1.1	_RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51			_RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
52			_RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
53			_RL myTime
54			INTEGER myIter, myThid
55
56			#ifdef INCLUDE_PHIHYD_CALCULATION_CODE
57
58			C !LOCAL VARIABLES:
59			C == Local variables ==
60			C i,j :: Loop counters
61			INTEGER i,j
62	jmc	1.2	_RL ddRloc, ratioRm, ratioRp
63	jmc	1.9	#ifdef NONLIN_FRSURF
64			_RL facP, dPhiRef
65			#endif /* NONLIN_FRSURF */
66	jmc	1.1	CEOP
67
68	jmc	1.6	IF ( usingZCoords ) THEN
69	jmc	1.1
70	jmc	1.2	C----- Compute bottom pressure deviation from gravityrho0H
71	jmc	1.7	C Start from phiHyd at the (bottom) tracer point and add Del_hgrho_prime
72	jmc	1.2	C with Del_h = distance from the bottom up to tracer point
73	jmc	1.1
74	jmc	1.6	C-- Initialise to zero (otherwise phi0surf accumulate over land)
75			IF ( k.EQ.1 ) THEN
76	jmc	1.10	DO j=1-OLy,sNy+OLy
77			DO i=1-OLx,sNx+OLx
78	jmc	1.6	phiHydLow(i,j,bi,bj) = 0. _d 0
79			ENDDO
80			ENDDO
81			ENDIF
82
83	jmc	1.1	IF (integr_GeoPot.EQ.1) THEN
84			C -- Finite Volume Form
85
86			DO j=jMin,jMax
87			DO i=iMin,iMax
88			IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
89	jmc	1.2	ddRloc = rC(k)-R_low(i,j,bi,bj)
90			phiHydLow(i,j,bi,bj) = phiHydC(i,j)
91	jmc	1.12	& + ddRlocgravFacC(k)gravityalphRho(i,j)recip_rhoConst
92	jmc	1.1	ENDIF
93			ENDDO
94			ENDDO
95
96			ELSE
97			C -- Finite Difference Form
98
99	jmc	1.10	ratioRm = oneRL
100			ratioRp = oneRL
101			IF (k.GT.1 ) ratioRm = halfRL*drC(k)/(rF(k)-rC(k))
102			IF (k.LT.Nr) ratioRp = halfRL*drC(k+1)/(rC(k)-rF(k+1))
103	jmc	1.12	ratioRm = ratioRm*gravFacF(k)
104			ratioRp = ratioRp*gravFacF(k+1)
105	jmc	1.1
106			DO j=jMin,jMax
107			DO i=iMin,iMax
108	jmc	1.2	IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
109			ddRloc = rC(k)-R_low(i,j,bi,bj)
110			phiHydLow(i,j,bi,bj) = phiHydC(i,j)
111	jmc	1.10	& +( MIN(zeroRL,ddRloc)*ratioRm
112			& +MAX(zeroRL,ddRloc)*ratioRp
113	jmc	1.2	& )gravityalphRho(i,j)*recip_rhoConst
114	jmc	1.1	ENDIF
115			ENDDO
116			ENDDO
117
118			C -- end if integr_GeoPot = ...
119			ENDIF
120	jmc	1.5
121	jmc	1.6	C -- end usingZCoords
122	jmc	1.2	ENDIF
123
124	jmc	1.6	IF ( k.EQ.Nr ) THEN
125	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
126	jmc	1.2	C -- last level (bottom): rescale (r*) and add surface contribution
127	jmc	1.1
128	jmc	1.6	IF ( usingPCoords ) THEN
129	jmc	1.2	C -- P coordinate : Phi(R_low) is simply at the top :
130			DO j=jMin,jMax
131			DO i=iMin,iMax
132			phiHydLow(i,j,bi,bj) = phiHydF(i,j)
133	jmc	1.1	ENDDO
134	jmc	1.2	ENDDO
135			ENDIF
136	jmc	1.1
137	jmc	1.5	#ifdef NONLIN_FRSURF
138			c IF ( select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
139	jmc	1.8	IF ( select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
140			C- Integral of b.dr = rStarFac * Integral of b.dr* :
141	jmc	1.9	IF ( fluidIsAir ) THEN
142	jmc	1.8	C- Consistent with Phi'= Integr[ theta'.dPi ] :
143			DO j=jMin,jMax
144			DO i=iMin,iMax
145	jmc	1.11	facP = pStarFacK(i,j,bi,bj)
146	jmc	1.9	dPhiRef = phiRef(2k+1) - gravitytopoZ(i,j,bi,bj)
147			& - phi0surf(i,j,bi,bj)
148			phiHydLow(i,j,bi,bj) =
149			& phiHydLow(i,j,bi,bj)*facP
150	jmc	1.10	& + MAX( dPhiRef, zeroRL )*( facP - 1. _d 0 )
151	jmc	1.9	& + phi0surf(i,j,bi,bj)
152			ENDDO
153			ENDDO
154			ELSEIF ( usingPCoords ) THEN
155			C- Note: dPhiRef(rStarFacC -1) = 1/rhoPSo((eta+PSo)/PSo -1) = Bo_surfetaN
156			C so that this expression is the same as before (same as in "ELSE" block below)
157			DO j=jMin,jMax
158			DO i=iMin,iMax
159			dPhiRef = ( Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) )
160			& *recip_rhoConst
161			phiHydLow(i,j,bi,bj) =
162			& phiHydLow(i,j,bi,bj)*rStarFacC(i,j,bi,bj)
163			& + dPhiRef*( rStarFacC(i,j,bi,bj) - 1. _d 0 )
164			& + phi0surf(i,j,bi,bj)
165	jmc	1.8	ENDDO
166			ENDDO
167			ELSE
168	jmc	1.9	C- Note: dPhiRef(rStarFacC -1) = gH((eta+H)/H -1) = Bo_surfetaN
169			C so that this expression is the same as before (same as in "ELSE" block below)
170	jmc	1.8	DO j=jMin,jMax
171			DO i=iMin,iMax
172	jmc	1.9	dPhiRef = ( Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj) )
173			& *gravity
174			phiHydLow(i,j,bi,bj) =
175			& phiHydLow(i,j,bi,bj)*rStarFacC(i,j,bi,bj)
176			& + dPhiRef*( rStarFacC(i,j,bi,bj) - 1. _d 0 )
177			& + phi0surf(i,j,bi,bj)
178	jmc	1.8	ENDDO
179	jmc	1.5	ENDDO
180	jmc	1.8	ENDIF
181	jmc	1.9	ELSE
182			#else /* NONLIN_FRSURF */
183			IF ( .TRUE. ) THEN
184	jmc	1.5	#endif /* NONLIN_FRSURF */
185
186	jmc	1.9	DO j=jMin,jMax
187			DO i=iMin,iMax
188	jmc	1.2	phiHydLow(i,j,bi,bj) = phiHydLow(i,j,bi,bj)
189			& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
190	jmc	1.5	& + phi0surf(i,j,bi,bj)
191	jmc	1.9	ENDDO
192			ENDDO
193			ENDIF
194	jmc	1.1
195			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
196	jmc	1.2	C -- end if k=Nr.
197	jmc	1.1	ENDIF
198
199			#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
200
201			RETURN
202			END