/[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.17 by adcroft, Thu Sep 27 18:14:20 2001 UTC revision 1.17.6.1 by heimbach, Fri Mar 7 23:10:20 2003 UTC
# Line 7  CBOP Line 7  CBOP
7  C     !ROUTINE: CALC_PHI_HYD  C     !ROUTINE: CALC_PHI_HYD
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE CALC_PHI_HYD(        SUBROUTINE CALC_PHI_HYD(
10       I                         bi, bj, iMin, iMax, jMin, jMax, K,       I                         bi, bj, iMin, iMax, jMin, jMax, k,
11       I                         theta, salt,       I                         tFld, sFld,
12       U                         phiHyd,       U                         phiHydF,
13       I                         myThid)       O                         phiHydC, dPhiHydX, dPhiHydY,
14         I                         myTime, myIter, myThid)
15  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
16  C     *==========================================================*  C     *==========================================================*
17  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
18  C     | o Integrate the hydrostatic relation to find the Hydros. |  C     | o Integrate the hydrostatic relation to find the Hydros. |
19  C     *==========================================================*  C     *==========================================================*
20  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)
21  C     | On entry:                                                |  C     | On entry:
22  C     |   theta,salt    are the current thermodynamics quantities|  C     |   tFld,sFld     are the current thermodynamics quantities
23  C     |                 (unchanged on exit)                      |  C     |                 (unchanged on exit)
24  C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
25  C     |                 at cell centers (tracer points)          |  C     |                at middle between tracer points k-1,k
26  C     |                 - 1:k-1 layers are valid                 |  C     | On exit:
27  C     |                 - k:Nr layers are invalid                |  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
28  C     |   phiHyd(i,j,k) is the hydrostatic Potential             |  C     |                at cell centers (tracer points), level k
29  C     |  (ocean only_^) at cell the interface k (w point above)  |  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
30  C     | On exit:                                                 |  C     |                at middle between tracer points k,k+1
31  C     |   phiHyd(i,j,1:k) is the hydrostatic Potential           |  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
32  C     |                 at cell centers (tracer points)          |  C     |                at cell centers (tracer points), level k
33  C     |                 - 1:k layers are valid                   |  C     | integr_GeoPot allows to select one integration method
34  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)                                  |  
35  C     *==========================================================*  C     *==========================================================*
36  C     \ev  C     \ev
37  C     !USES:  C     !USES:
# Line 49  C     == Global variables == Line 45  C     == Global variables ==
45  #include "tamc.h"  #include "tamc.h"
46  #include "tamc_keys.h"  #include "tamc_keys.h"
47  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
48    #include "SURFACE.h"
49    #include "DYNVARS.h"
50    
51  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
52  C     == Routine arguments ==  C     == Routine arguments ==
53        INTEGER bi,bj,iMin,iMax,jMin,jMax,K  C     bi, bj, k  :: tile and level indices
54        _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  C     iMin,iMax,jMin,jMax :: computational domain
55        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  C     tFld       :: potential temperature
56        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     sFld       :: salinity
57        INTEGER myThid  C     phiHydF    :: hydrostatic potential anomaly at middle between
58    C                   2 centers (entry: Interf_k ; output: Interf_k+1)
59    C     phiHydC    :: hydrostatic potential anomaly at cell center
60    C     dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
61    C     myTime     :: current time
62    C     myIter     :: current iteration number
63    C     myThid     :: thread number for this instance of the routine.
64          INTEGER bi,bj,iMin,iMax,jMin,jMax,k
65          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
66          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
67    c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
68          _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69          _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70          _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72          _RL myTime
73          INTEGER myIter, myThid
74                
75  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
76    
77  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
78  C     == Local variables ==  C     == Local variables ==
79        INTEGER i,j, Kp1        INTEGER i,j
80        _RL zero, one, half        _RL zero, one, half
81        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL dRloc,dRlocKp1        _RL dRlocM,dRlocP, ddRloc, locAlpha
83        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
84          _RL surfPhiFac
85          INTEGER iMnLoc,jMnLoc
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 100  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    C---    Calculate density
143    #ifdef ALLOW_AUTODIFF_TAMC
144            kkey = (ikey-1)*Nr + k
145    CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
146    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
147    #endif /* ALLOW_AUTODIFF_TAMC */
148            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
149         &                 tFld, sFld,
150         &                 alphaRho, myThid)
151    
152          dRloc=drC(k)  C Quasi-hydrostatic terms are added in as if they modify the buoyancy
153          IF (k.EQ.1) dRloc=drF(1)          IF (quasiHydrostatic) THEN
154          IF (k.EQ.Nr) THEN           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
           dRlocKp1=0.  
         ELSE  
           dRlocKp1=drC(k+1)  
155          ENDIF          ENDIF
156    
157  C--     If this is the top layer we impose the boundary condition  #ifdef NONLIN_FRSURF
158  C       P(z=eta) = P(atmospheric_loading)          IF (k.EQ.1 .AND. addSurfPhiAnom) THEN
         IF (k.EQ.1) THEN  
159            DO j=jMin,jMax            DO j=jMin,jMax
160              DO i=iMin,iMax              DO i=iMin,iMax
161  C             *NOTE* The loading should go here but has not been implemented yet                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
162                phiHyd(i,j,k)=0.       &                      *gravity*alphaRho(i,j)*recip_rhoConst
163              ENDDO              ENDDO
164            ENDDO            ENDDO
165          ENDIF          ENDIF
166    #endif /* NONLIN_FRSURF */
167    
168    C----  Hydrostatic pressure at cell centers
169    
170           IF (integr_GeoPot.EQ.1) THEN
171    C  --  Finite Volume Form
172    
173             DO j=jMin,jMax
174              DO i=iMin,iMax
175    
176    C---------- This discretization is the "finite volume" form
177    C           which has not been used to date since it does not
178    C           conserve KE+PE exactly even though it is more natural
179    C
180               phiHydC(i,j)=phiHydF(i,j)
181         &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
182               phiHydF(i,j)=phiHydF(i,j)
183         &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
184              ENDDO
185             ENDDO
186    
187           ELSE
188    C  --  Finite Difference Form
189    
190             dRlocM=half*drC(k)
191             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
192             IF (k.EQ.Nr) THEN
193               dRlocP=rC(k)-rF(k+1)
194             ELSE
195               dRlocP=half*drC(k+1)
196             ENDIF
197    
198  C       Calculate density           DO j=jMin,jMax
199              DO i=iMin,iMax
200    
201    C---------- This discretization is the "energy conserving" form
202    C           which has been used since at least Adcroft et al., MWR 1997
203    C
204                phiHydC(i,j)=phiHydF(i,j)
205         &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
206                phiHydF(i,j)=phiHydC(i,j)
207         &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
208              ENDDO
209             ENDDO
210    
211    C  --  end if integr_GeoPot = ...
212           ENDIF
213            
214    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
215          ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
216    C       This is the hydrostatic pressure calculation for the Ocean
217    C       which uses the FIND_RHO() routine to calculate density
218    C       before integrating (1/rho)'*dp over the current layer/interface
219    #ifdef      ALLOW_AUTODIFF_TAMC
220    CADJ GENERAL
221    #endif      /* ALLOW_AUTODIFF_TAMC */
222    
223    C--     Calculate density
224  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
225              kkey = (ikey-1)*Nr + k              kkey = (ikey-1)*Nr + k
226  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
227  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
228  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
229          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
230       &                 theta, salt,       &                 tFld, sFld,
231       &                 alphaRho, myThid)       &                 alphaRho, myThid)
232    #ifdef ALLOW_AUTODIFF_TAMC
233    CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
234    #endif /* ALLOW_AUTODIFF_TAMC */
235    
236  C       Hydrostatic pressure at cell centers  C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
237          DO j=jMin,jMax          DO j=jMin,jMax
238            DO i=iMin,iMax            DO i=iMin,iMax
239  #ifdef      ALLOW_AUTODIFF_TAMC              locAlpha=alphaRho(i,j)+rhoConst
240  c           Patrick, is this directive correct or even necessary in              alphaRho(i,j)=maskC(i,j,k,bi,bj)*
241  c           this new code?       &              (one/locAlpha - recip_rhoConst)
242  c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...            ENDDO
243  c           within the k-loop.          ENDDO
244  CADJ GENERAL  
245  #endif      /* ALLOW_AUTODIFF_TAMC */  C----  Hydrostatic pressure at cell centers
246    
247           IF (integr_GeoPot.EQ.1) THEN
248    C  --  Finite Volume Form
249    
250             DO j=jMin,jMax
251              DO i=iMin,iMax
252    
253  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
254  C           which has not been used to date since it does not  C           which has not been used to date since it does not
255  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
256  C  C
257  c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
258  c    &              drF(K)*gravity*alphaRho(i,j)*recip_rhoConst               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
259  c           phiHyd(i,j,k)=phiHyd(i,j,k)+  #ifdef NONLIN_FRSURF
260  c    &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
261  C-----------------------------------------------------------------------  #endif
262                 phiHydC(i,j) = ddRloc*alphaRho(i,j)
263    c--to reproduce results of c48d_post: uncomment those 4+1 lines
264    c            phiHydC(i,j)=phiHydF(i,j)
265    c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)
266    c            phiHydF(i,j)=phiHydF(i,j)
267    c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
268               ELSE
269                 phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)
270    c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)
271               ENDIF
272    c-- and comment this last one:
273                 phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)
274    c-----
275              ENDDO
276             ENDDO
277    
278           ELSE
279    C  --  Finite Difference Form, with Part-Cell Bathy
280    
281             dRlocM=half*drC(k)
282             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
283             IF (k.EQ.Nr) THEN
284               dRlocP=rC(k)-rF(k+1)
285             ELSE
286               dRlocP=half*drC(k+1)
287             ENDIF
288             rec_dRm = one/(rF(k)-rC(k))
289             rec_dRp = one/(rC(k)-rF(k+1))
290    
291             DO j=jMin,jMax
292              DO i=iMin,iMax
293    
294  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
295  C           which has been used since at least Adcroft et al., MWR 1997  
296  C             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
297              phiHyd(i,j,k)=phiHyd(i,j,k)+               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
298       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst  #ifdef NONLIN_FRSURF
299              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
300       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst  #endif
301  C-----------------------------------------------------------------------               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
302         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
303         &                     )*alphaRho(i,j)
304               ELSE
305                 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
306               ENDIF
307                 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
308            ENDDO            ENDDO
309          ENDDO           ENDDO
           
310    
311    C  --  end if integr_GeoPot = ...
312           ENDIF
313    
314        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
315  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
316  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
317  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
318  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
319    
320  C       Integrate d Phi / d pi  C---    Integrate d Phi / d pi
321    
322        IF (Integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
323  C  --  Energy Conserving Form, No hFac  --  C  --  Energy Conserving Form, accurate with Full cell topo  --
324  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
325  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
326  Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
327  C             condition is simply Phi-prime(Ro_surf)=0.  C             condition is simply Phi-prime(Ro_surf)=0.
328  C           o convention ddPI > 0 (same as drF & drC)  C           o convention ddPI > 0 (same as drF & drC)
329  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
330          IF (K.EQ.1) THEN           IF (k.EQ.1) THEN
331            ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
332       &                  -((rC(K)/atm_po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
333            DO j=jMin,jMax           ELSE
334             DO i=iMin,iMax             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
335               phiHyd(i,j,K)=       &                   -((rC( k )/atm_Po)**atm_kappa) )*half
336       &          ddPIp*maskC(i,j,K,bi,bj)           ENDIF
337       &               *(theta(I,J,K,bi,bj)-tRef(K))           IF (k.EQ.Nr) THEN
338             ENDDO             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
339            ENDDO       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
340          ELSE           ELSE
341               ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
342         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
343             ENDIF
344  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
345            ddPI=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)           DO j=jMin,jMax
346       &                 -((rC( K )/atm_po)**atm_kappa) )*0.5            DO i=iMin,iMax
347            DO j=jMin,jMax               phiHydC(i,j) = phiHydF(i,j)
348             DO i=iMin,iMax       &          +ddPIm*maskC(i,j,k,bi,bj)
349                phiHyd(i,j,K)=phiHyd(i,j,K-1)       &                *(tFld(i,j,k,bi,bj)-tRef(k))
350       &           +ddPI*maskC(i,j,K-1,bi,bj)               phiHydF(i,j) = phiHydC(i,j)
351       &                *(theta(I,J,K-1,bi,bj)-tRef(K-1))       &          +ddPIp*maskC(i,j,k,bi,bj)
352       &           +ddPI*maskC(i,j, K ,bi,bj)       &                *(tFld(i,j,k,bi,bj)-tRef(k))
      &                *(theta(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 &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))  
            ENDDO  
353            ENDDO            ENDDO
354          ENDIF           ENDDO
355  C end: Energy Conserving Form, No hFac  --  C end: Energy Conserving Form, No hFac  --
356  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
357    
358        ELSEIF (Integr_GeoPot.EQ.1) THEN         ELSEIF (integr_GeoPot.EQ.1) THEN
359  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
360  C---------  C---------
361  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
362  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 :
# Line 224  C   also: if Interface_W at the middle b Line 365  C   also: if Interface_W at the middle b
365  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
366  C     non-linearity in PI(p)  C     non-linearity in PI(p)
367  C---------  C---------
368          IF (K.EQ.1) THEN             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
369            ddPIp=atm_cp*( ((rF(K)/atm_po)**atm_kappa)       &                   -((rC( k )/atm_Po)**atm_kappa) )
370       &                  -((rC(K)/atm_po)**atm_kappa) )             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
371            DO j=jMin,jMax       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
372             DO i=iMin,iMax           DO j=jMin,jMax
373               phiHyd(i,j,K) =            DO i=iMin,iMax
374       &          ddPIp*hFacC(I,J, K ,bi,bj)             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
375       &               *(theta(I,J, K ,bi,bj)-tRef( K ))               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
376             ENDDO  #ifdef NONLIN_FRSURF
377            ENDDO               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
378          ELSE  #endif
379            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)               phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
380       &                  -((rF( K )/atm_po)**atm_kappa) )       &          *ddPIm*maskC(i,j,k,bi,bj)
381            ddPIp=atm_cp*( ((rF( K )/atm_po)**atm_kappa)       &                *(tFld(i,j,k,bi,bj)-tRef(k))
382       &                  -((rC( K )/atm_po)**atm_kappa) )             ELSE
383            DO j=jMin,jMax               phiHydC(i,j) = phiHydF(i,j)
384             DO i=iMin,iMax       &          +ddPIm*maskC(i,j,k,bi,bj)
385               phiHyd(i,j,K) = phiHyd(i,j,K-1)       &                *(tFld(i,j,k,bi,bj)-tRef(k))
386       &         +ddPIm*hFacC(I,J,K-1,bi,bj)             ENDIF
387       &               *(theta(I,J,K-1,bi,bj)-tRef(K-1))               phiHydF(i,j) = phiHydC(i,j)
388       &         +ddPIp*hFacC(I,J, K ,bi,bj)       &          +ddPIp*maskC(i,j,k,bi,bj)
389       &               *(theta(I,J, K ,bi,bj)-tRef( K ))       &                *(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) =  
      &        ( ddPIm*max(zero, hFacC(i,j,K,bi,bj)-half)  
      &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)-half) )  
      &               *(theta(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  
      &               *(theta(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) )  
      &               *(theta(i,j, K ,bi,bj)-tRef( K ))  
      &               * maskC(i,j, K ,bi,bj)  
            ENDDO  
390            ENDDO            ENDDO
391          ENDIF           ENDDO
392  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
393  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
394    
395        ELSEIF (Integr_GeoPot.EQ.3) THEN         ELSEIF ( integr_GeoPot.EQ.2
396  C  --  Finite Difference Form, with hFac, Interface_W = middle  --       &     .OR. integr_GeoPot.EQ.3 ) THEN
397    C  --  Finite Difference Form, with Part-Cell Topo,
398    C       works with Interface_W  at the middle between 2.Tracer_Level
399    C        and  with Tracer_Level at the middle between 2.Interface_W.
400  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
401  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
402  C   Valid & accurate if Interface_W at middle between tracer levels  C   Valid & accurate if Interface_W at middle between tracer levels
403  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
404  C---------  C---------
405          Kp1 = min(Nr,K+1)           IF (k.EQ.1) THEN
406          IF (K.EQ.1) THEN             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
407            ratioRm=0.5*drF(K)/(rF(k)-rC(K))       &                   -((rC( k )/atm_Po)**atm_kappa) )
408            ratioRp=drF(K)*recip_drC(Kp1)           ELSE
409            ddPIm=atm_cp*( ((rF( K )/atm_po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
410       &                  -((rC( K )/atm_po)**atm_kappa) ) * 2. _d 0       &                   -((rC( k )/atm_Po)**atm_kappa) )*half
411            ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)           ENDIF
412       &                  -((rC(Kp1)/atm_po)**atm_kappa) )             IF (k.EQ.Nr) THEN
413            DO j=jMin,jMax             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
414             DO i=iMin,iMax       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
415               phiHyd(i,j,K) =           ELSE
416       &        ( ddPIm*max(zero,(hFacC(i,j,K,bi,bj)-one)*ratioRm+half)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
417       &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
418       &               *(theta(i,j, K ,bi,bj)-tRef( K ))           ENDIF
419       &               * maskC(i,j, K ,bi,bj)           rec_dRm = one/(rF(k)-rC(k))
420             ENDDO           rec_dRp = one/(rC(k)-rF(k+1))
421            ENDDO           DO j=jMin,jMax
422          ELSE            DO i=iMin,iMax
423            ratioRm=drF(K)*recip_drC(K)             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
424            ratioRp=drF(K)*recip_drC(Kp1)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
425            ddPIm=atm_cp*( ((rC(K-1)/atm_po)**atm_kappa)  #ifdef NONLIN_FRSURF
426       &                  -((rC( K )/atm_po)**atm_kappa) )               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
427            ddPIp=atm_cp*( ((rC( K )/atm_po)**atm_kappa)  #endif
428       &                  -((rC(Kp1)/atm_po)**atm_kappa) )               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm
429            DO j=jMin,jMax       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )
430             DO i=iMin,iMax       &                *(tFld(i,j,k,bi,bj)-tRef(k))
431               phiHyd(i,j,K) = phiHyd(i,j,K-1)             ELSE
432       &        + ddPIm*0.5               phiHydC(i,j) = phiHydF(i,j)
433       &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))       &          +ddPIm*maskC(i,j,k,bi,bj)
434       &               * maskC(i,j,K-1,bi,bj)       &                *(tFld(I,J,k,bi,bj)-tRef(k))
435       &        +(ddPIm*max(zero,(hFacC(i,j,K,bi,bj)-one)*ratioRm+half)             ENDIF
436       &         +ddPIp*min(zero, hFacC(i,j,K,bi,bj)*ratioRp     -half) )               phiHydF(i,j) = phiHydC(i,j)
437       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &          +ddPIp*maskC(i,j,k,bi,bj)
438       &               * maskC(i,j, K ,bi,bj)       &                *(tFld(I,J,k,bi,bj)-tRef(k))
            ENDDO  
439            ENDDO            ENDDO
440          ENDIF           ENDDO
441  C end: Finite Difference Form, with hFac, Interface_W = middle  --  C end: Finite Difference Form, with Part-Cell Topo
442  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
443    
444        ELSE         ELSE
445          STOP 'CALC_PHI_HYD: Bad Integr_GeoPot option !'           STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
446        ENDIF         ENDIF
447    
448  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
449        ELSE        ELSE
450          STOP 'CALC_PHI_HYD: We should never reach this point!'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
451          ENDIF
452    
453    C---   Diagnose Phi at boundary r=R_low :
454    C       = Ocean bottom pressure (Ocean, Z-coord.)
455    C       = Sea-surface height    (Ocean, P-coord.)
456    C       = Top atmosphere height (Atmos, P-coord.)
457          IF (useDiagPhiRlow) THEN
458            CALL DIAGS_PHI_RLOW(
459         I                      k, bi, bj, iMin,iMax, jMin,jMax,
460         I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
461         I                      myTime, myIter, myThid)  
462          ENDIF
463    
464    C---   Diagnose Full Hydrostatic Potential at cell center level
465            CALL DIAGS_PHI_HYD(
466         I                      k, bi, bj, iMin,iMax, jMin,jMax,
467         I                      phiHydC,
468         I                      myTime, myIter, myThid)
469    
470          IF (momPressureForcing) THEN
471            iMnLoc = MAX(1-Olx+1,iMin)
472            jMnLoc = MAX(1-Oly+1,jMin)
473            CALL CALC_GRAD_PHI_HYD(
474         I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
475         I                         phiHydC, alphaRho, tFld, sFld,
476         O                         dPhiHydX, dPhiHydY,
477         I                         myTime, myIter, myThid)  
478        ENDIF        ENDIF
479    
480  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.17.6.1
  ViewVC Help
Powered by ViewVC 1.1.22