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