/[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.21 by mlosch, Wed Sep 25 19:36:50 2002 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                         phiHyd,
13         I                         myThid)
14    C     !DESCRIPTION: \bv
15    C     *==========================================================*
16  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
17  C     | o Integrate the hydrostatic relation to find phiHyd.     |  C     | o Integrate the hydrostatic relation to find the Hydros. |
18  C     |                                                          |  C     *==========================================================*
19  C     \==========================================================/  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|
20    C     | On entry:                                                |
21    C     |   tFld,sFld     are the current thermodynamics quantities|
22    C     |                 (unchanged on exit)                      |
23    C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |
24    C     |                 at cell centers (tracer points)          |
25    C     |                 - 1:k-1 layers are valid                 |
26    C     |                 - k:Nr layers are invalid                |
27    C     |   phiHyd(i,j,k) is the hydrostatic Potential             |
28    C     |  (ocean only_^) at cell the interface k (w point above)  |
29    C     | On exit:                                                 |
30    C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |
31    C     |                 at cell centers (tracer points)          |
32    C     |                 - 1:k layers are valid                   |
33    C     |                 - k+1:Nr layers are invalid              |
34    C     |   phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho)   |
35    C     |  (ocean only-^) at cell the interface k+1 (w point below)|
36    C     | Atmosphere:                                              |
37    C     |   Integr_GeoPot allows to select one integration method  |
38    C     |    (see the list below)                                  |
39    C     *==========================================================*
40    C     \ev
41    C     !USES:
42        IMPLICIT NONE        IMPLICIT NONE
43  C     == Global variables ==  C     == Global variables ==
44  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
45  #include "GRID.h"  #include "GRID.h"
46  #include "EEPARAMS.h"  #include "EEPARAMS.h"
47  #include "PARAMS.h"  #include "PARAMS.h"
48    #include "FFIELDS.h"
49    #ifdef ALLOW_AUTODIFF_TAMC
50    #include "tamc.h"
51    #include "tamc_keys.h"
52    #endif /* ALLOW_AUTODIFF_TAMC */
53    #include "SURFACE.h"
54    #include "DYNVARS.h"
55    
56    C     !INPUT/OUTPUT PARAMETERS:
57  C     == Routine arguments ==  C     == Routine arguments ==
58        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax,K
59        _RL buoyKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
60        _RL buoyKP1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
61        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
62        integer myThid        INTEGER myThid
63          
64    #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
65    
66    C     !LOCAL VARIABLES:
67  C     == Local variables ==  C     == Local variables ==
68        INTEGER i,j,Km1        INTEGER i,j, Kp1
69        _RL halfLayer        _RL zero, one, half
70          _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71          _RL dRloc,dRlocKp1,locAlpha
72          _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
73    CEOP
74    
75          zero = 0. _d 0
76          one  = 1. _d 0
77          half = .5 _d 0
78    
79    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
80    C  Atmosphere:  
81    C Integr_GeoPot => select one option for the integration of the Geopotential:
82    C   = 0 : Energy Conserving Form, No hFac ;
83    C   = 1 : Finite Volume Form, with hFac, linear in P by Half level;
84    C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
85    C     2 : case Tracer level at the middle of InterFace_W;
86    C     3 : case InterFace_W  at the middle of Tracer levels;
87    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
88    
89    #ifdef ALLOW_AUTODIFF_TAMC
90              act1 = bi - myBxLo(myThid)
91              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
92    
93              act2 = bj - myByLo(myThid)
94              max2 = myByHi(myThid) - myByLo(myThid) + 1
95    
96              act3 = myThid - 1
97              max3 = nTx*nTy
98    
99              act4 = ikey_dynamics - 1
100    
101              ikey = (act1 + 1) + act2*max1
102         &                      + act3*max1*max2
103         &                      + act4*max1*max2*max3
104    #endif /* ALLOW_AUTODIFF_TAMC */
105    
106          IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
107    C       This is the hydrostatic pressure calculation for the Ocean
108    C       which uses the FIND_RHO() routine to calculate density
109    C       before integrating g*rho over the current layer/interface
110    
111            dRloc=drC(k)
112            IF (k.EQ.1) dRloc=drF(1)
113            IF (k.EQ.Nr) THEN
114              dRlocKp1=0.
115            ELSE
116              dRlocKp1=drC(k+1)
117            ENDIF
118    
119    C--     If this is the top layer we impose the boundary condition
120    C       P(z=eta) = P(atmospheric_loading)
121            IF (k.EQ.1) THEN
122              DO j=jMin,jMax
123                DO i=iMin,iMax
124    #ifdef ATMOSPHERIC_LOADING
125                  phiHyd(i,j,k)=pload(i,j,bi,bj)*recip_rhoConst
126    #else
127                  phiHyd(i,j,k)=0. _d 0
128    #endif
129                ENDDO
130              ENDDO
131            ENDIF
132    
133    C       Calculate density
134    #ifdef ALLOW_AUTODIFF_TAMC
135                kkey = (ikey-1)*Nr + k
136    CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
137    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
138    #endif /* ALLOW_AUTODIFF_TAMC */
139            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
140         &                 tFld, sFld,
141         &                 alphaRho, myThid)
142    
143    C       Hydrostatic pressure at cell centers
144            DO j=jMin,jMax
145              DO i=iMin,iMax
146    #ifdef      ALLOW_AUTODIFF_TAMC
147    c           Patrick, is this directive correct or even necessary in
148    c           this new code?
149    c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...
150    c           within the k-loop.
151    CADJ GENERAL
152    #endif      /* ALLOW_AUTODIFF_TAMC */
153    
154    CmlC---------- This discretization is the "finite volume" form
155    CmlC           which has not been used to date since it does not
156    CmlC           conserve KE+PE exactly even though it is more natural
157    CmlC
158    Cml          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
159    Cml           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
160    Cml     &          + hFacC(i,j,k,bi,bj)
161    Cml     &            *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
162    Cml     &          + gravity*etaN(i,j,bi,bj)
163    Cml          ENDIF
164    Cml           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
165    Cml     &         drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
166    Cml           phiHyd(i,j,k)=phiHyd(i,j,k)+
167    Cml     &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
168    CmlC-----------------------------------------------------------------------
169    
170    C---------- This discretization is the "energy conserving" form
171    C           which has been used since at least Adcroft et al., MWR 1997
172    C
173                
174                phiHyd(i,j,k)=phiHyd(i,j,k)+
175         &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
176                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
177         &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
178    C-----------------------------------------------------------------------
179    
180    C---------- Compute bottom pressure deviation from gravity*rho0*H
181    C           This has to be done starting from phiHyd at the current
182    C           tracer point and .5 of the cell's thickness has to be
183    C           substracted from hFacC
184                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
185                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
186         &              + (hFacC(i,j,k,bi,bj)-.5)*drF(K)
187         &                   *gravity*alphaRho(i,j)*recip_rhoConst
188         &              + gravity*etaN(i,j,bi,bj)
189                ENDIF
190    C-----------------------------------------------------------------------
191    
192              ENDDO
193            ENDDO
194            
195          ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
196    C       This is the hydrostatic pressure calculation for the Ocean
197    C       which uses the FIND_RHO() routine to calculate density
198    C       before integrating g*rho over the current layer/interface
199    #ifdef      ALLOW_AUTODIFF_TAMC
200    CADJ GENERAL
201    #endif      /* ALLOW_AUTODIFF_TAMC */
202    
203            dRloc=drC(k)
204            IF (k.EQ.1) dRloc=drF(1)
205            IF (k.EQ.Nr) THEN
206              dRlocKp1=0.
207            ELSE
208              dRlocKp1=drC(k+1)
209            ENDIF
210    
211            IF (k.EQ.1) THEN
212              DO j=jMin,jMax
213                DO i=iMin,iMax
214                  phiHyd(i,j,k)=0.
215                  phiHyd(i,j,k)=pload(i,j,bi,bj)
216                ENDDO
217              ENDDO
218            ENDIF
219    
220    C       Calculate density
221    #ifdef ALLOW_AUTODIFF_TAMC
222                kkey = (ikey-1)*Nr + k
223    CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
224    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
225    #endif /* ALLOW_AUTODIFF_TAMC */
226            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
227         &                 tFld, sFld,
228         &                 alphaRho, myThid)
229    
230    C       Hydrostatic pressure at cell centers
231            DO j=jMin,jMax
232              DO i=iMin,iMax
233                locAlpha=alphaRho(i,j)+rhoConst
234                IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
235    
236    CmlC---------- This discretization is the "finite volume" form
237    CmlC           which has not been used to date since it does not
238    CmlC           conserve KE+PE exactly even though it is more natural
239    CmlC
240    Cml            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
241    Cml             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
242    Cml     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
243    Cml     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
244    Cml            ENDIF
245    Cml            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
246    Cml     &           drF(K)*locAlpha
247    Cml            phiHyd(i,j,k)=phiHyd(i,j,k)+
248    Cml     &           0.5*drF(K)*locAlpha
249    CmlC-----------------------------------------------------------------------
250    
251    C---------- This discretization is the "energy conserving" form
252    C           which has been used since at least Adcroft et al., MWR 1997
253    C
254    
255                phiHyd(i,j,k)=phiHyd(i,j,k)+
256         &          0.5*dRloc*locAlpha
257                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
258         &          0.5*dRlocKp1*locAlpha
259    
260    C-----------------------------------------------------------------------
261    
262    C---------- Compute gravity*(sea surface elevation) first
263    C           This has to be done starting from phiHyd at the current
264    C           tracer point and .5 of the cell's thickness has to be
265    C           substracted from hFacC
266                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
267                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
268         &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha
269         &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
270                ENDIF
271    C-----------------------------------------------------------------------
272    
273              ENDDO
274            ENDDO
275    
276          ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
277    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
278    C       This is the hydrostatic geopotential calculation for the Atmosphere
279    C       The ideal gas law is used implicitly here rather than calculating
280    C       the specific volume, analogous to the oceanic case.
281    
282    C       Integrate d Phi / d pi
283    
284          IF (Integr_GeoPot.EQ.0) THEN
285    C  --  Energy Conserving Form, No hFac  --
286    C------------ The integration for the first level phi(k=1) is the same
287    C             for both the "finite volume" and energy conserving methods.
288    Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
289    C             condition is simply Phi-prime(Ro_surf)=0.
290    C           o convention ddPI > 0 (same as drF & drC)
291    C-----------------------------------------------------------------------
292            IF (K.EQ.1) THEN
293              ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
294         &                  -((rC(K)/atm_po)**atm_kappa) )
295              DO j=jMin,jMax
296               DO i=iMin,iMax
297                 phiHyd(i,j,K)=
298         &          ddPIp*maskC(i,j,K,bi,bj)
299         &               *(tFld(I,J,K,bi,bj)-tRef(K))
300               ENDDO
301              ENDDO
302            ELSE
303    C-------- This discretization is the energy conserving form
304              ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
305         &                 -((rC( K )/atm_po)**atm_kappa) )*0.5
306              DO j=jMin,jMax
307               DO i=iMin,iMax
308                  phiHyd(i,j,K)=phiHyd(i,j,K-1)
309         &           +ddPI*maskC(i,j,K-1,bi,bj)
310         &                *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
311         &           +ddPI*maskC(i,j, K ,bi,bj)
312         &                *(tFld(I,J, K ,bi,bj)-tRef( K ))
313    C             Old code (atmos-exact) looked like this
314    Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
315    Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K))
316               ENDDO
317              ENDDO
318            ENDIF
319    C end: Energy Conserving Form, No hFac  --
320    C-----------------------------------------------------------------------
321    
322          ELSEIF (Integr_GeoPot.EQ.1) THEN
323    C  --  Finite Volume Form, with hFac, linear in P by Half level  --
324    C---------
325    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
326    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
327    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
328    C   also: if Interface_W at the middle between tracer levels, this form
329    C     is close to the Energy Cons. form in the Interior, except for the
330    C     non-linearity in PI(p)
331    C---------
332            IF (K.EQ.1) THEN
333              ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
334         &                  -((rC(K)/atm_po)**atm_kappa) )
335              DO j=jMin,jMax
336               DO i=iMin,iMax
337                 phiHyd(i,j,K) =
338         &          ddPIp*_hFacC(I,J, K ,bi,bj)
339         &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
340               ENDDO
341              ENDDO
342            ELSE
343              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
344         &                  -((rF( K )/atm_po)**atm_kappa) )
345              ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
346         &                  -((rC( K )/atm_po)**atm_kappa) )
347              DO j=jMin,jMax
348               DO i=iMin,iMax
349                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
350         &         +ddPIm*_hFacC(I,J,K-1,bi,bj)
351         &               *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
352         &         +ddPIp*_hFacC(I,J, K ,bi,bj)
353         &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
354               ENDDO
355              ENDDO
356            ENDIF
357    C end: Finite Volume Form, with hFac, linear in P by Half level  --
358    C-----------------------------------------------------------------------
359    
360          ELSEIF (Integr_GeoPot.EQ.2) THEN
361    C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
362    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
363    C  Finite Difference formulation consistent with Partial Cell,
364    C    case Tracer level at the middle of InterFace_W
365    C    linear between 2 Tracer levels ; conserve energy in the Interior
366    C---------
367            Kp1 = min(Nr,K+1)
368            IF (K.EQ.1) THEN
369              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
370         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
371              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
372         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
373              DO j=jMin,jMax
374               DO i=iMin,iMax
375                 phiHyd(i,j,K) =
376         &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
377         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
378         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
379         &               * maskC(i,j, K ,bi,bj)
380               ENDDO
381              ENDDO
382            ELSE
383              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
384         &                  -((rC( K )/atm_po)**atm_kappa) )
385              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
386         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
387              DO j=jMin,jMax
388               DO i=iMin,iMax
389                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
390         &        + ddPIm*0.5
391         &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
392         &               * maskC(i,j,K-1,bi,bj)
393         &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
394         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
395         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
396         &               * maskC(i,j, K ,bi,bj)
397               ENDDO
398              ENDDO
399            ENDIF
400    C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
401    C-----------------------------------------------------------------------
402    
403          ELSEIF (Integr_GeoPot.EQ.3) THEN
404    C  --  Finite Difference Form, with hFac, Interface_W = middle  --
405    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
406    C  Finite Difference formulation consistent with Partial Cell,
407    C   Valid & accurate if Interface_W at middle between tracer levels
408    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
409    C---------
410            Kp1 = min(Nr,K+1)
411            IF (K.EQ.1) THEN
412              ratioRm=0.5*drF(K)/(rF(k)-rC(K))
413              ratioRp=drF(K)*recip_drC(Kp1)
414              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
415         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
416              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
417         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
418              DO j=jMin,jMax
419               DO i=iMin,iMax
420                 phiHyd(i,j,K) =
421         &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
422         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
423         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
424         &               * maskC(i,j, K ,bi,bj)
425               ENDDO
426              ENDDO
427            ELSE
428              ratioRm=drF(K)*recip_drC(K)
429              ratioRp=drF(K)*recip_drC(Kp1)
430              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
431         &                  -((rC( K )/atm_po)**atm_kappa) )
432              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
433         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
434              DO j=jMin,jMax
435               DO i=iMin,iMax
436                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
437         &        + ddPIm*0.5
438         &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
439         &               * maskC(i,j,K-1,bi,bj)
440         &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
441         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
442         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
443         &               * maskC(i,j, K ,bi,bj)
444               ENDDO
445              ENDDO
446            ENDIF
447    C end: Finite Difference Form, with hFac, Interface_W = middle  --
448    C-----------------------------------------------------------------------
449    
450          ELSE
451            STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !'
452          ENDIF
453    
454    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
455          ELSE
456            STOP 'CALC_PHI_HYD: We should never reach this point!'
457          ENDIF
458    
459    #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
460    
461        if (K.eq.1) then        RETURN
462         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.21
  ViewVC Help
Powered by ViewVC 1.1.22