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