/[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.35 by jmc, Mon Feb 5 03:20:39 2007 UTC revision 1.44 by jmc, Fri Apr 4 20:54:11 2014 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_AUTODIFF
7    # include "AUTODIFF_OPTIONS.h"
8    #endif
9    
10  CBOP  CBOP
11  C     !ROUTINE: CALC_PHI_HYD  C     !ROUTINE: CALC_PHI_HYD
# Line 12  C     !INTERFACE: Line 15  C     !INTERFACE:
15       I                         tFld, sFld,       I                         tFld, sFld,
16       U                         phiHydF,       U                         phiHydF,
17       O                         phiHydC, dPhiHydX, dPhiHydY,       O                         phiHydC, dPhiHydX, dPhiHydY,
18       I                         myTime, myIter, myThid)       I                         myTime, myIter, myThid )
19  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
20  C     *==========================================================*  C     *==========================================================*
21  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
22  C     | o Integrate the hydrostatic relation to find the Hydros. |  C     | o Integrate the hydrostatic relation to find the Hydros. |
23  C     *==========================================================*  C     *==========================================================*
24  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)
25  C     | On entry:  C     | On entry:
26  C     |   tFld,sFld     are the current thermodynamics quantities  C     |   tFld,sFld     are the current thermodynamics quantities
27  C     |                 (unchanged on exit)  C     |                 (unchanged on exit)
28  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
29  C     |                at middle between tracer points k-1,k  C     |                at middle between tracer points k-1,k
30  C     | On exit:  C     | On exit:
31  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly  C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
32  C     |                at cell centers (tracer points), level k  C     |                at cell centers (tracer points), level k
33  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly  C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
34  C     |                at middle between tracer points k,k+1  C     |                at middle between tracer points k,k+1
35  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)  C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
36  C     |                at cell centers (tracer points), level k  C     |                at cell centers (tracer points), level k
37  C     | integr_GeoPot allows to select one integration method  C     | integr_GeoPot allows to select one integration method
# Line 42  C     == Global variables == Line 45  C     == Global variables ==
45  #include "GRID.h"  #include "GRID.h"
46  #include "EEPARAMS.h"  #include "EEPARAMS.h"
47  #include "PARAMS.h"  #include "PARAMS.h"
48  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
49  #include "tamc.h"  #include "tamc.h"
50  #include "tamc_keys.h"  #include "tamc_keys.h"
51  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
52  #include "SURFACE.h"  #include "SURFACE.h"
53  #include "DYNVARS.h"  #include "DYNVARS.h"
54    
55  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
56  C     == Routine arguments ==  C     == Routine arguments ==
57  C     bi, bj, k  :: tile and level indices  C     bi, bj, k  :: tile and level indices
58  C     iMin,iMax,jMin,jMax :: computational domain  C     iMin,iMax,jMin,jMax :: computational domain
59  C     tFld       :: potential temperature  C     tFld       :: potential temperature
60  C     sFld       :: salinity  C     sFld       :: salinity
# Line 65  C     myThid     :: thread number for th Line 68  C     myThid     :: thread number for th
68        INTEGER bi,bj,iMin,iMax,jMin,jMax,k        INTEGER bi,bj,iMin,iMax,jMin,jMax,k
69        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
70        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
 c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
71        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RL myTime        _RL myTime
76        INTEGER myIter, myThid        INTEGER myIter, myThid
77          
78  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
79    
80  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
81  C     == Local variables ==  C     == Local variables ==
82    C     phiHydU    :: hydrostatic potential anomaly at interface k+1 (Upper k)
83    C     pKappaF    :: (p/Po)^kappa at interface k
84    C     pKappaU    :: (p/Po)^kappa at interface k+1 (Upper k)
85    C     pKappaC    :: (p/Po)^kappa at cell center k
86        INTEGER i,j        INTEGER i,j
       _RL zero, one, half  
87        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88    #ifndef DISABLE_SIGMA_CODE
89          _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90          _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL pKappaC, locDepth, fullDepth
93    #endif /* DISABLE_SIGMA_CODE */
94        _RL dRlocM,dRlocP, ddRloc, locAlpha        _RL dRlocM,dRlocP, ddRloc, locAlpha
95        _RL ddPIm, ddPIp, rec_dRm, rec_dRp        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
96        _RL surfPhiFac        _RL surfPhiFac
       PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )  
97        LOGICAL useDiagPhiRlow, addSurfPhiAnom        LOGICAL useDiagPhiRlow, addSurfPhiAnom
98          LOGICAL useFVgradPhi
99  CEOP  CEOP
100        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
101        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GT.3        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
102          useFVgradPhi   = selectSigmaCoord.NE.0
103    
104        surfPhiFac = 0.        surfPhiFac = 0.
105        IF (addSurfPhiAnom) surfPhiFac = 1.        IF (addSurfPhiAnom) surfPhiFac = 1.
106    
107  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
108  C  Atmosphere:    C  Atmosphere:
109  C integr_GeoPot => select one option for the integration of the Geopotential:  C integr_GeoPot => select one option for the integration of the Geopotential:
110  C   = 0 : Energy Conserving Form, accurate with Topo full cell;  C   = 0 : Energy Conserving Form, accurate with Topo full cell;
111  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;
112  C   =2,3: Finite Difference Form, with Part-Cell,  C   =2,3: Finite Difference Form, with Part-Cell,
113  C         linear in P between 2 Tracer levels.  C         linear in P between 2 Tracer levels.
114  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
115  C                          and InterFace_W at the middle of Tracer lev;  C                          and InterFace_W at the middle of Tracer lev;
116  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
117    
# Line 119  C---+----1----+----2----+----3----+----4 Line 132  C---+----1----+----2----+----3----+----4
132       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
133  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
134    
135  C--   Initialize phiHydF to zero :  C--   Initialize phiHydF to zero :
136  C     note: atmospheric_loading or Phi_topo anomaly are incorporated  C     note: atmospheric_loading or Phi_topo anomaly are incorporated
137  C           later in S/R calc_grad_phi_hyd  C           later in S/R calc_grad_phi_hyd
138        IF (k.EQ.1) THEN        IF (k.EQ.1) THEN
139          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
140           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
141             phiHydF(i,j) = 0.             phiHydF(i,j) = 0.
142           ENDDO           ENDDO
143          ENDDO          ENDDO
# Line 135  C---+----1----+----2----+----3----+----4 Line 148  C---+----1----+----2----+----3----+----4
148  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
149  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
150  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
151  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
152  CADJ GENERAL  CADJ GENERAL
153  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
154    
155            IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
156  C---    Calculate density  C---    Calculate density
157  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
158          kkey = (ikey-1)*Nr + k            kkey = (ikey-1)*Nr + k
159  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,
160  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
161    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
162    CADJ &     kind = isbyte
163  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
164          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,            CALL FIND_RHO_2D(
165       &                 tFld, sFld,       I              iMin, iMax, jMin, jMax, k,
166       &                 alphaRho, myThid)       I              tFld(1-OLx,1-OLy,k,bi,bj),
167         I              sFld(1-OLx,1-OLy,k,bi,bj),
168         O              alphaRho,
169         I              k, bi, bj, myThid )
170            ELSE
171              DO j=jMin,jMax
172               DO i=iMin,iMax
173                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
174               ENDDO
175              ENDDO
176            ENDIF
177    
178  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
179  C     mask rho, so that there is no contribution of phiHyd from  C     mask rho, so that there is no contribution of phiHyd from
180  C     overlying shelfice (whose density we do not know)  C     overlying shelfice (whose density we do not know)
181          IF ( useShelfIce ) THEN          IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
182    C- note: does not work for down_slope pkg which needs rho below the bottom.
183    C    setting rho=0 above the ice-shelf base is enough (and works in both cases)
184    C    but might be slower (--> keep original masking if not using down_slope pkg)
185             DO j=jMin,jMax
186              DO i=iMin,iMax
187               IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
188              ENDDO
189             ENDDO
190            ELSEIF ( useShelfIce ) THEN
191           DO j=jMin,jMax           DO j=jMin,jMax
192            DO i=iMin,iMax            DO i=iMin,iMax
193             alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)             alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
# Line 160  C     overlying shelfice (whose density Line 196  C     overlying shelfice (whose density
196          ENDIF          ENDIF
197  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
198    
 #ifdef ALLOW_DIAGNOSTICS  
         IF ( useDiagnostics )  
      &   CALL DIAGNOSTICS_FILL(alphaRho,'RHOAnoma',k,1,2,bi,bj,myThid)  
 #endif  
   
199  #ifdef ALLOW_MOM_COMMON  #ifdef ALLOW_MOM_COMMON
200  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
201          IF (quasiHydrostatic) THEN          IF (quasiHydrostatic) THEN
202           CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)           CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
203          ENDIF          ENDIF
204  #endif /* ALLOW_MOM_COMMON */  #endif /* ALLOW_MOM_COMMON */
205    
206  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
207          IF (k.EQ.1 .AND. addSurfPhiAnom) THEN          IF ( addSurfPhiAnom .AND.
208         &       uniformFreeSurfLev .AND. k.EQ.1 ) THEN
209            DO j=jMin,jMax            DO j=jMin,jMax
210              DO i=iMin,iMax              DO i=iMin,iMax
211                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
# Line 188  C----  Hydrostatic pressure at cell cent Line 220  C----  Hydrostatic pressure at cell cent
220         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
221  C  --  Finite Volume Form  C  --  Finite Volume Form
222    
          DO j=jMin,jMax  
           DO i=iMin,iMax  
   
223  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
224  C           which has not been used to date since it does not  C           which has not been used to date since it does not
225  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
226  C  
227             phiHydC(i,j)=phiHydF(i,j)          IF ( uniformFreeSurfLev ) THEN
228       &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst           DO j=jMin,jMax
229             phiHydF(i,j)=phiHydF(i,j)            DO i=iMin,iMax
230       &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst              phiHydC(i,j) = phiHydF(i,j)
231         &          + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
232                phiHydF(i,j) = phiHydF(i,j)
233         &                 + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
234            ENDDO            ENDDO
235           ENDDO           ENDDO
236            ELSE
237             DO j=jMin,jMax
238              DO i=iMin,iMax
239               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
240                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
241    #ifdef NONLIN_FRSURF
242                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
243    #endif
244                phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst
245               ELSE
246                phiHydC(i,j) = phiHydF(i,j)
247         &            + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
248               ENDIF
249                phiHydF(i,j) = phiHydC(i,j)
250         &            + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
251              ENDDO
252             ENDDO
253            ENDIF
254    
255         ELSE         ELSE
256  C  --  Finite Difference Form  C  --  Finite Difference Form
257    
258           dRlocM=half*drC(k)  C---------- This discretization is the "energy conserving" form
259           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  
260    
261            dRlocM = halfRL*drC(k)
262            IF (k.EQ.1) dRlocM=rF(k)-rC(k)
263            IF (k.EQ.Nr) THEN
264              dRlocP=rC(k)-rF(k+1)
265            ELSE
266              dRlocP=halfRL*drC(k+1)
267            ENDIF
268            IF ( uniformFreeSurfLev ) THEN
269           DO j=jMin,jMax           DO j=jMin,jMax
270            DO i=iMin,iMax            DO i=iMin,iMax
271                phiHydC(i,j) = phiHydF(i,j)
272  C---------- This discretization is the "energy conserving" form       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
273  C           which has been used since at least Adcroft et al., MWR 1997              phiHydF(i,j) = phiHydC(i,j)
274  C       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
275              phiHydC(i,j)=phiHydF(i,j)            ENDDO
276             ENDDO
277            ELSE
278             rec_dRm = oneRL/(rF(k)-rC(k))
279             rec_dRp = oneRL/(rC(k)-rF(k+1))
280             DO j=jMin,jMax
281              DO i=iMin,iMax
282               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
283                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
284    #ifdef NONLIN_FRSURF
285                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
286    #endif
287                phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
288         &                     +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
289         &                    )*gravity*alphaRho(i,j)*recip_rhoConst
290               ELSE
291                phiHydC(i,j) = phiHydF(i,j)
292       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
293              phiHydF(i,j)=phiHydC(i,j)             ENDIF
294                phiHydF(i,j) = phiHydC(i,j)
295       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
296            ENDDO            ENDDO
297           ENDDO           ENDDO
298            ENDIF
299    
300  C  --  end if integr_GeoPot = ...  C  --  end if integr_GeoPot = ...
301         ENDIF         ENDIF
302            
303  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
304        ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN        ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
305  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
306  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density before
307  C       before integrating (1/rho)'*dp over the current layer/interface  C       integrating (1/rho)_prime*dp over the current layer/interface
308  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
309  CADJ GENERAL  CADJ GENERAL
310  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
311    
312            IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
313  C--     Calculate density  C--     Calculate density
314  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
315              kkey = (ikey-1)*Nr + k            kkey = (ikey-1)*Nr + k
316  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,
317  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
318    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte,
319    CADJ &     kind = isbyte
320  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
321          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,            CALL FIND_RHO_2D(
322       &                 tFld, sFld,       I              iMin, iMax, jMin, jMax, k,
323       &                 alphaRho, myThid)       I              tFld(1-OLx,1-OLy,k,bi,bj),
324         I              sFld(1-OLx,1-OLy,k,bi,bj),
325         O              alphaRho,
326         I              k, bi, bj, myThid )
327  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
328  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte,
329    CADJ &     kind = isbyte
330  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
331            ELSE
332              DO j=jMin,jMax
333               DO i=iMin,iMax
334                 alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
335               ENDDO
336              ENDDO
337            ENDIF
338    
339  #ifdef ALLOW_DIAGNOSTICS  C--     Calculate specific volume anomaly : alpha_prime = 1/rho - alpha_Cst
         IF ( useDiagnostics )  
      &   CALL DIAGNOSTICS_FILL(alphaRho,'RHOAnoma',k,1,2,bi,bj,myThid)  
 #endif  
   
 C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst  
340          DO j=jMin,jMax          DO j=jMin,jMax
341            DO i=iMin,iMax            DO i=iMin,iMax
342              locAlpha=alphaRho(i,j)+rhoConst              locAlpha=alphaRho(i,j)+rhoConst
343              alphaRho(i,j)=maskC(i,j,k,bi,bj)*              alphaRho(i,j)=maskC(i,j,k,bi,bj)*
344       &              (one/locAlpha - recip_rhoConst)       &              (oneRL/locAlpha - recip_rhoConst)
345            ENDDO            ENDDO
346          ENDDO          ENDDO
347    
348    #ifdef ALLOW_MOM_COMMON
349    C--     Quasi-hydrostatic terms are added as if they modify the specific-volume
350            IF (quasiHydrostatic) THEN
351             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
352            ENDIF
353    #endif /* ALLOW_MOM_COMMON */
354    
355  C----  Hydrostatic pressure at cell centers  C----  Hydrostatic pressure at cell centers
356    
357         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
# Line 276  C  --  Finite Volume Form Line 363  C  --  Finite Volume Form
363  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
364  C           which has not been used to date since it does not  C           which has not been used to date since it does not
365  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
366  C  
367             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
368               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
369  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
370               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
371  #endif  #endif
372               phiHydC(i,j) = ddRloc*alphaRho(i,j)               phiHydC(i,j) = ddRloc*alphaRho(i,j)
373  c--to reproduce results of c48d_post: uncomment those 4+1 lines  c--to reproduce results of c48d_post: uncomment those 4+1 lines
374  c            phiHydC(i,j)=phiHydF(i,j)  c            phiHydC(i,j)=phiHydF(i,j)
375  c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)  c    &          +(hFacC(i,j,k,bi,bj)-halfRL)*drF(k)*alphaRho(i,j)
376  c            phiHydF(i,j)=phiHydF(i,j)  c            phiHydF(i,j)=phiHydF(i,j)
377  c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)  c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
378             ELSE             ELSE
379               phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)               phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
380  c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)  c            phiHydF(i,j) = phiHydF(i,j) +        drF(k)*alphaRho(i,j)
381             ENDIF             ENDIF
382  c-- and comment this last one:  c-- and comment this last one:
383               phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)               phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
384  c-----  c-----
385            ENDDO            ENDDO
386           ENDDO           ENDDO
# Line 301  c----- Line 388  c-----
388         ELSE         ELSE
389  C  --  Finite Difference Form, with Part-Cell Bathy  C  --  Finite Difference Form, with Part-Cell Bathy
390    
391           dRlocM=half*drC(k)           dRlocM = halfRL*drC(k)
392           IF (k.EQ.1) dRlocM=rF(k)-rC(k)           IF (k.EQ.1) dRlocM=rF(k)-rC(k)
393           IF (k.EQ.Nr) THEN           IF (k.EQ.Nr) THEN
394             dRlocP=rC(k)-rF(k+1)             dRlocP=rC(k)-rF(k+1)
395           ELSE           ELSE
396             dRlocP=half*drC(k+1)             dRlocP=halfRL*drC(k+1)
397           ENDIF           ENDIF
398           rec_dRm = one/(rF(k)-rC(k))           rec_dRm = oneRL/(rF(k)-rC(k))
399           rec_dRp = one/(rC(k)-rF(k+1))           rec_dRp = oneRL/(rC(k)-rF(k+1))
400    
401           DO j=jMin,jMax           DO j=jMin,jMax
402            DO i=iMin,iMax            DO i=iMin,iMax
403    
404  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
405    
406             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
407               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
408  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
409               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
410  #endif  #endif
411               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM               phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
412       &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP       &                      +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
413       &                     )*alphaRho(i,j)       &                     )*alphaRho(i,j)
414             ELSE             ELSE
415               phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)               phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
# Line 343  C       the specific volume, analogous t Line 430  C       the specific volume, analogous t
430  C--     virtual potential temperature anomaly (including water vapour effect)  C--     virtual potential temperature anomaly (including water vapour effect)
431          DO j=jMin,jMax          DO j=jMin,jMax
432           DO i=iMin,iMax           DO i=iMin,iMax
433            alphaRho(i,j)=maskC(i,j,k,bi,bj)            alphaRho(i,j) = ( tFld(i,j,k,bi,bj)
434       &             *( tFld(i,j,k,bi,bj)*(sFld(i,j,k,bi,bj)*atm_Rq+one)       &                      *( sFld(i,j,k,bi,bj)*atm_Rq + oneRL )
435       &               -tRef(k) )       &                    - tRef(k) )*maskC(i,j,k,bi,bj)
436           ENDDO           ENDDO
437          ENDDO          ENDDO
438    
439    #ifdef ALLOW_MOM_COMMON
440    C--     Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
441            IF (quasiHydrostatic) THEN
442             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
443            ENDIF
444    #endif /* ALLOW_MOM_COMMON */
445    
446  C---    Integrate d Phi / d pi  C---    Integrate d Phi / d pi
447    
448    #ifndef DISABLE_SIGMA_CODE
449    C  --  Finite Volume Form, integrated to r-level (cell center & upper interface)
450           IF ( useFVgradPhi ) THEN
451    
452             fullDepth = rF(1)-rF(Nr+1)
453             DO j=jMin,jMax
454              DO i=iMin,iMax
455               locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
456    #ifdef NONLIN_FRSURF
457               locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
458    #endif
459               pKappaF(i,j) = (
460         &           ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
461         &                              + bHybSigmF( k )*locDepth
462         &           )/atm_Po )**atm_kappa
463               pKappaC      = (
464         &           ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
465         &                              + bHybSigmC( k )*locDepth
466         &           )/atm_Po )**atm_kappa
467               pKappaU(i,j) = (
468         &           ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
469         &                             + bHybSigmF(k+1)*locDepth
470         &           )/atm_Po )**atm_kappa
471               phiHydC(i,j) = phiHydF(i,j)
472         &        + atm_Cp*( pKappaF(i,j) - pKappaC      )*alphaRho(i,j)
473               phiHydU(i,j) = phiHydF(i,j)
474         &        + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
475              ENDDO
476             ENDDO
477    C end: Finite Volume Form, integrated to r-level.
478    
479           ELSEIF (integr_GeoPot.EQ.0) THEN
480    #else /* DISABLE_SIGMA_CODE */
481         IF (integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
482    #endif /* DISABLE_SIGMA_CODE */
483  C  --  Energy Conserving Form, accurate with Full cell topo  --  C  --  Energy Conserving Form, accurate with Full cell topo  --
484  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
485  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
486  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary  C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
487  C             condition is simply Phi-prime(Ro_surf)=0.  C             condition is simply Phi-prime(Ro_surf)=0.
488  C           o convention ddPI > 0 (same as drF & drC)  C           o convention ddPI > 0 (same as drF & drC)
489  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
# Line 364  C--------------------------------------- Line 492  C---------------------------------------
492       &                   -((rC( k )/atm_Po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
493           ELSE           ELSE
494             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
495       &                   -((rC( k )/atm_Po)**atm_kappa) )*half       &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
496           ENDIF           ENDIF
497           IF (k.EQ.Nr) THEN           IF (k.EQ.Nr) THEN
498             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
499       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
500           ELSE           ELSE
501             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
502       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
503           ENDIF           ENDIF
504  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
505           DO j=jMin,jMax           DO j=jMin,jMax
# Line 390  C  Finite Volume formulation consistent Line 518  C  Finite Volume formulation consistent
518  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 :
519  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)
520  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
521  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
522  C     non-linearity in PI(p)  C     non-linearity in PI(p)
523  C---------  C---------
524             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
# Line 399  C--------- Line 527  C---------
527       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
528           DO j=jMin,jMax           DO j=jMin,jMax
529            DO i=iMin,iMax            DO i=iMin,iMax
530             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
531               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
532  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
533               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
# Line 417  C--------------------------------------- Line 545  C---------------------------------------
545    
546         ELSEIF ( integr_GeoPot.EQ.2         ELSEIF ( integr_GeoPot.EQ.2
547       &     .OR. integr_GeoPot.EQ.3 ) THEN       &     .OR. integr_GeoPot.EQ.3 ) THEN
548  C  --  Finite Difference Form, with Part-Cell Topo,  C  --  Finite Difference Form, with Part-Cell Topo,
549  C       works with Interface_W  at the middle between 2.Tracer_Level  C       works with Interface_W  at the middle between 2.Tracer_Level
550  C        and  with Tracer_Level at the middle between 2.Interface_W.  C        and  with Tracer_Level at the middle between 2.Interface_W.
551  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
552  C  Finite Difference formulation consistent with Partial Cell,  C  Finite Difference formulation consistent with Partial Cell,
553  C   Valid & accurate if Interface_W at middle between tracer levels  C   Valid & accurate if Interface_W at middle between tracer levels
554  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
555  C---------  C---------
556           IF (k.EQ.1) THEN           IF (k.EQ.1) THEN
557             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
558       &                   -((rC( k )/atm_Po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
559           ELSE           ELSE
560             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
561       &                   -((rC( k )/atm_Po)**atm_kappa) )*half       &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
562           ENDIF           ENDIF
563           IF (k.EQ.Nr) THEN           IF (k.EQ.Nr) THEN
564             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
565       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
566           ELSE           ELSE
567             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
568       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
569           ENDIF           ENDIF
570           rec_dRm = one/(rF(k)-rC(k))           rec_dRm = oneRL/(rF(k)-rC(k))
571           rec_dRp = one/(rC(k)-rF(k+1))           rec_dRp = oneRL/(rC(k)-rF(k+1))
572           DO j=jMin,jMax           DO j=jMin,jMax
573            DO i=iMin,iMax            DO i=iMin,iMax
574             IF (k.EQ.ksurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
575               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
576  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
577               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
578  #endif  #endif
579               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm               phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
580       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )       &                      +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
581       &                    *alphaRho(i,j)       &                     )*alphaRho(i,j)
582             ELSE             ELSE
583               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
584             ENDIF             ENDIF
# Line 469  C---+----1----+----2----+----3----+----4 Line 597  C---+----1----+----2----+----3----+----4
597          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
598        ENDIF        ENDIF
599    
600          IF ( .NOT. useFVgradPhi ) THEN
601    C--   r-coordinate and r*-coordinate cases:
602    
603           IF ( momPressureForcing ) THEN
604            CALL CALC_GRAD_PHI_HYD(
605         I                         k, bi, bj, iMin,iMax, jMin,jMax,
606         I                         phiHydC, alphaRho, tFld, sFld,
607         O                         dPhiHydX, dPhiHydY,
608         I                         myTime, myIter, myThid)
609           ENDIF
610    
611    #ifndef DISABLE_SIGMA_CODE
612          ELSE
613    C--   else (SigmaCoords part)
614    
615           IF ( fluidIsWater ) THEN
616            STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
617           ENDIF
618           IF ( momPressureForcing ) THEN
619            CALL CALC_GRAD_PHI_FV(
620         I                         k, bi, bj, iMin,iMax, jMin,jMax,
621         I                         phiHydF, phiHydU, pKappaF, pKappaU,
622         O                         dPhiHydX, dPhiHydY,
623         I                         myTime, myIter, myThid)
624           ENDIF
625           DO j=jMin,jMax
626             DO i=iMin,iMax
627               phiHydF(i,j) = phiHydU(i,j)
628             ENDDO
629           ENDDO
630    
631    #endif /* DISABLE_SIGMA_CODE */
632    C--   end if-not/else useFVgradPhi
633          ENDIF
634    
635  C---   Diagnose Phi at boundary r=R_low :  C---   Diagnose Phi at boundary r=R_low :
636  C       = Ocean bottom pressure (Ocean, Z-coord.)  C       = Ocean bottom pressure (Ocean, Z-coord.)
637  C       = Sea-surface height    (Ocean, P-coord.)  C       = Sea-surface height    (Ocean, P-coord.)
# Line 477  C       = Top atmosphere height (Atmos, Line 640  C       = Top atmosphere height (Atmos,
640          CALL DIAGS_PHI_RLOW(          CALL DIAGS_PHI_RLOW(
641       I                      k, bi, bj, iMin,iMax, jMin,jMax,       I                      k, bi, bj, iMin,iMax, jMin,jMax,
642       I                      phiHydF, phiHydC, alphaRho, tFld, sFld,       I                      phiHydF, phiHydC, alphaRho, tFld, sFld,
643       I                      myTime, myIter, myThid)         I                      myTime, myIter, myThid)
644        ENDIF        ENDIF
645    
646  C---   Diagnose Full Hydrostatic Potential at cell center level  C---   Diagnose Full Hydrostatic Potential at cell center level
# Line 486  C---   Diagnose Full Hydrostatic Potenti Line 649  C---   Diagnose Full Hydrostatic Potenti
649       I                      phiHydC,       I                      phiHydC,
650       I                      myTime, myIter, myThid)       I                      myTime, myIter, myThid)
651    
       IF (momPressureForcing) THEN  
         CALL CALC_GRAD_PHI_HYD(  
      I                         k, bi, bj, iMin,iMax, jMin,jMax,  
      I                         phiHydC, alphaRho, tFld, sFld,  
      O                         dPhiHydX, dPhiHydY,  
      I                         myTime, myIter, myThid)    
       ENDIF  
   
652  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
653    
654        RETURN        RETURN
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22