/[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.11 by jmc, Thu Mar 8 20:21:34 2001 UTC revision 1.16 by cnh, Wed Sep 26 18:09:14 2001 UTC
# Line 3  C $Name$ Line 3  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                         theta, salt,
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 the Hydros. |  C     | o Integrate the hydrostatic relation to find the Hydros. |
18    C     *==========================================================*
19  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|
20  C     | On entry:                                                |  C     | On entry:                                                |
21  C     |   theta,salt    are the current thermodynamics quantities|  C     |   theta,salt    are the current thermodynamics quantities|
# Line 20  C     |                 at cell centers Line 25  C     |                 at cell centers
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 Potential             |  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 Potential           |  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 Potential (P/rho)   |  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    #ifdef ALLOW_AUTODIFF_TAMC
49    #include "tamc.h"
50    #include "tamc_keys.h"
51    #endif /* ALLOW_AUTODIFF_TAMC */
52    
53    C     !INPUT/OUTPUT PARAMETERS:
54  C     == Routine arguments ==  C     == Routine arguments ==
55        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax,K
56        _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
57        _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)
58        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59        INTEGER myThid        INTEGER myThid
60          
61  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
62    
63    C     !LOCAL VARIABLES:
64  C     == Local variables ==  C     == Local variables ==
65        INTEGER i,j        INTEGER i,j, Kp1
66          _RL zero, one, half
67        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68        _RL dRloc,dRlocKp1        _RL dRloc,dRlocKp1
69        _RL ddRm1, ddRp1, ddRm, ddRp        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
70        _RL atm_cp, atm_kappa, atm_po  CEOP
71    
72          zero = 0. _d 0
73          one  = 1. _d 0
74          half = .5 _d 0
75    
76    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
77    C  Atmosphere:  
78    C Integr_GeoPot => select one option for the integration of the Geopotential:
79    C   = 0 : Energy Conserving Form, No hFac ;
80    C   = 1 : Finite Volume Form, with hFac, linear in P by Half level;
81    C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
82    C     2 : case Tracer level at the middle of InterFace_W;
83    C     3 : case InterFace_W  at the middle of Tracer levels;
84    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
85    
86    #ifdef ALLOW_AUTODIFF_TAMC
87              act1 = bi - myBxLo(myThid)
88              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
89    
90              act2 = bj - myByLo(myThid)
91              max2 = myByHi(myThid) - myByLo(myThid) + 1
92    
93              act3 = myThid - 1
94              max3 = nTx*nTy
95    
96              act4 = ikey_dynamics - 1
97    
98              ikey = (act1 + 1) + act2*max1
99         &                      + act3*max1*max2
100         &                      + act4*max1*max2*max3
101    #endif /* ALLOW_AUTODIFF_TAMC */
102    
103        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
104  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
# Line 77  C             *NOTE* The loading should Line 125  C             *NOTE* The loading should
125          ENDIF          ENDIF
126    
127  C       Calculate density  C       Calculate density
128    #ifdef ALLOW_AUTODIFF_TAMC
129                kkey = (ikey-1)*Nr + k
130    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
131    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
132    #endif /* ALLOW_AUTODIFF_TAMC */
133          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
134       &                 theta, salt,       &                 theta, salt,
135       &                 alphaRho, myThid)       &                 alphaRho, myThid)
# Line 85  C       Hydrostatic pressure at cell cen Line 138  C       Hydrostatic pressure at cell cen
138          DO j=jMin,jMax          DO j=jMin,jMax
139            DO i=iMin,iMax            DO i=iMin,iMax
140  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
141              Is this directive correct or even necessary in this new code?  c           Patrick, is this directive correct or even necessary in
142    c           this new code?
143    c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...
144    c           within the k-loop.
145  CADJ GENERAL  CADJ GENERAL
146  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
147    
# Line 113  C--------------------------------------- Line 169  C---------------------------------------
169    
170    
171        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
172    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
173  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
174  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
175  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
176    
177  C       Integrate d Phi / d pi  C       Integrate d Phi / d pi
178    
179  C *NOTE* These constants should be in the data file and PARAMS.h        IF (Integr_GeoPot.EQ.0) THEN
180          atm_cp=1004. _d 0  C  --  Energy Conserving Form, No hFac  --
181          atm_kappa=2. _d 0/7. _d 0  C------------ The integration for the first level phi(k=1) is the same
182          atm_po=1. _d 5  C             for both the "finite volume" and energy conserving methods.
183    C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
184    C             condition is simply Phi'(Ro_surf)=0.
185    C           o convention ddPI > 0 (same as drF & drC)
186    C-----------------------------------------------------------------------
187          IF (K.EQ.1) THEN          IF (K.EQ.1) THEN
188            ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa)            ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
189       &                  -((rF(K)/atm_po)**atm_kappa) )       &                  -((rC(K)/atm_po)**atm_kappa) )
190            DO j=jMin,jMax            DO j=jMin,jMax
191              DO i=iMin,iMax             DO i=iMin,iMax
192                ddRp=ddRp1               phiHyd(i,j,K)=
193                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.       &          ddPIp*maskC(i,j,K,bi,bj)
194  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-----------------------------------------------------------------------  
195             ENDDO             ENDDO
196            ENDDO            ENDDO
197          ELSE          ELSE
198    C-------- This discretization is the energy conserving form
199  C-------- This discretization is the "finite volume" form which            ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
200  C         integrates the hydrostatic equation of each half/sub-layer.       &                 -((rC( K )/atm_po)**atm_kappa) )*0.5
201  C         This seems most natural and could easily allow for lopped cells            DO j=jMin,jMax
202  C         by replacing rF(K) with the height of the surface (not implemented).             DO i=iMin,iMax
203  C         in the lower layers (e.g. at k=1).                phiHyd(i,j,K)=phiHyd(i,j,K-1)
204  C       &           +ddPI*maskC(i,j,K-1,bi,bj)
205  c         ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa)       &                *(theta(I,J,K-1,bi,bj)-tRef(K-1))
206  c    &                  -((rC(K-1)/atm_po)**atm_kappa) )       &           +ddPI*maskC(i,j, K ,bi,bj)
207  c         ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)       &                *(theta(I,J, K ,bi,bj)-tRef( K ))
208  c    &                  -((rF( K )/atm_po)**atm_kappa) )  C             Old code (atmos-exact) looked like this
209    Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
210    Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))
211               ENDDO
212              ENDDO
213            ENDIF
214    C end: Energy Conserving Form, No hFac  --
215  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
216    
217          ELSEIF (Integr_GeoPot.EQ.1) THEN
218    C  --  Finite Volume Form, with hFac, linear in P by Half level  --
219    C---------
220    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
221    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
222    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
223    C   also: if Interface_W at the middle between tracer levels, this form
224    C     is close to the Energy Cons. form in the Interior, except for the
225    C     non-linearity in PI(p)
226    C---------
227            IF (K.EQ.1) THEN
228              ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
229         &                  -((rC(K)/atm_po)**atm_kappa) )
230              DO j=jMin,jMax
231               DO i=iMin,iMax
232                 phiHyd(i,j,K) =
233         &          ddPIp*hFacC(I,J, K ,bi,bj)
234         &               *(theta(I,J, K ,bi,bj)-tRef( K ))
235               ENDDO
236              ENDDO
237            ELSE
238              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
239         &                  -((rF( K )/atm_po)**atm_kappa) )
240              ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
241         &                  -((rC( K )/atm_po)**atm_kappa) )
242              DO j=jMin,jMax
243               DO i=iMin,iMax
244                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
245         &         +ddPIm*hFacC(I,J,K-1,bi,bj)
246         &               *(theta(I,J,K-1,bi,bj)-tRef(K-1))
247         &         +ddPIp*hFacC(I,J, K ,bi,bj)
248         &               *(theta(I,J, K ,bi,bj)-tRef( K ))
249               ENDDO
250              ENDDO
251            ENDIF
252    C end: Finite Volume Form, with hFac, linear in P by Half level  --
253    C-----------------------------------------------------------------------
254    
255  C-------- This discretization is the energy conserving form        ELSEIF (Integr_GeoPot.EQ.2) THEN
256            ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)  C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
257       &                  -((rC(K-1)/atm_po)**atm_kappa) )*0.5  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
258            ddRm1=ddRp1  C  Finite Difference formulation consistent with Partial Cell,
259    C    case Tracer level at the middle of InterFace_W
260    C    linear between 2 Tracer levels ; conserve energy in the Interior
261    C---------
262            Kp1 = min(Nr,K+1)
263            IF (K.EQ.1) THEN
264              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
265         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
266              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
267         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
268              DO j=jMin,jMax
269               DO i=iMin,iMax
270                 phiHyd(i,j,K) =
271         &        ( ddPIm*max(zero, hFacC(i,j,K,bi,bj)-half)
272         &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)-half) )
273         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
274         &               * maskC(i,j, K ,bi,bj)
275               ENDDO
276              ENDDO
277            ELSE
278              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
279         &                  -((rC( K )/atm_po)**atm_kappa) )
280              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
281         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
282              DO j=jMin,jMax
283               DO i=iMin,iMax
284                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
285         &        + ddPIm*0.5
286         &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))
287         &               * maskC(i,j,K-1,bi,bj)
288         &        +(ddPIm*max(zero, hFacC(i,j,K,bi,bj)-half)
289         &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)-half) )
290         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
291         &               * maskC(i,j, K ,bi,bj)
292               ENDDO
293              ENDDO
294            ENDIF
295    C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
296  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
297    
298          ELSEIF (Integr_GeoPot.EQ.3) THEN
299    C  --  Finite Difference Form, with hFac, Interface_W = middle  --
300    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
301    C  Finite Difference formulation consistent with Partial Cell,
302    C   Valid & accurate if Interface_W at middle between tracer levels
303    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
304    C---------
305            Kp1 = min(Nr,K+1)
306            IF (K.EQ.1) THEN
307              ratioRm=0.5*drF(K)/(rF(k)-rC(K))
308              ratioRp=drF(K)*recip_drC(Kp1)
309              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
310         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
311              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
312         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
313            DO j=jMin,jMax            DO j=jMin,jMax
314              DO i=iMin,iMax             DO i=iMin,iMax
315                ddRp=ddRp1               phiHyd(i,j,K) =
316                ddRm=ddRm1       &        ( ddPIm*max(zero,(hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
317                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.       &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)*ratioRp     -half) )
318                IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0.       &               *(theta(i,j, K ,bi,bj)-tRef( K ))
319                phiHyd(i,j,K)=phiHyd(i,j,K-1)       &               * maskC(i,j, K ,bi,bj)
320       &           -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1))             ENDDO
321       &             +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) )            ENDDO
322  C             Old code (atmos-exact) looked like this          ELSE
323  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddRm1*            ratioRm=drF(K)*recip_drC(K)
324  Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))            ratioRp=drF(K)*recip_drC(Kp1)
325              ENDDO            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
326         &                  -((rC( K )/atm_po)**atm_kappa) )
327              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
328         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
329              DO j=jMin,jMax
330               DO i=iMin,iMax
331                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
332         &        + ddPIm*0.5
333         &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))
334         &               * maskC(i,j,K-1,bi,bj)
335         &        +(ddPIm*max(zero,(hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
336         &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)*ratioRp     -half) )
337         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
338         &               * maskC(i,j, K ,bi,bj)
339               ENDDO
340            ENDDO            ENDDO
341          ENDIF          ENDIF
342    C end: Finite Difference Form, with hFac, Interface_W = middle  --
343    C-----------------------------------------------------------------------
344    
345          ELSE
346            STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !'
347          ENDIF
348    
349    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
350        ELSE        ELSE
351          STOP 'CALC_PHI_HYD: We should never reach this point!'          STOP 'CALC_PHI_HYD: We should never reach this point!'
352        ENDIF        ENDIF
353    
354  #endif  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
355    
356        RETURN        RETURN
357        END        END
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22