/[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.16 by jmc, Tue Oct 4 02:49:13 2005 UTC revision 1.17 by baylor, Mon Oct 10 17:23:07 2005 UTC
# Line 17  C     Calculate horizontal viscosities ( Line 17  C     Calculate horizontal viscosities (
17  C     harmonic viscosity=  C     harmonic viscosity=
18  C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)  C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
19  C       +0.25*L**2*viscAhGrid/deltaT  C       +0.25*L**2*viscAhGrid/deltaT
20  C       +sqrt(viscC2leith**2*grad(Vort3)**2  C       +sqrt((viscC2leith/pi)**6*grad(Vort3)**2
21  C             +viscC2leithD**2*grad(hDiv)**2)*L**3  C             +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3
22  C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)  C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
23  C  C
24  C     biharmonic viscosity=  C     biharmonic viscosity=
25  C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)  C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
26  C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)  C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
27  C       +0.125*L**5*sqrt(viscC4leith**2*grad(Vort3)**2  C       +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2
28  C                        +viscC4leithD**2*grad(hDiv)**2)  C                        +(viscC4leithD/pi)**6*grad(hDiv)**2)
29  C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)  C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
30  C  C
31  C     Note that often 0.125*L**2 is the scale between harmonic and  C     Note that often 0.125*L**2 is the scale between harmonic and
# Line 91  C     == Routine arguments == Line 91  C     == Routine arguments ==
91  C     == Local variables ==  C     == Local variables ==
92        INTEGER I,J        INTEGER I,J
93        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
94          _RL leith2fac, leith4fac
95          _RL leithD2fac, leithD4fac
96        _RL viscAhRe_max, viscA4Re_max        _RL viscAhRe_max, viscA4Re_max
97        _RL Alin,grdVrt,grdDiv, keZpt        _RL Alin,grdVrt,grdDiv, keZpt
98        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
# Line 183  C     == Local variables == Line 185  C     == Local variables ==
185          smag4fac=0. _d 0          smag4fac=0. _d 0
186        ENDIF        ENDIF
187    
188          IF (calcleith) THEN
189            IF (useFullLeith) THEN
190             leith2fac=(viscC2leith/pi)**6
191             leith4fac=0.015625 _d 0*(viscC4leith/pi)**6
192             leithD2fac=(viscC2leithD/pi)**6
193             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
194            ELSE
195             leith2fac=(viscC2leith/pi)**3
196             leithD2fac=(viscC2leithD/pi)**3
197             leith4fac=0.0125 _d 0*(viscC4leith/pi)**3
198             leithD4fac=0.0125 _d 0*(viscC4leithD/pi)**3
199            ENDIF
200          ELSE
201            leith2fac=0. _d 0
202            leith4fac=0. _d 0
203            leithD2fac=0. _d 0
204            leithD4fac=0. _d 0
205          ENDIF
206    
207  C     - Viscosity  C     - Viscosity
208        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
209    
# Line 276  C Using it in Leith serves to damp insta Line 297  C Using it in Leith serves to damp insta
297       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
298    
299            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
300       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
301            viscA4_DLth(i,j)=0.125 _d 0*            viscA4_DLth(i,j)=
302       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
303            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
304       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
305            viscA4_DLthd(i,j)=0.125 _d 0*            viscA4_DLthd(i,j)=
306       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
307           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
308  C but this approximation will work on cube  C but this approximation will work on cube
309  c (and differs by as much as 4X)  c (and differs by as much as 4X)
# Line 299  c (and differs by as much as 4X) Line 320  c (and differs by as much as 4X)
320            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
321    
322  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
323            viscAh_Dlth(i,j)=            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
324       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
325            viscA4_Dlth(i,j)=0.125 _d 0*            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
326       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
           viscAh_DlthD(i,j)=  
      &      ((viscC2leithD*grdDiv))*L3  
           viscA4_DlthD(i,j)=0.125 _d 0*  
      &      ((viscC4leithD*grdDiv))*L5  
327           ELSE           ELSE
328            viscAh_Dlth(i,j)=0. _d 0            viscAh_Dlth(i,j)=0. _d 0
329            viscA4_Dlth(i,j)=0. _d 0            viscA4_Dlth(i,j)=0. _d 0
# Line 394  C This is the vector magnitude of grad(d Line 411  C This is the vector magnitude of grad(d
411       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
412    
413            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
414       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
415            viscA4_ZLth(i,j)=0.125 _d 0*            viscA4_ZLth(i,j)=
416       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
417            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
418       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
419            viscA4_ZLthD(i,j)=0.125 _d 0*            viscA4_ZLthD(i,j)=
420       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
421    
422           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
423  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
# Line 416  C but this approximation will work on cu Line 433  C but this approximation will work on cu
433            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
434            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
435    
436            viscAh_ZLth(i,j)=(viscC2leith*grdVrt            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
437       &                     +(viscC2leithD*grdDiv))*L3            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
438            viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
439       &                     +(viscC4leithD*grdDiv))*L5            viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
           viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3  
           viscA4_ZLthD(i,j)=0.125 _d 0*((viscC4leithD*grdDiv))*L5  
440           ELSE           ELSE
441            viscAh_ZLth(i,j)=0. _d 0            viscAh_ZLth(i,j)=0. _d 0
442            viscA4_ZLth(i,j)=0. _d 0            viscA4_ZLth(i,j)=0. _d 0
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22