/[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.24 by jmc, Tue Dec 10 02:55:47 2002 UTC revision 1.27 by jmc, Sun Feb 9 02:58:39 2003 UTC
# Line 10  C     !INTERFACE: Line 10  C     !INTERFACE:
10       I                         bi, bj, iMin, iMax, jMin, jMax, K,       I                         bi, bj, iMin, iMax, jMin, jMax, K,
11       I                         tFld, sFld,       I                         tFld, sFld,
12       U                         phiHyd,       U                         phiHyd,
13       I                         myThid)       O                         dPhiHydX, dPhiHydY,
14         I                         myTime, myIter, myThid)
15  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
16  C     *==========================================================*  C     *==========================================================*
17  C     | SUBROUTINE CALC_PHI_HYD                                  |  C     | SUBROUTINE CALC_PHI_HYD                                  |
# Line 59  C     == Routine arguments == Line 60  C     == Routine arguments ==
60        _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)
61        _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)
62        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63        INTEGER myThid        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
64          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
65          _RL myTime
66          INTEGER myIter, myThid
67                
68  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE  #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
69    
# Line 70  C     == Local variables == Line 74  C     == Local variables ==
74        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RL dRloc,dRlocKp1,locAlpha        _RL dRloc,dRlocKp1,locAlpha
76        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm        _RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm
77          INTEGER iMnLoc,jMnLoc
78          PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
79          LOGICAL useDiagPhiRlow
80  CEOP  CEOP
81          useDiagPhiRlow = .TRUE.
       zero = 0. _d 0  
       one  = 1. _d 0  
       half = .5 _d 0  
82    
83  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
84  C  Atmosphere:    C  Atmosphere:  
# Line 103  C---+----1----+----2----+----3----+----4 Line 107  C---+----1----+----2----+----3----+----4
107       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
108  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
109    
110    
111    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
112        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
113  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
114  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
115  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
116    #ifdef      ALLOW_AUTODIFF_TAMC
117    CADJ GENERAL
118    #endif      /* ALLOW_AUTODIFF_TAMC */
119    
120          dRloc=drC(k)          dRloc=drC(k)
121          IF (k.EQ.1) dRloc=drF(1)          IF (k.EQ.1) dRloc=drF(1)
# Line 121  C       P(z=eta) = P(atmospheric_loading Line 130  C       P(z=eta) = P(atmospheric_loading
130          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
131            DO j=jMin,jMax            DO j=jMin,jMax
132              DO i=iMin,iMax              DO i=iMin,iMax
133                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)  c             phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
134                  phiHyd(i,j,k) = 0.
135              ENDDO              ENDDO
136            ENDDO            ENDDO
137          ENDIF          ENDIF
# Line 141  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       Hydrostatic pressure at cell centers  C---   Diagnose Hydrostatic pressure at the bottom:
155          DO j=jMin,jMax         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
163    
164           IF (integr_GeoPot.EQ.1) THEN
165    C  --  Finite Volume Form
166    
167             DO j=jMin,jMax
168            DO i=iMin,iMax            DO i=iMin,iMax
 #ifdef      ALLOW_AUTODIFF_TAMC  
 c           Patrick, is this directive correct or even necessary in  
 c           this new code?  
 c           Yes, because of phiHyd(i,j,k+1)=phiHyd(i,j,k)+...  
 c           within the k-loop.  
 CADJ GENERAL  
 #endif      /* ALLOW_AUTODIFF_TAMC */  
169    
170  CmlC---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
171  CmlC           which has not been used to date since it does not  C           which has not been used to date since it does not
172  CmlC           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
173  CmlC  C
174  Cml          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN             IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
175  Cml           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)       &            + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
176  Cml     &          + hFacC(i,j,k,bi,bj)             phiHyd(i,j,k)=phiHyd(i,j,k)+
177  Cml     &            *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &       + half*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
178  Cml     &          + gravity*etaN(i,j,bi,bj)  
179  Cml          ENDIF            ENDDO
180  Cml           IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+           ENDDO
181  Cml     &         drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  
182  Cml           phiHyd(i,j,k)=phiHyd(i,j,k)+         ELSE
183  Cml     &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst  C  --  Finite Difference Form
184  CmlC-----------------------------------------------------------------------  
185             DO j=jMin,jMax
186              DO i=iMin,iMax
187    
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)
192              phiHyd(i,j,k)=phiHyd(i,j,k)+       &        +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
193       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
194              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+       &        +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
      &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst  
 C-----------------------------------------------------------------------  
   
 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)-.5)*drF(K)  
      &                   *gravity*alphaRho(i,j)*recip_rhoConst  
      &              + gravity*etaN(i,j,bi,bj)  
             ENDIF  
 C-----------------------------------------------------------------------  
195    
196            ENDDO            ENDDO
197          ENDDO           ENDDO
198    
199    C  --  end if integr_GeoPot = ...
200           ENDIF
201                    
202    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
203        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
204  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
205  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
206  C       before integrating g*rho over the current layer/interface  C       before integrating (1/rho)'*dp over the current layer/interface
207  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
208  CADJ GENERAL  CADJ GENERAL
209  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
# Line 212  CADJ GENERAL Line 219  CADJ GENERAL
219          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
220            DO j=jMin,jMax            DO j=jMin,jMax
221              DO i=iMin,iMax              DO i=iMin,iMax
222                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)  c             phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
223                  phiHyd(i,j,k) = 0.
224              ENDDO              ENDDO
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 230  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
 C       Hydrostatic pressure at cell centers  
242          DO j=jMin,jMax          DO j=jMin,jMax
243            DO i=iMin,iMax            DO i=iMin,iMax
244              locAlpha=alphaRho(i,j)+rhoConst              locAlpha=alphaRho(i,j)+rhoConst
245              IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha              alphaRho(i,j)=maskC(i,j,k,bi,bj)*
246         &              (one/locAlpha - recip_rhoConst)
247              ENDDO
248            ENDDO
249    
250  CmlC---------- This discretization is the "finite volume" form  C---   Diagnose Sea-surface height (Hydrostatic geopotential at r=Rlow):
251  CmlC           which has not been used to date since it does not         IF (useDiagPhiRlow) THEN
252  CmlC           conserve KE+PE exactly even though it is more natural            CALL DIAGS_PHI_RLOW(
253  CmlC       I                        k, bi, bj, iMin,iMax, jMin,jMax,
254  Cml            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN       I                        phiHyd, alphaRho, tFld, sFld,
255  Cml             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)       I                        myTime, myIter, myThid)  
256  Cml     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha         ENDIF
 Cml     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
 Cml            ENDIF  
 Cml            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  
 Cml     &           drF(K)*locAlpha  
 Cml            phiHyd(i,j,k)=phiHyd(i,j,k)+  
 Cml     &           0.5*drF(K)*locAlpha  
 CmlC-----------------------------------------------------------------------  
257    
258  C---------- This discretization is the "energy conserving" form  C----  Hydrostatic pressure at cell centers
259  C           which has been used since at least Adcroft et al., MWR 1997  
260           IF (integr_GeoPot.EQ.1) THEN
261    C  --  Finite Volume Form
262    
263             DO j=jMin,jMax
264              DO i=iMin,iMax
265    
266    C---------- This discretization is the "finite volume" form
267    C           which has not been used to date since it does not
268    C           conserve KE+PE exactly even though it is more natural
269  C  C
270                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
271         &          + hFacC(i,j,k,bi,bj)*drF(K)*alphaRho(i,j)
272                phiHyd(i,j,k)=phiHyd(i,j,k)
273         &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*alphaRho(i,j)
274    
275              phiHyd(i,j,k)=phiHyd(i,j,k)+            ENDDO
276       &          0.5*dRloc*locAlpha           ENDDO
             IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+  
      &          0.5*dRlocKp1*locAlpha  
277    
278  C-----------------------------------------------------------------------         ELSE
279    C  --  Finite Difference Form
280    
281  C---------- Compute gravity*(sea surface elevation) first           DO j=jMin,jMax
282  C           This has to be done starting from phiHyd at the current            DO i=iMin,iMax
283  C           tracer point and .5 of the cell's thickness has to be  
284  C           substracted from hFacC  C---------- This discretization is the "energy conserving" form
285              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  
286               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)              phiHyd(i,j,k)=phiHyd(i,j,k)
287       &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha       &          + half*dRloc*alphaRho(i,j)
288       &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)              IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
289              ENDIF       &          + half*dRlocKp1*alphaRho(i,j)
 C-----------------------------------------------------------------------  
290    
291            ENDDO            ENDDO
292          ENDDO           ENDDO
293    
294    C  --  end if integr_GeoPot = ...
295           ENDIF
296    
297        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
298  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 298  C--------------------------------------- Line 315  C---------------------------------------
315       &                  -((rC(K)/atm_Po)**atm_kappa) )       &                  -((rC(K)/atm_Po)**atm_kappa) )
316            DO j=jMin,jMax            DO j=jMin,jMax
317             DO i=iMin,iMax             DO i=iMin,iMax
318               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
319       &         +ddPIp*maskC(i,j,K,bi,bj)               phiHyd(i,j,K)=
320         &          ddPIp*maskC(i,j,K,bi,bj)
321       &               *(tFld(I,J,K,bi,bj)-tRef(K))       &               *(tFld(I,J,K,bi,bj)-tRef(K))
322             ENDDO             ENDDO
323            ENDDO            ENDDO
# Line 338  C--------- Line 356  C---------
356       &                  -((rC(K)/atm_Po)**atm_kappa) )       &                  -((rC(K)/atm_Po)**atm_kappa) )
357            DO j=jMin,jMax            DO j=jMin,jMax
358             DO i=iMin,iMax             DO i=iMin,iMax
359               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
360       &         +ddPIp*_hFacC(I,J, K ,bi,bj)               phiHyd(i,j,K)=
361         &          ddPIp*_hFacC(I,J, K ,bi,bj)
362       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
363             ENDDO             ENDDO
364            ENDDO            ENDDO
# Line 376  C--------- Line 395  C---------
395       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
396            DO j=jMin,jMax            DO j=jMin,jMax
397             DO i=iMin,iMax             DO i=iMin,iMax
398               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
399       &       +( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)               phiHyd(i,j,K)=
400         &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
401       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
402       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
403       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
# Line 421  C--------- Line 441  C---------
441       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
442            DO j=jMin,jMax            DO j=jMin,jMax
443             DO i=iMin,iMax             DO i=iMin,iMax
444               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
445       &       +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)               phiHyd(i,j,K)=
446         &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
447       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
448       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
449       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
# Line 460  C---+----1----+----2----+----3----+----4 Line 481  C---+----1----+----2----+----3----+----4
481          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
482        ENDIF        ENDIF
483    
484          IF (momPressureForcing) THEN
485            iMnLoc = MAX(1-Olx+1,iMin)
486            jMnLoc = MAX(1-Oly+1,jMin)
487            CALL CALC_GRAD_PHI_HYD(
488         I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
489         I                         phiHyd, alphaRho, tFld, sFld,
490         O                         dPhiHydX, dPhiHydY,
491         I                         myTime, myIter, myThid)  
492          ENDIF
493    
494  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
495    
496        RETURN        RETURN
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22