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

Annotation of /MITgcm/model/src/calc_phi_hyd.F

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

Revision 1.43 - (hide annotations) (download)
Fri Jan 4 21:07:07 2013 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64c, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.42: +83 -10 lines 
implement Finite-Volume method for (hydrostatic) presure gradient
from S.-J. Lin (QJRMS 1997), for atmosphere using sigma-coordinate.
1 jmc 1.43 C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.42 2012/12/18 01:16:53 jmc Exp $
2 cnh 1.16 C $Name: $
3 cnh 1.1
4 jmc 1.32 #include "PACKAGES_CONFIG.h"
5 cnh 1.6 #include "CPP_OPTIONS.h"
6 cnh 1.1
7 cnh 1.16 CBOP
8     C !ROUTINE: CALC_PHI_HYD
9     C !INTERFACE:
10 adcroft 1.9 SUBROUTINE CALC_PHI_HYD(
11 jmc 1.29 I bi, bj, iMin, iMax, jMin, jMax, k,
12 mlosch 1.20 I tFld, sFld,
13 jmc 1.29 U phiHydF,
14     O phiHydC, dPhiHydX, dPhiHydY,
15 jmc 1.38 I myTime, myIter, myThid )
16 cnh 1.16 C !DESCRIPTION: \bv
17     C *==========================================================*
18 cnh 1.1 C | SUBROUTINE CALC_PHI_HYD |
19 jmc 1.36 C | o Integrate the hydrostatic relation to find the Hydros. |
20 cnh 1.16 C *==========================================================*
21 jmc 1.29 C | Potential (ocean: Pressure/rho ; atmos = geopotential)
22     C | On entry:
23     C | tFld,sFld are the current thermodynamics quantities
24     C | (unchanged on exit)
25     C | phiHydF(i,j) is the hydrostatic Potential anomaly
26 jmc 1.36 C | at middle between tracer points k-1,k
27 jmc 1.29 C | On exit:
28     C | phiHydC(i,j) is the hydrostatic Potential anomaly
29     C | at cell centers (tracer points), level k
30     C | phiHydF(i,j) is the hydrostatic Potential anomaly
31 jmc 1.36 C | at middle between tracer points k,k+1
32 jmc 1.29 C | dPhiHydX,Y hydrostatic Potential gradient (X&Y dir)
33     C | at cell centers (tracer points), level k
34     C | integr_GeoPot allows to select one integration method
35     C | 1= Finite volume form ; else= Finite difference form
36 cnh 1.16 C *==========================================================*
37     C \ev
38     C !USES:
39 cnh 1.1 IMPLICIT NONE
40     C == Global variables ==
41     #include "SIZE.h"
42     #include "GRID.h"
43     #include "EEPARAMS.h"
44     #include "PARAMS.h"
45 heimbach 1.13 #ifdef ALLOW_AUTODIFF_TAMC
46     #include "tamc.h"
47     #include "tamc_keys.h"
48     #endif /* ALLOW_AUTODIFF_TAMC */
49 adcroft 1.19 #include "SURFACE.h"
50 mlosch 1.20 #include "DYNVARS.h"
51 heimbach 1.13
52 cnh 1.16 C !INPUT/OUTPUT PARAMETERS:
53 cnh 1.1 C == Routine arguments ==
54 jmc 1.36 C bi, bj, k :: tile and level indices
55 jmc 1.29 C iMin,iMax,jMin,jMax :: computational domain
56     C tFld :: potential temperature
57     C sFld :: salinity
58     C phiHydF :: hydrostatic potential anomaly at middle between
59     C 2 centers (entry: Interf_k ; output: Interf_k+1)
60     C phiHydC :: hydrostatic potential anomaly at cell center
61     C dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
62     C myTime :: current time
63     C myIter :: current iteration number
64     C myThid :: thread number for this instance of the routine.
65     INTEGER bi,bj,iMin,iMax,jMin,jMax,k
66 mlosch 1.20 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
67     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
68 jmc 1.29 _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69     _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 jmc 1.41 _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 jmc 1.25 _RL myTime
73     INTEGER myIter, myThid
74 jmc 1.36
75 adcroft 1.9 #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
76    
77 cnh 1.16 C !LOCAL VARIABLES:
78 cnh 1.1 C == Local variables ==
79 jmc 1.43 C phiHydU :: hydrostatic potential anomaly at interface k+1 (Upper k)
80     C pKappaF :: (p/Po)^kappa at interface k
81     C pKappaU :: (p/Po)^kappa at interface k+1 (Upper k)
82     C pKappaC :: (p/Po)^kappa at cell center k
83 jmc 1.29 INTEGER i,j
84 adcroft 1.9 _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85 jmc 1.43 #ifndef DISABLE_SIGMA_CODE
86     _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87     _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88     _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89     _RL pKappaC, locDepth, fullDepth
90     #endif /* DISABLE_SIGMA_CODE */
91 jmc 1.29 _RL dRlocM,dRlocP, ddRloc, locAlpha
92     _RL ddPIm, ddPIp, rec_dRm, rec_dRp
93     _RL surfPhiFac
94     LOGICAL useDiagPhiRlow, addSurfPhiAnom
95 jmc 1.43 LOGICAL useFVgradPhi
96 cnh 1.16 CEOP
97 jmc 1.27 useDiagPhiRlow = .TRUE.
98 jmc 1.41 addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
99 jmc 1.43 useFVgradPhi = selectSigmaCoord.NE.0
100    
101 jmc 1.29 surfPhiFac = 0.
102     IF (addSurfPhiAnom) surfPhiFac = 1.
103 jmc 1.14
104     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105 jmc 1.36 C Atmosphere:
106 jmc 1.24 C integr_GeoPot => select one option for the integration of the Geopotential:
107 jmc 1.29 C = 0 : Energy Conserving Form, accurate with Topo full cell;
108     C = 1 : Finite Volume Form, with Part-Cell, linear in P by Half level;
109 jmc 1.36 C =2,3: Finite Difference Form, with Part-Cell,
110 jmc 1.29 C linear in P between 2 Tracer levels.
111 jmc 1.36 C can handle both cases: Tracer lev at the middle of InterFace_W
112 jmc 1.29 C and InterFace_W at the middle of Tracer lev;
113 jmc 1.14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
114 adcroft 1.9
115 heimbach 1.13 #ifdef ALLOW_AUTODIFF_TAMC
116     act1 = bi - myBxLo(myThid)
117     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
118    
119     act2 = bj - myByLo(myThid)
120     max2 = myByHi(myThid) - myByLo(myThid) + 1
121    
122     act3 = myThid - 1
123     max3 = nTx*nTy
124    
125     act4 = ikey_dynamics - 1
126    
127     ikey = (act1 + 1) + act2*max1
128     & + act3*max1*max2
129     & + act4*max1*max2*max3
130     #endif /* ALLOW_AUTODIFF_TAMC */
131    
132 jmc 1.36 C-- Initialize phiHydF to zero :
133 jmc 1.29 C note: atmospheric_loading or Phi_topo anomaly are incorporated
134     C later in S/R calc_grad_phi_hyd
135     IF (k.EQ.1) THEN
136 jmc 1.41 DO j=1-OLy,sNy+OLy
137     DO i=1-OLx,sNx+OLx
138 jmc 1.29 phiHydF(i,j) = 0.
139     ENDDO
140     ENDDO
141     ENDIF
142 jmc 1.25
143     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
144 jmc 1.29 IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
145 adcroft 1.9 C This is the hydrostatic pressure calculation for the Ocean
146     C which uses the FIND_RHO() routine to calculate density
147     C before integrating g*rho over the current layer/interface
148 jmc 1.38 #ifdef ALLOW_AUTODIFF_TAMC
149 jmc 1.25 CADJ GENERAL
150 jmc 1.38 #endif /* ALLOW_AUTODIFF_TAMC */
151 adcroft 1.9
152 jmc 1.38 IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
153 jmc 1.29 C--- Calculate density
154 heimbach 1.13 #ifdef ALLOW_AUTODIFF_TAMC
155 jmc 1.38 kkey = (ikey-1)*Nr + k
156 heimbach 1.39 CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
157     CADJ & kind = isbyte
158     CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
159     CADJ & kind = isbyte
160 heimbach 1.13 #endif /* ALLOW_AUTODIFF_TAMC */
161 jmc 1.38 CALL FIND_RHO_2D(
162     I iMin, iMax, jMin, jMax, k,
163     I tFld(1-OLx,1-OLy,k,bi,bj),
164     I sFld(1-OLx,1-OLy,k,bi,bj),
165     O alphaRho,
166     I k, bi, bj, myThid )
167     ELSE
168     DO j=jMin,jMax
169     DO i=iMin,iMax
170     alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
171     ENDDO
172     ENDDO
173     ENDIF
174 jmc 1.36
175 mlosch 1.33 #ifdef ALLOW_SHELFICE
176 jmc 1.36 C mask rho, so that there is no contribution of phiHyd from
177 mlosch 1.33 C overlying shelfice (whose density we do not know)
178 jmc 1.37 IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
179     C- note: does not work for down_slope pkg which needs rho below the bottom.
180     C setting rho=0 above the ice-shelf base is enough (and works in both cases)
181     C but might be slower (--> keep original masking if not using down_slope pkg)
182     DO j=jMin,jMax
183     DO i=iMin,iMax
184     IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
185     ENDDO
186     ENDDO
187     ELSEIF ( useShelfIce ) THEN
188 mlosch 1.33 DO j=jMin,jMax
189     DO i=iMin,iMax
190     alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
191     ENDDO
192     ENDDO
193     ENDIF
194     #endif /* ALLOW_SHELFICE */
195 adcroft 1.22
196 jmc 1.34 #ifdef ALLOW_MOM_COMMON
197 jmc 1.42 C-- Quasi-hydrostatic terms are added in as if they modify the buoyancy
198 adcroft 1.22 IF (quasiHydrostatic) THEN
199 jmc 1.34 CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
200 adcroft 1.22 ENDIF
201 jmc 1.34 #endif /* ALLOW_MOM_COMMON */
202 adcroft 1.9
203 jmc 1.29 #ifdef NONLIN_FRSURF
204 jmc 1.41 IF ( addSurfPhiAnom .AND.
205     & uniformFreeSurfLev .AND. k.EQ.1 ) THEN
206 jmc 1.29 DO j=jMin,jMax
207     DO i=iMin,iMax
208     phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
209     & *gravity*alphaRho(i,j)*recip_rhoConst
210     ENDDO
211     ENDDO
212     ENDIF
213     #endif /* NONLIN_FRSURF */
214 jmc 1.27
215 jmc 1.29 C---- Hydrostatic pressure at cell centers
216 jmc 1.25
217     IF (integr_GeoPot.EQ.1) THEN
218     C -- Finite Volume Form
219    
220     C---------- This discretization is the "finite volume" form
221     C which has not been used to date since it does not
222     C conserve KE+PE exactly even though it is more natural
223 jmc 1.41
224     IF ( uniformFreeSurfLev ) THEN
225     DO j=jMin,jMax
226     DO i=iMin,iMax
227     phiHydC(i,j) = phiHydF(i,j)
228 jmc 1.42 & + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
229 jmc 1.41 phiHydF(i,j) = phiHydF(i,j)
230     & + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
231     ENDDO
232     ENDDO
233     ELSE
234     DO j=jMin,jMax
235     DO i=iMin,iMax
236     IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
237     ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
238     #ifdef NONLIN_FRSURF
239     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
240     #endif
241     phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst
242     ELSE
243     phiHydC(i,j) = phiHydF(i,j)
244 jmc 1.42 & + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
245 jmc 1.41 ENDIF
246     phiHydF(i,j) = phiHydC(i,j)
247 jmc 1.42 & + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
248 jmc 1.25 ENDDO
249     ENDDO
250 jmc 1.41 ENDIF
251 jmc 1.25
252     ELSE
253     C -- Finite Difference Form
254    
255 jmc 1.41 C---------- This discretization is the "energy conserving" form
256     C which has been used since at least Adcroft et al., MWR 1997
257 jmc 1.29
258 jmc 1.42 dRlocM = halfRL*drC(k)
259 jmc 1.41 IF (k.EQ.1) dRlocM=rF(k)-rC(k)
260     IF (k.EQ.Nr) THEN
261     dRlocP=rC(k)-rF(k+1)
262     ELSE
263 jmc 1.42 dRlocP=halfRL*drC(k+1)
264 jmc 1.41 ENDIF
265     IF ( uniformFreeSurfLev ) THEN
266     DO j=jMin,jMax
267     DO i=iMin,iMax
268     phiHydC(i,j) = phiHydF(i,j)
269     & +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
270     phiHydF(i,j) = phiHydC(i,j)
271     & +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
272     ENDDO
273     ENDDO
274     ELSE
275 jmc 1.42 rec_dRm = oneRL/(rF(k)-rC(k))
276     rec_dRp = oneRL/(rC(k)-rF(k+1))
277 jmc 1.25 DO j=jMin,jMax
278     DO i=iMin,iMax
279 jmc 1.41 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
280     ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
281     #ifdef NONLIN_FRSURF
282     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
283     #endif
284 jmc 1.42 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
285     & +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
286     & )*gravity*alphaRho(i,j)*recip_rhoConst
287 jmc 1.41 ELSE
288     phiHydC(i,j) = phiHydF(i,j)
289 jmc 1.29 & +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
290 jmc 1.41 ENDIF
291     phiHydF(i,j) = phiHydC(i,j)
292 jmc 1.29 & +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
293 adcroft 1.9 ENDDO
294 jmc 1.25 ENDDO
295 jmc 1.41 ENDIF
296 jmc 1.25
297     C -- end if integr_GeoPot = ...
298     ENDIF
299 jmc 1.36
300 jmc 1.25 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
301 jmc 1.29 ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
302 adcroft 1.19 C This is the hydrostatic pressure calculation for the Ocean
303 jmc 1.40 C which uses the FIND_RHO() routine to calculate density before
304     C integrating (1/rho)_prime*dp over the current layer/interface
305 mlosch 1.21 #ifdef ALLOW_AUTODIFF_TAMC
306     CADJ GENERAL
307     #endif /* ALLOW_AUTODIFF_TAMC */
308 adcroft 1.19
309 jmc 1.38 IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
310 jmc 1.27 C-- Calculate density
311 adcroft 1.19 #ifdef ALLOW_AUTODIFF_TAMC
312 jmc 1.38 kkey = (ikey-1)*Nr + k
313 heimbach 1.39 CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
314     CADJ & kind = isbyte
315     CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
316     CADJ & kind = isbyte
317 adcroft 1.19 #endif /* ALLOW_AUTODIFF_TAMC */
318 jmc 1.38 CALL FIND_RHO_2D(
319     I iMin, iMax, jMin, jMax, k,
320     I tFld(1-OLx,1-OLy,k,bi,bj),
321     I sFld(1-OLx,1-OLy,k,bi,bj),
322     O alphaRho,
323     I k, bi, bj, myThid )
324 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
325 heimbach 1.39 CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte,
326     CADJ & kind = isbyte
327 heimbach 1.23 #endif /* ALLOW_AUTODIFF_TAMC */
328 jmc 1.38 ELSE
329     DO j=jMin,jMax
330     DO i=iMin,iMax
331     alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
332     ENDDO
333     ENDDO
334     ENDIF
335 jmc 1.31
336 jmc 1.40 C-- Calculate specific volume anomaly : alpha_prime = 1/rho - alpha_Cst
337 jmc 1.27 DO j=jMin,jMax
338     DO i=iMin,iMax
339     locAlpha=alphaRho(i,j)+rhoConst
340     alphaRho(i,j)=maskC(i,j,k,bi,bj)*
341 jmc 1.42 & (oneRL/locAlpha - recip_rhoConst)
342 jmc 1.27 ENDDO
343     ENDDO
344    
345 jmc 1.42 #ifdef ALLOW_MOM_COMMON
346     C-- Quasi-hydrostatic terms are added as if they modify the specific-volume
347     IF (quasiHydrostatic) THEN
348     CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
349     ENDIF
350     #endif /* ALLOW_MOM_COMMON */
351    
352 jmc 1.25 C---- Hydrostatic pressure at cell centers
353    
354     IF (integr_GeoPot.EQ.1) THEN
355     C -- Finite Volume Form
356    
357     DO j=jMin,jMax
358 adcroft 1.19 DO i=iMin,iMax
359 jmc 1.25
360     C---------- This discretization is the "finite volume" form
361     C which has not been used to date since it does not
362     C conserve KE+PE exactly even though it is more natural
363 jmc 1.41
364 jmc 1.37 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
365 jmc 1.29 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
366     #ifdef NONLIN_FRSURF
367     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
368     #endif
369     phiHydC(i,j) = ddRloc*alphaRho(i,j)
370 jmc 1.36 c--to reproduce results of c48d_post: uncomment those 4+1 lines
371 jmc 1.29 c phiHydC(i,j)=phiHydF(i,j)
372 jmc 1.42 c & +(hFacC(i,j,k,bi,bj)-halfRL)*drF(k)*alphaRho(i,j)
373 jmc 1.29 c phiHydF(i,j)=phiHydF(i,j)
374     c & + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
375     ELSE
376 jmc 1.42 phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
377     c phiHydF(i,j) = phiHydF(i,j) + drF(k)*alphaRho(i,j)
378 jmc 1.29 ENDIF
379     c-- and comment this last one:
380 jmc 1.42 phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
381 jmc 1.29 c-----
382 jmc 1.25 ENDDO
383     ENDDO
384    
385     ELSE
386 jmc 1.29 C -- Finite Difference Form, with Part-Cell Bathy
387    
388 jmc 1.42 dRlocM = halfRL*drC(k)
389 jmc 1.29 IF (k.EQ.1) dRlocM=rF(k)-rC(k)
390     IF (k.EQ.Nr) THEN
391     dRlocP=rC(k)-rF(k+1)
392     ELSE
393 jmc 1.42 dRlocP=halfRL*drC(k+1)
394 jmc 1.29 ENDIF
395 jmc 1.42 rec_dRm = oneRL/(rF(k)-rC(k))
396     rec_dRp = oneRL/(rC(k)-rF(k+1))
397 jmc 1.25
398     DO j=jMin,jMax
399     DO i=iMin,iMax
400 adcroft 1.9
401 adcroft 1.19 C---------- This discretization is the "energy conserving" form
402 mlosch 1.21
403 jmc 1.37 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
404 jmc 1.29 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
405     #ifdef NONLIN_FRSURF
406     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
407     #endif
408 jmc 1.42 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
409     & +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
410 jmc 1.29 & )*alphaRho(i,j)
411     ELSE
412     phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
413     ENDIF
414     phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
415 adcroft 1.19 ENDDO
416 jmc 1.25 ENDDO
417    
418     C -- end if integr_GeoPot = ...
419     ENDIF
420 adcroft 1.9
421 jmc 1.29 ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
422 jmc 1.14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
423 adcroft 1.9 C This is the hydrostatic geopotential calculation for the Atmosphere
424     C The ideal gas law is used implicitly here rather than calculating
425     C the specific volume, analogous to the oceanic case.
426    
427 jmc 1.30 C-- virtual potential temperature anomaly (including water vapour effect)
428     DO j=jMin,jMax
429     DO i=iMin,iMax
430 jmc 1.42 alphaRho(i,j) = ( tFld(i,j,k,bi,bj)
431     & *( sFld(i,j,k,bi,bj)*atm_Rq + oneRL )
432     & - tRef(k) )*maskC(i,j,k,bi,bj)
433 jmc 1.30 ENDDO
434     ENDDO
435    
436 jmc 1.42 #ifdef ALLOW_MOM_COMMON
437     C-- Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
438     IF (quasiHydrostatic) THEN
439     CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
440     ENDIF
441     #endif /* ALLOW_MOM_COMMON */
442    
443 jmc 1.29 C--- Integrate d Phi / d pi
444 adcroft 1.9
445 jmc 1.43 #ifndef DISABLE_SIGMA_CODE
446     C -- Finite Volume Form, integrated to r-level (cell center & upper interface)
447     IF ( useFVgradPhi ) THEN
448    
449     fullDepth = rF(1)-rF(Nr+1)
450     DO j=jMin,jMax
451     DO i=iMin,iMax
452     locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
453     #ifdef NONLIN_FRSURF
454     locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
455     #endif
456     pKappaF(i,j) = (
457     & ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
458     & + bHybSigmF( k )*locDepth
459     & )/atm_Po )**atm_kappa
460     pKappaC = (
461     & ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
462     & + bHybSigmC( k )*locDepth
463     & )/atm_Po )**atm_kappa
464     pKappaU(i,j) = (
465     & ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
466     & + bHybSigmF(k+1)*locDepth
467     & )/atm_Po )**atm_kappa
468     phiHydC(i,j) = phiHydF(i,j)
469     & + atm_Cp*( pKappaF(i,j) - pKappaC )*alphaRho(i,j)
470     phiHydU(i,j) = phiHydF(i,j)
471     & + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
472     ENDDO
473     ENDDO
474     C end: Finite Volume Form, integrated to r-level.
475    
476     ELSEIF (integr_GeoPot.EQ.0) THEN
477     #else /* DISABLE_SIGMA_CODE */
478 jmc 1.29 IF (integr_GeoPot.EQ.0) THEN
479 jmc 1.43 #endif /* DISABLE_SIGMA_CODE */
480 jmc 1.29 C -- Energy Conserving Form, accurate with Full cell topo --
481 jmc 1.14 C------------ The integration for the first level phi(k=1) is the same
482     C for both the "finite volume" and energy conserving methods.
483 jmc 1.36 C *NOTE* o Working with geopotential Anomaly, the geopotential boundary
484 adcroft 1.17 C condition is simply Phi-prime(Ro_surf)=0.
485 jmc 1.14 C o convention ddPI > 0 (same as drF & drC)
486     C-----------------------------------------------------------------------
487 jmc 1.29 IF (k.EQ.1) THEN
488     ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
489     & -((rC( k )/atm_Po)**atm_kappa) )
490     ELSE
491     ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
492 jmc 1.42 & -((rC( k )/atm_Po)**atm_kappa) )*halfRL
493 jmc 1.29 ENDIF
494     IF (k.EQ.Nr) THEN
495     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
496     & -((rF(k+1)/atm_Po)**atm_kappa) )
497     ELSE
498     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
499 jmc 1.42 & -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
500 jmc 1.29 ENDIF
501 jmc 1.14 C-------- This discretization is the energy conserving form
502 jmc 1.29 DO j=jMin,jMax
503     DO i=iMin,iMax
504 jmc 1.30 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
505     phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
506 jmc 1.14 ENDDO
507 jmc 1.29 ENDDO
508 jmc 1.14 C end: Energy Conserving Form, No hFac --
509 adcroft 1.9 C-----------------------------------------------------------------------
510 jmc 1.14
511 jmc 1.29 ELSEIF (integr_GeoPot.EQ.1) THEN
512     C -- Finite Volume Form, with Part-Cell Topo, linear in P by Half level
513 jmc 1.14 C---------
514     C Finite Volume formulation consistent with Partial Cell, linear in p by piece
515     C Note: a true Finite Volume form should be linear between 2 Interf_W :
516     C phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
517     C also: if Interface_W at the middle between tracer levels, this form
518 jmc 1.36 C is close to the Energy Cons. form in the Interior, except for the
519 jmc 1.14 C non-linearity in PI(p)
520     C---------
521 jmc 1.29 ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
522     & -((rC( k )/atm_Po)**atm_kappa) )
523     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
524     & -((rF(k+1)/atm_Po)**atm_kappa) )
525     DO j=jMin,jMax
526     DO i=iMin,iMax
527 jmc 1.37 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
528 jmc 1.29 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
529     #ifdef NONLIN_FRSURF
530     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
531     #endif
532     phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
533 jmc 1.30 & *ddPIm*alphaRho(i,j)
534 jmc 1.29 ELSE
535 jmc 1.30 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
536 jmc 1.29 ENDIF
537 jmc 1.30 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
538 adcroft 1.9 ENDDO
539 jmc 1.29 ENDDO
540     C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
541 adcroft 1.9 C-----------------------------------------------------------------------
542    
543 jmc 1.29 ELSEIF ( integr_GeoPot.EQ.2
544     & .OR. integr_GeoPot.EQ.3 ) THEN
545 jmc 1.36 C -- Finite Difference Form, with Part-Cell Topo,
546 jmc 1.29 C works with Interface_W at the middle between 2.Tracer_Level
547     C and with Tracer_Level at the middle between 2.Interface_W.
548 jmc 1.14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
549     C Finite Difference formulation consistent with Partial Cell,
550     C Valid & accurate if Interface_W at middle between tracer levels
551 jmc 1.36 C linear in p between 2 Tracer levels ; conserve energy in the Interior
552 jmc 1.14 C---------
553 jmc 1.29 IF (k.EQ.1) THEN
554     ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
555     & -((rC( k )/atm_Po)**atm_kappa) )
556     ELSE
557     ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
558 jmc 1.42 & -((rC( k )/atm_Po)**atm_kappa) )*halfRL
559 jmc 1.29 ENDIF
560     IF (k.EQ.Nr) THEN
561     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
562     & -((rF(k+1)/atm_Po)**atm_kappa) )
563     ELSE
564     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
565 jmc 1.42 & -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
566 jmc 1.29 ENDIF
567 jmc 1.42 rec_dRm = oneRL/(rF(k)-rC(k))
568     rec_dRp = oneRL/(rC(k)-rF(k+1))
569 jmc 1.29 DO j=jMin,jMax
570     DO i=iMin,iMax
571 jmc 1.37 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
572 jmc 1.29 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
573     #ifdef NONLIN_FRSURF
574     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
575     #endif
576 jmc 1.42 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
577     & +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
578     & )*alphaRho(i,j)
579 jmc 1.29 ELSE
580 jmc 1.30 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
581 jmc 1.29 ENDIF
582 jmc 1.30 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
583 jmc 1.14 ENDDO
584 jmc 1.29 ENDDO
585     C end: Finite Difference Form, with Part-Cell Topo
586 jmc 1.14 C-----------------------------------------------------------------------
587 cnh 1.1
588 jmc 1.29 ELSE
589     STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
590     ENDIF
591 cnh 1.6
592 jmc 1.14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
593 adcroft 1.9 ELSE
594 jmc 1.24 STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
595 jmc 1.25 ENDIF
596    
597 jmc 1.43 IF ( .NOT. useFVgradPhi ) THEN
598     C-- r-coordinate and r*-coordinate cases:
599    
600     IF ( momPressureForcing ) THEN
601     CALL CALC_GRAD_PHI_HYD(
602     I k, bi, bj, iMin,iMax, jMin,jMax,
603     I phiHydC, alphaRho, tFld, sFld,
604     O dPhiHydX, dPhiHydY,
605     I myTime, myIter, myThid)
606     ENDIF
607    
608     #ifndef DISABLE_SIGMA_CODE
609     ELSE
610     C-- else (SigmaCoords part)
611    
612     IF ( fluidIsWater ) THEN
613     STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
614     ENDIF
615     IF ( momPressureForcing ) THEN
616     CALL CALC_GRAD_PHI_FV(
617     I k, bi, bj, iMin,iMax, jMin,jMax,
618     I phiHydF, phiHydU, pKappaF, pKappaU,
619     O dPhiHydX, dPhiHydY,
620     I myTime, myIter, myThid)
621     ENDIF
622     DO j=jMin,jMax
623     DO i=iMin,iMax
624     phiHydF(i,j) = phiHydU(i,j)
625     ENDDO
626     ENDDO
627    
628     #endif /* DISABLE_SIGMA_CODE */
629     C-- end if-not/else useFVgradPhi
630     ENDIF
631    
632 jmc 1.29 C--- Diagnose Phi at boundary r=R_low :
633     C = Ocean bottom pressure (Ocean, Z-coord.)
634     C = Sea-surface height (Ocean, P-coord.)
635     C = Top atmosphere height (Atmos, P-coord.)
636     IF (useDiagPhiRlow) THEN
637     CALL DIAGS_PHI_RLOW(
638     I k, bi, bj, iMin,iMax, jMin,jMax,
639     I phiHydF, phiHydC, alphaRho, tFld, sFld,
640 jmc 1.36 I myTime, myIter, myThid)
641 jmc 1.29 ENDIF
642    
643     C--- Diagnose Full Hydrostatic Potential at cell center level
644     CALL DIAGS_PHI_HYD(
645     I k, bi, bj, iMin,iMax, jMin,jMax,
646     I phiHydC,
647     I myTime, myIter, myThid)
648    
649 jmc 1.14 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
650 cnh 1.6
651 jmc 1.11 RETURN
652     END
  ViewVC Help
Powered by ViewVC 1.1.22