/[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.13 by heimbach, Mon May 14 21:51:24 2001 UTC revision 1.14 by jmc, Fri Jul 6 21:47:00 2001 UTC
# Line 20  C     |                 at cell centers Line 20  C     |                 at cell centers
20  C     |                 - 1:k-1 layers are valid                 |  C     |                 - 1:k-1 layers are valid                 |
21  C     |                 - k:Nr layers are invalid                |  C     |                 - k:Nr layers are invalid                |
22  C     |   phiHyd(i,j,k) is the hydrostatic Potential             |  C     |   phiHyd(i,j,k) is the hydrostatic Potential             |
23  C     |                 at cell the interface k (w point above)  |  C     |  (ocean only_^) at cell the interface k (w point above)  |
24  C     | On exit:                                                 |  C     | On exit:                                                 |
25  C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |  C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |
26  C     |                 at cell centers (tracer points)          |  C     |                 at cell centers (tracer points)          |
27  C     |                 - 1:k layers are valid                   |  C     |                 - 1:k layers are valid                   |
28  C     |                 - k+1:Nr layers are invalid              |  C     |                 - k+1:Nr layers are invalid              |
29  C     |   phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho)   |  C     |   phiHyd(i,j,k+1) is the hydrostatic Potential (P/rho)   |
30  C     |                 at cell the interface k+1 (w point below)|  C     |  (ocean only-^) at cell the interface k+1 (w point below)|
31  C     |                                                          |  C     | Atmosphere:                                              |
32    C     |   Integr_GeoPot allows to select one integration method  |
33    C     |    (see the list below)                                  |
34  C     \==========================================================/  C     \==========================================================/
35        IMPLICIT NONE        IMPLICIT NONE
36  C     == Global variables ==  C     == Global variables ==
# Line 47  C     == Routine arguments == Line 49  C     == Routine arguments ==
49        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
50        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        INTEGER myThid        INTEGER myThid
52          
53  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
54    
55  C     == Local variables ==  C     == Local variables ==
56        INTEGER i,j        INTEGER i,j, Kp1
57          _RL zero, one, half
58        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59        _RL dRloc,dRlocKp1        _RL dRloc,dRlocKp1
60        _RL ddRm1, ddRp1, ddRm, ddRp        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
61        _RL atm_cp, atm_kappa, atm_po  
62          zero = 0. _d 0
63          one  = 1. _d 0
64          half = .5 _d 0
65    
66    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
67    C  Atmosphere:  
68    C Integr_GeoPot => select one option for the integration of the Geopotential:
69    C   = 0 : Energy Conserving Form, No hFac ;
70    C   = 1 : Finite Volume Form, with hFac, linear in P by Half level;
71    C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
72    C     2 : case Tracer level at the middle of InterFace_W;
73    C     3 : case InterFace_W  at the middle of Tracer levels;
74    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
75    
76  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
77            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
# Line 112  C       Hydrostatic pressure at cell cen Line 128  C       Hydrostatic pressure at cell cen
128          DO j=jMin,jMax          DO j=jMin,jMax
129            DO i=iMin,iMax            DO i=iMin,iMax
130  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
131  c           Patrick, is this directive correct or even necessary in  c           Patrick, is this directive correct or even necessary in
132  c           this new code?  c           this new code?
133  c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...  c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...
134  c           within the k-loop.  c           within the k-loop.
# Line 143  C--------------------------------------- Line 159  C---------------------------------------
159    
160    
161        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
162    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
163  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
164  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
165  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
166    
167  C       Integrate d Phi / d pi  C       Integrate d Phi / d pi
168    
169  C *NOTE* These constants should be in the data file and PARAMS.h        IF (Integr_GeoPot.EQ.0) THEN
170          atm_cp=1004. _d 0  C  --  Energy Conserving Form, No hFac  --
171          atm_kappa=2. _d 0/7. _d 0  C------------ The integration for the first level phi(k=1) is the same
172          atm_po=1. _d 5  C             for both the "finite volume" and energy conserving methods.
173    C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
174    C             condition is simply Phi'(Ro_surf)=0.
175    C           o convention ddPI > 0 (same as drF & drC)
176    C-----------------------------------------------------------------------
177          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
178            ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa)            ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
179       &                  -((rF(K)/atm_po)**atm_kappa) )       &                  -((rC(K)/atm_po)**atm_kappa) )
180            DO j=jMin,jMax            DO j=jMin,jMax
181              DO i=iMin,iMax             DO i=iMin,iMax
182                ddRp=ddRp1               phiHyd(i,j,K)=
183                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.       &          ddPIp*maskC(i,j,K,bi,bj)
184  C------------ The integration for the first level phi(k=1) is the       &               *(theta(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-----------------------------------------------------------------------  
185             ENDDO             ENDDO
186            ENDDO            ENDDO
187          ELSE          ELSE
188    C-------- This discretization is the energy conserving form
189  C-------- This discretization is the "finite volume" form which            ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
190  C         integrates the hydrostatic equation of each half/sub-layer.       &                 -((rC( K )/atm_po)**atm_kappa) )*0.5
191  C         This seems most natural and could easily allow for lopped cells            DO j=jMin,jMax
192  C         by replacing rF(K) with the height of the surface (not implemented).             DO i=iMin,iMax
193  C         in the lower layers (e.g. at k=1).                phiHyd(i,j,K)=phiHyd(i,j,K-1)
194  C       &           +ddPI*maskC(i,j,K-1,bi,bj)
195  c         ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa)       &                *(theta(I,J,K-1,bi,bj)-tRef(K-1))
196  c    &                  -((rC(K-1)/atm_po)**atm_kappa) )       &           +ddPI*maskC(i,j, K ,bi,bj)
197  c         ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)       &                *(theta(I,J, K ,bi,bj)-tRef( K ))
198  c    &                  -((rF( K )/atm_po)**atm_kappa) )  C             Old code (atmos-exact) looked like this
199    Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
200    Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))
201               ENDDO
202              ENDDO
203            ENDIF
204    C end: Energy Conserving Form, No hFac  --
205  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
206    
207          ELSEIF (Integr_GeoPot.EQ.1) THEN
208    C  --  Finite Volume Form, with hFac, linear in P by Half level  --
209    C---------
210    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
211    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
212    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
213    C   also: if Interface_W at the middle between tracer levels, this form
214    C     is close to the Energy Cons. form in the Interior, except for the
215    C     non-linearity in PI(p)
216    C---------
217            IF (K.EQ.1) THEN
218              ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
219         &                  -((rC(K)/atm_po)**atm_kappa) )
220              DO j=jMin,jMax
221               DO i=iMin,iMax
222                 phiHyd(i,j,K) =
223         &          ddPIp*hFacC(I,J, K ,bi,bj)
224         &               *(theta(I,J, K ,bi,bj)-tRef( K ))
225               ENDDO
226              ENDDO
227            ELSE
228              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
229         &                  -((rF( K )/atm_po)**atm_kappa) )
230              ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
231         &                  -((rC( K )/atm_po)**atm_kappa) )
232              DO j=jMin,jMax
233               DO i=iMin,iMax
234                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
235         &         +ddPIm*hFacC(I,J,K-1,bi,bj)
236         &               *(theta(I,J,K-1,bi,bj)-tRef(K-1))
237         &         +ddPIp*hFacC(I,J, K ,bi,bj)
238         &               *(theta(I,J, K ,bi,bj)-tRef( K ))
239               ENDDO
240              ENDDO
241            ENDIF
242    C end: Finite Volume Form, with hFac, linear in P by Half level  --
243    C-----------------------------------------------------------------------
244    
245  C-------- This discretization is the energy conserving form        ELSEIF (Integr_GeoPot.EQ.2) THEN
246            ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)  C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
247       &                  -((rC(K-1)/atm_po)**atm_kappa) )*0.5  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
248            ddRm1=ddRp1  C  Finite Difference formulation consistent with Partial Cell,
249    C    case Tracer level at the middle of InterFace_W
250    C    linear between 2 Tracer levels ; conserve energy in the Interior
251    C---------
252            Kp1 = min(Nr,K+1)
253            IF (K.EQ.1) THEN
254              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
255         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
256              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
257         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
258              DO j=jMin,jMax
259               DO i=iMin,iMax
260                 phiHyd(i,j,K) =
261         &        ( ddPIm*max(zero, hFacC(i,j,K,bi,bj)-half)
262         &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)-half) )
263         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
264         &               * maskC(i,j, K ,bi,bj)
265               ENDDO
266              ENDDO
267            ELSE
268              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
269         &                  -((rC( K )/atm_po)**atm_kappa) )
270              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
271         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
272              DO j=jMin,jMax
273               DO i=iMin,iMax
274                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
275         &        + ddPIm*0.5
276         &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))
277         &               * maskC(i,j,K-1,bi,bj)
278         &        +(ddPIm*max(zero, hFacC(i,j,K,bi,bj)-half)
279         &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)-half) )
280         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
281         &               * maskC(i,j, K ,bi,bj)
282               ENDDO
283              ENDDO
284            ENDIF
285    C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
286  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
287    
288          ELSEIF (Integr_GeoPot.EQ.3) THEN
289    C  --  Finite Difference Form, with hFac, Interface_W = middle  --
290    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
291    C  Finite Difference formulation consistent with Partial Cell,
292    C   Valid & accurate if Interface_W at middle between tracer levels
293    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
294    C---------
295            Kp1 = min(Nr,K+1)
296            IF (K.EQ.1) THEN
297              ratioRm=0.5*drF(K)/(rF(k)-rC(K))
298              ratioRp=drF(K)*recip_drC(Kp1)
299              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
300         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
301              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
302         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
303            DO j=jMin,jMax            DO j=jMin,jMax
304              DO i=iMin,iMax             DO i=iMin,iMax
305                ddRp=ddRp1               phiHyd(i,j,K) =
306                ddRm=ddRm1       &        ( ddPIm*max(zero,(hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
307                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.       &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)*ratioRp     -half) )
308                IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0.       &               *(theta(i,j, K ,bi,bj)-tRef( K ))
309                phiHyd(i,j,K)=phiHyd(i,j,K-1)       &               * maskC(i,j, K ,bi,bj)
310       &           -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1))             ENDDO
311       &             +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) )            ENDDO
312  C             Old code (atmos-exact) looked like this          ELSE
313  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddRm1*            ratioRm=drF(K)*recip_drC(K)
314  Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))            ratioRp=drF(K)*recip_drC(Kp1)
315              ENDDO            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
316         &                  -((rC( K )/atm_po)**atm_kappa) )
317              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
318         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
319              DO j=jMin,jMax
320               DO i=iMin,iMax
321                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
322         &        + ddPIm*0.5
323         &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))
324         &               * maskC(i,j,K-1,bi,bj)
325         &        +(ddPIm*max(zero,(hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
326         &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)*ratioRp     -half) )
327         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
328         &               * maskC(i,j, K ,bi,bj)
329               ENDDO
330            ENDDO            ENDDO
331          ENDIF          ENDIF
332    C end: Finite Difference Form, with hFac, Interface_W = middle  --
333    C-----------------------------------------------------------------------
334    
335          ELSE
336            STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !'
337          ENDIF
338    
339    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
340        ELSE        ELSE
341          STOP 'CALC_PHI_HYD: We should never reach this point!'          STOP 'CALC_PHI_HYD: We should never reach this point!'
342        ENDIF        ENDIF
343    
344  #endif  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
345    
346        RETURN        RETURN
347        END        END
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22