/[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.8.2.4 by adcroft, Thu Jan 25 19:43:32 2001 UTC revision 1.24 by jmc, Tue Dec 10 02:55:47 2002 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    CBOP
7    C     !ROUTINE: CALC_PHI_HYD
8    C     !INTERFACE:
9        SUBROUTINE CALC_PHI_HYD(        SUBROUTINE CALC_PHI_HYD(
10       I                         bi, bj, iMin, iMax, jMin, jMax, K,       I                         bi, bj, iMin, iMax, jMin, jMax, K,
11       I                         theta, salt,       I                         tFld, sFld,
12       U                         phiHyd,       U                         phiHyd,
13       I                         myThid)       I                         myThid)
14  C     /==========================================================\  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     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|
20  C     | On entry:                                                |  C     | On entry:                                                |
21  C     |   theta,salt    are the current thermodynamics quantities|  C     |   tFld,sFld     are the current thermodynamics quantities|
22  C     |                 (unchanged on exit)                      |  C     |                 (unchanged on exit)                      |
23  C     |   phiHyd(i,j,1:k-1) is the hydrostatic pressure/geopot.  |  C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |
24  C     |                 at cell centers (tracer points)          |  C     |                 at cell centers (tracer points)          |
25  C     |                 - 1:k-1 layers are valid                 |  C     |                 - 1:k-1 layers are valid                 |
26  C     |                 - k:Nr layers are invalid                |  C     |                 - k:Nr layers are invalid                |
27  C     |   phiHyd(i,j,k) is the hydrostatic pressure/geop.        |  C     |   phiHyd(i,j,k) is the hydrostatic Potential             |
28  C     |                 at cell the interface k (w point above)  |  C     |  (ocean only_^) at cell the interface k (w point above)  |
29  C     | On exit:                                                 |  C     | On exit:                                                 |
30  C     |   phiHyd(i,j,1:k) is the hydrostatic pressure/geopot.    |  C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |
31  C     |                 at cell centers (tracer points)          |  C     |                 at cell centers (tracer points)          |
32  C     |                 - 1:k layers are valid                   |  C     |                 - 1:k layers are valid                   |
33  C     |                 - k+1:Nr layers are invalid              |  C     |                 - k+1:Nr layers are invalid              |
34  C     |   phiHyd(i,j,k+1) is the hydrostatic pressure/geop.      |  C     |   phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho)   |
35  C     |                 at cell the interface k+1 (w point below)|  C     |  (ocean only-^) at cell the interface k+1 (w point below)|
36  C     |                                                          |  C     | Atmosphere:                                              |
37  C     \==========================================================/  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"
45  #include "GRID.h"  #include "GRID.h"
46  #include "EEPARAMS.h"  #include "EEPARAMS.h"
47  #include "PARAMS.h"  #include "PARAMS.h"
48    c #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 theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
60        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _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  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
65    
66    C     !LOCAL VARIABLES:
67  C     == Local variables ==  C     == Local variables ==
68        INTEGER i,j        INTEGER i,j, Kp1
69          _RL zero, one, half
70        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RL dRloc,dRlocKp1        _RL dRloc,dRlocKp1,locAlpha
72        _RL ddRm1, ddRp1, ddRm, ddRp        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
73        _RL atm_cp, atm_kappa, atm_po  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        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
107  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
# Line 69  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  C             *NOTE* The loading should go here but has not been implemented yet                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
               phiHyd(i,j,k)=0.  
125              ENDDO              ENDDO
126            ENDDO            ENDDO
127          ENDIF          ENDIF
128    
129  C       Calculate density  C       Calculate density
130          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  #ifdef ALLOW_AUTODIFF_TAMC
131       &                 theta, salt,          kkey = (ikey-1)*Nr + k
132    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
134    #endif /* ALLOW_AUTODIFF_TAMC */
135            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
136         &                 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
147  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
148              Is this directive correct or even necessary in this new code?  c           Patrick, is this directive correct or even necessary in
149    c           this new code?
150    c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...
151    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)  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)  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
173  C  C
174                
175              phiHyd(i,j,k)=phiHyd(i,j,k)+              phiHyd(i,j,k)=phiHyd(i,j,k)+
176       &          0.5*dRloc*gravity*alphaRho(i,j)       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
177              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)+
178       &          0.5*dRlocKp1*gravity*alphaRho(i,j)       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
179  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
180    
181    C---------- Compute bottom pressure deviation from gravity*rho0*H
182    C           This has to be done starting from phiHyd at the current
183    C           tracer point and .5 of the cell's thickness has to be
184    C           substracted from hFacC
185                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
186                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
187         &              + (hFacC(i,j,k,bi,bj)-.5)*drF(K)
188         &                   *gravity*alphaRho(i,j)*recip_rhoConst
189         &              + gravity*etaN(i,j,bi,bj)
190                ENDIF
191    C-----------------------------------------------------------------------
192    
193            ENDDO            ENDDO
194          ENDDO          ENDDO
195                    
196          ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
197    C       This is the hydrostatic pressure calculation for the Ocean
198    C       which uses the FIND_RHO() routine to calculate density
199    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)
205            IF (k.EQ.1) dRloc=drF(1)
206            IF (k.EQ.Nr) THEN
207              dRlocKp1=0.
208            ELSE
209              dRlocKp1=drC(k+1)
210            ENDIF
211    
212            IF (k.EQ.1) THEN
213              DO j=jMin,jMax
214                DO i=iMin,iMax
215                  phiHyd(i,j,k) = phi0surf(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    #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
235            DO j=jMin,jMax
236              DO i=iMin,iMax
237                locAlpha=alphaRho(i,j)+rhoConst
238                IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
239    
240    CmlC---------- This discretization is the "finite volume" form
241    CmlC           which has not been used to date since it does not
242    CmlC           conserve KE+PE exactly even though it is more natural
243    CmlC
244    Cml            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
245    Cml             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
246    Cml     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
247    Cml     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
248    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
256    C           which has been used since at least Adcroft et al., MWR 1997
257    C
258    
259                phiHyd(i,j,k)=phiHyd(i,j,k)+
260         &          0.5*dRloc*locAlpha
261                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
262         &          0.5*dRlocKp1*locAlpha
263    
264    C-----------------------------------------------------------------------
265    
266    C---------- Compute gravity*(sea surface elevation) first
267    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
269    C           substracted from hFacC
270                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
271                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
272         &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha
273         &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
274                ENDIF
275    C-----------------------------------------------------------------------
276    
277              ENDDO
278            ENDDO
279    
280        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
281    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
282  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
283  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
284  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
285    
286  C       Integrate d Phi / d pi  C       Integrate d Phi / d pi
287    
288  C *NOTE* These constants should be in the data file and PARAMS.h        IF (integr_GeoPot.EQ.0) THEN
289          atm_cp=1004. _d 0  C  --  Energy Conserving Form, No hFac  --
290          atm_kappa=2. _d 0/7. _d 0  C------------ The integration for the first level phi(k=1) is the same
291          atm_po=1. _d 5  C             for both the "finite volume" and energy conserving methods.
292    Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
293    C             condition is simply Phi-prime(Ro_surf)=0.
294    C           o convention ddPI > 0 (same as drF & drC)
295    C-----------------------------------------------------------------------
296          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
297            ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa)            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
298       &                  -((rF(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                ddRp=ddRp1               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
302                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.       &         +ddPIp*maskC(i,j,K,bi,bj)
303  C------------ The integration for the first level phi(k=1) is the       &               *(tFld(I,J,K,bi,bj)-tRef(K))
 C             same for both the "finite volume" and energy conserving  
 C             methods.  
 C             *NOTE* The geopotential boundary condition should go  
 C                    here but has not been implemented yet  
               phiHyd(i,j,K)=0.  
      &          -ddRp*(theta(I,J,K,bi,bj)-tRef(K))  
 C-----------------------------------------------------------------------  
304             ENDDO             ENDDO
305            ENDDO            ENDDO
306          ELSE          ELSE
307    C-------- This discretization is the energy conserving form
308  C-------- This discretization is the "finite volume" form which            ddPI=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
309  C         integrates the hydrostatic equation of each half/sub-layer.       &                 -((rC( K )/atm_Po)**atm_kappa) )*0.5
310  C         This seems most natural and could easily allow for lopped cells            DO j=jMin,jMax
311  C         by replacing rF(K) with the height of the surface (not implemented).             DO i=iMin,iMax
312  C         in the lower layers (e.g. at k=1).                phiHyd(i,j,K)=phiHyd(i,j,K-1)
313  C       &           +ddPI*maskC(i,j,K-1,bi,bj)
314  c         ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa)       &                *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
315  c    &                  -((rC(K-1)/atm_po)**atm_kappa) )       &           +ddPI*maskC(i,j, K ,bi,bj)
316  c         ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)       &                *(tFld(I,J, K ,bi,bj)-tRef( K ))
317  c    &                  -((rF( K )/atm_po)**atm_kappa) )  C             Old code (atmos-exact) looked like this
318    Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
319    Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K))
320               ENDDO
321              ENDDO
322            ENDIF
323    C end: Energy Conserving Form, No hFac  --
324  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
325    
326          ELSEIF (integr_GeoPot.EQ.1) THEN
327    C  --  Finite Volume Form, with hFac, linear in P by Half level  --
328    C---------
329    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
330    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
331    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
332    C   also: if Interface_W at the middle between tracer levels, this form
333    C     is close to the Energy Cons. form in the Interior, except for the
334    C     non-linearity in PI(p)
335    C---------
336            IF (K.EQ.1) THEN
337              ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)
338         &                  -((rC(K)/atm_Po)**atm_kappa) )
339              DO j=jMin,jMax
340               DO i=iMin,iMax
341                 phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
342         &         +ddPIp*_hFacC(I,J, K ,bi,bj)
343         &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
344               ENDDO
345              ENDDO
346            ELSE
347              ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
348         &                  -((rF( K )/atm_Po)**atm_kappa) )
349              ddPIp=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
350         &                  -((rC( K )/atm_Po)**atm_kappa) )
351              DO j=jMin,jMax
352               DO i=iMin,iMax
353                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
354         &         +ddPIm*_hFacC(I,J,K-1,bi,bj)
355         &               *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
356         &         +ddPIp*_hFacC(I,J, K ,bi,bj)
357         &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
358               ENDDO
359              ENDDO
360            ENDIF
361    C end: Finite Volume Form, with hFac, linear in P by Half level  --
362    C-----------------------------------------------------------------------
363    
364  C-------- This discretization is the energy conserving form        ELSEIF (integr_GeoPot.EQ.2) THEN
365            ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)  C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
366       &                  -((rC(K-1)/atm_po)**atm_kappa) )*0.5  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
367            ddRm1=ddRp1  C  Finite Difference formulation consistent with Partial Cell,
368    C    case Tracer level at the middle of InterFace_W
369    C    linear between 2 Tracer levels ; conserve energy in the Interior
370    C---------
371            Kp1 = min(Nr,K+1)
372            IF (K.EQ.1) THEN
373              ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
374         &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0
375              ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
376         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
377              DO j=jMin,jMax
378               DO i=iMin,iMax
379                 phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
380         &       +( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
381         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
382         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
383         &               * maskC(i,j, K ,bi,bj)
384               ENDDO
385              ENDDO
386            ELSE
387              ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
388         &                  -((rC( K )/atm_Po)**atm_kappa) )
389              ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
390         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )
391              DO j=jMin,jMax
392               DO i=iMin,iMax
393                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
394         &        + ddPIm*0.5
395         &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
396         &               * maskC(i,j,K-1,bi,bj)
397         &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
398         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
399         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
400         &               * maskC(i,j, K ,bi,bj)
401               ENDDO
402              ENDDO
403            ENDIF
404    C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
405  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
406    
407          ELSEIF (integr_GeoPot.EQ.3) THEN
408    C  --  Finite Difference Form, with hFac, Interface_W = middle  --
409    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
410    C  Finite Difference formulation consistent with Partial Cell,
411    C   Valid & accurate if Interface_W at middle between tracer levels
412    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
413    C---------
414            Kp1 = min(Nr,K+1)
415            IF (K.EQ.1) THEN
416              ratioRm=0.5*drF(K)/(rF(k)-rC(K))
417              ratioRp=drF(K)*recip_drC(Kp1)
418              ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)
419         &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0
420              ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
421         &                  -((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                ddRp=ddRp1               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)
425                ddRm=ddRm1       &       +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
426                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
427                IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0.       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
428                phiHyd(i,j,K)=phiHyd(i,j,K-1)       &               * maskC(i,j, K ,bi,bj)
429       &           -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1))             ENDDO
430       &             +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) )            ENDDO
431  C             Old code bug looked like this          ELSE
432  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1)-(ddRm1*            ratioRm=drF(K)*recip_drC(K)
433  Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj))-tRef(K))            ratioRp=drF(K)*recip_drC(Kp1)
434              ENDDO            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)
435         &                  -((rC( K )/atm_Po)**atm_kappa) )
436              ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)
437         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )
438              DO j=jMin,jMax
439               DO i=iMin,iMax
440                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
441         &        + ddPIm*0.5
442         &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
443         &               * maskC(i,j,K-1,bi,bj)
444         &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
445         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
446         &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
447         &               * maskC(i,j, K ,bi,bj)
448               ENDDO
449            ENDDO            ENDDO
450          ENDIF          ENDIF
451    C end: Finite Difference Form, with hFac, Interface_W = middle  --
452    C-----------------------------------------------------------------------
453    
454          ELSE
455            STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
456          ENDIF
457    
458    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  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
464    
465        return        RETURN
466        end        END
Legend:
Removed from v.1.8.2.4  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22