/[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.19 by adcroft, Thu Aug 15 17:25:31 2002 UTC revision 1.20 by mlosch, Wed Sep 18 16:38:01 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 51  C     == Global variables == Line 51  C     == Global variables ==
51  #include "tamc_keys.h"  #include "tamc_keys.h"
52  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
53  #include "SURFACE.h"  #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 132  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       Hydrostatic pressure at cell centers  C       Hydrostatic pressure at cell centers
# Line 163  C--------------------------------------- Line 164  C---------------------------------------
164  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
165  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
166  C  C
167                
168              phiHyd(i,j,k)=phiHyd(i,j,k)+              phiHyd(i,j,k)=phiHyd(i,j,k)+
169       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
170              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)+
171       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
172  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
173    
174    C---------- Compute bottom pressure deviation from gravity*rho0*H
175    C           This has to be done starting from phiHyd at the current
176    C           tracer point and .5 of the cell's thickness has to be
177    C           substracted from hFacC
178                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
179                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
180         &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(K)
181         &                   *gravity*alphaRho(i,j)*recip_rhoConst
182         &              + gravity*etaN(i,j,bi,bj)
183                ENDIF
184    C-----------------------------------------------------------------------
185    
186            ENDDO            ENDDO
187          ENDDO          ENDDO
188                    
# Line 199  c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSu Line 214  c    &  -(Ro_surf(i,j,bi,bj)-.5*drF( kSu
214  C       Calculate density  C       Calculate density
215  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
216              kkey = (ikey-1)*Nr + k              kkey = (ikey-1)*Nr + k
217  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
218  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
219  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
220          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,          CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k,
221       &                 theta, salt,       &                 tFld, sFld,
222       &                 alphaRho, myThid)       &                 alphaRho, myThid)
223    
224  C       Hydrostatic pressure at cell centers  C       Hydrostatic pressure at cell centers
# Line 230  C Line 245  C
245              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)+
246       &          0.5*dRlocKp1*locAlpha       &          0.5*dRlocKp1*locAlpha
247  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
248    
249    C---------- Compute gravity*(sea surface elevation)
250    C           This has to be done starting from phiHyd at the current
251    C           tracer point and .5 of the cell's thickness has to be
252    C           substracted from hFacC
253                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
254                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
255         &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha
256         &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
257                ENDIF
258    C-----------------------------------------------------------------------
259    
260            ENDDO            ENDDO
261          ENDDO          ENDDO
262    
# Line 256  C--------------------------------------- Line 283  C---------------------------------------
283             DO i=iMin,iMax             DO i=iMin,iMax
284               phiHyd(i,j,K)=               phiHyd(i,j,K)=
285       &          ddPIp*maskC(i,j,K,bi,bj)       &          ddPIp*maskC(i,j,K,bi,bj)
286       &               *(theta(I,J,K,bi,bj)-tRef(K))       &               *(tFld(I,J,K,bi,bj)-tRef(K))
287             ENDDO             ENDDO
288            ENDDO            ENDDO
289          ELSE          ELSE
# Line 267  C-------- This discretization is the ene Line 294  C-------- This discretization is the ene
294             DO i=iMin,iMax             DO i=iMin,iMax
295                phiHyd(i,j,K)=phiHyd(i,j,K-1)                phiHyd(i,j,K)=phiHyd(i,j,K-1)
296       &           +ddPI*maskC(i,j,K-1,bi,bj)       &           +ddPI*maskC(i,j,K-1,bi,bj)
297       &                *(theta(I,J,K-1,bi,bj)-tRef(K-1))       &                *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
298       &           +ddPI*maskC(i,j, K ,bi,bj)       &           +ddPI*maskC(i,j, K ,bi,bj)
299       &                *(theta(I,J, K ,bi,bj)-tRef( K ))       &                *(tFld(I,J, K ,bi,bj)-tRef( K ))
300  C             Old code (atmos-exact) looked like this  C             Old code (atmos-exact) looked like this
301  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*  Cold          phiHyd(i,j,K)=phiHyd(i,j,K-1) - ddPI*
302  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))
303             ENDDO             ENDDO
304            ENDDO            ENDDO
305          ENDIF          ENDIF
# Line 296  C--------- Line 323  C---------
323             DO i=iMin,iMax             DO i=iMin,iMax
324               phiHyd(i,j,K) =               phiHyd(i,j,K) =
325       &          ddPIp*_hFacC(I,J, K ,bi,bj)       &          ddPIp*_hFacC(I,J, K ,bi,bj)
326       &               *(theta(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
327             ENDDO             ENDDO
328            ENDDO            ENDDO
329          ELSE          ELSE
# Line 308  C--------- Line 335  C---------
335             DO i=iMin,iMax             DO i=iMin,iMax
336               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
337       &         +ddPIm*_hFacC(I,J,K-1,bi,bj)       &         +ddPIm*_hFacC(I,J,K-1,bi,bj)
338       &               *(theta(I,J,K-1,bi,bj)-tRef(K-1))       &               *(tFld(I,J,K-1,bi,bj)-tRef(K-1))
339       &         +ddPIp*_hFacC(I,J, K ,bi,bj)       &         +ddPIp*_hFacC(I,J, K ,bi,bj)
340       &               *(theta(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
341             ENDDO             ENDDO
342            ENDDO            ENDDO
343          ENDIF          ENDIF
# Line 335  C--------- Line 362  C---------
362               phiHyd(i,j,K) =               phiHyd(i,j,K) =
363       &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)       &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
364       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
365       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
366       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
367             ENDDO             ENDDO
368            ENDDO            ENDDO
# Line 348  C--------- Line 375  C---------
375             DO i=iMin,iMax             DO i=iMin,iMax
376               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
377       &        + ddPIm*0.5       &        + ddPIm*0.5
378       &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))       &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
379       &               * maskC(i,j,K-1,bi,bj)       &               * maskC(i,j,K-1,bi,bj)
380       &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)       &        +(ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
381       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
382       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
383       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
384             ENDDO             ENDDO
385            ENDDO            ENDDO
# Line 380  C--------- Line 407  C---------
407               phiHyd(i,j,K) =               phiHyd(i,j,K) =
408       &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)       &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
409       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
410       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
411       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
412             ENDDO             ENDDO
413            ENDDO            ENDDO
# Line 395  C--------- Line 422  C---------
422             DO i=iMin,iMax             DO i=iMin,iMax
423               phiHyd(i,j,K) = phiHyd(i,j,K-1)               phiHyd(i,j,K) = phiHyd(i,j,K-1)
424       &        + ddPIm*0.5       &        + ddPIm*0.5
425       &               *(theta(i,j,K-1,bi,bj)-tRef(K-1))       &               *(tFld(i,j,K-1,bi,bj)-tRef(K-1))
426       &               * maskC(i,j,K-1,bi,bj)       &               * maskC(i,j,K-1,bi,bj)
427       &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)       &        +(ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
428       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
429       &               *(theta(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
430       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
431             ENDDO             ENDDO
432            ENDDO            ENDDO
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22