/[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.25 by jmc, Sat Feb 8 02:09:20 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  CEOP  CEOP
80    
       zero = 0. _d 0  
       one  = 1. _d 0  
       half = .5 _d 0  
   
81  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
82  C  Atmosphere:    C  Atmosphere:  
83  C integr_GeoPot => select one option for the integration of the Geopotential:  C integr_GeoPot => select one option for the integration of the Geopotential:
# Line 103  C---+----1----+----2----+----3----+----4 Line 105  C---+----1----+----2----+----3----+----4
105       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
106  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
107    
108    
109    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
110        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
111  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
112  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
113  C       before integrating g*rho over the current layer/interface  C       before integrating g*rho over the current layer/interface
114    #ifdef      ALLOW_AUTODIFF_TAMC
115    CADJ GENERAL
116    #endif      /* ALLOW_AUTODIFF_TAMC */
117    
118          dRloc=drC(k)          dRloc=drC(k)
119          IF (k.EQ.1) dRloc=drF(1)          IF (k.EQ.1) dRloc=drF(1)
# Line 122  C       P(z=eta) = P(atmospheric_loading Line 129  C       P(z=eta) = P(atmospheric_loading
129            DO j=jMin,jMax            DO j=jMin,jMax
130              DO i=iMin,iMax              DO i=iMin,iMax
131                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
132    c             phiHyd(i,j,k) = 0.
133              ENDDO              ENDDO
134            ENDDO            ENDDO
135          ENDIF          ENDIF
# Line 141  C Quasi-hydrostatic terms are added in a Line 149  C Quasi-hydrostatic terms are added in a
149           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)           CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid)
150          ENDIF          ENDIF
151    
152  C       Hydrostatic pressure at cell centers  C---   Hydrostatic pressure at cell centers
153          DO j=jMin,jMax  
154           IF (integr_GeoPot.EQ.1) THEN
155    C  --  Finite Volume Form
156    
157             DO j=jMin,jMax
158            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 */  
159    
160  CmlC---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
161  CmlC           which has not been used to date since it does not  C           which has not been used to date since it does not
162  CmlC           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
163  CmlC  C
164  Cml          IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
165  Cml           phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
166  Cml     &          + hFacC(i,j,k,bi,bj)       &          + hFacC(i,j,k,bi,bj)
167  Cml     &            *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &             *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
168  Cml     &          + gravity*etaN(i,j,bi,bj)       &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
169  Cml          ENDIF            ENDIF
170  Cml           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  Cml     &         drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &            + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
172  Cml           phiHyd(i,j,k)=phiHyd(i,j,k)+             phiHyd(i,j,k)=phiHyd(i,j,k)+
173  Cml     &          0.5*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst       &       + half*drF(K)*gravity*alphaRho(i,j)*recip_rhoConst
174  CmlC-----------------------------------------------------------------------  
175              ENDDO
176             ENDDO
177    
178           ELSE
179    C  --  Finite Difference Form
180    
181             DO j=jMin,jMax
182              DO i=iMin,iMax
183    
184  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
185  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
186  C  C
187                            
188              phiHyd(i,j,k)=phiHyd(i,j,k)+              phiHyd(i,j,k)=phiHyd(i,j,k)
189       &          0.5*dRloc*gravity*alphaRho(i,j)*recip_rhoConst       &        +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst
190              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)
191       &          0.5*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst       &        +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst
 C-----------------------------------------------------------------------  
192    
193  C---------- Compute bottom pressure deviation from gravity*rho0*H  C---------- Compute bottom pressure deviation from gravity*rho0*H
194  C           This has to be done starting from phiHyd at the current  C           This has to be done starting from phiHyd at the current
# Line 184  C           tracer point and .5 of the c Line 196  C           tracer point and .5 of the c
196  C           substracted from hFacC  C           substracted from hFacC
197              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
198               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
199       &              + (hFacC(i,j,k,bi,bj)-.5)*drF(K)       &              + (hFacC(i,j,k,bi,bj)-half)*drF(K)
200       &                   *gravity*alphaRho(i,j)*recip_rhoConst       &                   *gravity*alphaRho(i,j)*recip_rhoConst
201       &              + gravity*etaN(i,j,bi,bj)       &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
202              ENDIF              ENDIF
 C-----------------------------------------------------------------------  
203    
204            ENDDO            ENDDO
205          ENDDO           ENDDO
206    
207    C  --  end if integr_GeoPot = ...
208           ENDIF
209                    
210    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN
212  C       This is the hydrostatic pressure calculation for the Ocean  C       This is the hydrostatic pressure calculation for the Ocean
213  C       which uses the FIND_RHO() routine to calculate density  C       which uses the FIND_RHO() routine to calculate density
214  C       before integrating g*rho over the current layer/interface  C       before integrating (1/rho)'*dp over the current layer/interface
215  #ifdef      ALLOW_AUTODIFF_TAMC  #ifdef      ALLOW_AUTODIFF_TAMC
216  CADJ GENERAL  CADJ GENERAL
217  #endif      /* ALLOW_AUTODIFF_TAMC */  #endif      /* ALLOW_AUTODIFF_TAMC */
# Line 213  CADJ GENERAL Line 228  CADJ GENERAL
228            DO j=jMin,jMax            DO j=jMin,jMax
229              DO i=iMin,iMax              DO i=iMin,iMax
230                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)                phiHyd(i,j,k) = phi0surf(i,j,bi,bj)
231    c             phiHyd(i,j,k) = 0.
232              ENDDO              ENDDO
233            ENDDO            ENDDO
234          ENDIF          ENDIF
# Line 231  CADJ STORE alphaRho (:,:) = comlev1_bibj Line 247  CADJ STORE alphaRho (:,:) = comlev1_bibj
247  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
248    
249    
250  C       Hydrostatic pressure at cell centers  C----  Hydrostatic pressure at cell centers
251          DO j=jMin,jMax  
252           IF (integr_GeoPot.EQ.1) THEN
253    C  --  Finite Volume Form
254    
255             DO j=jMin,jMax
256            DO i=iMin,iMax            DO i=iMin,iMax
257              locAlpha=alphaRho(i,j)+rhoConst              locAlpha=alphaRho(i,j)+rhoConst
258              IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha              locAlpha=maskC(i,j,k,bi,bj)*
259         &              (one/locAlpha - recip_rhoConst)
260    c           IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
261    
262    C---------- This discretization is the "finite volume" form
263    C           which has not been used to date since it does not
264    C           conserve KE+PE exactly even though it is more natural
265    C
266                IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
267                 phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
268         &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
269         &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
270                ENDIF
271                IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)
272         &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha
273                phiHyd(i,j,k)=phiHyd(i,j,k)
274         &          +(hFacC(i,j,k,bi,bj)-half)*drF(K)*locAlpha
275    
276  CmlC---------- This discretization is the "finite volume" form            ENDDO
277  CmlC           which has not been used to date since it does not           ENDDO
278  CmlC           conserve KE+PE exactly even though it is more natural  
279  CmlC         ELSE
280  Cml            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  C  --  Finite Difference Form
281  Cml             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)  
282  Cml     &          + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha           DO j=jMin,jMax
283  Cml     &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)            DO i=iMin,iMax
284  Cml            ENDIF              locAlpha=alphaRho(i,j)+rhoConst
285  Cml            IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+              locAlpha=maskC(i,j,k,bi,bj)*
286  Cml     &           drF(K)*locAlpha       &              (one/locAlpha - recip_rhoConst)
287  Cml            phiHyd(i,j,k)=phiHyd(i,j,k)+  c           IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha
 Cml     &           0.5*drF(K)*locAlpha  
 CmlC-----------------------------------------------------------------------  
288    
289  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
 C           which has been used since at least Adcroft et al., MWR 1997  
 C  
290    
291              phiHyd(i,j,k)=phiHyd(i,j,k)+              phiHyd(i,j,k)=phiHyd(i,j,k)
292       &          0.5*dRloc*locAlpha       &          + half*dRloc*locAlpha
293              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)
294       &          0.5*dRlocKp1*locAlpha       &          + half*dRlocKp1*locAlpha
295    
 C-----------------------------------------------------------------------  
296    
297  C---------- Compute gravity*(sea surface elevation) first  C---------- Compute gravity*(sea surface elevation) first
298  C           This has to be done starting from phiHyd at the current  C           This has to be done starting from phiHyd at the current
# Line 269  C           tracer point and .5 of the c Line 300  C           tracer point and .5 of the c
300  C           substracted from hFacC  C           substracted from hFacC
301              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN              IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN
302               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)               phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)
303       &              + (hFacC(i,j,k,bi,bj)-0.5)*drF(k)*locAlpha       &              + (hFacC(i,j,k,bi,bj)-half)*drF(k)*locAlpha
304       &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)       &              + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
305              ENDIF              ENDIF
 C-----------------------------------------------------------------------  
306    
307            ENDDO            ENDDO
308          ENDDO           ENDDO
309    
310    C  --  end if integr_GeoPot = ...
311           ENDIF
312    
313        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN        ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
314  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 298  C--------------------------------------- Line 331  C---------------------------------------
331       &                  -((rC(K)/atm_Po)**atm_kappa) )       &                  -((rC(K)/atm_Po)**atm_kappa) )
332            DO j=jMin,jMax            DO j=jMin,jMax
333             DO i=iMin,iMax             DO i=iMin,iMax
334               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)=
335       &         +ddPIp*maskC(i,j,K,bi,bj)               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
336         &          ddPIp*maskC(i,j,K,bi,bj)
337       &               *(tFld(I,J,K,bi,bj)-tRef(K))       &               *(tFld(I,J,K,bi,bj)-tRef(K))
338             ENDDO             ENDDO
339            ENDDO            ENDDO
# Line 338  C--------- Line 372  C---------
372       &                  -((rC(K)/atm_Po)**atm_kappa) )       &                  -((rC(K)/atm_Po)**atm_kappa) )
373            DO j=jMin,jMax            DO j=jMin,jMax
374             DO i=iMin,iMax             DO i=iMin,iMax
375               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)=
376       &         +ddPIp*_hFacC(I,J, K ,bi,bj)               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
377         &          ddPIp*_hFacC(I,J, K ,bi,bj)
378       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))       &               *(tFld(I,J, K ,bi,bj)-tRef( K ))
379             ENDDO             ENDDO
380            ENDDO            ENDDO
# Line 376  C--------- Line 411  C---------
411       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
412            DO j=jMin,jMax            DO j=jMin,jMax
413             DO i=iMin,iMax             DO i=iMin,iMax
414               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)=
415       &       +( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
416         &        ( ddPIm*max(zero, _hFacC(i,j,K,bi,bj)-half)
417       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)-half) )
418       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
419       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
# Line 421  C--------- Line 457  C---------
457       &                  -((rC(Kp1)/atm_Po)**atm_kappa) )         &                  -((rC(Kp1)/atm_Po)**atm_kappa) )  
458            DO j=jMin,jMax            DO j=jMin,jMax
459             DO i=iMin,iMax             DO i=iMin,iMax
460               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)  c            phiHyd(i,j,K)=
461       &       +( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)               phiHyd(i,j,K)= phi0surf(i,j,bi,bj)+
462         &        ( ddPIm*max(zero,(_hFacC(i,j,K,bi,bj)-one)*ratioRm+half)
463       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )       &         +ddPIp*min(zero, _hFacC(i,j,K,bi,bj)*ratioRp     -half) )
464       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))       &               *(tFld(i,j, K ,bi,bj)-tRef( K ))
465       &               * maskC(i,j, K ,bi,bj)       &               * maskC(i,j, K ,bi,bj)
# Line 460  C---+----1----+----2----+----3----+----4 Line 497  C---+----1----+----2----+----3----+----4
497          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
498        ENDIF        ENDIF
499    
500          IF (momPressureForcing) THEN
501            iMnLoc = MAX(1-Olx+1,iMin)
502            jMnLoc = MAX(1-Oly+1,jMin)
503            CALL CALC_GRAD_PHI_HYD(
504         I                         k, bi, bj, iMnLoc,iMax, jMnLoc,jMax,
505         I                         phiHyd, alphaRho, tFld, sFld,
506         O                         dPhiHydX, dPhiHydY,
507         I                         myTime, myIter, myThid)  
508          ENDIF
509    
510  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
511    
512        RETURN        RETURN
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22