/[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.45 - (hide annotations) (download)
Mon May 12 01:27:47 2014 UTC (11 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65a, checkpoint65, checkpoint64y, checkpoint64x, checkpoint64z
Changes since 1.44: +11 -3 lines 
- With p* or Sigma-P, use constant reference Pot.Temp (thetaConst) instead
  of vertical profile tRef in geopotential background and anomaly.
1 jmc 1.45 C $Header: /u/gcmpack/MITgcm/model/src/calc_phi_hyd.F,v 1.44 2014/04/04 20:54:11 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.30 C-- virtual potential temperature anomaly (including water vapour effect)
432 jmc 1.45 IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
433     C- isothermal (theta=const) reference state
434     thetaRef = thetaConst
435     ELSE
436     C- horizontally uniform (tRef) reference state
437     thetaRef = tRef(k)
438     ENDIF
439 jmc 1.30 DO j=jMin,jMax
440     DO i=iMin,iMax
441 jmc 1.42 alphaRho(i,j) = ( tFld(i,j,k,bi,bj)
442     & *( sFld(i,j,k,bi,bj)*atm_Rq + oneRL )
443 jmc 1.45 & - thetaRef )*maskC(i,j,k,bi,bj)
444 jmc 1.30 ENDDO
445     ENDDO
446    
447 jmc 1.42 #ifdef ALLOW_MOM_COMMON
448     C-- Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
449     IF (quasiHydrostatic) THEN
450     CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
451     ENDIF
452     #endif /* ALLOW_MOM_COMMON */
453    
454 jmc 1.29 C--- Integrate d Phi / d pi
455 adcroft 1.9
456 jmc 1.43 #ifndef DISABLE_SIGMA_CODE
457     C -- Finite Volume Form, integrated to r-level (cell center & upper interface)
458     IF ( useFVgradPhi ) THEN
459    
460     fullDepth = rF(1)-rF(Nr+1)
461     DO j=jMin,jMax
462     DO i=iMin,iMax
463     locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
464     #ifdef NONLIN_FRSURF
465     locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
466     #endif
467     pKappaF(i,j) = (
468     & ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
469     & + bHybSigmF( k )*locDepth
470     & )/atm_Po )**atm_kappa
471     pKappaC = (
472     & ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
473     & + bHybSigmC( k )*locDepth
474     & )/atm_Po )**atm_kappa
475     pKappaU(i,j) = (
476     & ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
477     & + bHybSigmF(k+1)*locDepth
478     & )/atm_Po )**atm_kappa
479     phiHydC(i,j) = phiHydF(i,j)
480     & + atm_Cp*( pKappaF(i,j) - pKappaC )*alphaRho(i,j)
481     phiHydU(i,j) = phiHydF(i,j)
482     & + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
483     ENDDO
484     ENDDO
485     C end: Finite Volume Form, integrated to r-level.
486    
487     ELSEIF (integr_GeoPot.EQ.0) THEN
488     #else /* DISABLE_SIGMA_CODE */
489 jmc 1.29 IF (integr_GeoPot.EQ.0) THEN
490 jmc 1.43 #endif /* DISABLE_SIGMA_CODE */
491 jmc 1.29 C -- Energy Conserving Form, accurate with Full cell topo --
492 jmc 1.14 C------------ The integration for the first level phi(k=1) is the same
493     C for both the "finite volume" and energy conserving methods.
494 jmc 1.36 C *NOTE* o Working with geopotential Anomaly, the geopotential boundary
495 adcroft 1.17 C condition is simply Phi-prime(Ro_surf)=0.
496 jmc 1.14 C o convention ddPI > 0 (same as drF & drC)
497     C-----------------------------------------------------------------------
498 jmc 1.29 IF (k.EQ.1) THEN
499     ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
500     & -((rC( k )/atm_Po)**atm_kappa) )
501     ELSE
502     ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
503 jmc 1.42 & -((rC( k )/atm_Po)**atm_kappa) )*halfRL
504 jmc 1.29 ENDIF
505     IF (k.EQ.Nr) THEN
506     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
507     & -((rF(k+1)/atm_Po)**atm_kappa) )
508     ELSE
509     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
510 jmc 1.42 & -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
511 jmc 1.29 ENDIF
512 jmc 1.14 C-------- This discretization is the energy conserving form
513 jmc 1.29 DO j=jMin,jMax
514     DO i=iMin,iMax
515 jmc 1.30 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
516     phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
517 jmc 1.14 ENDDO
518 jmc 1.29 ENDDO
519 jmc 1.14 C end: Energy Conserving Form, No hFac --
520 adcroft 1.9 C-----------------------------------------------------------------------
521 jmc 1.14
522 jmc 1.29 ELSEIF (integr_GeoPot.EQ.1) THEN
523     C -- Finite Volume Form, with Part-Cell Topo, linear in P by Half level
524 jmc 1.14 C---------
525     C Finite Volume formulation consistent with Partial Cell, linear in p by piece
526     C Note: a true Finite Volume form should be linear between 2 Interf_W :
527     C phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
528     C also: if Interface_W at the middle between tracer levels, this form
529 jmc 1.36 C is close to the Energy Cons. form in the Interior, except for the
530 jmc 1.14 C non-linearity in PI(p)
531     C---------
532 jmc 1.29 ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
533     & -((rC( k )/atm_Po)**atm_kappa) )
534     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
535     & -((rF(k+1)/atm_Po)**atm_kappa) )
536     DO j=jMin,jMax
537     DO i=iMin,iMax
538 jmc 1.37 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
539 jmc 1.29 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
540     #ifdef NONLIN_FRSURF
541     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
542     #endif
543     phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
544 jmc 1.30 & *ddPIm*alphaRho(i,j)
545 jmc 1.29 ELSE
546 jmc 1.30 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
547 jmc 1.29 ENDIF
548 jmc 1.30 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
549 adcroft 1.9 ENDDO
550 jmc 1.29 ENDDO
551     C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
552 adcroft 1.9 C-----------------------------------------------------------------------
553    
554 jmc 1.29 ELSEIF ( integr_GeoPot.EQ.2
555     & .OR. integr_GeoPot.EQ.3 ) THEN
556 jmc 1.36 C -- Finite Difference Form, with Part-Cell Topo,
557 jmc 1.29 C works with Interface_W at the middle between 2.Tracer_Level
558     C and with Tracer_Level at the middle between 2.Interface_W.
559 jmc 1.14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
560     C Finite Difference formulation consistent with Partial Cell,
561     C Valid & accurate if Interface_W at middle between tracer levels
562 jmc 1.36 C linear in p between 2 Tracer levels ; conserve energy in the Interior
563 jmc 1.14 C---------
564 jmc 1.29 IF (k.EQ.1) THEN
565     ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
566     & -((rC( k )/atm_Po)**atm_kappa) )
567     ELSE
568     ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
569 jmc 1.42 & -((rC( k )/atm_Po)**atm_kappa) )*halfRL
570 jmc 1.29 ENDIF
571     IF (k.EQ.Nr) THEN
572     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
573     & -((rF(k+1)/atm_Po)**atm_kappa) )
574     ELSE
575     ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
576 jmc 1.42 & -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
577 jmc 1.29 ENDIF
578 jmc 1.42 rec_dRm = oneRL/(rF(k)-rC(k))
579     rec_dRp = oneRL/(rC(k)-rF(k+1))
580 jmc 1.29 DO j=jMin,jMax
581     DO i=iMin,iMax
582 jmc 1.37 IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
583 jmc 1.29 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
584     #ifdef NONLIN_FRSURF
585     ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
586     #endif
587 jmc 1.42 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
588     & +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
589     & )*alphaRho(i,j)
590 jmc 1.29 ELSE
591 jmc 1.30 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
592 jmc 1.29 ENDIF
593 jmc 1.30 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
594 jmc 1.14 ENDDO
595 jmc 1.29 ENDDO
596     C end: Finite Difference Form, with Part-Cell Topo
597 jmc 1.14 C-----------------------------------------------------------------------
598 cnh 1.1
599 jmc 1.29 ELSE
600     STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
601     ENDIF
602 cnh 1.6
603 jmc 1.14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
604 adcroft 1.9 ELSE
605 jmc 1.24 STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
606 jmc 1.25 ENDIF
607    
608 jmc 1.43 IF ( .NOT. useFVgradPhi ) THEN
609     C-- r-coordinate and r*-coordinate cases:
610    
611     IF ( momPressureForcing ) THEN
612     CALL CALC_GRAD_PHI_HYD(
613     I k, bi, bj, iMin,iMax, jMin,jMax,
614     I phiHydC, alphaRho, tFld, sFld,
615     O dPhiHydX, dPhiHydY,
616     I myTime, myIter, myThid)
617     ENDIF
618    
619     #ifndef DISABLE_SIGMA_CODE
620     ELSE
621     C-- else (SigmaCoords part)
622    
623     IF ( fluidIsWater ) THEN
624     STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
625     ENDIF
626     IF ( momPressureForcing ) THEN
627     CALL CALC_GRAD_PHI_FV(
628     I k, bi, bj, iMin,iMax, jMin,jMax,
629     I phiHydF, phiHydU, pKappaF, pKappaU,
630     O dPhiHydX, dPhiHydY,
631     I myTime, myIter, myThid)
632     ENDIF
633     DO j=jMin,jMax
634     DO i=iMin,iMax
635     phiHydF(i,j) = phiHydU(i,j)
636     ENDDO
637     ENDDO
638    
639     #endif /* DISABLE_SIGMA_CODE */
640     C-- end if-not/else useFVgradPhi
641     ENDIF
642    
643 jmc 1.29 C--- Diagnose Phi at boundary r=R_low :
644     C = Ocean bottom pressure (Ocean, Z-coord.)
645     C = Sea-surface height (Ocean, P-coord.)
646     C = Top atmosphere height (Atmos, P-coord.)
647     IF (useDiagPhiRlow) THEN
648     CALL DIAGS_PHI_RLOW(
649     I k, bi, bj, iMin,iMax, jMin,jMax,
650     I phiHydF, phiHydC, alphaRho, tFld, sFld,
651 jmc 1.36 I myTime, myIter, myThid)
652 jmc 1.29 ENDIF
653    
654     C--- Diagnose Full Hydrostatic Potential at cell center level
655     CALL DIAGS_PHI_HYD(
656     I k, bi, bj, iMin,iMax, jMin,jMax,
657     I phiHydC,
658     I myTime, myIter, myThid)
659    
660 jmc 1.14 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
661 cnh 1.6
662 jmc 1.11 RETURN
663     END
  ViewVC Help
Powered by ViewVC 1.1.22