/[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.20 by mlosch, Wed Sep 18 16:38:01 2002 UTC revision 1.24 by jmc, Tue Dec 10 02:55:47 2002 UTC
# Line 34  C     |                 - k+1:Nr layers Line 34  C     |                 - k+1:Nr layers
34  C     |   phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho)   |  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)|  C     |  (ocean only-^) at cell the interface k+1 (w point below)|
36  C     | Atmosphere:                                              |  C     | Atmosphere:                                              |
37  C     |   Integr_GeoPot allows to select one integration method  |  C     |   integr_GeoPot allows to select one integration method  |
38  C     |    (see the list below)                                  |  C     |    (see the list below)                                  |
39  C     *==========================================================*  C     *==========================================================*
40  C     \ev  C     \ev
# Line 45  C     == Global variables == Line 45  C     == Global variables ==
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"  c #include "FFIELDS.h"
49  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
50  #include "tamc.h"  #include "tamc.h"
51  #include "tamc_keys.h"  #include "tamc_keys.h"
# Line 78  CEOP Line 78  CEOP
78    
79  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
80  C  Atmosphere:    C  Atmosphere:  
81  C Integr_GeoPot => select one option for the integration of the Geopotential:  C integr_GeoPot => select one option for the integration of the Geopotential:
82  C   = 0 : Energy Conserving Form, No hFac ;  C   = 0 : Energy Conserving Form, No hFac ;
83  C   = 1 : Finite Volume Form, with hFac, linear in P by Half level;  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  C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
# Line 121  C       P(z=eta) = P(atmospheric_loading Line 121  C       P(z=eta) = P(atmospheric_loading
121          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
122            DO j=jMin,jMax            DO j=jMin,jMax
123              DO i=iMin,iMax              DO i=iMin,iMax
124  #ifdef ATMOSPHERIC_LOADING                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
               phiHyd(i,j,k)=pload(i,j,bi,bj)*recip_rhoConst  
 #else  
               phiHyd(i,j,k)=0. _d 0  
 #endif  
125              ENDDO              ENDDO
126            ENDDO            ENDDO
127          ENDIF          ENDIF
128    
129  C       Calculate density  C       Calculate density
130  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
131              kkey = (ikey-1)*Nr + k          kkey = (ikey-1)*Nr + k
132  CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
133  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
134  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
135          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
136       &                 tFld, sFld,       &                 tFld, sFld,
137       &                 alphaRho, myThid)       &                 alphaRho, myThid)
138    
139    C Quasi-hydrostatic terms are added in as if they modify the buoyancy
140            IF (quasiHydrostatic) THEN
141             CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
142            ENDIF
143    
144  C       Hydrostatic pressure at cell centers  C       Hydrostatic pressure at cell centers
145          DO j=jMin,jMax          DO j=jMin,jMax
146            DO i=iMin,iMax            DO i=iMin,iMax
# Line 151  c           within the k-loop. Line 152  c           within the k-loop.
152  CADJ GENERAL  CADJ GENERAL
153  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
154    
155  C---------- This discretization is the "finite volume" form  CmlC---------- This discretization is the "finite volume" form
156  C           which has not been used to date since it does not  CmlC           which has not been used to date since it does not
157  C           conserve KE+PE exactly even though it is more natural  CmlC           conserve KE+PE exactly even though it is more natural
158  C  CmlC
159  c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  Cml          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
160  c    &              drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  Cml           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
161  c           phiHyd(i,j,k)=phiHyd(i,j,k)+  Cml     &          + hFacC(i,j,k,bi,bj)
162  c    &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  Cml     &            *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
163  C-----------------------------------------------------------------------  Cml     &          + gravity*etaN(i,j,bi,bj)
164    Cml          ENDIF
165    Cml           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
166    Cml     &         drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
167    Cml           phiHyd(i,j,k)=phiHyd(i,j,k)+
168    Cml     &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
169    CmlC-----------------------------------------------------------------------
170    
171  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
172  C           which has been used since at least Adcroft et al., MWR 1997  C           which has been used since at least Adcroft et al., MWR 1997
# Line 177  C           tracer point and .5 of the c Line 184  C           tracer point and .5 of the c
184  C           substracted from hFacC  C           substracted from hFacC
185              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
186               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
187       &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(K)       &              + (hFacC(i,j,k,bi,bj)-.5)*drF(K)
188       &                   *gravity*alphaRho(i,j)*recip_rhoConst       &                   *gravity*alphaRho(i,j)*recip_rhoConst
189       &              + gravity*etaN(i,j,bi,bj)       &              + gravity*etaN(i,j,bi,bj)
190              ENDIF              ENDIF
# Line 190  C--------------------------------------- Line 197  C---------------------------------------
197  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
198  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
199  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
200    #ifdef      ALLOW_AUTODIFF_TAMC
201    CADJ GENERAL
202    #endif      /* ALLOW_AUTODIFF_TAMC */
203    
204          dRloc=drC(k)          dRloc=drC(k)
205          IF (k.EQ.1) dRloc=drF(1)          IF (k.EQ.1) dRloc=drF(1)
# Line 202  C       before integrating g*rho over th Line 212  C       before integrating g*rho over th
212          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
213            DO j=jMin,jMax            DO j=jMin,jMax
214              DO i=iMin,iMax              DO i=iMin,iMax
215                phiHyd(i,j,k)=0.                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
               phiHyd(i,j,k)=pload(i,j,bi,bj)  
 c    &  -Ro_surf(i,j,bi,bj)*recip_rhoNil  
 c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSurfC(i,j,bi,bj) ))/1000.  
 c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSurfC(i,j,bi,bj) ))*recip_rhoNil  
216              ENDDO              ENDDO
217            ENDDO            ENDDO
218          ENDIF          ENDIF
# Line 214  c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSu Line 220  c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSu
220  C       Calculate density  C       Calculate density
221  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
222              kkey = (ikey-1)*Nr + k              kkey = (ikey-1)*Nr + k
223  CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  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  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
225  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
226          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
227       &                 tFld, sFld,       &                 tFld, sFld,
228       &                 alphaRho, myThid)       &                 alphaRho, myThid)
229    #ifdef ALLOW_AUTODIFF_TAMC
230    CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
231    #endif /* ALLOW_AUTODIFF_TAMC */
232    
233    
234  C       Hydrostatic pressure at cell centers  C       Hydrostatic pressure at cell centers
235          DO j=jMin,jMax          DO j=jMin,jMax
236            DO i=iMin,iMax            DO i=iMin,iMax
237              locAlpha=alphaRho(i,j)+rhoNil              locAlpha=alphaRho(i,j)+rhoConst
238              IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha              IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
239    
240  C---------- This discretization is the "finite volume" form  CmlC---------- This discretization is the "finite volume" form
241  C           which has not been used to date since it does not  CmlC           which has not been used to date since it does not
242  C           conserve KE+PE exactly even though it is more natural  CmlC           conserve KE+PE exactly even though it is more natural
243  C  CmlC
244  c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  Cml            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
245  c    &              drF(K)*locAlpha  Cml             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
246  c           phiHyd(i,j,k)=phiHyd(i,j,k)+  Cml     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
247  c    &          0.5*drF(K)*locAlpha  Cml     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
248  C-----------------------------------------------------------------------  Cml            ENDIF
249    Cml            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
250    Cml     &           drF(K)*locAlpha
251    Cml            phiHyd(i,j,k)=phiHyd(i,j,k)+
252    Cml     &           0.5*drF(K)*locAlpha
253    CmlC-----------------------------------------------------------------------
254    
255  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
256  C           which has been used since at least Adcroft et al., MWR 1997  C           which has been used since at least Adcroft et al., MWR 1997
257  C  C
258    
259              phiHyd(i,j,k)=phiHyd(i,j,k)+              phiHyd(i,j,k)=phiHyd(i,j,k)+
260       &          0.5*dRloc*locAlpha       &          0.5*dRloc*locAlpha
261              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
262       &          0.5*dRlocKp1*locAlpha       &          0.5*dRlocKp1*locAlpha
263    
264  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
265    
266  C---------- Compute gravity*(sea surface elevation)  C---------- Compute gravity*(sea surface elevation) first
267  C           This has to be done starting from phiHyd at the current  C           This has to be done starting from phiHyd at the current
268  C           tracer point and .5 of the cell's thickness has to be  C           tracer point and .5 of the cell's thickness has to be
269  C           substracted from hFacC  C           substracted from hFacC
# Line 268  C       the specific volume, analogous t Line 285  C       the specific volume, analogous t
285    
286  C       Integrate d Phi / d pi  C       Integrate d Phi / d pi
287    
288        IF (Integr_GeoPot.EQ.0) THEN        IF (integr_GeoPot.EQ.0) THEN
289  C  --  Energy Conserving Form, No hFac  --  C  --  Energy Conserving Form, No hFac  --
290  C------------ The integration for the first level phi(k=1) is the same  C------------ The integration for the first level phi(k=1) is the same
291  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
# Line 277  C             condition is simply Phi-pr Line 294  C             condition is simply Phi-pr
294  C           o convention ddPI > 0 (same as drF & drC)  C           o convention ddPI > 0 (same as drF & drC)
295  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
296          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
297            ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
298       &                  -((rC(K)/atm_po)**atm_kappa) )       &                  -((rC(K)/atm_Po)**atm_kappa) )
299            DO j=jMin,jMax            DO j=jMin,jMax
300             DO i=iMin,iMax             DO i=iMin,iMax
301               phiHyd(i,j,K)=               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
302       &          ddPIp*maskC(i,j,K,bi,bj)       &         +ddPIp*maskC(i,j,K,bi,bj)
303       &               *(tFld(I,J,K,bi,bj)-tRef(K))       &               *(tFld(I,J,K,bi,bj)-tRef(K))
304             ENDDO             ENDDO
305            ENDDO            ENDDO
306          ELSE          ELSE
307  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
308            ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)            ddPI=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
309       &                 -((rC( K )/atm_po)**atm_kappa) )*0.5       &                 -((rC( K )/atm_Po)**atm_kappa) )*0.5
310            DO j=jMin,jMax            DO j=jMin,jMax
311             DO i=iMin,iMax             DO i=iMin,iMax
312                phiHyd(i,j,K)=phiHyd(i,j,K-1)                phiHyd(i,j,K)=phiHyd(i,j,K-1)
# Line 306  Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I, Line 323  Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I,
323  C end: Energy Conserving Form, No hFac  --  C end: Energy Conserving Form, No hFac  --
324  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
325    
326        ELSEIF (Integr_GeoPot.EQ.1) THEN        ELSEIF (integr_GeoPot.EQ.1) THEN
327  C  --  Finite Volume Form, with hFac, linear in P by Half level  --  C  --  Finite Volume Form, with hFac, linear in P by Half level  --
328  C---------  C---------
329  C  Finite Volume formulation consistent with Partial Cell, linear in p by piece  C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
# Line 317  C     is close to the Energy Cons. form Line 334  C     is close to the Energy Cons. form
334  C     non-linearity in PI(p)  C     non-linearity in PI(p)
335  C---------  C---------
336          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
337            ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
338       &                  -((rC(K)/atm_po)**atm_kappa) )       &                  -((rC(K)/atm_Po)**atm_kappa) )
339            DO j=jMin,jMax            DO j=jMin,jMax
340             DO i=iMin,iMax             DO i=iMin,iMax
341               phiHyd(i,j,K) =               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
342       &          ddPIp*_hFacC(I,J, K ,bi,bj)       &         +ddPIp*_hFacC(I,J, K ,bi,bj)
343       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
344             ENDDO             ENDDO
345            ENDDO            ENDDO
346          ELSE          ELSE
347            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
348       &                  -((rF( K )/atm_po)**atm_kappa) )       &                  -((rF( K )/atm_Po)**atm_kappa) )
349            ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
350       &                  -((rC( K )/atm_po)**atm_kappa) )       &                  -((rC( K )/atm_Po)**atm_kappa) )
351            DO j=jMin,jMax            DO j=jMin,jMax
352             DO i=iMin,iMax             DO i=iMin,iMax
353               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
# Line 344  C--------- Line 361  C---------
361  C end: Finite Volume Form, with hFac, linear in P by Half level  --  C end: Finite Volume Form, with hFac, linear in P by Half level  --
362  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
363    
364        ELSEIF (Integr_GeoPot.EQ.2) THEN        ELSEIF (integr_GeoPot.EQ.2) THEN
365  C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --  C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
366  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
367  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
# Line 353  C    linear between 2 Tracer levels ; co Line 370  C    linear between 2 Tracer levels ; co
370  C---------  C---------
371          Kp1 = min(Nr,K+1)          Kp1 = min(Nr,K+1)
372          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
373            ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)            ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
374       &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0       &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0
375            ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
376       &                  -((rC(Kp1)/atm_po)**atm_kappa) )         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
377            DO j=jMin,jMax            DO j=jMin,jMax
378             DO i=iMin,iMax             DO i=iMin,iMax
379               phiHyd(i,j,K) =               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
380       &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)       &       +( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
381       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
382       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
383       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
384             ENDDO             ENDDO
385            ENDDO            ENDDO
386          ELSE          ELSE
387            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
388       &                  -((rC( K )/atm_po)**atm_kappa) )       &                  -((rC( K )/atm_Po)**atm_kappa) )
389            ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
390       &                  -((rC(Kp1)/atm_po)**atm_kappa) )       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )
391            DO j=jMin,jMax            DO j=jMin,jMax
392             DO i=iMin,iMax             DO i=iMin,iMax
393               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
# Line 387  C--------- Line 404  C---------
404  C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --  C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
405  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
406    
407        ELSEIF (Integr_GeoPot.EQ.3) THEN        ELSEIF (integr_GeoPot.EQ.3) THEN
408  C  --  Finite Difference Form, with hFac, Interface_W = middle  --  C  --  Finite Difference Form, with hFac, Interface_W = middle  --
409  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
410  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
# Line 398  C--------- Line 415  C---------
415          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
416            ratioRm=0.5*drF(K)/(rF(k)-rC(K))            ratioRm=0.5*drF(K)/(rF(k)-rC(K))
417            ratioRp=drF(K)*recip_drC(Kp1)            ratioRp=drF(K)*recip_drC(Kp1)
418            ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)            ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
419       &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0       &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0
420            ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
421       &                  -((rC(Kp1)/atm_po)**atm_kappa) )         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
422            DO j=jMin,jMax            DO j=jMin,jMax
423             DO i=iMin,iMax             DO i=iMin,iMax
424               phiHyd(i,j,K) =               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
425       &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)       &       +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
426       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
427       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
428       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
# Line 414  C--------- Line 431  C---------
431          ELSE          ELSE
432            ratioRm=drF(K)*recip_drC(K)            ratioRm=drF(K)*recip_drC(K)
433            ratioRp=drF(K)*recip_drC(Kp1)            ratioRp=drF(K)*recip_drC(Kp1)
434            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
435       &                  -((rC( K )/atm_po)**atm_kappa) )       &                  -((rC( K )/atm_Po)**atm_kappa) )
436            ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
437       &                  -((rC(Kp1)/atm_po)**atm_kappa) )       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )
438            DO j=jMin,jMax            DO j=jMin,jMax
439             DO i=iMin,iMax             DO i=iMin,iMax
440               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
# Line 435  C end: Finite Difference Form, with hFac Line 452  C end: Finite Difference Form, with hFac
452  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
453    
454        ELSE        ELSE
455          STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !'          STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
456        ENDIF        ENDIF
457    
458  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
459        ELSE        ELSE
460          STOP 'CALC_PHI_HYD: We should never reach this point!'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
461        ENDIF        ENDIF
462    
463  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22