/[MITgcm]/MITgcm/model/src/diags_phi_rlow.F
ViewVC logotype

Contents of /MITgcm/model/src/diags_phi_rlow.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.9 - (show annotations) (download)
Tue Jul 19 23:49:34 2011 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63g, checkpoint64, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64a, checkpoint64b
Changes since 1.8: +41 -11 lines 
- include r* effect on reference-state (PhiRef) in totPhiHyd
  (so that it's really the potential anomaly at the cell center);
  affect solution only if z* with full pressure in EOS.
- keep the previous diagnostics as "PHIHYDcR" (closer to potential anomaly
  @ constant r) until interpolation at constant r is implemented.
1 C $Header: /u/gcmpack/MITgcm/model/src/diags_phi_rlow.F,v 1.8 2011/06/10 20:55:51 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: DIAGS_PHI_RLOW
8 C !INTERFACE:
9 SUBROUTINE DIAGS_PHI_RLOW(
10 I k, bi, bj, iMin,iMax, jMin,jMax,
11 I phiHydF, phiHydC, alphRho, tFld, sFld,
12 I myTime, myIter, myThid)
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | S/R DIAGS_PHI_RLOW
16 C | o Diagnose Phi-Hydrostatic at r-lower boundary
17 C | = bottom pressure (ocean in z-coord) ;
18 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 C bi,bj :: tile index
36 C iMin,iMax,jMin,jMax :: Loop counters
37 C phiHydF :: hydrostatic potential anomaly at middle between
38 C 2 centers k & k+1 (interface k+1)
39 C phiHydC :: hydrostatic potential anomaly at cell center
40 C (atmos: =Geopotential ; ocean-z: =Pressure/rho)
41 C alphRho :: Density (z-coord) or specific volume (p-coord)
42 C tFld :: Potential temp.
43 C sFld :: Salinity
44 C myTime :: Current time
45 C myIter :: Current iteration number
46 C myThid :: my Thread Id number
47 INTEGER k, bi,bj, iMin,iMax, jMin,jMax
48 _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49 _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50 _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 _RL zero, one, half
63 _RL ddRloc, ratioRm, ratioRp
64 PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
65 #ifdef NONLIN_FRSURF
66 _RL facP, dPhiRef
67 #endif /* NONLIN_FRSURF */
68 CEOP
69
70 IF ( usingZCoords ) THEN
71
72 C----- Compute bottom pressure deviation from gravity*rho0*H
73 C Start from phiHyd at the (bottom) tracer point and add Del_h*g*rho_prime
74 C with Del_h = distance from the bottom up to tracer point
75
76 C-- Initialise to zero (otherwise phi0surf accumulate over land)
77 IF ( k.EQ.1 ) THEN
78 DO j=1-Oly,sNy+Oly
79 DO i=1-Olx,sNx+Olx
80 phiHydLow(i,j,bi,bj) = 0. _d 0
81 ENDDO
82 ENDDO
83 ENDIF
84
85 IF (integr_GeoPot.EQ.1) THEN
86 C -- Finite Volume Form
87
88 DO j=jMin,jMax
89 DO i=iMin,iMax
90 IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
91 ddRloc = rC(k)-R_low(i,j,bi,bj)
92 phiHydLow(i,j,bi,bj) = phiHydC(i,j)
93 & + ddRloc*gravity*alphRho(i,j)*recip_rhoConst
94 ENDIF
95 ENDDO
96 ENDDO
97
98 ELSE
99 C -- Finite Difference Form
100
101 ratioRm = one
102 ratioRp = one
103 IF (k.GT.1 ) ratioRm = half*drC(k)/(rF(k)-rC(k))
104 IF (k.LT.Nr) ratioRp = half*drC(k+1)/(rC(k)-rF(k+1))
105
106 DO j=jMin,jMax
107 DO i=iMin,iMax
108 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 & +( MIN(zero,ddRloc)*ratioRm
112 & +MAX(zero,ddRloc)*ratioRp
113 & )*gravity*alphRho(i,j)*recip_rhoConst
114 ENDIF
115 ENDDO
116 ENDDO
117
118 C -- end if integr_GeoPot = ...
119 ENDIF
120
121 C -- end usingZCoords
122 ENDIF
123
124 IF ( k.EQ.Nr ) THEN
125 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
126 C -- last level (bottom): rescale (r*) and add surface contribution
127
128 IF ( usingPCoords ) THEN
129 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 ENDDO
134 ENDDO
135 ENDIF
136
137 #ifdef NONLIN_FRSURF
138 c IF ( select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
139 IF ( select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
140 C- Integral of b.dr = rStarFac * Integral of b.dr* :
141 IF ( fluidIsAir ) THEN
142 C- Consistent with Phi'= Integr[ theta'.dPi ] :
143 DO j=jMin,jMax
144 DO i=iMin,iMax
145 facP = rStarFacC(i,j,bi,bj)**atm_kappa
146 dPhiRef = phiRef(2*k+1) - gravity*topoZ(i,j,bi,bj)
147 & - phi0surf(i,j,bi,bj)
148 phiHydLow(i,j,bi,bj) =
149 & phiHydLow(i,j,bi,bj)*facP
150 & + MAX( dPhiRef, zero )*( facP - 1. _d 0 )
151 & + phi0surf(i,j,bi,bj)
152 ENDDO
153 ENDDO
154 ELSEIF ( usingPCoords ) THEN
155 C- Note: dPhiRef*(rStarFacC -1) = 1/rho*PSo*((eta+PSo)/PSo -1) = Bo_surf*etaN
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 ENDDO
166 ENDDO
167 ELSE
168 C- Note: dPhiRef*(rStarFacC -1) = g*H*((eta+H)/H -1) = Bo_surf*etaN
169 C so that this expression is the same as before (same as in "ELSE" block below)
170 DO j=jMin,jMax
171 DO i=iMin,iMax
172 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 ENDDO
179 ENDDO
180 ENDIF
181 ELSE
182 #else /* NONLIN_FRSURF */
183 IF ( .TRUE. ) THEN
184 #endif /* NONLIN_FRSURF */
185
186 DO j=jMin,jMax
187 DO i=iMin,iMax
188 phiHydLow(i,j,bi,bj) = phiHydLow(i,j,bi,bj)
189 & + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
190 & + phi0surf(i,j,bi,bj)
191 ENDDO
192 ENDDO
193 ENDIF
194
195 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
196 C -- end if k=Nr.
197 ENDIF
198
199 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
200
201 RETURN
202 END
  ViewVC Help
Powered by ViewVC 1.1.22