/[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.29 by jmc, Tue Feb 18 15:30:47 2003 UTC revision 1.41 by jmc, Wed Apr 11 04:02:05 2012 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
# Line 11  C     !INTERFACE: Line 12  C     !INTERFACE:
12       I                         tFld, sFld,       I                         tFld, sFld,
13       U                         phiHydF,       U                         phiHydF,
14       O                         phiHydC, dPhiHydX, dPhiHydY,       O                         phiHydC, dPhiHydX, dPhiHydY,
15       I                         myTime, myIter, myThid)       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     |   phiHydF(i,j) is the hydrostatic Potential anomaly  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
26  C     |                at middle between tracer points k-1,k  C     |                at middle between tracer points k-1,k
27  C     | On exit:  C     | On exit:
28  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
29  C     |                at cell centers (tracer points), level k  C     |                at cell centers (tracer points), level k
30  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
31  C     |                at middle between tracer points k,k+1  C     |                at middle between tracer points k,k+1
32  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
33  C     |                at cell centers (tracer points), level k  C     |                at cell centers (tracer points), level k
34  C     | integr_GeoPot allows to select one integration method  C     | integr_GeoPot allows to select one integration method
# Line 50  C     == Global variables == Line 51  C     == Global variables ==
51    
52  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
53  C     == Routine arguments ==  C     == Routine arguments ==
54  C     bi, bj, k  :: tile and level indices  C     bi, bj, k  :: tile and level indices
55  C     iMin,iMax,jMin,jMax :: computational domain  C     iMin,iMax,jMin,jMax :: computational domain
56  C     tFld       :: potential temperature  C     tFld       :: potential temperature
57  C     sFld       :: salinity  C     sFld       :: salinity
# Line 67  C     myThid     :: thread number for th Line 68  C     myThid     :: thread number for th
68  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
69        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _RL myTime        _RL myTime
74        INTEGER myIter, myThid        INTEGER myIter, myThid
75          
76  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
77    
78  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
# Line 82  C     == Local variables == Line 83  C     == Local variables ==
83        _RL dRlocM,dRlocP, ddRloc, locAlpha        _RL dRlocM,dRlocP, ddRloc, locAlpha
84        _RL ddPIm, ddPIp, rec_dRm, rec_dRp        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
85        _RL surfPhiFac        _RL surfPhiFac
       INTEGER iMnLoc,jMnLoc  
86        PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )        PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
87        LOGICAL useDiagPhiRlow, addSurfPhiAnom        LOGICAL useDiagPhiRlow, addSurfPhiAnom
88  CEOP  CEOP
89        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
90        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GT.3        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
91        surfPhiFac = 0.        surfPhiFac = 0.
92        IF (addSurfPhiAnom) surfPhiFac = 1.        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, accurate with Topo full cell;  C   = 0 : Energy Conserving Form, accurate with Topo full cell;
98  C   = 1 : Finite Volume Form, with Part-Cell, 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 Part-Cell,  C   =2,3: Finite Difference Form, with Part-Cell,
100  C         linear in P between 2 Tracer levels.  C         linear in P between 2 Tracer levels.
101  C       can handle both cases: Tracer lev at the middle of InterFace_W  C       can handle both cases: Tracer lev at the middle of InterFace_W
102  C                          and InterFace_W at the middle of Tracer lev;  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    
# Line 119  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  C--   Initialize phiHydF to zero :  C--   Initialize phiHydF to zero :
123  C     note: atmospheric_loading or Phi_topo anomaly are incorporated  C     note: atmospheric_loading or Phi_topo anomaly are incorporated
124  C           later in S/R calc_grad_phi_hyd  C           later in S/R calc_grad_phi_hyd
125        IF (k.EQ.1) THEN        IF (k.EQ.1) THEN
126          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
127           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
128             phiHydF(i,j) = 0.             phiHydF(i,j) = 0.
129           ENDDO           ENDDO
130          ENDDO          ENDDO
# Line 135  C---+----1----+----2----+----3----+----4 Line 135  C---+----1----+----2----+----3----+----4
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  #ifdef ALLOW_AUTODIFF_TAMC
139  CADJ GENERAL  CADJ GENERAL
140  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
141    
142            IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
143  C---    Calculate density  C---    Calculate density
144  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
145          kkey = (ikey-1)*Nr + k            kkey = (ikey-1)*Nr + k
146  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,
147  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  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 */  #endif /* ALLOW_AUTODIFF_TAMC */
151          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,            CALL FIND_RHO_2D(
152       &                 tFld, sFld,       I              iMin, iMax, jMin, jMax, k,
153       &                 alphaRho, myThid)       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
158              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
164    
165    #ifdef ALLOW_SHELFICE
166    C     mask rho, so that there is no contribution of phiHyd from
167    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  C Quasi-hydrostatic terms are added in as if they modify the buoyancy
188          IF (quasiHydrostatic) THEN          IF (quasiHydrostatic) THEN
189           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)           CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
190          ENDIF          ENDIF
191    #endif /* ALLOW_MOM_COMMON */
192    
193  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
194          IF (k.EQ.1 .AND. addSurfPhiAnom) THEN          IF ( addSurfPhiAnom .AND.
195         &       uniformFreeSurfLev .AND. k.EQ.1 ) THEN
196            DO j=jMin,jMax            DO j=jMin,jMax
197              DO i=iMin,iMax              DO i=iMin,iMax
198                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
# Line 170  C----  Hydrostatic pressure at cell cent Line 207  C----  Hydrostatic pressure at cell cent
207         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
208  C  --  Finite Volume Form  C  --  Finite Volume Form
209    
          DO j=jMin,jMax  
           DO i=iMin,iMax  
   
210  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
211  C           which has not been used to date since it does not  C           which has not been used to date since it does not
212  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
213  C  
214             phiHydC(i,j)=phiHydF(i,j)          IF ( uniformFreeSurfLev ) THEN
215       &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst           DO j=jMin,jMax
216             phiHydF(i,j)=phiHydF(i,j)            DO i=iMin,iMax
217       &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst              phiHydC(i,j) = phiHydF(i,j)
218         &            + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
219                phiHydF(i,j) = phiHydF(i,j)
220         &                 + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
221              ENDDO
222             ENDDO
223            ELSE
224             DO j=jMin,jMax
225              DO i=iMin,iMax
226               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
227                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
228    #ifdef NONLIN_FRSURF
229                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
230    #endif
231                phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst
232               ELSE
233                phiHydC(i,j) = phiHydF(i,j)
234         &              + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
235               ENDIF
236                phiHydF(i,j) = phiHydC(i,j)
237         &              + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
238            ENDDO            ENDDO
239           ENDDO           ENDDO
240            ENDIF
241    
242         ELSE         ELSE
243  C  --  Finite Difference Form  C  --  Finite Difference Form
244    
245           dRlocM=half*drC(k)  C---------- This discretization is the "energy conserving" form
246           IF (k.EQ.1) dRlocM=rF(k)-rC(k)  C           which has been used since at least Adcroft et al., MWR 1997
          IF (k.EQ.Nr) THEN  
            dRlocP=rC(k)-rF(k+1)  
          ELSE  
            dRlocP=half*drC(k+1)  
          ENDIF  
247    
248            dRlocM=half*drC(k)
249            IF (k.EQ.1) dRlocM=rF(k)-rC(k)
250            IF (k.EQ.Nr) THEN
251              dRlocP=rC(k)-rF(k+1)
252            ELSE
253              dRlocP=half*drC(k+1)
254            ENDIF
255            IF ( uniformFreeSurfLev ) THEN
256           DO j=jMin,jMax           DO j=jMin,jMax
257            DO i=iMin,iMax            DO i=iMin,iMax
258                phiHydC(i,j) = phiHydF(i,j)
259  C---------- This discretization is the "energy conserving" form       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
260  C           which has been used since at least Adcroft et al., MWR 1997              phiHydF(i,j) = phiHydC(i,j)
261  C       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
262              phiHydC(i,j)=phiHydF(i,j)            ENDDO
263             ENDDO
264            ELSE
265             rec_dRm = one/(rF(k)-rC(k))
266             rec_dRp = one/(rC(k)-rF(k+1))
267             DO j=jMin,jMax
268              DO i=iMin,iMax
269               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
270                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
271    #ifdef NONLIN_FRSURF
272                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
273    #endif
274                phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
275         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
276         &                     )*gravity*alphaRho(i,j)*recip_rhoConst
277               ELSE
278                phiHydC(i,j) = phiHydF(i,j)
279       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
280              phiHydF(i,j)=phiHydC(i,j)             ENDIF
281                phiHydF(i,j) = phiHydC(i,j)
282       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
283            ENDDO            ENDDO
284           ENDDO           ENDDO
285            ENDIF
286    
287  C  --  end if integr_GeoPot = ...  C  --  end if integr_GeoPot = ...
288         ENDIF         ENDIF
289            
290  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
291        ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
292  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
293  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density before
294  C       before integrating (1/rho)'*dp over the current layer/interface  C       integrating (1/rho)_prime*dp over the current layer/interface
295  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
296  CADJ GENERAL  CADJ GENERAL
297  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
298    
299            IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
300  C--     Calculate density  C--     Calculate density
301  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
302              kkey = (ikey-1)*Nr + k            kkey = (ikey-1)*Nr + k
303  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,
304  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
305    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
306    CADJ &     kind = isbyte
307  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
308          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,            CALL FIND_RHO_2D(
309       &                 tFld, sFld,       I              iMin, iMax, jMin, jMax, k,
310       &                 alphaRho, myThid)       I              tFld(1-OLx,1-OLy,k,bi,bj),
311         I              sFld(1-OLx,1-OLy,k,bi,bj),
312         O              alphaRho,
313         I              k, bi, bj, myThid )
314  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
315  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte,
316    CADJ &     kind = isbyte
317  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
318            ELSE
319              DO j=jMin,jMax
320               DO i=iMin,iMax
321                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
322               ENDDO
323              ENDDO
324            ENDIF
325    
326  C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst  C--     Calculate specific volume anomaly : alpha_prime = 1/rho - alpha_Cst
327          DO j=jMin,jMax          DO j=jMin,jMax
328            DO i=iMin,iMax            DO i=iMin,iMax
329              locAlpha=alphaRho(i,j)+rhoConst              locAlpha=alphaRho(i,j)+rhoConst
# Line 253  C  --  Finite Volume Form Line 343  C  --  Finite Volume Form
343  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
344  C           which has not been used to date since it does not  C           which has not been used to date since it does not
345  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
346  C  
347             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
348               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
349  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
350               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
351  #endif  #endif
352               phiHydC(i,j) = ddRloc*alphaRho(i,j)               phiHydC(i,j) = ddRloc*alphaRho(i,j)
353  c--to reproduce results of c48d_post: uncomment those 4+1 lines  c--to reproduce results of c48d_post: uncomment those 4+1 lines
354  c            phiHydC(i,j)=phiHydF(i,j)  c            phiHydC(i,j)=phiHydF(i,j)
355  c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)  c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)
356  c            phiHydF(i,j)=phiHydF(i,j)  c            phiHydF(i,j)=phiHydF(i,j)
# Line 293  C  --  Finite Difference Form, with Part Line 383  C  --  Finite Difference Form, with Part
383    
384  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
385    
386             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
387               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
388  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
389               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
# Line 317  C       This is the hydrostatic geopoten Line 407  C       This is the hydrostatic geopoten
407  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
408  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
409    
410    C--     virtual potential temperature anomaly (including water vapour effect)
411            DO j=jMin,jMax
412             DO i=iMin,iMax
413              alphaRho(i,j)=maskC(i,j,k,bi,bj)
414         &             *( tFld(i,j,k,bi,bj)*(sFld(i,j,k,bi,bj)*atm_Rq+one)
415         &               -tRef(k) )
416             ENDDO
417            ENDDO
418    
419  C---    Integrate d Phi / d pi  C---    Integrate d Phi / d pi
420    
421         IF (integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
422  C  --  Energy Conserving Form, accurate with Full cell topo  --  C  --  Energy Conserving Form, accurate with Full cell topo  --
423  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
424  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
425  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
426  C             condition is simply Phi-prime(Ro_surf)=0.  C             condition is simply Phi-prime(Ro_surf)=0.
427  C           o convention ddPI > 0 (same as drF & drC)  C           o convention ddPI > 0 (same as drF & drC)
428  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
# Line 339  C--------------------------------------- Line 438  C---------------------------------------
438       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
439           ELSE           ELSE
440             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
441       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
442           ENDIF           ENDIF
443  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
444           DO j=jMin,jMax           DO j=jMin,jMax
445            DO i=iMin,iMax            DO i=iMin,iMax
446               phiHydC(i,j) = phiHydF(i,j)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
447       &          +ddPIm*maskC(i,j,k,bi,bj)               phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
      &                *(tFld(i,j,k,bi,bj)-tRef(k))  
              phiHydF(i,j) = phiHydC(i,j)  
      &          +ddPIp*maskC(i,j,k,bi,bj)  
      &                *(tFld(i,j,k,bi,bj)-tRef(k))  
448            ENDDO            ENDDO
449           ENDDO           ENDDO
450  C end: Energy Conserving Form, No hFac  --  C end: Energy Conserving Form, No hFac  --
# Line 362  C  Finite Volume formulation consistent Line 457  C  Finite Volume formulation consistent
457  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 :
458  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)
459  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
460  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
461  C     non-linearity in PI(p)  C     non-linearity in PI(p)
462  C---------  C---------
463             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
# Line 371  C--------- Line 466  C---------
466       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
467           DO j=jMin,jMax           DO j=jMin,jMax
468            DO i=iMin,iMax            DO i=iMin,iMax
469             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
470               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
471  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
472               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
473  #endif  #endif
474               phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0               phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
475       &          *ddPIm*maskC(i,j,k,bi,bj)       &          *ddPIm*alphaRho(i,j)
      &                *(tFld(i,j,k,bi,bj)-tRef(k))  
476             ELSE             ELSE
477               phiHydC(i,j) = phiHydF(i,j)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
      &          +ddPIm*maskC(i,j,k,bi,bj)  
      &                *(tFld(i,j,k,bi,bj)-tRef(k))  
478             ENDIF             ENDIF
479               phiHydF(i,j) = phiHydC(i,j)               phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
      &          +ddPIp*maskC(i,j,k,bi,bj)  
      &                *(tFld(i,j,k,bi,bj)-tRef(k))  
480            ENDDO            ENDDO
481           ENDDO           ENDDO
482  C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level  C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
# Line 394  C--------------------------------------- Line 484  C---------------------------------------
484    
485         ELSEIF ( integr_GeoPot.EQ.2         ELSEIF ( integr_GeoPot.EQ.2
486       &     .OR. integr_GeoPot.EQ.3 ) THEN       &     .OR. integr_GeoPot.EQ.3 ) THEN
487  C  --  Finite Difference Form, with Part-Cell Topo,  C  --  Finite Difference Form, with Part-Cell Topo,
488  C       works with Interface_W  at the middle between 2.Tracer_Level  C       works with Interface_W  at the middle between 2.Tracer_Level
489  C        and  with Tracer_Level at the middle between 2.Interface_W.  C        and  with Tracer_Level at the middle between 2.Interface_W.
490  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
491  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
492  C   Valid & accurate if Interface_W at middle between tracer levels  C   Valid & accurate if Interface_W at middle between tracer levels
493  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
494  C---------  C---------
495           IF (k.EQ.1) THEN           IF (k.EQ.1) THEN
496             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
# Line 414  C--------- Line 504  C---------
504       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
505           ELSE           ELSE
506             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
507       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
508           ENDIF           ENDIF
509           rec_dRm = one/(rF(k)-rC(k))           rec_dRm = one/(rF(k)-rC(k))
510           rec_dRp = one/(rC(k)-rF(k+1))           rec_dRp = one/(rC(k)-rF(k+1))
511           DO j=jMin,jMax           DO j=jMin,jMax
512            DO i=iMin,iMax            DO i=iMin,iMax
513             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
514               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
515  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
516               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
517  #endif  #endif
518               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm
519       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )
520       &                *(tFld(i,j,k,bi,bj)-tRef(k))       &                    *alphaRho(i,j)
521             ELSE             ELSE
522               phiHydC(i,j) = phiHydF(i,j)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
      &          +ddPIm*maskC(i,j,k,bi,bj)  
      &                *(tFld(I,J,k,bi,bj)-tRef(k))  
523             ENDIF             ENDIF
524               phiHydF(i,j) = phiHydC(i,j)               phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
      &          +ddPIp*maskC(i,j,k,bi,bj)  
      &                *(tFld(I,J,k,bi,bj)-tRef(k))  
525            ENDDO            ENDDO
526           ENDDO           ENDDO
527  C end: Finite Difference Form, with Part-Cell Topo  C end: Finite Difference Form, with Part-Cell Topo
# Line 458  C       = Top atmosphere height (Atmos, Line 544  C       = Top atmosphere height (Atmos,
544          CALL DIAGS_PHI_RLOW(          CALL DIAGS_PHI_RLOW(
545       I                      k, bi, bj, iMin,iMax, jMin,jMax,       I                      k, bi, bj, iMin,iMax, jMin,jMax,
546       I                      phiHydF, phiHydC, alphaRho, tFld, sFld,       I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
547       I                      myTime, myIter, myThid)         I                      myTime, myIter, myThid)
548        ENDIF        ENDIF
549    
550  C---   Diagnose Full Hydrostatic Potential at cell center level  C---   Diagnose Full Hydrostatic Potential at cell center level
# Line 467  C---   Diagnose Full Hydrostatic Potenti Line 553  C---   Diagnose Full Hydrostatic Potenti
553       I                      phiHydC,       I                      phiHydC,
554       I                      myTime, myIter, myThid)       I                      myTime, myIter, myThid)
555    
556        IF (momPressureForcing) THEN        IF (momPressureForcing) THEN
         iMnLoc = MAX(1-Olx+1,iMin)  
         jMnLoc = MAX(1-Oly+1,jMin)  
557          CALL CALC_GRAD_PHI_HYD(          CALL CALC_GRAD_PHI_HYD(
558       I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,       I                         k, bi, bj, iMin,iMax, jMin,jMax,
559       I                         phiHydC, alphaRho, tFld, sFld,       I                         phiHydC, alphaRho, tFld, sFld,
560       O                         dPhiHydX, dPhiHydY,       O                         dPhiHydX, dPhiHydY,
561       I                         myTime, myIter, myThid)         I                         myTime, myIter, myThid)
562        ENDIF        ENDIF
563    
564  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22