/[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.4 by cnh, Wed Sep 9 15:04:44 1998 UTC revision 1.18 by mlosch, Wed Jul 31 16:38:30 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                         theta, salt,
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     |   theta,salt    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    
54    C     !INPUT/OUTPUT PARAMETERS:
55  C     == Routine arguments ==  C     == Routine arguments ==
56        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax,K
57        _RL buoyKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
58        _RL buoyKP1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
59        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
60        integer myThid        INTEGER myThid
61          
62    #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
63    
64    C     !LOCAL VARIABLES:
65  C     == Local variables ==  C     == Local variables ==
66        INTEGER i,j,Km1        INTEGER i,j, Kp1
67        _RL halfLayer        _RL zero, one, half
68          _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69          _RL dRloc,dRlocKp1
70          _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
71    CEOP
72    
73          zero = 0. _d 0
74          one  = 1. _d 0
75          half = .5 _d 0
76    
77    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
78    C  Atmosphere:  
79    C Integr_GeoPot => select one option for the integration of the Geopotential:
80    C   = 0 : Energy Conserving Form, No hFac ;
81    C   = 1 : Finite Volume Form, with hFac, linear in P by Half level;
82    C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels
83    C     2 : case Tracer level at the middle of InterFace_W;
84    C     3 : case InterFace_W  at the middle of Tracer levels;
85    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
86    
87    #ifdef ALLOW_AUTODIFF_TAMC
88              act1 = bi - myBxLo(myThid)
89              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
90    
91              act2 = bj - myByLo(myThid)
92              max2 = myByHi(myThid) - myByLo(myThid) + 1
93    
94              act3 = myThid - 1
95              max3 = nTx*nTy
96    
97              act4 = ikey_dynamics - 1
98    
99              ikey = (act1 + 1) + act2*max1
100         &                      + act3*max1*max2
101         &                      + act4*max1*max2*max3
102    #endif /* ALLOW_AUTODIFF_TAMC */
103    
104          IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
105    C       This is the hydrostatic pressure calculation for the Ocean
106    C       which uses the FIND_RHO() routine to calculate density
107    C       before integrating g*rho over the current layer/interface
108    
109            dRloc=drC(k)
110            IF (k.EQ.1) dRloc=drF(1)
111            IF (k.EQ.Nr) THEN
112              dRlocKp1=0.
113            ELSE
114              dRlocKp1=drC(k+1)
115            ENDIF
116    
117    C--     If this is the top layer we impose the boundary condition
118    C       P(z=eta) = P(atmospheric_loading)
119            IF (k.EQ.1) THEN
120              DO j=jMin,jMax
121                DO i=iMin,iMax
122    #ifdef ATMOSPHERIC_LOADING
123                  phiHyd(i,j,k)=pload(i,j,bi,bj)*recip_rhoConst
124    #else
125                  phiHyd(i,j,k)=0. _d 0
126    #endif
127                ENDDO
128              ENDDO
129            ENDIF
130    
131    C       Calculate density
132    #ifdef ALLOW_AUTODIFF_TAMC
133                kkey = (ikey-1)*Nr + k
134    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
135    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
136    #endif /* ALLOW_AUTODIFF_TAMC */
137            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
138         &                 theta, salt,
139         &                 alphaRho, myThid)
140    
141    C       Hydrostatic pressure at cell centers
142            DO j=jMin,jMax
143              DO i=iMin,iMax
144    #ifdef      ALLOW_AUTODIFF_TAMC
145    c           Patrick, is this directive correct or even necessary in
146    c           this new code?
147    c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...
148    c           within the k-loop.
149    CADJ GENERAL
150    #endif      /* ALLOW_AUTODIFF_TAMC */
151    
152    C---------- This discretization is the "finite volume" form
153    C           which has not been used to date since it does not
154    C           conserve KE+PE exactly even though it is more natural
155    C
156    c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
157    c    &              drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
158    c           phiHyd(i,j,k)=phiHyd(i,j,k)+
159    c    &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
160    C-----------------------------------------------------------------------
161    
162    C---------- This discretization is the "energy conserving" form
163    C           which has been used since at least Adcroft et al., MWR 1997
164    C
165                phiHyd(i,j,k)=phiHyd(i,j,k)+
166         &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
167                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
168         &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
169    C-----------------------------------------------------------------------
170              ENDDO
171            ENDDO
172            
173    
174    
175          ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
176    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
177    C       This is the hydrostatic geopotential calculation for the Atmosphere
178    C       The ideal gas law is used implicitly here rather than calculating
179    C       the specific volume, analogous to the oceanic case.
180    
181    C       Integrate d Phi / d pi
182    
183          IF (Integr_GeoPot.EQ.0) THEN
184    C  --  Energy Conserving Form, No hFac  --
185    C------------ The integration for the first level phi(k=1) is the same
186    C             for both the "finite volume" and energy conserving methods.
187    Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
188    C             condition is simply Phi-prime(Ro_surf)=0.
189    C           o convention ddPI > 0 (same as drF & drC)
190    C-----------------------------------------------------------------------
191            IF (K.EQ.1) THEN
192              ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
193         &                  -((rC(K)/atm_po)**atm_kappa) )
194              DO j=jMin,jMax
195               DO i=iMin,iMax
196                 phiHyd(i,j,K)=
197         &          ddPIp*maskC(i,j,K,bi,bj)
198         &               *(theta(I,J,K,bi,bj)-tRef(K))
199               ENDDO
200              ENDDO
201            ELSE
202    C-------- This discretization is the energy conserving form
203              ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
204         &                 -((rC( K )/atm_po)**atm_kappa) )*0.5
205              DO j=jMin,jMax
206               DO i=iMin,iMax
207                  phiHyd(i,j,K)=phiHyd(i,j,K-1)
208         &           +ddPI*maskC(i,j,K-1,bi,bj)
209         &                *(theta(I,J,K-1,bi,bj)-tRef(K-1))
210         &           +ddPI*maskC(i,j, K ,bi,bj)
211         &                *(theta(I,J, K ,bi,bj)-tRef( K ))
212    C             Old code (atmos-exact) looked like this
213    Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
214    Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))
215               ENDDO
216              ENDDO
217            ENDIF
218    C end: Energy Conserving Form, No hFac  --
219    C-----------------------------------------------------------------------
220    
221          ELSEIF (Integr_GeoPot.EQ.1) THEN
222    C  --  Finite Volume Form, with hFac, linear in P by Half level  --
223    C---------
224    C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
225    C   Note: a true Finite Volume form should be linear between 2 Interf_W :
226    C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
227    C   also: if Interface_W at the middle between tracer levels, this form
228    C     is close to the Energy Cons. form in the Interior, except for the
229    C     non-linearity in PI(p)
230    C---------
231            IF (K.EQ.1) THEN
232              ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)
233         &                  -((rC(K)/atm_po)**atm_kappa) )
234              DO j=jMin,jMax
235               DO i=iMin,iMax
236                 phiHyd(i,j,K) =
237         &          ddPIp*_hFacC(I,J, K ,bi,bj)
238         &               *(theta(I,J, K ,bi,bj)-tRef( K ))
239               ENDDO
240              ENDDO
241            ELSE
242              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
243         &                  -((rF( K )/atm_po)**atm_kappa) )
244              ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
245         &                  -((rC( K )/atm_po)**atm_kappa) )
246              DO j=jMin,jMax
247               DO i=iMin,iMax
248                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
249         &         +ddPIm*_hFacC(I,J,K-1,bi,bj)
250         &               *(theta(I,J,K-1,bi,bj)-tRef(K-1))
251         &         +ddPIp*_hFacC(I,J, K ,bi,bj)
252         &               *(theta(I,J, K ,bi,bj)-tRef( K ))
253               ENDDO
254              ENDDO
255            ENDIF
256    C end: Finite Volume Form, with hFac, linear in P by Half level  --
257    C-----------------------------------------------------------------------
258    
259          ELSEIF (Integr_GeoPot.EQ.2) THEN
260    C  --  Finite Difference Form, with hFac, Tracer Lev. = middle  --
261    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
262    C  Finite Difference formulation consistent with Partial Cell,
263    C    case Tracer level at the middle of InterFace_W
264    C    linear between 2 Tracer levels ; conserve energy in the Interior
265    C---------
266            Kp1 = min(Nr,K+1)
267            IF (K.EQ.1) THEN
268              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
269         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
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) =
275         &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
276         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
277         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
278         &               * maskC(i,j, K ,bi,bj)
279               ENDDO
280              ENDDO
281            ELSE
282              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
283         &                  -((rC( K )/atm_po)**atm_kappa) )
284              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
285         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
286              DO j=jMin,jMax
287               DO i=iMin,iMax
288                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
289         &        + ddPIm*0.5
290         &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))
291         &               * maskC(i,j,K-1,bi,bj)
292         &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
293         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
294         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
295         &               * maskC(i,j, K ,bi,bj)
296               ENDDO
297              ENDDO
298            ENDIF
299    C end: Finite Difference Form, with hFac, Tracer Lev. = middle  --
300    C-----------------------------------------------------------------------
301    
302          ELSEIF (Integr_GeoPot.EQ.3) THEN
303    C  --  Finite Difference Form, with hFac, Interface_W = middle  --
304    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
305    C  Finite Difference formulation consistent with Partial Cell,
306    C   Valid & accurate if Interface_W at middle between tracer levels
307    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
308    C---------
309            Kp1 = min(Nr,K+1)
310            IF (K.EQ.1) THEN
311              ratioRm=0.5*drF(K)/(rF(k)-rC(K))
312              ratioRp=drF(K)*recip_drC(Kp1)
313              ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
314         &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0
315              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
316         &                  -((rC(Kp1)/atm_po)**atm_kappa) )  
317              DO j=jMin,jMax
318               DO i=iMin,iMax
319                 phiHyd(i,j,K) =
320         &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
321         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
322         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
323         &               * maskC(i,j, K ,bi,bj)
324               ENDDO
325              ENDDO
326            ELSE
327              ratioRm=drF(K)*recip_drC(K)
328              ratioRp=drF(K)*recip_drC(Kp1)
329              ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)
330         &                  -((rC( K )/atm_po)**atm_kappa) )
331              ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
332         &                  -((rC(Kp1)/atm_po)**atm_kappa) )
333              DO j=jMin,jMax
334               DO i=iMin,iMax
335                 phiHyd(i,j,K) = phiHyd(i,j,K-1)
336         &        + ddPIm*0.5
337         &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))
338         &               * maskC(i,j,K-1,bi,bj)
339         &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
340         &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
341         &               *(theta(i,j, K ,bi,bj)-tRef( K ))
342         &               * maskC(i,j, K ,bi,bj)
343               ENDDO
344              ENDDO
345            ENDIF
346    C end: Finite Difference Form, with hFac, Interface_W = middle  --
347    C-----------------------------------------------------------------------
348    
349          ELSE
350            STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !'
351          ENDIF
352    
353    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
354          ELSE
355            STOP 'CALC_PHI_HYD: We should never reach this point!'
356          ENDIF
357    
358    #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
359    
360        if (K.eq.1) then        RETURN
361         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) )*recip_HoriVertRatio  
      &  *0.5 _d 0*( buoyKM1(i,j)+buoyKP1(i,j) )  
 C    &  *rkFac  
        ENDDO  
       ENDDO  
   
 ! ------------------------------------------------------------------------------  
       return  
       end  
 ! ==============================================================================  
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22