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

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

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

revision 1.1 by cnh, Thu Aug 20 18:47:27 1998 UTC revision 1.33 by mlosch, Tue Feb 7 11:47:49 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #include "CPP_OPTIONS.h"
6    
7        SUBROUTINE CALC_PHI_HYD( bi, bj, iMin, iMax, jMin, jMax, K,  CBOP
8       I    buoyKM1, buoyKP1k, phiHyd, myThid)  C     !ROUTINE: CALC_PHI_HYD
9  C     /==========================================================\  C     !INTERFACE:
10          SUBROUTINE CALC_PHI_HYD(
11         I                         bi, bj, iMin, iMax, jMin, jMax, k,
12         I                         tFld, sFld,
13         U                         phiHydF,
14         O                         phiHydC, dPhiHydX, dPhiHydY,
15         I                         myTime, myIter, myThid)
16    C     !DESCRIPTION: \bv
17    C     *==========================================================*
18  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
19  C     | o Integrate the hydrostatic relation to find phiHyd.     |  C     | o Integrate the hydrostatic relation to find the Hydros. |
20  C     |                                                          |  C     *==========================================================*
21  C     \==========================================================/  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    C     |                at middle between tracer points k-1,k
27    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    C     |                at middle between tracer points k,k+1
32    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    C     *==========================================================*
37    C     \ev
38    C     !USES:
39        IMPLICIT NONE        IMPLICIT NONE
40  C     == Global variables ==  C     == Global variables ==
41  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
42  #include "GRID.h"  #include "GRID.h"
43  #include "EEPARAMS.h"  #include "EEPARAMS.h"
44  #include "PARAMS.h"  #include "PARAMS.h"
45    #ifdef ALLOW_AUTODIFF_TAMC
46    #include "tamc.h"
47    #include "tamc_keys.h"
48    #endif /* ALLOW_AUTODIFF_TAMC */
49    #include "SURFACE.h"
50    #include "DYNVARS.h"
51    
52    C     !INPUT/OUTPUT PARAMETERS:
53  C     == Routine arguments ==  C     == Routine arguments ==
54        INTEGER bi,bj,iMin,iMax,jMin,jMax,K  C     bi, bj, k  :: tile and level indices
55        _RL buoyKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     iMin,iMax,jMin,jMax :: computational domain
56        _RL buoyKP1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     tFld       :: potential temperature
57        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)  C     sFld       :: salinity
58        integer myThid  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          _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    c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
69          _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70          _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71          _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73          _RL myTime
74          INTEGER myIter, myThid
75          
76    #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
77    
78    C     !LOCAL VARIABLES:
79  C     == Local variables ==  C     == Local variables ==
80        INTEGER i,j,Km1        INTEGER i,j
81        _RL halfLayer        _RL zero, one, half
82          _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL dRlocM,dRlocP, ddRloc, locAlpha
84          _RL ddPIm, ddPIp, rec_dRm, rec_dRp
85          _RL surfPhiFac
86          INTEGER iMnLoc,jMnLoc
87          PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
88          LOGICAL useDiagPhiRlow, addSurfPhiAnom
89    CEOP
90          useDiagPhiRlow = .TRUE.
91          addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GT.3
92          surfPhiFac = 0.
93          IF (addSurfPhiAnom) surfPhiFac = 1.
94    
95    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
96    C  Atmosphere:  
97    C integr_GeoPot => select one option for the integration of the Geopotential:
98    C   = 0 : Energy Conserving Form, accurate with Topo full cell;
99    C   = 1 : Finite Volume Form, with Part-Cell, linear in P by Half level;
100    C   =2,3: Finite Difference Form, with Part-Cell,
101    C         linear in P between 2 Tracer levels.
102    C       can handle both cases: Tracer lev at the middle of InterFace_W
103    C                          and InterFace_W at the middle of Tracer lev;
104    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105    
106    #ifdef ALLOW_AUTODIFF_TAMC
107              act1 = bi - myBxLo(myThid)
108              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
109    
110              act2 = bj - myByLo(myThid)
111              max2 = myByHi(myThid) - myByLo(myThid) + 1
112    
113              act3 = myThid - 1
114              max3 = nTx*nTy
115    
116              act4 = ikey_dynamics - 1
117    
118              ikey = (act1 + 1) + act2*max1
119         &                      + act3*max1*max2
120         &                      + act4*max1*max2*max3
121    #endif /* ALLOW_AUTODIFF_TAMC */
122    
123    C--   Initialize phiHydF to zero :
124    C     note: atmospheric_loading or Phi_topo anomaly are incorporated
125    C           later in S/R calc_grad_phi_hyd
126          IF (k.EQ.1) THEN
127            DO j=1-Oly,sNy+Oly
128             DO i=1-Olx,sNx+Olx
129               phiHydF(i,j) = 0.
130             ENDDO
131            ENDDO
132          ENDIF
133    
134    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
135          IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
136    C       This is the hydrostatic pressure calculation for the Ocean
137    C       which uses the FIND_RHO() routine to calculate density
138    C       before integrating g*rho over the current layer/interface
139    #ifdef      ALLOW_AUTODIFF_TAMC
140    CADJ GENERAL
141    #endif      /* ALLOW_AUTODIFF_TAMC */
142    
143    C---    Calculate density
144    #ifdef ALLOW_AUTODIFF_TAMC
145            kkey = (ikey-1)*Nr + k
146    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
147    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
148    #endif /* ALLOW_AUTODIFF_TAMC */
149            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
150         &                 tFld, sFld,
151         &                 alphaRho, myThid)
152    #ifdef ALLOW_SHELFICE
153    C     mask rho, so that there is no contribution of phiHyd from
154    C     overlying shelfice (whose density we do not know)
155            IF ( useShelfIce ) THEN
156             DO j=jMin,jMax
157              DO i=iMin,iMax
158               alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
159              ENDDO
160             ENDDO
161            ENDIF
162    #endif /* ALLOW_SHELFICE */
163    
164    #ifdef ALLOW_DIAGNOSTICS
165            IF ( useDiagnostics )
166         &   CALL DIAGNOSTICS_FILL(alphaRho,'RHOAnoma',k,1,2,bi,bj,myThid)
167    #endif
168    
169    C Quasi-hydrostatic terms are added in as if they modify the buoyancy
170            IF (quasiHydrostatic) THEN
171             CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
172            ENDIF
173    
174    #ifdef NONLIN_FRSURF
175            IF (k.EQ.1 .AND. addSurfPhiAnom) THEN
176              DO j=jMin,jMax
177                DO i=iMin,iMax
178                  phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
179         &                      *gravity*alphaRho(i,j)*recip_rhoConst
180                ENDDO
181              ENDDO
182            ENDIF
183    #endif /* NONLIN_FRSURF */
184    
185    C----  Hydrostatic pressure at cell centers
186    
187           IF (integr_GeoPot.EQ.1) THEN
188    C  --  Finite Volume Form
189    
190             DO j=jMin,jMax
191              DO i=iMin,iMax
192    
193    C---------- This discretization is the "finite volume" form
194    C           which has not been used to date since it does not
195    C           conserve KE+PE exactly even though it is more natural
196    C
197               phiHydC(i,j)=phiHydF(i,j)
198         &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
199               phiHydF(i,j)=phiHydF(i,j)
200         &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
201              ENDDO
202             ENDDO
203    
204           ELSE
205    C  --  Finite Difference Form
206    
207             dRlocM=half*drC(k)
208             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
209             IF (k.EQ.Nr) THEN
210               dRlocP=rC(k)-rF(k+1)
211             ELSE
212               dRlocP=half*drC(k+1)
213             ENDIF
214    
215             DO j=jMin,jMax
216              DO i=iMin,iMax
217    
218    C---------- This discretization is the "energy conserving" form
219    C           which has been used since at least Adcroft et al., MWR 1997
220    C
221                phiHydC(i,j)=phiHydF(i,j)
222         &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
223                phiHydF(i,j)=phiHydC(i,j)
224         &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
225              ENDDO
226             ENDDO
227    
228    C  --  end if integr_GeoPot = ...
229           ENDIF
230            
231    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
232          ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
233    C       This is the hydrostatic pressure calculation for the Ocean
234    C       which uses the FIND_RHO() routine to calculate density
235    C       before integrating (1/rho)'*dp over the current layer/interface
236    #ifdef      ALLOW_AUTODIFF_TAMC
237    CADJ GENERAL
238    #endif      /* ALLOW_AUTODIFF_TAMC */
239    
240    C--     Calculate density
241    #ifdef ALLOW_AUTODIFF_TAMC
242                kkey = (ikey-1)*Nr + k
243    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
244    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
245    #endif /* ALLOW_AUTODIFF_TAMC */
246            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
247         &                 tFld, sFld,
248         &                 alphaRho, myThid)
249    #ifdef ALLOW_AUTODIFF_TAMC
250    CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
251    #endif /* ALLOW_AUTODIFF_TAMC */
252    
253    #ifdef ALLOW_DIAGNOSTICS
254            IF ( useDiagnostics )
255         &   CALL DIAGNOSTICS_FILL(alphaRho,'RHOAnoma',k,1,2,bi,bj,myThid)
256    #endif
257    
258    C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
259            DO j=jMin,jMax
260              DO i=iMin,iMax
261                locAlpha=alphaRho(i,j)+rhoConst
262                alphaRho(i,j)=maskC(i,j,k,bi,bj)*
263         &              (one/locAlpha - recip_rhoConst)
264              ENDDO
265            ENDDO
266    
267    C----  Hydrostatic pressure at cell centers
268    
269           IF (integr_GeoPot.EQ.1) THEN
270    C  --  Finite Volume Form
271    
272             DO j=jMin,jMax
273              DO i=iMin,iMax
274    
275    C---------- This discretization is the "finite volume" form
276    C           which has not been used to date since it does not
277    C           conserve KE+PE exactly even though it is more natural
278    C
279               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                 phiHydC(i,j) = ddRloc*alphaRho(i,j)
285    c--to reproduce results of c48d_post: uncomment those 4+1 lines
286    c            phiHydC(i,j)=phiHydF(i,j)
287    c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)
288    c            phiHydF(i,j)=phiHydF(i,j)
289    c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
290               ELSE
291                 phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)
292    c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)
293               ENDIF
294    c-- and comment this last one:
295                 phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)
296    c-----
297              ENDDO
298             ENDDO
299    
300           ELSE
301    C  --  Finite Difference Form, with Part-Cell Bathy
302    
303             dRlocM=half*drC(k)
304             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
305             IF (k.EQ.Nr) THEN
306               dRlocP=rC(k)-rF(k+1)
307             ELSE
308               dRlocP=half*drC(k+1)
309             ENDIF
310             rec_dRm = one/(rF(k)-rC(k))
311             rec_dRp = one/(rC(k)-rF(k+1))
312    
313             DO j=jMin,jMax
314              DO i=iMin,iMax
315    
316    C---------- This discretization is the "energy conserving" form
317    
318               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
319                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
320    #ifdef NONLIN_FRSURF
321                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
322    #endif
323                 phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
324         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
325         &                     )*alphaRho(i,j)
326               ELSE
327                 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
328               ENDIF
329                 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
330              ENDDO
331             ENDDO
332    
333    C  --  end if integr_GeoPot = ...
334           ENDIF
335    
336          ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
337    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
338    C       This is the hydrostatic geopotential calculation for the Atmosphere
339    C       The ideal gas law is used implicitly here rather than calculating
340    C       the specific volume, analogous to the oceanic case.
341    
342    C--     virtual potential temperature anomaly (including water vapour effect)
343            DO j=jMin,jMax
344             DO i=iMin,iMax
345              alphaRho(i,j)=maskC(i,j,k,bi,bj)
346         &             *( tFld(i,j,k,bi,bj)*(sFld(i,j,k,bi,bj)*atm_Rq+one)
347         &               -tRef(k) )
348             ENDDO
349            ENDDO
350    
351    C---    Integrate d Phi / d pi
352    
353           IF (integr_GeoPot.EQ.0) THEN
354    C  --  Energy Conserving Form, accurate with Full cell topo  --
355    C------------ The integration for the first level phi(k=1) is the same
356    C             for both the "finite volume" and energy conserving methods.
357    C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
358    C             condition is simply Phi-prime(Ro_surf)=0.
359    C           o convention ddPI > 0 (same as drF & drC)
360    C-----------------------------------------------------------------------
361             IF (k.EQ.1) THEN
362               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
363         &                   -((rC( k )/atm_Po)**atm_kappa) )
364             ELSE
365               ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
366         &                   -((rC( k )/atm_Po)**atm_kappa) )*half
367             ENDIF
368             IF (k.EQ.Nr) THEN
369               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
370         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
371             ELSE
372               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
373         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
374             ENDIF
375    C-------- This discretization is the energy conserving form
376             DO j=jMin,jMax
377              DO i=iMin,iMax
378                 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
379                 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
380              ENDDO
381             ENDDO
382    C end: Energy Conserving Form, No hFac  --
383    C-----------------------------------------------------------------------
384    
385           ELSEIF (integr_GeoPot.EQ.1) THEN
386    C  --  Finite Volume Form, with Part-Cell Topo, linear in P by Half level
387    C---------
388    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
389    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
390    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
391    C   also: if Interface_W at the middle between tracer levels, this form
392    C     is close to the Energy Cons. form in the Interior, except for the
393    C     non-linearity in PI(p)
394    C---------
395               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
396         &                   -((rC( k )/atm_Po)**atm_kappa) )
397               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
398         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
399             DO j=jMin,jMax
400              DO i=iMin,iMax
401               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
402                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
403    #ifdef NONLIN_FRSURF
404                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
405    #endif
406                 phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
407         &          *ddPIm*alphaRho(i,j)
408               ELSE
409                 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
410               ENDIF
411                 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
412              ENDDO
413             ENDDO
414    C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
415    C-----------------------------------------------------------------------
416    
417           ELSEIF ( integr_GeoPot.EQ.2
418         &     .OR. integr_GeoPot.EQ.3 ) THEN
419    C  --  Finite Difference Form, with Part-Cell Topo,
420    C       works with Interface_W  at the middle between 2.Tracer_Level
421    C        and  with Tracer_Level at the middle between 2.Interface_W.
422    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
423    C  Finite Difference formulation consistent with Partial Cell,
424    C   Valid & accurate if Interface_W at middle between tracer levels
425    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
426    C---------
427             IF (k.EQ.1) THEN
428               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
429         &                   -((rC( k )/atm_Po)**atm_kappa) )
430             ELSE
431               ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
432         &                   -((rC( k )/atm_Po)**atm_kappa) )*half
433             ENDIF
434             IF (k.EQ.Nr) THEN
435               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
436         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
437             ELSE
438               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
439         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
440             ENDIF
441             rec_dRm = one/(rF(k)-rC(k))
442             rec_dRp = one/(rC(k)-rF(k+1))
443             DO j=jMin,jMax
444              DO i=iMin,iMax
445               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
446                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
447    #ifdef NONLIN_FRSURF
448                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
449    #endif
450                 phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm
451         &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )
452         &                    *alphaRho(i,j)
453               ELSE
454                 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
455               ENDIF
456                 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
457              ENDDO
458             ENDDO
459    C end: Finite Difference Form, with Part-Cell Topo
460    C-----------------------------------------------------------------------
461    
462           ELSE
463             STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
464           ENDIF
465    
466    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
467          ELSE
468            STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
469          ENDIF
470    
471    C---   Diagnose Phi at boundary r=R_low :
472    C       = Ocean bottom pressure (Ocean, Z-coord.)
473    C       = Sea-surface height    (Ocean, P-coord.)
474    C       = Top atmosphere height (Atmos, P-coord.)
475          IF (useDiagPhiRlow) THEN
476            CALL DIAGS_PHI_RLOW(
477         I                      k, bi, bj, iMin,iMax, jMin,jMax,
478         I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
479         I                      myTime, myIter, myThid)  
480          ENDIF
481    
482    C---   Diagnose Full Hydrostatic Potential at cell center level
483            CALL DIAGS_PHI_HYD(
484         I                      k, bi, bj, iMin,iMax, jMin,jMax,
485         I                      phiHydC,
486         I                      myTime, myIter, myThid)
487    
488          IF (momPressureForcing) THEN
489            iMnLoc = MAX(1-Olx+1,iMin)
490            jMnLoc = MAX(1-Oly+1,jMin)
491            CALL CALC_GRAD_PHI_HYD(
492         I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
493         I                         phiHydC, alphaRho, tFld, sFld,
494         O                         dPhiHydX, dPhiHydY,
495         I                         myTime, myIter, myThid)  
496          ENDIF
497    
498    #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
499    
500        if (K.eq.1) then        RETURN
501         Km1=1        END
        halfLayer=0.5 _d 0  
       else  
        Km1=K-1  
        halfLayer=1.0 _d 0  
       endif  
   
 C--   Contribution to phiHyd(:,:,K) from buoy(:,:,K-1) + buoy(:,:,K)  
 C     (This is now the actual hydrostatic pressure|height at the T/S points)  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         phiHyd(i,j,K)=phiHyd(i,j,Km1)+rhoConst*halfLayer  
      &  *0.5 _d 0*(drF(Km1)+drF(K))  
      &  *0.5 _d 0*( buoyKM1(i,j)+buoyKP1(i,j) )  
        ENDDO  
       ENDDO  
   
 ! ------------------------------------------------------------------------------  
       return  
       end  
 ! ==============================================================================  
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22