/[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.18 by mlosch, Wed Jul 31 16:38:30 2002 UTC revision 1.22 by adcroft, Thu Nov 7 21:51:15 2002 UTC
# Line 8  C     !ROUTINE: CALC_PHI_HYD Line 8  C     !ROUTINE: CALC_PHI_HYD
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE CALC_PHI_HYD(        SUBROUTINE CALC_PHI_HYD(
10       I                         bi, bj, iMin, iMax, jMin, jMax, K,       I                         bi, bj, iMin, iMax, jMin, jMax, K,
11       I                         theta, salt,       I                         tFld, sFld,
12       U                         phiHyd,       U                         phiHyd,
13       I                         myThid)       I                         myThid)
14  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
# Line 18  C     | o Integrate the hydrostatic rela Line 18  C     | o Integrate the hydrostatic rela
18  C     *==========================================================*  C     *==========================================================*
19  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|  C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)|
20  C     | On entry:                                                |  C     | On entry:                                                |
21  C     |   theta,salt    are the current thermodynamics quantities|  C     |   tFld,sFld     are the current thermodynamics quantities|
22  C     |                 (unchanged on exit)                      |  C     |                 (unchanged on exit)                      |
23  C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |  C     |   phiHyd(i,j,1:k-1) is the hydrostatic Potential         |
24  C     |                 at cell centers (tracer points)          |  C     |                 at cell centers (tracer points)          |
# Line 50  C     == Global variables == Line 50  C     == Global variables ==
50  #include "tamc.h"  #include "tamc.h"
51  #include "tamc_keys.h"  #include "tamc_keys.h"
52  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
53    #include "SURFACE.h"
54    #include "DYNVARS.h"
55    
56  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
57  C     == Routine arguments ==  C     == Routine arguments ==
58        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax,K
59        _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
60        _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
61        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
62        INTEGER myThid        INTEGER myThid
63                
# Line 66  C     == Local variables == Line 68  C     == Local variables ==
68        INTEGER i,j, Kp1        INTEGER i,j, Kp1
69        _RL zero, one, half        _RL zero, one, half
70        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RL dRloc,dRlocKp1        _RL dRloc,dRlocKp1,locAlpha
72        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
73  CEOP  CEOP
74    
# Line 131  C       P(z=eta) = P(atmospheric_loading Line 133  C       P(z=eta) = P(atmospheric_loading
133  C       Calculate density  C       Calculate density
134  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
135              kkey = (ikey-1)*Nr + k              kkey = (ikey-1)*Nr + k
136  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
137  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
138  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
139          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
140       &                 theta, salt,       &                 tFld, sFld,
141       &                 alphaRho, myThid)       &                 alphaRho, myThid)
142    
143    C Quasi-hydrostatic terms are added in as if they modify the buoyancy
144            IF (quasiHydrostatic) THEN
145             CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
146            ENDIF
147    
148  C       Hydrostatic pressure at cell centers  C       Hydrostatic pressure at cell centers
149          DO j=jMin,jMax          DO j=jMin,jMax
150            DO i=iMin,iMax            DO i=iMin,iMax
# Line 149  c           within the k-loop. Line 156  c           within the k-loop.
156  CADJ GENERAL  CADJ GENERAL
157  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
158    
159  C---------- This discretization is the "finite volume" form  CmlC---------- This discretization is the "finite volume" form
160  C           which has not been used to date since it does not  CmlC           which has not been used to date since it does not
161  C           conserve KE+PE exactly even though it is more natural  CmlC           conserve KE+PE exactly even though it is more natural
162  C  CmlC
163  c           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  Cml          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
164  c    &              drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  Cml           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
165  c           phiHyd(i,j,k)=phiHyd(i,j,k)+  Cml     &          + hFacC(i,j,k,bi,bj)
166  c    &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  Cml     &            *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
167  C-----------------------------------------------------------------------  Cml     &          + gravity*etaN(i,j,bi,bj)
168    Cml          ENDIF
169    Cml           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
170    Cml     &         drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
171    Cml           phiHyd(i,j,k)=phiHyd(i,j,k)+
172    Cml     &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
173    CmlC-----------------------------------------------------------------------
174    
175  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
176  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
177  C  C
178                
179              phiHyd(i,j,k)=phiHyd(i,j,k)+              phiHyd(i,j,k)=phiHyd(i,j,k)+
180       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
181              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
182       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
183  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
184    
185    C---------- Compute bottom pressure deviation from gravity*rho0*H
186    C           This has to be done starting from phiHyd at the current
187    C           tracer point and .5 of the cell's thickness has to be
188    C           substracted from hFacC
189                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
190                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
191         &              + (hFacC(i,j,k,bi,bj)-.5)*drF(K)
192         &                   *gravity*alphaRho(i,j)*recip_rhoConst
193         &              + gravity*etaN(i,j,bi,bj)
194                ENDIF
195    C-----------------------------------------------------------------------
196    
197            ENDDO            ENDDO
198          ENDDO          ENDDO
199                    
200          ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
201    C       This is the hydrostatic pressure calculation for the Ocean
202    C       which uses the FIND_RHO() routine to calculate density
203    C       before integrating g*rho over the current layer/interface
204    #ifdef      ALLOW_AUTODIFF_TAMC
205    CADJ GENERAL
206    #endif      /* ALLOW_AUTODIFF_TAMC */
207    
208            dRloc=drC(k)
209            IF (k.EQ.1) dRloc=drF(1)
210            IF (k.EQ.Nr) THEN
211              dRlocKp1=0.
212            ELSE
213              dRlocKp1=drC(k+1)
214            ENDIF
215    
216            IF (k.EQ.1) THEN
217              DO j=jMin,jMax
218                DO i=iMin,iMax
219                  phiHyd(i,j,k)=0.
220                  phiHyd(i,j,k)=pload(i,j,bi,bj)
221                ENDDO
222              ENDDO
223            ENDIF
224    
225    C       Calculate density
226    #ifdef ALLOW_AUTODIFF_TAMC
227                kkey = (ikey-1)*Nr + k
228    CADJ STORE tFld(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
229    CADJ STORE sFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
230    #endif /* ALLOW_AUTODIFF_TAMC */
231            CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
232         &                 tFld, sFld,
233         &                 alphaRho, myThid)
234    
235    C       Hydrostatic pressure at cell centers
236            DO j=jMin,jMax
237              DO i=iMin,iMax
238                locAlpha=alphaRho(i,j)+rhoConst
239                IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
240    
241    CmlC---------- This discretization is the "finite volume" form
242    CmlC           which has not been used to date since it does not
243    CmlC           conserve KE+PE exactly even though it is more natural
244    CmlC
245    Cml            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
246    Cml             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
247    Cml     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
248    Cml     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
249    Cml            ENDIF
250    Cml            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
251    Cml     &           drF(K)*locAlpha
252    Cml            phiHyd(i,j,k)=phiHyd(i,j,k)+
253    Cml     &           0.5*drF(K)*locAlpha
254    CmlC-----------------------------------------------------------------------
255    
256    C---------- This discretization is the "energy conserving" form
257    C           which has been used since at least Adcroft et al., MWR 1997
258    C
259    
260                phiHyd(i,j,k)=phiHyd(i,j,k)+
261         &          0.5*dRloc*locAlpha
262                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+
263         &          0.5*dRlocKp1*locAlpha
264    
265    C-----------------------------------------------------------------------
266    
267    C---------- Compute gravity*(sea surface elevation) first
268    C           This has to be done starting from phiHyd at the current
269    C           tracer point and .5 of the cell's thickness has to be
270    C           substracted from hFacC
271                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
272                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
273         &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha
274         &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
275                ENDIF
276    C-----------------------------------------------------------------------
277    
278              ENDDO
279            ENDDO
280    
281        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
282  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 195  C--------------------------------------- Line 301  C---------------------------------------
301             DO i=iMin,iMax             DO i=iMin,iMax
302               phiHyd(i,j,K)=               phiHyd(i,j,K)=
303       &          ddPIp*maskC(i,j,K,bi,bj)       &          ddPIp*maskC(i,j,K,bi,bj)
304       &               *(theta(I,J,K,bi,bj)-tRef(K))       &               *(tFld(I,J,K,bi,bj)-tRef(K))
305             ENDDO             ENDDO
306            ENDDO            ENDDO
307          ELSE          ELSE
# Line 206  C-------- This discretization is the ene Line 312  C-------- This discretization is the ene
312             DO i=iMin,iMax             DO i=iMin,iMax
313                phiHyd(i,j,K)=phiHyd(i,j,K-1)                phiHyd(i,j,K)=phiHyd(i,j,K-1)
314       &           +ddPI*maskC(i,j,K-1,bi,bj)       &           +ddPI*maskC(i,j,K-1,bi,bj)
315       &                *(theta(I,J,K-1,bi,bj)-tRef(K-1))       &                *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
316       &           +ddPI*maskC(i,j, K ,bi,bj)       &           +ddPI*maskC(i,j, K ,bi,bj)
317       &                *(theta(I,J, K ,bi,bj)-tRef( K ))       &                *(tFld(I,J, K ,bi,bj)-tRef( K ))
318  C             Old code (atmos-exact) looked like this  C             Old code (atmos-exact) looked like this
319  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
320  Cold &      (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj)-2.*tRef(K))  Cold &      (tFld(I,J,K-1,bi,bj)+tFld(I,J,K,bi,bj)-2.*tRef(K))
321             ENDDO             ENDDO
322            ENDDO            ENDDO
323          ENDIF          ENDIF
# Line 235  C--------- Line 341  C---------
341             DO i=iMin,iMax             DO i=iMin,iMax
342               phiHyd(i,j,K) =               phiHyd(i,j,K) =
343       &          ddPIp*_hFacC(I,J, K ,bi,bj)       &          ddPIp*_hFacC(I,J, K ,bi,bj)
344       &               *(theta(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
345             ENDDO             ENDDO
346            ENDDO            ENDDO
347          ELSE          ELSE
# Line 247  C--------- Line 353  C---------
353             DO i=iMin,iMax             DO i=iMin,iMax
354               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
355       &         +ddPIm*_hFacC(I,J,K-1,bi,bj)       &         +ddPIm*_hFacC(I,J,K-1,bi,bj)
356       &               *(theta(I,J,K-1,bi,bj)-tRef(K-1))       &               *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
357       &         +ddPIp*_hFacC(I,J, K ,bi,bj)       &         +ddPIp*_hFacC(I,J, K ,bi,bj)
358       &               *(theta(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
359             ENDDO             ENDDO
360            ENDDO            ENDDO
361          ENDIF          ENDIF
# Line 274  C--------- Line 380  C---------
380               phiHyd(i,j,K) =               phiHyd(i,j,K) =
381       &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)       &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
382       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
383       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
384       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
385             ENDDO             ENDDO
386            ENDDO            ENDDO
# Line 287  C--------- Line 393  C---------
393             DO i=iMin,iMax             DO i=iMin,iMax
394               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
395       &        + ddPIm*0.5       &        + ddPIm*0.5
396       &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))       &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
397       &               * maskC(i,j,K-1,bi,bj)       &               * maskC(i,j,K-1,bi,bj)
398       &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)       &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
399       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
400       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
401       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
402             ENDDO             ENDDO
403            ENDDO            ENDDO
# Line 319  C--------- Line 425  C---------
425               phiHyd(i,j,K) =               phiHyd(i,j,K) =
426       &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)       &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
427       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
428       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
429       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
430             ENDDO             ENDDO
431            ENDDO            ENDDO
# Line 334  C--------- Line 440  C---------
440             DO i=iMin,iMax             DO i=iMin,iMax
441               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
442       &        + ddPIm*0.5       &        + ddPIm*0.5
443       &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))       &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
444       &               * maskC(i,j,K-1,bi,bj)       &               * maskC(i,j,K-1,bi,bj)
445       &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)       &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
446       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
447       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
448       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
449             ENDDO             ENDDO
450            ENDDO            ENDDO
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22