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