/[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.3 by cnh, Sat Sep 5 17:52:13 1998 UTC revision 1.29 by jmc, Tue Feb 18 15:30:47 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6        SUBROUTINE CALC_PHI_HYD( bi, bj, iMin, iMax, jMin, jMax, K,  CBOP
7       I    buoyKM1, buoyKP1, phiHyd, myThid)  C     !ROUTINE: CALC_PHI_HYD
8  C     /==========================================================\  C     !INTERFACE:
9          SUBROUTINE CALC_PHI_HYD(
10         I                         bi, bj, iMin, iMax, jMin, jMax, k,
11         I                         tFld, sFld,
12         U                         phiHydF,
13         O                         phiHydC, dPhiHydX, dPhiHydY,
14         I                         myTime, myIter, myThid)
15    C     !DESCRIPTION: \bv
16    C     *==========================================================*
17  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
18  C     | o Integrate the hydrostatic relation to find phiHyd.     |  C     | o Integrate the hydrostatic relation to find the Hydros. |
19  C     |                                                          |  C     *==========================================================*
20  C     \==========================================================/  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)
21    C     | On entry:
22    C     |   tFld,sFld     are the current thermodynamics quantities
23    C     |                 (unchanged on exit)
24    C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
25    C     |                at middle between tracer points k-1,k
26    C     | On exit:
27    C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
28    C     |                at cell centers (tracer points), level k
29    C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
30    C     |                at middle between tracer points k,k+1
31    C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
32    C     |                at cell centers (tracer points), level k
33    C     | integr_GeoPot allows to select one integration method
34    C     |    1= Finite volume form ; else= Finite difference form
35    C     *==========================================================*
36    C     \ev
37    C     !USES:
38        IMPLICIT NONE        IMPLICIT NONE
39  C     == Global variables ==  C     == Global variables ==
40  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
41  #include "GRID.h"  #include "GRID.h"
42  #include "EEPARAMS.h"  #include "EEPARAMS.h"
43  #include "PARAMS.h"  #include "PARAMS.h"
44    #ifdef ALLOW_AUTODIFF_TAMC
45    #include "tamc.h"
46    #include "tamc_keys.h"
47    #endif /* ALLOW_AUTODIFF_TAMC */
48    #include "SURFACE.h"
49    #include "DYNVARS.h"
50    
51    C     !INPUT/OUTPUT PARAMETERS:
52  C     == Routine arguments ==  C     == Routine arguments ==
53        INTEGER bi,bj,iMin,iMax,jMin,jMax,K  C     bi, bj, k  :: tile and level indices
54        _RL buoyKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     iMin,iMax,jMin,jMax :: computational domain
55        _RL buoyKP1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     tFld       :: potential temperature
56        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     sFld       :: salinity
57        integer myThid  C     phiHydF    :: hydrostatic potential anomaly at middle between
58    C                   2 centers (entry: Interf_k ; output: Interf_k+1)
59    C     phiHydC    :: hydrostatic potential anomaly at cell center
60    C     dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
61    C     myTime     :: current time
62    C     myIter     :: current iteration number
63    C     myThid     :: thread number for this instance of the routine.
64          INTEGER bi,bj,iMin,iMax,jMin,jMax,k
65          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
66          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
67    c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
68          _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69          _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70          _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72          _RL myTime
73          INTEGER myIter, myThid
74          
75    #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
76    
77    C     !LOCAL VARIABLES:
78  C     == Local variables ==  C     == Local variables ==
79        INTEGER i,j,Km1        INTEGER i,j
80        _RL halfLayer        _RL zero, one, half
81          _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL dRlocM,dRlocP, ddRloc, locAlpha
83          _RL ddPIm, ddPIp, rec_dRm, rec_dRp
84          _RL surfPhiFac
85          INTEGER iMnLoc,jMnLoc
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    
152    C Quasi-hydrostatic terms are added in as if they modify the buoyancy
153            IF (quasiHydrostatic) THEN
154             CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
155            ENDIF
156    
157    #ifdef NONLIN_FRSURF
158            IF (k.EQ.1 .AND. addSurfPhiAnom) THEN
159              DO j=jMin,jMax
160                DO i=iMin,iMax
161                  phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
162         &                      *gravity*alphaRho(i,j)*recip_rhoConst
163                ENDDO
164              ENDDO
165            ENDIF
166    #endif /* NONLIN_FRSURF */
167    
168    C----  Hydrostatic pressure at cell centers
169    
170           IF (integr_GeoPot.EQ.1) THEN
171    C  --  Finite Volume Form
172    
173             DO j=jMin,jMax
174              DO i=iMin,iMax
175    
176    C---------- This discretization is the "finite volume" form
177    C           which has not been used to date since it does not
178    C           conserve KE+PE exactly even though it is more natural
179    C
180               phiHydC(i,j)=phiHydF(i,j)
181         &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
182               phiHydF(i,j)=phiHydF(i,j)
183         &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
184              ENDDO
185             ENDDO
186    
187           ELSE
188    C  --  Finite Difference Form
189    
190             dRlocM=half*drC(k)
191             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
192             IF (k.EQ.Nr) THEN
193               dRlocP=rC(k)-rF(k+1)
194             ELSE
195               dRlocP=half*drC(k+1)
196             ENDIF
197    
198             DO j=jMin,jMax
199              DO i=iMin,iMax
200    
201    C---------- This discretization is the "energy conserving" form
202    C           which has been used since at least Adcroft et al., MWR 1997
203    C
204                phiHydC(i,j)=phiHydF(i,j)
205         &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
206                phiHydF(i,j)=phiHydC(i,j)
207         &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
208              ENDDO
209             ENDDO
210    
211    C  --  end if integr_GeoPot = ...
212           ENDIF
213            
214    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
215          ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
216    C       This is the hydrostatic pressure calculation for the Ocean
217    C       which uses the FIND_RHO() routine to calculate density
218    C       before integrating (1/rho)'*dp over the current layer/interface
219    #ifdef      ALLOW_AUTODIFF_TAMC
220    CADJ GENERAL
221    #endif      /* ALLOW_AUTODIFF_TAMC */
222    
223    C--     Calculate density
224    #ifdef ALLOW_AUTODIFF_TAMC
225                kkey = (ikey-1)*Nr + k
226    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
227    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
228    #endif /* ALLOW_AUTODIFF_TAMC */
229            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
230         &                 tFld, sFld,
231         &                 alphaRho, myThid)
232    #ifdef ALLOW_AUTODIFF_TAMC
233    CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
234    #endif /* ALLOW_AUTODIFF_TAMC */
235    
236    C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
237            DO j=jMin,jMax
238              DO i=iMin,iMax
239                locAlpha=alphaRho(i,j)+rhoConst
240                alphaRho(i,j)=maskC(i,j,k,bi,bj)*
241         &              (one/locAlpha - recip_rhoConst)
242              ENDDO
243            ENDDO
244    
245    C----  Hydrostatic pressure at cell centers
246    
247           IF (integr_GeoPot.EQ.1) THEN
248    C  --  Finite Volume Form
249    
250             DO j=jMin,jMax
251              DO i=iMin,iMax
252    
253    C---------- This discretization is the "finite volume" form
254    C           which has not been used to date since it does not
255    C           conserve KE+PE exactly even though it is more natural
256    C
257               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
258                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
259    #ifdef NONLIN_FRSURF
260                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
261    #endif
262                 phiHydC(i,j) = ddRloc*alphaRho(i,j)
263    c--to reproduce results of c48d_post: uncomment those 4+1 lines
264    c            phiHydC(i,j)=phiHydF(i,j)
265    c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)
266    c            phiHydF(i,j)=phiHydF(i,j)
267    c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
268               ELSE
269                 phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)
270    c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)
271               ENDIF
272    c-- and comment this last one:
273                 phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)
274    c-----
275              ENDDO
276             ENDDO
277    
278           ELSE
279    C  --  Finite Difference Form, with Part-Cell Bathy
280    
281             dRlocM=half*drC(k)
282             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
283             IF (k.EQ.Nr) THEN
284               dRlocP=rC(k)-rF(k+1)
285             ELSE
286               dRlocP=half*drC(k+1)
287             ENDIF
288             rec_dRm = one/(rF(k)-rC(k))
289             rec_dRp = one/(rC(k)-rF(k+1))
290    
291             DO j=jMin,jMax
292              DO i=iMin,iMax
293    
294    C---------- This discretization is the "energy conserving" form
295    
296               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
297                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
298    #ifdef NONLIN_FRSURF
299                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
300    #endif
301                 phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
302         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
303         &                     )*alphaRho(i,j)
304               ELSE
305                 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
306               ENDIF
307                 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
308              ENDDO
309             ENDDO
310    
311    C  --  end if integr_GeoPot = ...
312           ENDIF
313    
314          ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
315    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
316    C       This is the hydrostatic geopotential calculation for the Atmosphere
317    C       The ideal gas law is used implicitly here rather than calculating
318    C       the specific volume, analogous to the oceanic case.
319    
320    C---    Integrate d Phi / d pi
321    
322           IF (integr_GeoPot.EQ.0) THEN
323    C  --  Energy Conserving Form, accurate with Full cell topo  --
324    C------------ The integration for the first level phi(k=1) is the same
325    C             for both the "finite volume" and energy conserving methods.
326    C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
327    C             condition is simply Phi-prime(Ro_surf)=0.
328    C           o convention ddPI > 0 (same as drF & drC)
329    C-----------------------------------------------------------------------
330             IF (k.EQ.1) THEN
331               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
332         &                   -((rC( k )/atm_Po)**atm_kappa) )
333             ELSE
334               ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
335         &                   -((rC( k )/atm_Po)**atm_kappa) )*half
336             ENDIF
337             IF (k.EQ.Nr) THEN
338               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
339         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
340             ELSE
341               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
342         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
343             ENDIF
344    C-------- This discretization is the energy conserving form
345             DO j=jMin,jMax
346              DO i=iMin,iMax
347                 phiHydC(i,j) = phiHydF(i,j)
348         &          +ddPIm*maskC(i,j,k,bi,bj)
349         &                *(tFld(i,j,k,bi,bj)-tRef(k))
350                 phiHydF(i,j) = phiHydC(i,j)
351         &          +ddPIp*maskC(i,j,k,bi,bj)
352         &                *(tFld(i,j,k,bi,bj)-tRef(k))
353              ENDDO
354             ENDDO
355    C end: Energy Conserving Form, No hFac  --
356    C-----------------------------------------------------------------------
357    
358           ELSEIF (integr_GeoPot.EQ.1) THEN
359    C  --  Finite Volume Form, with Part-Cell Topo, linear in P by Half level
360    C---------
361    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
362    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
363    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
364    C   also: if Interface_W at the middle between tracer levels, this form
365    C     is close to the Energy Cons. form in the Interior, except for the
366    C     non-linearity in PI(p)
367    C---------
368               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
369         &                   -((rC( k )/atm_Po)**atm_kappa) )
370               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
371         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
372             DO j=jMin,jMax
373              DO i=iMin,iMax
374               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
375                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
376    #ifdef NONLIN_FRSURF
377                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
378    #endif
379                 phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
380         &          *ddPIm*maskC(i,j,k,bi,bj)
381         &                *(tFld(i,j,k,bi,bj)-tRef(k))
382               ELSE
383                 phiHydC(i,j) = phiHydF(i,j)
384         &          +ddPIm*maskC(i,j,k,bi,bj)
385         &                *(tFld(i,j,k,bi,bj)-tRef(k))
386               ENDIF
387                 phiHydF(i,j) = phiHydC(i,j)
388         &          +ddPIp*maskC(i,j,k,bi,bj)
389         &                *(tFld(i,j,k,bi,bj)-tRef(k))
390              ENDDO
391             ENDDO
392    C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
393    C-----------------------------------------------------------------------
394    
395           ELSEIF ( integr_GeoPot.EQ.2
396         &     .OR. integr_GeoPot.EQ.3 ) THEN
397    C  --  Finite Difference Form, with Part-Cell Topo,
398    C       works with Interface_W  at the middle between 2.Tracer_Level
399    C        and  with Tracer_Level at the middle between 2.Interface_W.
400    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
401    C  Finite Difference formulation consistent with Partial Cell,
402    C   Valid & accurate if Interface_W at middle between tracer levels
403    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
404    C---------
405             IF (k.EQ.1) THEN
406               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
407         &                   -((rC( k )/atm_Po)**atm_kappa) )
408             ELSE
409               ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
410         &                   -((rC( k )/atm_Po)**atm_kappa) )*half
411             ENDIF
412             IF (k.EQ.Nr) THEN
413               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
414         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
415             ELSE
416               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
417         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
418             ENDIF
419             rec_dRm = one/(rF(k)-rC(k))
420             rec_dRp = one/(rC(k)-rF(k+1))
421             DO j=jMin,jMax
422              DO i=iMin,iMax
423               IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
424                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
425    #ifdef NONLIN_FRSURF
426                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
427    #endif
428                 phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm
429         &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )
430         &                *(tFld(i,j,k,bi,bj)-tRef(k))
431               ELSE
432                 phiHydC(i,j) = phiHydF(i,j)
433         &          +ddPIm*maskC(i,j,k,bi,bj)
434         &                *(tFld(I,J,k,bi,bj)-tRef(k))
435               ENDIF
436                 phiHydF(i,j) = phiHydC(i,j)
437         &          +ddPIp*maskC(i,j,k,bi,bj)
438         &                *(tFld(I,J,k,bi,bj)-tRef(k))
439              ENDDO
440             ENDDO
441    C end: Finite Difference Form, with Part-Cell Topo
442    C-----------------------------------------------------------------------
443    
444           ELSE
445             STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
446           ENDIF
447    
448    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
449          ELSE
450            STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
451          ENDIF
452    
453    C---   Diagnose Phi at boundary r=R_low :
454    C       = Ocean bottom pressure (Ocean, Z-coord.)
455    C       = Sea-surface height    (Ocean, P-coord.)
456    C       = Top atmosphere height (Atmos, P-coord.)
457          IF (useDiagPhiRlow) THEN
458            CALL DIAGS_PHI_RLOW(
459         I                      k, bi, bj, iMin,iMax, jMin,jMax,
460         I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
461         I                      myTime, myIter, myThid)  
462          ENDIF
463    
464    C---   Diagnose Full Hydrostatic Potential at cell center level
465            CALL DIAGS_PHI_HYD(
466         I                      k, bi, bj, iMin,iMax, jMin,jMax,
467         I                      phiHydC,
468         I                      myTime, myIter, myThid)
469    
470          IF (momPressureForcing) THEN
471            iMnLoc = MAX(1-Olx+1,iMin)
472            jMnLoc = MAX(1-Oly+1,jMin)
473            CALL CALC_GRAD_PHI_HYD(
474         I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
475         I                         phiHydC, alphaRho, tFld, sFld,
476         O                         dPhiHydX, dPhiHydY,
477         I                         myTime, myIter, myThid)  
478          ENDIF
479    
480    #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
481    
482        if (K.eq.1) then        RETURN
483         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) )  
      &  *rkFac  
        ENDDO  
       ENDDO  
   
 ! ------------------------------------------------------------------------------  
       return  
       end  
 ! ==============================================================================  
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22