/[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.27 by jmc, Sun Feb 9 02:58:39 2003 UTC revision 1.33 by mlosch, Tue Feb 7 11:47:49 2006 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       O                         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     |   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 46  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 56  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          _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)        _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
# Line 69  C     == Routine arguments == Line 77  C     == Routine arguments ==
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        INTEGER iMnLoc,jMnLoc        INTEGER iMnLoc,jMnLoc
87        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 )
88        LOGICAL useDiagPhiRlow        LOGICAL useDiagPhiRlow, addSurfPhiAnom
89  CEOP  CEOP
90        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
91          addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GT.3
92          surfPhiFac = 0.
93          IF (addSurfPhiAnom) surfPhiFac = 1.
94    
95  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
96  C  Atmosphere:    C  Atmosphere:  
97  C integr_GeoPot => select one option for the integration of the Geopotential:  C integr_GeoPot => select one option for the integration of the Geopotential:
98  C   = 0 : Energy Conserving Form, No hFac ;  C   = 0 : Energy Conserving Form, accurate with Topo full cell;
99  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;
100  C   =2,3: Finite Difference Form, with hFac, linear in P between 2 Tracer levels  C   =2,3: Finite Difference Form, with Part-Cell,
101  C     2 : case Tracer level at the middle of InterFace_W;  C         linear in P between 2 Tracer levels.
102  C     3 : case InterFace_W  at the middle of Tracer levels;  C       can handle both cases: Tracer lev at the middle of InterFace_W
103    C                          and InterFace_W at the middle of Tracer lev;
104  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105    
106  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 107  C---+----1----+----2----+----3----+----4 Line 120  C---+----1----+----2----+----3----+----4
120       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
121  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
122    
123    C--   Initialize phiHydF to zero :
124    C     note: atmospheric_loading or Phi_topo anomaly are incorporated
125    C           later in S/R calc_grad_phi_hyd
126          IF (k.EQ.1) THEN
127            DO j=1-Oly,sNy+Oly
128             DO i=1-Olx,sNx+Olx
129               phiHydF(i,j) = 0.
130             ENDDO
131            ENDDO
132          ENDIF
133    
134  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
135        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN        IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
136  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
137  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
138  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
# Line 117  C       before integrating g*rho over th Line 140  C       before integrating g*rho over th
140  CADJ GENERAL  CADJ GENERAL
141  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
142    
143          dRloc=drC(k)  C---    Calculate density
         IF (k.EQ.1) dRloc=drF(1)  
         IF (k.EQ.Nr) THEN  
           dRlocKp1=0.  
         ELSE  
           dRlocKp1=drC(k+1)  
         ENDIF  
   
 C--     If this is the top layer we impose the boundary condition  
 C       P(z=eta) = P(atmospheric_loading)  
         IF (k.EQ.1) THEN  
           DO j=jMin,jMax  
             DO i=iMin,iMax  
 c             phiHyd(i,j,k) = phi0surf(i,j,bi,bj)  
               phiHyd(i,j,k) = 0.  
             ENDDO  
           ENDDO  
         ENDIF  
   
 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
# Line 145  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_ Line 149  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_
149          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
150       &                 tFld, sFld,       &                 tFld, sFld,
151       &                 alphaRho, myThid)       &                 alphaRho, myThid)
152    #ifdef ALLOW_SHELFICE
153    C     mask rho, so that there is no contribution of phiHyd from
154    C     overlying shelfice (whose density we do not know)
155            IF ( useShelfIce ) THEN
156             DO j=jMin,jMax
157              DO i=iMin,iMax
158               alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
159              ENDDO
160             ENDDO
161            ENDIF
162    #endif /* ALLOW_SHELFICE */
163    
164    #ifdef ALLOW_DIAGNOSTICS
165            IF ( useDiagnostics )
166         &   CALL DIAGNOSTICS_FILL(alphaRho,'RHOAnoma',k,1,2,bi,bj,myThid)
167    #endif
168    
169  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
170          IF (quasiHydrostatic) THEN          IF (quasiHydrostatic) THEN
171           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
172          ENDIF          ENDIF
173    
174  C---   Diagnose Hydrostatic pressure at the bottom:  #ifdef NONLIN_FRSURF
175         IF (useDiagPhiRlow) THEN          IF (k.EQ.1 .AND. addSurfPhiAnom) THEN
176            CALL DIAGS_PHI_RLOW(            DO j=jMin,jMax
177       I                        k, bi, bj, iMin,iMax, jMin,jMax,              DO i=iMin,iMax
178       I                        phiHyd, alphaRho, tFld, sFld,                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
179       I                        myTime, myIter, myThid)         &                      *gravity*alphaRho(i,j)*recip_rhoConst
180         ENDIF              ENDDO
181              ENDDO
182            ENDIF
183    #endif /* NONLIN_FRSURF */
184    
185  C---   Hydrostatic pressure at cell centers  C----  Hydrostatic pressure at cell centers
186    
187         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
188  C  --  Finite Volume Form  C  --  Finite Volume Form
# Line 171  C---------- This discretization is the " Line 194  C---------- This discretization is the "
194  C           which has not been used to date since it does not  C           which has not been used to date since it does not
195  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
196  C  C
197             IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)             phiHydC(i,j)=phiHydF(i,j)
198       &            + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
199             phiHyd(i,j,k)=phiHyd(i,j,k)+             phiHydF(i,j)=phiHydF(i,j)
200       &       + half*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
   
201            ENDDO            ENDDO
202           ENDDO           ENDDO
203    
204         ELSE         ELSE
205  C  --  Finite Difference Form  C  --  Finite Difference Form
206    
207             dRlocM=half*drC(k)
208             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
209             IF (k.EQ.Nr) THEN
210               dRlocP=rC(k)-rF(k+1)
211             ELSE
212               dRlocP=half*drC(k+1)
213             ENDIF
214    
215           DO j=jMin,jMax           DO j=jMin,jMax
216            DO i=iMin,iMax            DO i=iMin,iMax
217    
218  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
219  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
220  C  C
221              phiHyd(i,j,k)=phiHyd(i,j,k)              phiHydC(i,j)=phiHydF(i,j)
222       &        +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
223              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)              phiHydF(i,j)=phiHydC(i,j)
224       &        +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
   
225            ENDDO            ENDDO
226           ENDDO           ENDDO
227    
# Line 200  C  --  end if integr_GeoPot = ... Line 229  C  --  end if integr_GeoPot = ...
229         ENDIF         ENDIF
230                    
231  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
232        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
233  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
234  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
235  C       before integrating (1/rho)'*dp over the current layer/interface  C       before integrating (1/rho)'*dp over the current layer/interface
# Line 208  C       before integrating (1/rho)'*dp o Line 237  C       before integrating (1/rho)'*dp o
237  CADJ GENERAL  CADJ GENERAL
238  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
239    
         dRloc=drC(k)  
         IF (k.EQ.1) dRloc=drF(1)  
         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  
 c             phiHyd(i,j,k) = phi0surf(i,j,bi,bj)  
               phiHyd(i,j,k) = 0.  
             ENDDO  
           ENDDO  
         ENDIF  
   
240  C--     Calculate density  C--     Calculate density
241  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
242              kkey = (ikey-1)*Nr + k              kkey = (ikey-1)*Nr + k
# Line 238  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_ Line 250  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_
250  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
251  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
252    
253    #ifdef ALLOW_DIAGNOSTICS
254            IF ( useDiagnostics )
255         &   CALL DIAGNOSTICS_FILL(alphaRho,'RHOAnoma',k,1,2,bi,bj,myThid)
256    #endif
257    
258  C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst  C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
259          DO j=jMin,jMax          DO j=jMin,jMax
260            DO i=iMin,iMax            DO i=iMin,iMax
# Line 247  C--     Calculate specific volume anomal Line 264  C--     Calculate specific volume anomal
264            ENDDO            ENDDO
265          ENDDO          ENDDO
266    
 C---   Diagnose Sea-surface height (Hydrostatic geopotential at r=Rlow):  
        IF (useDiagPhiRlow) THEN  
           CALL DIAGS_PHI_RLOW(  
      I                        k, bi, bj, iMin,iMax, jMin,jMax,  
      I                        phiHyd, alphaRho, tFld, sFld,  
      I                        myTime, myIter, myThid)    
        ENDIF  
   
267  C----  Hydrostatic pressure at cell centers  C----  Hydrostatic pressure at cell centers
268    
269         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
# Line 267  C---------- This discretization is the " Line 276  C---------- This discretization is the "
276  C           which has not been used to date since it does not  C           which has not been used to date since it does not
277  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
278  C  C
279              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
280       &          + hFacC(i,j,k,bi,bj)*drF(K)*alphaRho(i,j)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
281              phiHyd(i,j,k)=phiHyd(i,j,k)  #ifdef NONLIN_FRSURF
282       &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*alphaRho(i,j)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
283    #endif
284                 phiHydC(i,j) = ddRloc*alphaRho(i,j)
285    c--to reproduce results of c48d_post: uncomment those 4+1 lines
286    c            phiHydC(i,j)=phiHydF(i,j)
287    c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)
288    c            phiHydF(i,j)=phiHydF(i,j)
289    c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
290               ELSE
291                 phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)
292    c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)
293               ENDIF
294    c-- and comment this last one:
295                 phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)
296    c-----
297            ENDDO            ENDDO
298           ENDDO           ENDDO
299    
300         ELSE         ELSE
301  C  --  Finite Difference Form  C  --  Finite Difference Form, with Part-Cell Bathy
302    
303             dRlocM=half*drC(k)
304             IF (k.EQ.1) dRlocM=rF(k)-rC(k)
305             IF (k.EQ.Nr) THEN
306               dRlocP=rC(k)-rF(k+1)
307             ELSE
308               dRlocP=half*drC(k+1)
309             ENDIF
310             rec_dRm = one/(rF(k)-rC(k))
311             rec_dRp = one/(rC(k)-rF(k+1))
312    
313           DO j=jMin,jMax           DO j=jMin,jMax
314            DO i=iMin,iMax            DO i=iMin,iMax
315    
316  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
317    
318              phiHyd(i,j,k)=phiHyd(i,j,k)             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
319       &          + half*dRloc*alphaRho(i,j)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
320              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)  #ifdef NONLIN_FRSURF
321       &          + half*dRlocKp1*alphaRho(i,j)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
322    #endif
323                 phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
324         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
325         &                     )*alphaRho(i,j)
326               ELSE
327                 phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
328               ENDIF
329                 phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
330            ENDDO            ENDDO
331           ENDDO           ENDDO
332    
333  C  --  end if integr_GeoPot = ...  C  --  end if integr_GeoPot = ...
334         ENDIF         ENDIF
335    
336        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
337  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
338  C       This is the hydrostatic geopotential calculation for the Atmosphere  C       This is the hydrostatic geopotential calculation for the Atmosphere
339  C       The ideal gas law is used implicitly here rather than calculating  C       The ideal gas law is used implicitly here rather than calculating
340  C       the specific volume, analogous to the oceanic case.  C       the specific volume, analogous to the oceanic case.
341    
342  C       Integrate d Phi / d pi  C--     virtual potential temperature anomaly (including water vapour effect)
343            DO j=jMin,jMax
344             DO i=iMin,iMax
345              alphaRho(i,j)=maskC(i,j,k,bi,bj)
346         &             *( tFld(i,j,k,bi,bj)*(sFld(i,j,k,bi,bj)*atm_Rq+one)
347         &               -tRef(k) )
348             ENDDO
349            ENDDO
350    
351        IF (integr_GeoPot.EQ.0) THEN  C---    Integrate d Phi / d pi
352  C  --  Energy Conserving Form, No hFac  --  
353           IF (integr_GeoPot.EQ.0) THEN
354    C  --  Energy Conserving Form, accurate with Full cell topo  --
355  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
356  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
357  Ci    *NOTE* o Working with geopotential Anomaly, the geopotential boundary  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
358  C             condition is simply Phi-prime(Ro_surf)=0.  C             condition is simply Phi-prime(Ro_surf)=0.
359  C           o convention ddPI > 0 (same as drF & drC)  C           o convention ddPI > 0 (same as drF & drC)
360  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
361          IF (K.EQ.1) THEN           IF (k.EQ.1) THEN
362            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
363       &                  -((rC(K)/atm_Po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
364            DO j=jMin,jMax           ELSE
365             DO i=iMin,iMax             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
366  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+       &                   -((rC( k )/atm_Po)**atm_kappa) )*half
367               phiHyd(i,j,K)=           ENDIF
368       &          ddPIp*maskC(i,j,K,bi,bj)           IF (k.EQ.Nr) THEN
369       &               *(tFld(I,J,K,bi,bj)-tRef(K))             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
370             ENDDO       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
371            ENDDO           ELSE
372          ELSE             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
373         &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
374             ENDIF
375  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
376            ddPI=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)           DO j=jMin,jMax
377       &                 -((rC( K )/atm_Po)**atm_kappa) )*0.5            DO i=iMin,iMax
378            DO j=jMin,jMax               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
379             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  
380            ENDDO            ENDDO
381          ENDIF           ENDDO
382  C end: Energy Conserving Form, No hFac  --  C end: Energy Conserving Form, No hFac  --
383  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
384    
385        ELSEIF (integr_GeoPot.EQ.1) THEN         ELSEIF (integr_GeoPot.EQ.1) THEN
386  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
387  C---------  C---------
388  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
389  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 351  C   also: if Interface_W at the middle b Line 392  C   also: if Interface_W at the middle b
392  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
393  C     non-linearity in PI(p)  C     non-linearity in PI(p)
394  C---------  C---------
395          IF (K.EQ.1) THEN             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
396            ddPIp=atm_Cp*( ((rF(K)/atm_Po)**atm_kappa)       &                   -((rC( k )/atm_Po)**atm_kappa) )
397       &                  -((rC(K)/atm_Po)**atm_kappa) )             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
398            DO j=jMin,jMax       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
399             DO i=iMin,iMax           DO j=jMin,jMax
400  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+            DO i=iMin,iMax
401               phiHyd(i,j,K)=             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
402       &          ddPIp*_hFacC(I,J, K ,bi,bj)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
403       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))  #ifdef NONLIN_FRSURF
404             ENDDO               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
405            ENDDO  #endif
406          ELSE               phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
407            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)       &          *ddPIm*alphaRho(i,j)
408       &                  -((rF( K )/atm_Po)**atm_kappa) )             ELSE
409            ddPIp=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
410       &                  -((rC( K )/atm_Po)**atm_kappa) )             ENDIF
411            DO j=jMin,jMax               phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
            DO i=iMin,iMax  
              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  
 c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+  
              phiHyd(i,j,K)=  
      &        ( 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  
412            ENDDO            ENDDO
413          ENDIF           ENDDO
414  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
415  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
416    
417        ELSEIF (integr_GeoPot.EQ.3) THEN         ELSEIF ( integr_GeoPot.EQ.2
418  C  --  Finite Difference Form, with hFac, Interface_W = middle  --       &     .OR. integr_GeoPot.EQ.3 ) THEN
419    C  --  Finite Difference Form, with Part-Cell Topo,
420    C       works with Interface_W  at the middle between 2.Tracer_Level
421    C        and  with Tracer_Level at the middle between 2.Interface_W.
422  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
423  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
424  C   Valid & accurate if Interface_W at middle between tracer levels  C   Valid & accurate if Interface_W at middle between tracer levels
425  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
426  C---------  C---------
427          Kp1 = min(Nr,K+1)           IF (k.EQ.1) THEN
428          IF (K.EQ.1) THEN             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
429            ratioRm=0.5*drF(K)/(rF(k)-rC(K))       &                   -((rC( k )/atm_Po)**atm_kappa) )
430            ratioRp=drF(K)*recip_drC(Kp1)           ELSE
431            ddPIm=atm_Cp*( ((rF( K )/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
432       &                  -((rC( K )/atm_Po)**atm_kappa) ) * 2. _d 0       &                   -((rC( k )/atm_Po)**atm_kappa) )*half
433            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)           ENDIF
434       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )             IF (k.EQ.Nr) THEN
435            DO j=jMin,jMax             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
436             DO i=iMin,iMax       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
437  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+           ELSE
438               phiHyd(i,j,K)=             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
439       &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half
440       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )           ENDIF
441       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))           rec_dRm = one/(rF(k)-rC(k))
442       &               * maskC(i,j, K ,bi,bj)           rec_dRp = one/(rC(k)-rF(k+1))
443             ENDDO           DO j=jMin,jMax
444            ENDDO            DO i=iMin,iMax
445          ELSE             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN
446            ratioRm=drF(K)*recip_drC(K)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
447            ratioRp=drF(K)*recip_drC(Kp1)  #ifdef NONLIN_FRSURF
448            ddPIm=atm_Cp*( ((rC(K-1)/atm_Po)**atm_kappa)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
449       &                  -((rC( K )/atm_Po)**atm_kappa) )  #endif
450            ddPIp=atm_Cp*( ((rC( K )/atm_Po)**atm_kappa)               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm
451       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )
452            DO j=jMin,jMax       &                    *alphaRho(i,j)
453             DO i=iMin,iMax             ELSE
454               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
455       &        + ddPIm*0.5             ENDIF
456       &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))               phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
      &               * maskC(i,j,K-1,bi,bj)  
      &        +(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  
457            ENDDO            ENDDO
458          ENDIF           ENDDO
459  C end: Finite Difference Form, with hFac, Interface_W = middle  --  C end: Finite Difference Form, with Part-Cell Topo
460  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
461    
462        ELSE         ELSE
463          STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'           STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
464        ENDIF         ENDIF
465    
466  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
467        ELSE        ELSE
468          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
469        ENDIF        ENDIF
470    
471    C---   Diagnose Phi at boundary r=R_low :
472    C       = Ocean bottom pressure (Ocean, Z-coord.)
473    C       = Sea-surface height    (Ocean, P-coord.)
474    C       = Top atmosphere height (Atmos, P-coord.)
475          IF (useDiagPhiRlow) THEN
476            CALL DIAGS_PHI_RLOW(
477         I                      k, bi, bj, iMin,iMax, jMin,jMax,
478         I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
479         I                      myTime, myIter, myThid)  
480          ENDIF
481    
482    C---   Diagnose Full Hydrostatic Potential at cell center level
483            CALL DIAGS_PHI_HYD(
484         I                      k, bi, bj, iMin,iMax, jMin,jMax,
485         I                      phiHydC,
486         I                      myTime, myIter, myThid)
487    
488        IF (momPressureForcing) THEN        IF (momPressureForcing) THEN
489          iMnLoc = MAX(1-Olx+1,iMin)          iMnLoc = MAX(1-Olx+1,iMin)
490          jMnLoc = MAX(1-Oly+1,jMin)          jMnLoc = MAX(1-Oly+1,jMin)
491          CALL CALC_GRAD_PHI_HYD(          CALL CALC_GRAD_PHI_HYD(
492       I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,       I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
493       I                         phiHyd, alphaRho, tFld, sFld,       I                         phiHydC, alphaRho, tFld, sFld,
494       O                         dPhiHydX, dPhiHydY,       O                         dPhiHydX, dPhiHydY,
495       I                         myTime, myIter, myThid)         I                         myTime, myIter, myThid)  
496        ENDIF        ENDIF
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22