/[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.15 by jmc, Mon Oct 3 21:43:03 2005 UTC revision 1.18 by baylor, Mon Oct 10 18:45:30 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 51  C      harmonic viscosity>0.25*viscAhgri Line 51  C      harmonic viscosity>0.25*viscAhgri
51  C      biharmonic viscosity>viscA4gridmax*L**4/32/deltaT (approx)  C      biharmonic viscosity>viscA4gridmax*L**4/32/deltaT (approx)
52  C  C
53  C     RECOMMENDED VALUES  C     RECOMMENDED VALUES
54  C     viscC2Leith=?  C     viscC2Leith=1-3
55  C     viscC2LeithD=?  C     viscC2LeithD=1-3
56  C     viscC4Leith=?  C     viscC4Leith=1-3
57  C     viscC4LeithD=?  C     viscC4LeithD=1.5-3
58  C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)  C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
59  C               0.2-0.9 (Smagorinsky,1993)  C               0.2-0.9 (Smagorinsky,1993)
60  C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)  C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
# 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
99        _RL Uscl,U4scl        _RL Uscl,U4scl
100          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 181  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    
210    C      horizontal gradient of horizontal divergence:
211           DO j=1-Oly,sNy+Oly
212            DO i=1-Olx,sNx+Olx
213              divDx(i,j) = 0.
214              divDy(i,j) = 0.
215            ENDDO
216           ENDDO
217           IF (calcleith) THEN
218    C-       gradient in x direction:
219    #ifndef ALLOW_AUTODIFF_TAMC
220             IF (useCubedSphereExchange) THEN
221    C        to compute d/dx(hDiv), fill corners with appropriate values:
222               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
223             ENDIF
224    #endif
225             DO j=2-Oly,sNy+Oly-1
226              DO i=2-Olx,sNx+Olx-1
227                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
228              ENDDO
229             ENDDO
230    
231    C-       gradient in y direction:
232    #ifndef ALLOW_AUTODIFF_TAMC
233             IF (useCubedSphereExchange) THEN
234    C        to compute d/dy(hDiv), fill corners with appropriate values:
235               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
236             ENDIF
237    #endif
238             DO j=2-Oly,sNy+Oly-1
239              DO i=2-Olx,sNx+Olx-1
240                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
241              ENDDO
242             ENDDO
243           ENDIF
244    
245         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
246          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
247  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
# Line 232  C This is the vector magnitude of the vo Line 291  C This is the vector magnitude of the vo
291    
292  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
293  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
294            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
295       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
296       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
297       &     +((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
      &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2)  
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 257  c (and differs by as much as 4X) Line 315  c (and differs by as much as 4X)
315            grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
316       &     abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))       &     abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))
317    
318            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
319            grdDiv=max(grdDiv,            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
320       &     abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i,j))   )
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)))  
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))  
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 355  C This is the vector magnitude of the vo Line 405  C This is the vector magnitude of the vo
405       &     +((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)
406    
407  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
408            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
409       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
410       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
411       &     +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
      &     +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)  
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 380  C but this approximation will work on cu Line 429  C but this approximation will work on cu
429            grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
430       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
431    
432            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
433            grdDiv=max(grdDiv,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
434       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
435            grdDiv=max(grdDiv,  
436       &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj)))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
437            grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
438       &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
439              viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
           viscAh_ZLth(i,j)=(viscC2leith*grdVrt  
      &                     +(viscC2leithD*grdDiv))*L3  
           viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt  
      &                     +(viscC4leithD*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.15  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22