/[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.40 by jmc, Tue Mar 16 00:08:27 2010 UTC revision 1.41 by jmc, Wed Apr 11 04:02:05 2012 UTC
# Line 68  C     myThid     :: thread number for th Line 68  C     myThid     :: thread number for th
68  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
69        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _RL myTime        _RL myTime
74        INTEGER myIter, myThid        INTEGER myIter, myThid
75    
# Line 87  C     == Local variables == Line 87  C     == Local variables ==
87        LOGICAL useDiagPhiRlow, addSurfPhiAnom        LOGICAL useDiagPhiRlow, addSurfPhiAnom
88  CEOP  CEOP
89        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
90        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GT.3        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
91        surfPhiFac = 0.        surfPhiFac = 0.
92        IF (addSurfPhiAnom) surfPhiFac = 1.        IF (addSurfPhiAnom) surfPhiFac = 1.
93    
# Line 123  C--   Initialize phiHydF to zero : Line 123  C--   Initialize phiHydF to zero :
123  C     note: atmospheric_loading or Phi_topo anomaly are incorporated  C     note: atmospheric_loading or Phi_topo anomaly are incorporated
124  C           later in S/R calc_grad_phi_hyd  C           later in S/R calc_grad_phi_hyd
125        IF (k.EQ.1) THEN        IF (k.EQ.1) THEN
126          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
127           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
128             phiHydF(i,j) = 0.             phiHydF(i,j) = 0.
129           ENDDO           ENDDO
130          ENDDO          ENDDO
# Line 191  C Quasi-hydrostatic terms are added in a Line 191  C Quasi-hydrostatic terms are added in a
191  #endif /* ALLOW_MOM_COMMON */  #endif /* ALLOW_MOM_COMMON */
192    
193  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
194          IF (k.EQ.1 .AND. addSurfPhiAnom) THEN          IF ( addSurfPhiAnom .AND.
195         &       uniformFreeSurfLev .AND. k.EQ.1 ) THEN
196            DO j=jMin,jMax            DO j=jMin,jMax
197              DO i=iMin,iMax              DO i=iMin,iMax
198                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)                phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
# Line 206  C----  Hydrostatic pressure at cell cent Line 207  C----  Hydrostatic pressure at cell cent
207         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
208  C  --  Finite Volume Form  C  --  Finite Volume Form
209    
          DO j=jMin,jMax  
           DO i=iMin,iMax  
   
210  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
211  C           which has not been used to date since it does not  C           which has not been used to date since it does not
212  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
213  C  
214             phiHydC(i,j)=phiHydF(i,j)          IF ( uniformFreeSurfLev ) THEN
215       &       + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst           DO j=jMin,jMax
216             phiHydF(i,j)=phiHydF(i,j)            DO i=iMin,iMax
217       &            + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst              phiHydC(i,j) = phiHydF(i,j)
218         &            + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
219                phiHydF(i,j) = phiHydF(i,j)
220         &                 + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
221              ENDDO
222             ENDDO
223            ELSE
224             DO j=jMin,jMax
225              DO i=iMin,iMax
226               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
227                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
228    #ifdef NONLIN_FRSURF
229                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
230    #endif
231                phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst
232               ELSE
233                phiHydC(i,j) = phiHydF(i,j)
234         &              + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
235               ENDIF
236                phiHydF(i,j) = phiHydC(i,j)
237         &              + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
238            ENDDO            ENDDO
239           ENDDO           ENDDO
240            ENDIF
241    
242         ELSE         ELSE
243  C  --  Finite Difference Form  C  --  Finite Difference Form
244    
245           dRlocM=half*drC(k)  C---------- This discretization is the "energy conserving" form
246           IF (k.EQ.1) dRlocM=rF(k)-rC(k)  C           which has been used since at least Adcroft et al., MWR 1997
          IF (k.EQ.Nr) THEN  
            dRlocP=rC(k)-rF(k+1)  
          ELSE  
            dRlocP=half*drC(k+1)  
          ENDIF  
247    
248            dRlocM=half*drC(k)
249            IF (k.EQ.1) dRlocM=rF(k)-rC(k)
250            IF (k.EQ.Nr) THEN
251              dRlocP=rC(k)-rF(k+1)
252            ELSE
253              dRlocP=half*drC(k+1)
254            ENDIF
255            IF ( uniformFreeSurfLev ) THEN
256           DO j=jMin,jMax           DO j=jMin,jMax
257            DO i=iMin,iMax            DO i=iMin,iMax
258                phiHydC(i,j) = phiHydF(i,j)
259  C---------- This discretization is the "energy conserving" form       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
260  C           which has been used since at least Adcroft et al., MWR 1997              phiHydF(i,j) = phiHydC(i,j)
261  C       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
262              phiHydC(i,j)=phiHydF(i,j)            ENDDO
263             ENDDO
264            ELSE
265             rec_dRm = one/(rF(k)-rC(k))
266             rec_dRp = one/(rC(k)-rF(k+1))
267             DO j=jMin,jMax
268              DO i=iMin,iMax
269               IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
270                ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
271    #ifdef NONLIN_FRSURF
272                ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
273    #endif
274                phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM
275         &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP
276         &                     )*gravity*alphaRho(i,j)*recip_rhoConst
277               ELSE
278                phiHydC(i,j) = phiHydF(i,j)
279       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
280              phiHydF(i,j)=phiHydC(i,j)             ENDIF
281                phiHydF(i,j) = phiHydC(i,j)
282       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
283            ENDDO            ENDDO
284           ENDDO           ENDDO
285            ENDIF
286    
287  C  --  end if integr_GeoPot = ...  C  --  end if integr_GeoPot = ...
288         ENDIF         ENDIF
# Line 303  C  --  Finite Volume Form Line 343  C  --  Finite Volume Form
343  C---------- This discretization is the "finite volume" form  C---------- This discretization is the "finite volume" form
344  C           which has not been used to date since it does not  C           which has not been used to date since it does not
345  C           conserve KE+PE exactly even though it is more natural  C           conserve KE+PE exactly even though it is more natural
346  C  
347             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
348               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)               ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
349  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
Legend:
Removed from v.1.40  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22