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