/[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.40 by jmc, Tue Mar 16 00:08:27 2010 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    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_PHI_HYD  C     !ROUTINE: CALC_PHI_HYD
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE CALC_PHI_HYD(        SUBROUTINE CALC_PHI_HYD(
11       I                         bi, bj, iMin, iMax, jMin, jMax, K,       I                         bi, bj, iMin, iMax, jMin, jMax, k,
12       I                         tFld, sFld,       I                         tFld, sFld,
13       U                         phiHyd,       U                         phiHydF,
14       I                         myThid)       O                         phiHydC, dPhiHydX, dPhiHydY,
15         I                         myTime, myIter, myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
19  C     | o Integrate the hydrostatic relation to find the Hydros. |  C     | o Integrate the hydrostatic relation to find the Hydros. |
20  C     *==========================================================*  C     *==========================================================*
21  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)
22  C     | On entry:                                                |  C     | On entry:
23  C     |   tFld,sFld     are the current thermodynamics quantities|  C     |   tFld,sFld     are the current thermodynamics quantities
24  C     |                 (unchanged on exit)                      |  C     |                 (unchanged on exit)
25  C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
26  C     |                 at cell centers (tracer points)          |  C     |                at middle between tracer points k-1,k
27  C     |                 - 1:k-1 layers are valid                 |  C     | On exit:
28  C     |                 - k:Nr layers are invalid                |  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
29  C     |   phiHyd(i,j,k) is the hydrostatic Potential             |  C     |                at cell centers (tracer points), level k
30  C     |  (ocean only_^) at cell the interface k (w point above)  |  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
31  C     | On exit:                                                 |  C     |                at middle between tracer points k,k+1
32  C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
33  C     |                 at cell centers (tracer points)          |  C     |                at cell centers (tracer points), level k
34  C     |                 - 1:k layers are valid                   |  C     | integr_GeoPot allows to select one integration method
35  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)                                  |  
36  C     *==========================================================*  C     *==========================================================*
37  C     \ev  C     \ev
38  C     !USES:  C     !USES:
# Line 45  C     == Global variables == Line 42  C     == Global variables ==
42  #include "GRID.h"  #include "GRID.h"
43  #include "EEPARAMS.h"  #include "EEPARAMS.h"
44  #include "PARAMS.h"  #include "PARAMS.h"
 c #include "FFIELDS.h"  
45  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
46  #include "tamc.h"  #include "tamc.h"
47  #include "tamc_keys.h"  #include "tamc_keys.h"
# Line 55  c #include "FFIELDS.h" Line 51  c #include "FFIELDS.h"
51    
52  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
53  C     == Routine arguments ==  C     == Routine arguments ==
54        INTEGER bi,bj,iMin,iMax,jMin,jMax,K  C     bi, bj, k  :: tile and level indices
55    C     iMin,iMax,jMin,jMax :: computational domain
56    C     tFld       :: potential temperature
57    C     sFld       :: salinity
58    C     phiHydF    :: hydrostatic potential anomaly at middle between
59    C                   2 centers (entry: Interf_k ; output: Interf_k+1)
60    C     phiHydC    :: hydrostatic potential anomaly at cell center
61    C     dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
62    C     myTime     :: current time
63    C     myIter     :: current iteration number
64    C     myThid     :: thread number for this instance of the routine.
65          INTEGER bi,bj,iMin,iMax,jMin,jMax,k
66        _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)
67        _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)
68        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
69        INTEGER myThid        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70                _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71          _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73          _RL myTime
74          INTEGER myIter, myThid
75    
76  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
77    
78  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
79  C     == Local variables ==  C     == Local variables ==
80        INTEGER i,j, Kp1        INTEGER i,j
81        _RL zero, one, half        _RL zero, one, half
82        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL dRloc,dRlocKp1,locAlpha        _RL dRlocM,dRlocP, ddRloc, locAlpha
84        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
85          _RL surfPhiFac
86          PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
87          LOGICAL useDiagPhiRlow, addSurfPhiAnom
88  CEOP  CEOP
89          useDiagPhiRlow = .TRUE.
90        zero = 0. _d 0        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GT.3
91        one  = 1. _d 0        surfPhiFac = 0.
92        half = .5 _d 0        IF (addSurfPhiAnom) surfPhiFac = 1.
93    
94  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
95  C  Atmosphere:    C  Atmosphere:
96  C integr_GeoPot => select one option for the integration of the Geopotential:  C integr_GeoPot => select one option for the integration of the Geopotential:
97  C   = 0 : Energy Conserving Form, No hFac ;  C   = 0 : Energy Conserving Form, accurate with Topo full cell;
98  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;
99  C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels  C   =2,3: Finite Difference Form, with Part-Cell,
100  C     2 : case Tracer level at the middle of InterFace_W;  C         linear in P between 2 Tracer levels.
101  C     3 : case InterFace_W  at the middle of Tracer levels;  C       can handle both cases: Tracer lev at the middle of InterFace_W
102    C                          and InterFace_W at the middle of Tracer lev;
103  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
104    
105  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 103  C---+----1----+----2----+----3----+----4 Line 119  C---+----1----+----2----+----3----+----4
119       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
120  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
121    
122        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN  C--   Initialize phiHydF to zero :
123    C     note: atmospheric_loading or Phi_topo anomaly are incorporated
124    C           later in S/R calc_grad_phi_hyd
125          IF (k.EQ.1) THEN
126            DO j=1-Oly,sNy+Oly
127             DO i=1-Olx,sNx+Olx
128               phiHydF(i,j) = 0.
129             ENDDO
130            ENDDO
131          ENDIF
132    
133    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
134          IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
135  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
136  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
137  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
138    #ifdef ALLOW_AUTODIFF_TAMC
139    CADJ GENERAL
140    #endif /* ALLOW_AUTODIFF_TAMC */
141    
142          dRloc=drC(k)          IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
143          IF (k.EQ.1) dRloc=drF(1)  C---    Calculate density
144          IF (k.EQ.Nr) THEN  #ifdef ALLOW_AUTODIFF_TAMC
145            dRlocKp1=0.            kkey = (ikey-1)*Nr + k
146    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
147    CADJ &     kind = isbyte
148    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
149    CADJ &     kind = isbyte
150    #endif /* ALLOW_AUTODIFF_TAMC */
151              CALL FIND_RHO_2D(
152         I              iMin, iMax, jMin, jMax, k,
153         I              tFld(1-OLx,1-OLy,k,bi,bj),
154         I              sFld(1-OLx,1-OLy,k,bi,bj),
155         O              alphaRho,
156         I              k, bi, bj, myThid )
157          ELSE          ELSE
158            dRlocKp1=drC(k+1)            DO j=jMin,jMax
159               DO i=iMin,iMax
160                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
161               ENDDO
162              ENDDO
163          ENDIF          ENDIF
164    
165  C--     If this is the top layer we impose the boundary condition  #ifdef ALLOW_SHELFICE
166  C       P(z=eta) = P(atmospheric_loading)  C     mask rho, so that there is no contribution of phiHyd from
167          IF (k.EQ.1) THEN  C     overlying shelfice (whose density we do not know)
168            IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
169    C- note: does not work for down_slope pkg which needs rho below the bottom.
170    C    setting rho=0 above the ice-shelf base is enough (and works in both cases)
171    C    but might be slower (--> keep original masking if not using down_slope pkg)
172             DO j=jMin,jMax
173              DO i=iMin,iMax
174               IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
175              ENDDO
176             ENDDO
177            ELSEIF ( useShelfIce ) THEN
178             DO j=jMin,jMax
179              DO i=iMin,iMax
180               alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
181              ENDDO
182             ENDDO
183            ENDIF
184    #endif /* ALLOW_SHELFICE */
185    
186    #ifdef ALLOW_MOM_COMMON
187    C Quasi-hydrostatic terms are added in as if they modify the buoyancy
188            IF (quasiHydrostatic) THEN
189             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
190            ENDIF
191    #endif /* ALLOW_MOM_COMMON */
192    
193    #ifdef NONLIN_FRSURF
194            IF (k.EQ.1 .AND. addSurfPhiAnom) THEN
195            DO j=jMin,jMax            DO j=jMin,jMax
196              DO i=iMin,iMax              DO i=iMin,iMax
197                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
198         &                      *gravity*alphaRho(i,j)*recip_rhoConst
199              ENDDO              ENDDO
200            ENDDO            ENDDO
201          ENDIF          ENDIF
202    #endif /* NONLIN_FRSURF */
203    
204  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)  
205    
206  C Quasi-hydrostatic terms are added in as if they modify the buoyancy         IF (integr_GeoPot.EQ.1) THEN
207          IF (quasiHydrostatic) THEN  C  --  Finite Volume Form
          CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)  
         ENDIF  
208    
209  C       Hydrostatic pressure at cell centers           DO j=jMin,jMax
         DO j=jMin,jMax  
210            DO i=iMin,iMax            DO i=iMin,iMax
 #ifdef      ALLOW_AUTODIFF_TAMC  
 c           Patrick, is this directive correct or even necessary in  
 c           this new code?  
 c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...  
 c           within the k-loop.  
 CADJ GENERAL  
 #endif      /* ALLOW_AUTODIFF_TAMC */  
211    
212  CmlC---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
213  CmlC           which has not been used to date since it does not  C           which has not been used to date since it does not
214  CmlC           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
215  CmlC  C
216  Cml          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN             phiHydC(i,j)=phiHydF(i,j)
217  Cml           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)       &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
218  Cml     &          + hFacC(i,j,k,bi,bj)             phiHydF(i,j)=phiHydF(i,j)
219  Cml     &            *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
220  Cml     &          + gravity*etaN(i,j,bi,bj)            ENDDO
221  Cml          ENDIF           ENDDO
222  Cml           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  
223  Cml     &         drF(K)*gravity*alphaRho(i,j)*recip_rhoConst         ELSE
224  Cml           phiHyd(i,j,k)=phiHyd(i,j,k)+  C  --  Finite Difference Form
225  Cml     &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  
226  CmlC-----------------------------------------------------------------------           dRlocM=half*drC(k)
227             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
228             IF (k.EQ.Nr) THEN
229               dRlocP=rC(k)-rF(k+1)
230             ELSE
231               dRlocP=half*drC(k+1)
232             ENDIF
233    
234             DO j=jMin,jMax
235              DO i=iMin,iMax
236    
237  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
238  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
239  C  C
240                            phiHydC(i,j)=phiHydF(i,j)
241              phiHyd(i,j,k)=phiHyd(i,j,k)+       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
242       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst              phiHydF(i,j)=phiHydC(i,j)
243              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
244       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst            ENDDO
245  C-----------------------------------------------------------------------           ENDDO
246    
247  C---------- Compute bottom pressure deviation from gravity*rho0*H  C  --  end if integr_GeoPot = ...
248  C           This has to be done starting from phiHyd at the current         ENDIF
 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-----------------------------------------------------------------------  
249    
250            ENDDO  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
251          ENDDO        ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
           
       ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN  
252  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
253  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density before
254  C       before integrating g*rho over the current layer/interface  C       integrating (1/rho)_prime*dp over the current layer/interface
255  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
256  CADJ GENERAL  CADJ GENERAL
257  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
258    
259          dRloc=drC(k)          IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
260          IF (k.EQ.1) dRloc=drF(1)  C--     Calculate density
         IF (k.EQ.Nr) THEN  
           dRlocKp1=0.  
         ELSE  
           dRlocKp1=drC(k+1)  
         ENDIF  
   
         IF (k.EQ.1) THEN  
           DO j=jMin,jMax  
             DO i=iMin,iMax  
               phiHyd(i,j,k) = phi0surf(i,j,bi,bj)  
             ENDDO  
           ENDDO  
         ENDIF  
   
 C       Calculate density  
261  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
262              kkey = (ikey-1)*Nr + k            kkey = (ikey-1)*Nr + k
263  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,
264  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
265    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
266    CADJ &     kind = isbyte
267  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
268          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,            CALL FIND_RHO_2D(
269       &                 tFld, sFld,       I              iMin, iMax, jMin, jMax, k,
270       &                 alphaRho, myThid)       I              tFld(1-OLx,1-OLy,k,bi,bj),
271         I              sFld(1-OLx,1-OLy,k,bi,bj),
272         O              alphaRho,
273         I              k, bi, bj, myThid )
274  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
275  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte,
276    CADJ &     kind = isbyte
277  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
278            ELSE
279              DO j=jMin,jMax
280               DO i=iMin,iMax
281                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
282               ENDDO
283              ENDDO
284            ENDIF
285    
286    C--     Calculate specific volume anomaly : alpha_prime = 1/rho - alpha_Cst
 C       Hydrostatic pressure at cell centers  
287          DO j=jMin,jMax          DO j=jMin,jMax
288            DO i=iMin,iMax            DO i=iMin,iMax
289              locAlpha=alphaRho(i,j)+rhoConst              locAlpha=alphaRho(i,j)+rhoConst
290              IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha              alphaRho(i,j)=maskC(i,j,k,bi,bj)*
291         &              (one/locAlpha - recip_rhoConst)
292              ENDDO
293            ENDDO
294    
295  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-----------------------------------------------------------------------  
296    
297  C---------- This discretization is the "energy conserving" form         IF (integr_GeoPot.EQ.1) THEN
298  C           which has been used since at least Adcroft et al., MWR 1997  C  --  Finite Volume Form
 C  
299    
300              phiHyd(i,j,k)=phiHyd(i,j,k)+           DO j=jMin,jMax
301       &          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  
302    
303  C-----------------------------------------------------------------------  C---------- This discretization is the "finite volume" form
304    C           which has not been used to date since it does not
305    C           conserve KE+PE exactly even though it is more natural
306    C
307               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
308                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
309    #ifdef NONLIN_FRSURF
310                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
311    #endif
312                 phiHydC(i,j) = ddRloc*alphaRho(i,j)
313    c--to reproduce results of c48d_post: uncomment those 4+1 lines
314    c            phiHydC(i,j)=phiHydF(i,j)
315    c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)
316    c            phiHydF(i,j)=phiHydF(i,j)
317    c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
318               ELSE
319                 phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)
320    c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)
321               ENDIF
322    c-- and comment this last one:
323                 phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)
324    c-----
325              ENDDO
326             ENDDO
327    
328           ELSE
329    C  --  Finite Difference Form, with Part-Cell Bathy
330    
331             dRlocM=half*drC(k)
332             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
333             IF (k.EQ.Nr) THEN
334               dRlocP=rC(k)-rF(k+1)
335             ELSE
336               dRlocP=half*drC(k+1)
337             ENDIF
338             rec_dRm = one/(rF(k)-rC(k))
339             rec_dRp = one/(rC(k)-rF(k+1))
340    
341  C---------- Compute gravity*(sea surface elevation) first           DO j=jMin,jMax
342  C           This has to be done starting from phiHyd at the current            DO i=iMin,iMax
343  C           tracer point and .5 of the cell's thickness has to be  
344  C           substracted from hFacC  C---------- This discretization is the "energy conserving" form
             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-----------------------------------------------------------------------  
345    
346               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
347                 ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
348    #ifdef NONLIN_FRSURF
349                 ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
350    #endif
351                 phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
352         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
353         &                     )*alphaRho(i,j)
354               ELSE
355                 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
356               ENDIF
357                 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
358            ENDDO            ENDDO
359          ENDDO           ENDDO
360    
361    C  --  end if integr_GeoPot = ...
362           ENDIF
363    
364        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
365  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
366  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
367  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
368  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
369    
370  C       Integrate d Phi / d pi  C--     virtual potential temperature anomaly (including water vapour effect)
371            DO j=jMin,jMax
372             DO i=iMin,iMax
373              alphaRho(i,j)=maskC(i,j,k,bi,bj)
374         &             *( tFld(i,j,k,bi,bj)*(sFld(i,j,k,bi,bj)*atm_Rq+one)
375         &               -tRef(k) )
376             ENDDO
377            ENDDO
378    
379    C---    Integrate d Phi / d pi
380    
381        IF (integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
382  C  --  Energy Conserving Form, No hFac  --  C  --  Energy Conserving Form, accurate with Full cell topo  --
383  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
384  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
385  Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
386  C             condition is simply Phi-prime(Ro_surf)=0.  C             condition is simply Phi-prime(Ro_surf)=0.
387  C           o convention ddPI > 0 (same as drF & drC)  C           o convention ddPI > 0 (same as drF & drC)
388  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
389          IF (K.EQ.1) THEN           IF (k.EQ.1) THEN
390            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
391       &                  -((rC(K)/atm_Po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
392            DO j=jMin,jMax           ELSE
393             DO i=iMin,iMax             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
394               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)       &                   -((rC( k )/atm_Po)**atm_kappa) )*half
395       &         +ddPIp*maskC(i,j,K,bi,bj)           ENDIF
396       &               *(tFld(I,J,K,bi,bj)-tRef(K))           IF (k.EQ.Nr) THEN
397             ENDDO             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
398            ENDDO       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
399          ELSE           ELSE
400               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
401         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
402             ENDIF
403  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
404            ddPI=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)           DO j=jMin,jMax
405       &                 -((rC( K )/atm_Po)**atm_kappa) )*0.5            DO i=iMin,iMax
406            DO j=jMin,jMax               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
407             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  
408            ENDDO            ENDDO
409          ENDIF           ENDDO
410  C end: Energy Conserving Form, No hFac  --  C end: Energy Conserving Form, No hFac  --
411  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
412    
413        ELSEIF (integr_GeoPot.EQ.1) THEN         ELSEIF (integr_GeoPot.EQ.1) THEN
414  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
415  C---------  C---------
416  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
417  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 :
418  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)
419  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
420  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
421  C     non-linearity in PI(p)  C     non-linearity in PI(p)
422  C---------  C---------
423          IF (K.EQ.1) THEN             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
424            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)       &                   -((rC( k )/atm_Po)**atm_kappa) )
425       &                  -((rC(K)/atm_Po)**atm_kappa) )             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
426            DO j=jMin,jMax       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
427             DO i=iMin,iMax           DO j=jMin,jMax
428               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)            DO i=iMin,iMax
429       &         +ddPIp*_hFacC(I,J, K ,bi,bj)             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
430       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
431             ENDDO  #ifdef NONLIN_FRSURF
432            ENDDO               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
433          ELSE  #endif
434            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)               phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
435       &                  -((rF( K )/atm_Po)**atm_kappa) )       &          *ddPIm*alphaRho(i,j)
436            ddPIp=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)             ELSE
437       &                  -((rC( K )/atm_Po)**atm_kappa) )               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
438            DO j=jMin,jMax             ENDIF
439             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  
440            ENDDO            ENDDO
441          ENDIF           ENDDO
442  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
443  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
444    
445        ELSEIF (integr_GeoPot.EQ.3) THEN         ELSEIF ( integr_GeoPot.EQ.2
446  C  --  Finite Difference Form, with hFac, Interface_W = middle  --       &     .OR. integr_GeoPot.EQ.3 ) THEN
447    C  --  Finite Difference Form, with Part-Cell Topo,
448    C       works with Interface_W  at the middle between 2.Tracer_Level
449    C        and  with Tracer_Level at the middle between 2.Interface_W.
450  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
451  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
452  C   Valid & accurate if Interface_W at middle between tracer levels  C   Valid & accurate if Interface_W at middle between tracer levels
453  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
454  C---------  C---------
455          Kp1 = min(Nr,K+1)           IF (k.EQ.1) THEN
456          IF (K.EQ.1) THEN             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
457            ratioRm=0.5*drF(K)/(rF(k)-rC(K))       &                   -((rC( k )/atm_Po)**atm_kappa) )
458            ratioRp=drF(K)*recip_drC(Kp1)           ELSE
459            ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
460       &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0       &                   -((rC( k )/atm_Po)**atm_kappa) )*half
461            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)           ENDIF
462       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )             IF (k.EQ.Nr) THEN
463            DO j=jMin,jMax             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
464             DO i=iMin,iMax       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
465               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)           ELSE
466       &       +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
467       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
468       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))           ENDIF
469       &               * maskC(i,j, K ,bi,bj)           rec_dRm = one/(rF(k)-rC(k))
470             ENDDO           rec_dRp = one/(rC(k)-rF(k+1))
471            ENDDO           DO j=jMin,jMax
472          ELSE            DO i=iMin,iMax
473            ratioRm=drF(K)*recip_drC(K)             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
474            ratioRp=drF(K)*recip_drC(Kp1)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
475            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)  #ifdef NONLIN_FRSURF
476       &                  -((rC( K )/atm_Po)**atm_kappa) )               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
477            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)  #endif
478       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm
479            DO j=jMin,jMax       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )
480             DO i=iMin,iMax       &                    *alphaRho(i,j)
481               phiHyd(i,j,K) = phiHyd(i,j,K-1)             ELSE
482       &        + ddPIm*0.5               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
483       &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))             ENDIF
484       &               * 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  
485            ENDDO            ENDDO
486          ENDIF           ENDDO
487  C end: Finite Difference Form, with hFac, Interface_W = middle  --  C end: Finite Difference Form, with Part-Cell Topo
488  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
489    
490        ELSE         ELSE
491          STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'           STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
492        ENDIF         ENDIF
493    
494  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
495        ELSE        ELSE
496          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
497        ENDIF        ENDIF
498    
499    C---   Diagnose Phi at boundary r=R_low :
500    C       = Ocean bottom pressure (Ocean, Z-coord.)
501    C       = Sea-surface height    (Ocean, P-coord.)
502    C       = Top atmosphere height (Atmos, P-coord.)
503          IF (useDiagPhiRlow) THEN
504            CALL DIAGS_PHI_RLOW(
505         I                      k, bi, bj, iMin,iMax, jMin,jMax,
506         I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
507         I                      myTime, myIter, myThid)
508          ENDIF
509    
510    C---   Diagnose Full Hydrostatic Potential at cell center level
511            CALL DIAGS_PHI_HYD(
512         I                      k, bi, bj, iMin,iMax, jMin,jMax,
513         I                      phiHydC,
514         I                      myTime, myIter, myThid)
515    
516          IF (momPressureForcing) THEN
517            CALL CALC_GRAD_PHI_HYD(
518         I                         k, bi, bj, iMin,iMax, jMin,jMax,
519         I                         phiHydC, alphaRho, tFld, sFld,
520         O                         dPhiHydX, dPhiHydY,
521         I                         myTime, myIter, myThid)
522          ENDIF
523    
524  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
525    
526        RETURN        RETURN
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22