/[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.5 by jmc, Wed Jan 31 17:11:16 2001 UTC revision 1.44 by jmc, Fri Apr 4 20:54:11 2014 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_AUTODIFF
7    # include "AUTODIFF_OPTIONS.h"
8    #endif
9    
10    CBOP
11    C     !ROUTINE: CALC_PHI_HYD
12    C     !INTERFACE:
13        SUBROUTINE CALC_PHI_HYD(        SUBROUTINE CALC_PHI_HYD(
14       I                         bi, bj, iMin, iMax, jMin, jMax, K,       I                         bi, bj, iMin, iMax, jMin, jMax, k,
15       I                         theta, salt,       I                         tFld, sFld,
16       U                         phiHyd,       U                         phiHydF,
17       I                         myThid)       O                         phiHydC, dPhiHydX, dPhiHydY,
18  C     /==========================================================\       I                         myTime, myIter, myThid )
19    C     !DESCRIPTION: \bv
20    C     *==========================================================*
21  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
22  C     | o Integrate the hydrostatic relation to find phiHyd.     |  C     | o Integrate the hydrostatic relation to find the Hydros. |
23  C     |                                                          |  C     *==========================================================*
24  C     | On entry:                                                |  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)
25  C     |   theta,salt    are the current thermodynamics quantities|  C     | On entry:
26  C     |                 (unchanged on exit)                      |  C     |   tFld,sFld     are the current thermodynamics quantities
27  C     |   phiHyd(i,j,1:k-1) is the hydrostatic pressure/geopot.  |  C     |                 (unchanged on exit)
28  C     |                 at cell centers (tracer points)          |  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
29  C     |                 - 1:k-1 layers are valid                 |  C     |                at middle between tracer points k-1,k
30  C     |                 - k:Nr layers are invalid                |  C     | On exit:
31  C     |   phiHyd(i,j,k) is the hydrostatic pressure/geop.        |  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
32  C     |                 at cell the interface k (w point above)  |  C     |                at cell centers (tracer points), level k
33  C     | On exit:                                                 |  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
34  C     |   phiHyd(i,j,1:k) is the hydrostatic pressure/geopot.    |  C     |                at middle between tracer points k,k+1
35  C     |                 at cell centers (tracer points)          |  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
36  C     |                 - 1:k layers are valid                   |  C     |                at cell centers (tracer points), level k
37  C     |                 - k+1:Nr layers are invalid              |  C     | integr_GeoPot allows to select one integration method
38  C     |   phiHyd(i,j,k+1) is the hydrostatic pressure/geop.      |  C     |    1= Finite volume form ; else= Finite difference form
39  C     |                 at cell the interface k+1 (w point below)|  C     *==========================================================*
40  C     |                                                          |  C     \ev
41  C     \==========================================================/  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
49    #include "tamc.h"
50    #include "tamc_keys.h"
51    #endif /* ALLOW_AUTODIFF */
52    #include "SURFACE.h"
53    #include "DYNVARS.h"
54    
55    C     !INPUT/OUTPUT PARAMETERS:
56  C     == Routine arguments ==  C     == Routine arguments ==
57        INTEGER bi,bj,iMin,iMax,jMin,jMax,K  C     bi, bj, k  :: tile and level indices
58        _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  C     iMin,iMax,jMin,jMax :: computational domain
59        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  C     tFld       :: potential temperature
60        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     sFld       :: salinity
61        INTEGER myThid  C     phiHydF    :: hydrostatic potential anomaly at middle between
62    C                   2 centers (entry: Interf_k ; output: Interf_k+1)
63    C     phiHydC    :: hydrostatic potential anomaly at cell center
64    C     dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
65    C     myTime     :: current time
66    C     myIter     :: current iteration number
67    C     myThid     :: thread number for this instance of the routine.
68          INTEGER bi,bj,iMin,iMax,jMin,jMax,k
69          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
70          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
71          _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72          _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73          _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74          _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75          _RL myTime
76          INTEGER myIter, myThid
77    
78  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
79    
80    C     !LOCAL VARIABLES:
81  C     == Local variables ==  C     == Local variables ==
82    C     phiHydU    :: hydrostatic potential anomaly at interface k+1 (Upper k)
83    C     pKappaF    :: (p/Po)^kappa at interface k
84    C     pKappaU    :: (p/Po)^kappa at interface k+1 (Upper k)
85    C     pKappaC    :: (p/Po)^kappa at cell center k
86        INTEGER i,j        INTEGER i,j
87        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _RL dRloc,dRlocKp1  #ifndef DISABLE_SIGMA_CODE
89        _RL ddRm1, ddRp1, ddRm, ddRp        _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL atm_cp, atm_kappa, atm_po        _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL pKappaC, locDepth, fullDepth
93    #endif /* DISABLE_SIGMA_CODE */
94          _RL dRlocM,dRlocP, ddRloc, locAlpha
95          _RL ddPIm, ddPIp, rec_dRm, rec_dRp
96          _RL surfPhiFac
97          LOGICAL useDiagPhiRlow, addSurfPhiAnom
98          LOGICAL useFVgradPhi
99    CEOP
100          useDiagPhiRlow = .TRUE.
101          addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
102          useFVgradPhi   = selectSigmaCoord.NE.0
103    
104          surfPhiFac = 0.
105          IF (addSurfPhiAnom) surfPhiFac = 1.
106    
107    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
108    C  Atmosphere:
109    C integr_GeoPot => select one option for the integration of the Geopotential:
110    C   = 0 : Energy Conserving Form, accurate with Topo full cell;
111    C   = 1 : Finite Volume Form, with Part-Cell, linear in P by Half level;
112    C   =2,3: Finite Difference Form, with Part-Cell,
113    C         linear in P between 2 Tracer levels.
114    C       can handle both cases: Tracer lev at the middle of InterFace_W
115    C                          and InterFace_W at the middle of Tracer lev;
116    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
117    
118    #ifdef ALLOW_AUTODIFF_TAMC
119              act1 = bi - myBxLo(myThid)
120              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
121    
122              act2 = bj - myByLo(myThid)
123              max2 = myByHi(myThid) - myByLo(myThid) + 1
124    
125              act3 = myThid - 1
126              max3 = nTx*nTy
127    
128              act4 = ikey_dynamics - 1
129    
130              ikey = (act1 + 1) + act2*max1
131         &                      + act3*max1*max2
132         &                      + act4*max1*max2*max3
133    #endif /* ALLOW_AUTODIFF_TAMC */
134    
135    C--   Initialize phiHydF to zero :
136    C     note: atmospheric_loading or Phi_topo anomaly are incorporated
137    C           later in S/R calc_grad_phi_hyd
138          IF (k.EQ.1) THEN
139            DO j=1-OLy,sNy+OLy
140             DO i=1-OLx,sNx+OLx
141               phiHydF(i,j) = 0.
142             ENDDO
143            ENDDO
144          ENDIF
145    
146        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
147          IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
148  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
149  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
150  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
151    #ifdef ALLOW_AUTODIFF_TAMC
152    CADJ GENERAL
153    #endif /* ALLOW_AUTODIFF_TAMC */
154    
155          dRloc=drC(k)          IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
156          IF (k.EQ.1) dRloc=drF(1)  C---    Calculate density
157          IF (k.EQ.Nr) THEN  #ifdef ALLOW_AUTODIFF_TAMC
158            dRlocKp1=0.            kkey = (ikey-1)*Nr + k
159    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
160    CADJ &     kind = isbyte
161    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
162    CADJ &     kind = isbyte
163    #endif /* ALLOW_AUTODIFF_TAMC */
164              CALL FIND_RHO_2D(
165         I              iMin, iMax, jMin, jMax, k,
166         I              tFld(1-OLx,1-OLy,k,bi,bj),
167         I              sFld(1-OLx,1-OLy,k,bi,bj),
168         O              alphaRho,
169         I              k, bi, bj, myThid )
170          ELSE          ELSE
171            dRlocKp1=drC(k+1)            DO j=jMin,jMax
172               DO i=iMin,iMax
173                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
174               ENDDO
175              ENDDO
176            ENDIF
177    
178    #ifdef ALLOW_SHELFICE
179    C     mask rho, so that there is no contribution of phiHyd from
180    C     overlying shelfice (whose density we do not know)
181            IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
182    C- note: does not work for down_slope pkg which needs rho below the bottom.
183    C    setting rho=0 above the ice-shelf base is enough (and works in both cases)
184    C    but might be slower (--> keep original masking if not using down_slope pkg)
185             DO j=jMin,jMax
186              DO i=iMin,iMax
187               IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
188              ENDDO
189             ENDDO
190            ELSEIF ( useShelfIce ) THEN
191             DO j=jMin,jMax
192              DO i=iMin,iMax
193               alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
194              ENDDO
195             ENDDO
196          ENDIF          ENDIF
197    #endif /* ALLOW_SHELFICE */
198    
199  C--     If this is the top layer we impose the boundary condition  #ifdef ALLOW_MOM_COMMON
200  C       P(z=eta) = P(atmospheric_loading)  C--     Quasi-hydrostatic terms are added in as if they modify the buoyancy
201          IF (k.EQ.1) THEN          IF (quasiHydrostatic) THEN
202             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
203            ENDIF
204    #endif /* ALLOW_MOM_COMMON */
205    
206    #ifdef NONLIN_FRSURF
207            IF ( addSurfPhiAnom .AND.
208         &       uniformFreeSurfLev .AND. k.EQ.1 ) THEN
209            DO j=jMin,jMax            DO j=jMin,jMax
210              DO i=iMin,iMax              DO i=iMin,iMax
211  C             *NOTE* The loading should go here but has not been implemented yet                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
212                phiHyd(i,j,k)=0.       &                      *gravity*alphaRho(i,j)*recip_rhoConst
213              ENDDO              ENDDO
214            ENDDO            ENDDO
215          ENDIF          ENDIF
216    #endif /* NONLIN_FRSURF */
217    
218  C       Calculate density  C----  Hydrostatic pressure at cell centers
         CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      &                 theta, salt,  
      &                 alphaRho, myThid)  
219    
220  C       Hydrostatic pressure at cell centers         IF (integr_GeoPot.EQ.1) THEN
221          DO j=jMin,jMax  C  --  Finite Volume Form
222    
223    C---------- This discretization is the "finite volume" form
224    C           which has not been used to date since it does not
225    C           conserve KE+PE exactly even though it is more natural
226    
227            IF ( uniformFreeSurfLev ) THEN
228             DO j=jMin,jMax
229            DO i=iMin,iMax            DO i=iMin,iMax
230                phiHydC(i,j) = phiHydF(i,j)
231         &          + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
232                phiHydF(i,j) = phiHydF(i,j)
233         &                 + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
234              ENDDO
235             ENDDO
236            ELSE
237             DO j=jMin,jMax
238              DO i=iMin,iMax
239               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
240                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
241    #ifdef NONLIN_FRSURF
242                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
243    #endif
244                phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst
245               ELSE
246                phiHydC(i,j) = phiHydF(i,j)
247         &            + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
248               ENDIF
249                phiHydF(i,j) = phiHydC(i,j)
250         &            + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
251              ENDDO
252             ENDDO
253            ENDIF
254    
255           ELSE
256    C  --  Finite Difference Form
257    
258    C---------- This discretization is the "energy conserving" form
259    C           which has been used since at least Adcroft et al., MWR 1997
260    
261            dRlocM = halfRL*drC(k)
262            IF (k.EQ.1) dRlocM=rF(k)-rC(k)
263            IF (k.EQ.Nr) THEN
264              dRlocP=rC(k)-rF(k+1)
265            ELSE
266              dRlocP=halfRL*drC(k+1)
267            ENDIF
268            IF ( uniformFreeSurfLev ) THEN
269             DO j=jMin,jMax
270              DO i=iMin,iMax
271                phiHydC(i,j) = phiHydF(i,j)
272         &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
273                phiHydF(i,j) = phiHydC(i,j)
274         &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
275              ENDDO
276             ENDDO
277            ELSE
278             rec_dRm = oneRL/(rF(k)-rC(k))
279             rec_dRp = oneRL/(rC(k)-rF(k+1))
280             DO j=jMin,jMax
281              DO i=iMin,iMax
282               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
283                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
284    #ifdef NONLIN_FRSURF
285                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
286    #endif
287                phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
288         &                     +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
289         &                    )*gravity*alphaRho(i,j)*recip_rhoConst
290               ELSE
291                phiHydC(i,j) = phiHydF(i,j)
292         &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
293               ENDIF
294                phiHydF(i,j) = phiHydC(i,j)
295         &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
296              ENDDO
297             ENDDO
298            ENDIF
299    
300    C  --  end if integr_GeoPot = ...
301           ENDIF
302    
303    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
304          ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
305    C       This is the hydrostatic pressure calculation for the Ocean
306    C       which uses the FIND_RHO() routine to calculate density before
307    C       integrating (1/rho)_prime*dp over the current layer/interface
308  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
             Is this directive correct or even necessary in this new code?  
309  CADJ GENERAL  CADJ GENERAL
310  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
311    
312            IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
313    C--     Calculate density
314    #ifdef ALLOW_AUTODIFF_TAMC
315              kkey = (ikey-1)*Nr + k
316    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
317    CADJ &     kind = isbyte
318    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
319    CADJ &     kind = isbyte
320    #endif /* ALLOW_AUTODIFF_TAMC */
321              CALL FIND_RHO_2D(
322         I              iMin, iMax, jMin, jMax, k,
323         I              tFld(1-OLx,1-OLy,k,bi,bj),
324         I              sFld(1-OLx,1-OLy,k,bi,bj),
325         O              alphaRho,
326         I              k, bi, bj, myThid )
327    #ifdef ALLOW_AUTODIFF_TAMC
328    CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte,
329    CADJ &     kind = isbyte
330    #endif /* ALLOW_AUTODIFF_TAMC */
331            ELSE
332              DO j=jMin,jMax
333               DO i=iMin,iMax
334                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
335               ENDDO
336              ENDDO
337            ENDIF
338    
339    C--     Calculate specific volume anomaly : alpha_prime = 1/rho - alpha_Cst
340            DO j=jMin,jMax
341              DO i=iMin,iMax
342                locAlpha=alphaRho(i,j)+rhoConst
343                alphaRho(i,j)=maskC(i,j,k,bi,bj)*
344         &              (oneRL/locAlpha - recip_rhoConst)
345              ENDDO
346            ENDDO
347    
348    #ifdef ALLOW_MOM_COMMON
349    C--     Quasi-hydrostatic terms are added as if they modify the specific-volume
350            IF (quasiHydrostatic) THEN
351             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
352            ENDIF
353    #endif /* ALLOW_MOM_COMMON */
354    
355    C----  Hydrostatic pressure at cell centers
356    
357           IF (integr_GeoPot.EQ.1) THEN
358    C  --  Finite Volume Form
359    
360             DO j=jMin,jMax
361              DO i=iMin,iMax
362    
363  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
364  C           which has not been used to date since it does not  C           which has not been used to date since it does not
365  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
 C  
 c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  
 c    &              drF(K)*gravity*alphaRho(i,j)  
 c           phiHyd(i,j,k)=phiHyd(i,j,k)+  
 c    &          0.5*drF(K)*gravity*alphaRho(i,j)  
 C-----------------------------------------------------------------------  
366    
367  C---------- This discretization is the "energy conserving" form             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
368  C           which has been used since at least Adcroft et al., MWR 1997               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
369  C  #ifdef NONLIN_FRSURF
370              phiHyd(i,j,k)=phiHyd(i,j,k)+               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
371       &          0.5*dRloc*gravity*alphaRho(i,j)  #endif
372              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+               phiHydC(i,j) = ddRloc*alphaRho(i,j)
373       &          0.5*dRlocKp1*gravity*alphaRho(i,j)  c--to reproduce results of c48d_post: uncomment those 4+1 lines
374  C-----------------------------------------------------------------------  c            phiHydC(i,j)=phiHydF(i,j)
375    c    &          +(hFacC(i,j,k,bi,bj)-halfRL)*drF(k)*alphaRho(i,j)
376    c            phiHydF(i,j)=phiHydF(i,j)
377    c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
378               ELSE
379                 phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
380    c            phiHydF(i,j) = phiHydF(i,j) +        drF(k)*alphaRho(i,j)
381               ENDIF
382    c-- and comment this last one:
383                 phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
384    c-----
385            ENDDO            ENDDO
386          ENDDO           ENDDO
387            
388           ELSE
389    C  --  Finite Difference Form, with Part-Cell Bathy
390    
391             dRlocM = halfRL*drC(k)
392             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
393             IF (k.EQ.Nr) THEN
394               dRlocP=rC(k)-rF(k+1)
395             ELSE
396               dRlocP=halfRL*drC(k+1)
397             ENDIF
398             rec_dRm = oneRL/(rF(k)-rC(k))
399             rec_dRp = oneRL/(rC(k)-rF(k+1))
400    
401             DO j=jMin,jMax
402              DO i=iMin,iMax
403    
404    C---------- This discretization is the "energy conserving" form
405    
406               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
407                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
408    #ifdef NONLIN_FRSURF
409                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
410    #endif
411                 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
412         &                      +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
413         &                     )*alphaRho(i,j)
414               ELSE
415                 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
416               ENDIF
417                 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
418              ENDDO
419             ENDDO
420    
421    C  --  end if integr_GeoPot = ...
422           ENDIF
423    
424        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
425    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
426  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
427  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
428  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
429    
430  C       Integrate d Phi / d pi  C--     virtual potential temperature anomaly (including water vapour effect)
431            DO j=jMin,jMax
432             DO i=iMin,iMax
433              alphaRho(i,j) = ( tFld(i,j,k,bi,bj)
434         &                      *( sFld(i,j,k,bi,bj)*atm_Rq + oneRL )
435         &                    - tRef(k) )*maskC(i,j,k,bi,bj)
436             ENDDO
437            ENDDO
438    
439  C *NOTE* These constants should be in the data file and PARAMS.h  #ifdef ALLOW_MOM_COMMON
440          atm_cp=1004. _d 0  C--     Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
441          atm_kappa=2. _d 0/7. _d 0          IF (quasiHydrostatic) THEN
442          atm_po=1. _d 5           CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
443          IF (K.EQ.1) THEN          ENDIF
444            ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa)  #endif /* ALLOW_MOM_COMMON */
      &                  -((rF(K)/atm_po)**atm_kappa) )  
           DO j=jMin,jMax  
             DO i=iMin,iMax  
               ddRp=ddRp1  
               IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.  
 C------------ The integration for the first level phi(k=1) is the  
 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-----------------------------------------------------------------------  
            ENDDO  
           ENDDO  
         ELSE  
445    
446  C-------- This discretization is the "finite volume" form which  C---    Integrate d Phi / d pi
 C         integrates the hydrostatic equation of each half/sub-layer.  
 C         This seems most natural and could easily allow for lopped cells  
 C         by replacing rF(K) with the height of the surface (not implemented).  
 C         in the lower layers (e.g. at k=1).  
 C  
 c         ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa)  
 c    &                  -((rC(K-1)/atm_po)**atm_kappa) )  
 c         ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)  
 c    &                  -((rF( K )/atm_po)**atm_kappa) )  
 C-----------------------------------------------------------------------  
447    
448    #ifndef DISABLE_SIGMA_CODE
449    C  --  Finite Volume Form, integrated to r-level (cell center & upper interface)
450           IF ( useFVgradPhi ) THEN
451    
452             fullDepth = rF(1)-rF(Nr+1)
453             DO j=jMin,jMax
454              DO i=iMin,iMax
455               locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
456    #ifdef NONLIN_FRSURF
457               locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
458    #endif
459               pKappaF(i,j) = (
460         &           ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
461         &                              + bHybSigmF( k )*locDepth
462         &           )/atm_Po )**atm_kappa
463               pKappaC      = (
464         &           ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
465         &                              + bHybSigmC( k )*locDepth
466         &           )/atm_Po )**atm_kappa
467               pKappaU(i,j) = (
468         &           ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
469         &                             + bHybSigmF(k+1)*locDepth
470         &           )/atm_Po )**atm_kappa
471               phiHydC(i,j) = phiHydF(i,j)
472         &        + atm_Cp*( pKappaF(i,j) - pKappaC      )*alphaRho(i,j)
473               phiHydU(i,j) = phiHydF(i,j)
474         &        + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
475              ENDDO
476             ENDDO
477    C end: Finite Volume Form, integrated to r-level.
478    
479           ELSEIF (integr_GeoPot.EQ.0) THEN
480    #else /* DISABLE_SIGMA_CODE */
481           IF (integr_GeoPot.EQ.0) THEN
482    #endif /* DISABLE_SIGMA_CODE */
483    C  --  Energy Conserving Form, accurate with Full cell topo  --
484    C------------ The integration for the first level phi(k=1) is the same
485    C             for both the "finite volume" and energy conserving methods.
486    C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
487    C             condition is simply Phi-prime(Ro_surf)=0.
488    C           o convention ddPI > 0 (same as drF & drC)
489    C-----------------------------------------------------------------------
490             IF (k.EQ.1) THEN
491               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
492         &                   -((rC( k )/atm_Po)**atm_kappa) )
493             ELSE
494               ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
495         &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
496             ENDIF
497             IF (k.EQ.Nr) THEN
498               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
499         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
500             ELSE
501               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
502         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
503             ENDIF
504  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
505            ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)           DO j=jMin,jMax
506       &                  -((rC(K-1)/atm_po)**atm_kappa) )*0.5            DO i=iMin,iMax
507            ddRm1=ddRp1               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
508                 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
509              ENDDO
510             ENDDO
511    C end: Energy Conserving Form, No hFac  --
512  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
513    
514            DO j=jMin,jMax         ELSEIF (integr_GeoPot.EQ.1) THEN
515              DO i=iMin,iMax  C  --  Finite Volume Form, with Part-Cell Topo, linear in P by Half level
516                ddRp=ddRp1  C---------
517                ddRm=ddRm1  C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
518                IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.  C   Note: a true Finite Volume form should be linear between 2 Interf_W :
519                IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0.  C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
520                phiHyd(i,j,K)=phiHyd(i,j,K-1)  C   also: if Interface_W at the middle between tracer levels, this form
521       &           -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1))  C     is close to the Energy Cons. form in the Interior, except for the
522       &             +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) )  C     non-linearity in PI(p)
523  C             Old code (atmos-exact) looked like this  C---------
524  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddRm1*             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
525  Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))       &                   -((rC( k )/atm_Po)**atm_kappa) )
526              ENDDO             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
527         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
528             DO j=jMin,jMax
529              DO i=iMin,iMax
530               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
531                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
532    #ifdef NONLIN_FRSURF
533                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
534    #endif
535                 phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
536         &          *ddPIm*alphaRho(i,j)
537               ELSE
538                 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
539               ENDIF
540                 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
541            ENDDO            ENDDO
542          ENDIF           ENDDO
543    C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
544    C-----------------------------------------------------------------------
545    
546           ELSEIF ( integr_GeoPot.EQ.2
547         &     .OR. integr_GeoPot.EQ.3 ) THEN
548    C  --  Finite Difference Form, with Part-Cell Topo,
549    C       works with Interface_W  at the middle between 2.Tracer_Level
550    C        and  with Tracer_Level at the middle between 2.Interface_W.
551    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
552    C  Finite Difference formulation consistent with Partial Cell,
553    C   Valid & accurate if Interface_W at middle between tracer levels
554    C   linear in p between 2 Tracer levels ; conserve energy in the Interior
555    C---------
556             IF (k.EQ.1) THEN
557               ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
558         &                   -((rC( k )/atm_Po)**atm_kappa) )
559             ELSE
560               ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
561         &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
562             ENDIF
563             IF (k.EQ.Nr) THEN
564               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
565         &                   -((rF(k+1)/atm_Po)**atm_kappa) )
566             ELSE
567               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
568         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
569             ENDIF
570             rec_dRm = oneRL/(rF(k)-rC(k))
571             rec_dRp = oneRL/(rC(k)-rF(k+1))
572             DO j=jMin,jMax
573              DO i=iMin,iMax
574               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
575                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
576    #ifdef NONLIN_FRSURF
577                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
578    #endif
579                 phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
580         &                      +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
581         &                     )*alphaRho(i,j)
582               ELSE
583                 phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
584               ENDIF
585                 phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
586              ENDDO
587             ENDDO
588    C end: Finite Difference Form, with Part-Cell Topo
589    C-----------------------------------------------------------------------
590    
591           ELSE
592             STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
593           ENDIF
594    
595    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
596        ELSE        ELSE
597          STOP 'CALC_PHI_HYD: We should never reach this point!'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
598        ENDIF        ENDIF
599    
600  #endif        IF ( .NOT. useFVgradPhi ) THEN
601    C--   r-coordinate and r*-coordinate cases:
602    
603           IF ( momPressureForcing ) THEN
604            CALL CALC_GRAD_PHI_HYD(
605         I                         k, bi, bj, iMin,iMax, jMin,jMax,
606         I                         phiHydC, alphaRho, tFld, sFld,
607         O                         dPhiHydX, dPhiHydY,
608         I                         myTime, myIter, myThid)
609           ENDIF
610    
611    #ifndef DISABLE_SIGMA_CODE
612          ELSE
613    C--   else (SigmaCoords part)
614    
615           IF ( fluidIsWater ) THEN
616            STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
617           ENDIF
618           IF ( momPressureForcing ) THEN
619            CALL CALC_GRAD_PHI_FV(
620         I                         k, bi, bj, iMin,iMax, jMin,jMax,
621         I                         phiHydF, phiHydU, pKappaF, pKappaU,
622         O                         dPhiHydX, dPhiHydY,
623         I                         myTime, myIter, myThid)
624           ENDIF
625           DO j=jMin,jMax
626             DO i=iMin,iMax
627               phiHydF(i,j) = phiHydU(i,j)
628             ENDDO
629           ENDDO
630    
631    #endif /* DISABLE_SIGMA_CODE */
632    C--   end if-not/else useFVgradPhi
633          ENDIF
634    
635    C---   Diagnose Phi at boundary r=R_low :
636    C       = Ocean bottom pressure (Ocean, Z-coord.)
637    C       = Sea-surface height    (Ocean, P-coord.)
638    C       = Top atmosphere height (Atmos, P-coord.)
639          IF (useDiagPhiRlow) THEN
640            CALL DIAGS_PHI_RLOW(
641         I                      k, bi, bj, iMin,iMax, jMin,jMax,
642         I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
643         I                      myTime, myIter, myThid)
644          ENDIF
645    
646    C---   Diagnose Full Hydrostatic Potential at cell center level
647            CALL DIAGS_PHI_HYD(
648         I                      k, bi, bj, iMin,iMax, jMin,jMax,
649         I                      phiHydC,
650         I                      myTime, myIter, myThid)
651    
652    #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
653    
654        return        RETURN
655        end        END
Legend:
Removed from v.1.8.2.5  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22