/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_visc.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_common/mom_calc_visc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.7 by baylor, Wed Sep 21 14:32:59 2005 UTC revision 1.8 by baylor, Wed Sep 21 15:11:22 2005 UTC
# Line 160  C     == Local variables == Line 160  C     == Local variables ==
160       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
161    
162        IF (deltaTmom.NE.0.) THEN        IF (deltaTmom.NE.0.) THEN
163         recip_dt=1./deltaTmom         recip_dt=1d0/deltaTmom
164        ELSE        ELSE
165         recip_dt=0.         recip_dt=0d0
166        ENDIF        ENDIF
167    
168        IF (calcsmag) THEN        IF (calcsmag) THEN
169          smag2fac=(viscC2smag/pi)**2          smag2fac=(viscC2smag/pi)**2
170          smag4fac=0.125*(viscC4smag/pi)**2          smag4fac=0.125d0*(viscC4smag/pi)**2
171        ENDIF        ENDIF
172    
173  C     - Viscosity  C     - Viscosity
# Line 177  C     - Viscosity Line 177  C     - Viscosity
177  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
178    
179  C These are (powers of) length scales  C These are (powers of) length scales
180           L2=2./((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))           L2=2d0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
181           L3=(L2**1.5)           L3=(L2**1.5)
182           L4=(L2**2)           L4=(L2**2)
183           L5=(L2**2.5)           L5=(L2**2.5)
184    
185           L2rdt=0.25*recip_dt*L2           L2rdt=0.25d0*recip_dt*L2
186    
187           L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4           L4rdt=recip_dt/( 6d0*(recip_DXF(I,J,bi,bj)**4
188       &                       +recip_DYF(I,J,bi,bj)**4)       &                       +recip_DYF(I,J,bi,bj)**4)
189       &                   +8.*((recip_DXF(I,J,bi,bj)       &                   +8d0*((recip_DXF(I,J,bi,bj)
190       &                        *recip_DYF(I,J,bi,bj))**2) )       &                        *recip_DYF(I,J,bi,bj))**2) )
191    
192  C Velocity Reynolds Scale  C Velocity Reynolds Scale
193           Uscl=sqrt(KE(i,j)*L2*0.5)/viscAhRe_max           Uscl=sqrt(KE(i,j)*L2*0.5d0)/viscAhRe_max
194           U4scl=0.125*L2*Uscl/viscA4Re_max           U4scl=0.125d0*L2*Uscl/viscA4Re_max
195    
196           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
197  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
198            grdVrt=0.25*(            grdVrt=0.25d0*(
199       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
200       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
201       &     +((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj))**2       &     +((vort3(i+1,j+1)-vort3(i,j+1))
202       &     +((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)       &               *recip_DXG(i,j+1,bi,bj))**2
203         &     +((vort3(i+1,j+1)-vort3(i+1,j))
204         &               *recip_DYG(i+1,j,bi,bj))**2)
205    
206  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
207  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
208            grdDiv=0.25*(            grdDiv=0.25d0*(
209       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2
210       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2
211       &     +((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &     +((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2
# Line 239  c (and differs by as much as 4X) Line 241  c (and differs by as much as 4X)
241  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
242            viscAh_Dlth(i,j)=            viscAh_Dlth(i,j)=
243       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3
244            viscA4_Dlth(i,j)=0.125*            viscA4_Dlth(i,j)=0.125d0*
245       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5
246            viscAh_DlthD(i,j)=            viscAh_DlthD(i,j)=
247       &      ((viscC2leithD*grdDiv))*L3       &      ((viscC2leithD*grdDiv))*L3
248            viscA4_DlthD(i,j)=0.125*            viscA4_DlthD(i,j)=0.125d0*
249       &      ((viscC4leithD*grdDiv))*L5       &      ((viscC4leithD*grdDiv))*L5
250           ELSE           ELSE
251            viscAh_Dlth(i,j)=0d0            viscAh_Dlth(i,j)=0d0
# Line 255  c This approximation is good to the same Line 257  c This approximation is good to the same
257           IF (calcsmag) THEN           IF (calcsmag) THEN
258            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
259       &       *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
260       &       +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2       &       +0.25d0*(strain(i+1, j )**2+strain( i ,j+1)**2
261       &              +strain(i  , j )**2+strain(i+1,j+1)**2))       &              +strain(i  , j )**2+strain(i+1,j+1)**2))
262            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
263            viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)            viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
# Line 282  C  BiHarmonic on Div.u points Line 284  C  BiHarmonic on Div.u points
284    
285  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
286  C These are (powers of) length scales  C These are (powers of) length scales
287           L2=2./((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))           L2=2d0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
288           L3=(L2**1.5)           L3=(L2**1.5)
289           L4=(L2**2)           L4=(L2**2)
290           L5=(L2**2.5)           L5=(L2**2.5)
291    
292           L2rdt=0.25*recip_dt*L2           L2rdt=0.25d0*recip_dt*L2
293           L4rdt=recip_dt/           L4rdt=recip_dt/
294       &     ( 6.*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)       &     ( 6d0*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)
295       &      +8.*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))       &      +8d0*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))
296    
297  C Velocity Reynolds Scale  C Velocity Reynolds Scale
298           Uscl=sqrt((KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))*L2*0.125)/           Uscl=sqrt((KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))
299       &         viscAhRe_max       &         *L2*0.125d0)/viscAhRe_max
300           U4scl=0.125*L2*Uscl/viscA4Re_max           U4scl=0.125d0*L2*Uscl/viscA4Re_max
301    
302  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
303           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
304            grdVrt=0.25*(            grdVrt=0.25d0*(
305       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
306       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
307       &     +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2       &     +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2
308       &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)       &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
309    
310  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
311            grdDiv=0.25*(            grdDiv=0.25d0*(
312       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2
313       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2
314       &     +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2       &     +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2
# Line 355  C but this approximation will work on cu Line 357  C but this approximation will work on cu
357           IF (calcsmag) THEN           IF (calcsmag) THEN
358            viscAh_ZSmg(i,j)=L2            viscAh_ZSmg(i,j)=L2
359       &      *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
360       &        +0.25*(tension( i , j )**2+tension( i ,j-1)**2       &        +0.25d0*(tension( i , j )**2+tension( i ,j-1)**2
361       &              +tension(i-1, j )**2+tension(i-1,j-1)**2))       &              +tension(i-1, j )**2+tension(i-1,j-1)**2))
362            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
363            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
# Line 412  C  BiHarmonic on Zeta points Line 414  C  BiHarmonic on Zeta points
414         CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
415    
416         CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'         CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
417       &    ,k,1,2,bi,bj,myThid)       &   ,k,1,2,bi,bj,myThid)
418         CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'         CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
419       &    ,k,1,2,bi,bj,myThid)       &   ,k,1,2,bi,bj,myThid)
420         CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'         CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
421       &    ,k,1,2,bi,bj,myThid)       &   ,k,1,2,bi,bj,myThid)
422         CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'         CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
423       &    ,k,1,2,bi,bj,myThid)       &   ,k,1,2,bi,bj,myThid)
424    
425         CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
426         CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22