/[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.26 by jmc, Sun Feb 9 02:00:50 2003 UTC revision 1.27 by jmc, Sun Feb 9 02:58:39 2003 UTC
# Line 76  C     == Local variables == Line 76  C     == Local variables ==
76        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
77        INTEGER iMnLoc,jMnLoc        INTEGER iMnLoc,jMnLoc
78        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 )
79          LOGICAL useDiagPhiRlow
80  CEOP  CEOP
81          useDiagPhiRlow = .TRUE.
82    
83  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
84  C  Atmosphere:    C  Atmosphere:  
# Line 149  C Quasi-hydrostatic terms are added in a Line 151  C Quasi-hydrostatic terms are added in a
151           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
152          ENDIF          ENDIF
153    
154    C---   Diagnose Hydrostatic pressure at the bottom:
155           IF (useDiagPhiRlow) THEN
156              CALL DIAGS_PHI_RLOW(
157         I                        k, bi, bj, iMin,iMax, jMin,jMax,
158         I                        phiHyd, alphaRho, tFld, sFld,
159         I                        myTime, myIter, myThid)  
160           ENDIF
161    
162  C---   Hydrostatic pressure at cell centers  C---   Hydrostatic pressure at cell centers
163    
164         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
# Line 161  C---------- This discretization is the " Line 171  C---------- This discretization is the "
171  C           which has not been used to date since it does not  C           which has not been used to date since it does not
172  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
173  C  C
           IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  
            phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)  
      &          + hFacC(i,j,k,bi,bj)  
      &             *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  
      &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
           ENDIF  
174             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)
175       &            + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &            + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
176             phiHyd(i,j,k)=phiHyd(i,j,k)+             phiHyd(i,j,k)=phiHyd(i,j,k)+
# Line 184  C  --  Finite Difference Form Line 188  C  --  Finite Difference Form
188  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
189  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
190  C  C
               
191              phiHyd(i,j,k)=phiHyd(i,j,k)              phiHyd(i,j,k)=phiHyd(i,j,k)
192       &        +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst       &        +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
193              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)
194       &        +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst       &        +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
195    
 C---------- Compute bottom pressure deviation from gravity*rho0*H  
 C           This has to be done starting from phiHyd at the current  
 C           tracer point and .5 of the cell's thickness has to be  
 C           substracted from hFacC  
             IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  
              phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)  
      &              + (hFacC(i,j,k,bi,bj)-half)*drF(K)  
      &                   *gravity*alphaRho(i,j)*recip_rhoConst  
      &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
             ENDIF  
   
196            ENDDO            ENDDO
197           ENDDO           ENDDO
198    
# Line 233  c             phiHyd(i,j,k) = phi0surf(i Line 225  c             phiHyd(i,j,k) = phi0surf(i
225            ENDDO            ENDDO
226          ENDIF          ENDIF
227    
228  C       Calculate density  C--     Calculate density
229  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
230              kkey = (ikey-1)*Nr + k              kkey = (ikey-1)*Nr + k
231  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 246  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_ Line 238  CADJ STORE sFld (:,:,k,bi,bj) = comlev1_
238  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
239  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
240    
241    C--     Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst
242            DO j=jMin,jMax
243              DO i=iMin,iMax
244                locAlpha=alphaRho(i,j)+rhoConst
245                alphaRho(i,j)=maskC(i,j,k,bi,bj)*
246         &              (one/locAlpha - recip_rhoConst)
247              ENDDO
248            ENDDO
249    
250    C---   Diagnose Sea-surface height (Hydrostatic geopotential at r=Rlow):
251           IF (useDiagPhiRlow) THEN
252              CALL DIAGS_PHI_RLOW(
253         I                        k, bi, bj, iMin,iMax, jMin,jMax,
254         I                        phiHyd, alphaRho, tFld, sFld,
255         I                        myTime, myIter, myThid)  
256           ENDIF
257    
258  C----  Hydrostatic pressure at cell centers  C----  Hydrostatic pressure at cell centers
259    
# Line 254  C  --  Finite Volume Form Line 262  C  --  Finite Volume Form
262    
263           DO j=jMin,jMax           DO j=jMin,jMax
264            DO i=iMin,iMax            DO i=iMin,iMax
             locAlpha=alphaRho(i,j)+rhoConst  
             locAlpha=maskC(i,j,k,bi,bj)*  
      &              (one/locAlpha - recip_rhoConst)  
 c           IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha  
265    
266  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
267  C           which has not been used to date since it does not  C           which has not been used to date since it does not
268  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
269  C  C
             IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  
              phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)  
      &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha  
      &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
             ENDIF  
270              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)
271       &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha       &          + hFacC(i,j,k,bi,bj)*drF(K)*alphaRho(i,j)
272              phiHyd(i,j,k)=phiHyd(i,j,k)              phiHyd(i,j,k)=phiHyd(i,j,k)
273       &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*locAlpha       &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*alphaRho(i,j)
274    
275            ENDDO            ENDDO
276           ENDDO           ENDDO
# Line 281  C  --  Finite Difference Form Line 280  C  --  Finite Difference Form
280    
281           DO j=jMin,jMax           DO j=jMin,jMax
282            DO i=iMin,iMax            DO i=iMin,iMax
             locAlpha=alphaRho(i,j)+rhoConst  
             locAlpha=maskC(i,j,k,bi,bj)*  
      &              (one/locAlpha - recip_rhoConst)  
 c           IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha  
283    
284  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
285    
286              phiHyd(i,j,k)=phiHyd(i,j,k)              phiHyd(i,j,k)=phiHyd(i,j,k)
287       &          + half*dRloc*locAlpha       &          + half*dRloc*alphaRho(i,j)
288              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)
289       &          + half*dRlocKp1*locAlpha       &          + half*dRlocKp1*alphaRho(i,j)
   
   
 C---------- Compute gravity*(sea surface elevation) first  
 C           This has to be done starting from phiHyd at the current  
 C           tracer point and .5 of the cell's thickness has to be  
 C           substracted from hFacC  
             IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  
              phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)  
      &              + (hFacC(i,j,k,bi,bj)-half)*drF(k)*locAlpha  
      &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
             ENDIF  
290    
291            ENDDO            ENDDO
292           ENDDO           ENDDO
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22