/[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.8 by baylor, Wed Sep 21 15:11:22 2005 UTC revision 1.10 by jmc, Thu Sep 22 00:21:23 2005 UTC
# Line 35  C     for roughly similar results with b Line 35  C     for roughly similar results with b
35  C  C
36  C     LIMITERS -- limit min and max values of viscosities  C     LIMITERS -- limit min and max values of viscosities
37  C     viscAhRemax is min value for grid point harmonic Reynolds num  C     viscAhRemax is min value for grid point harmonic Reynolds num
38  C      harmonic viscosity>sqrt(2*KE)*L/2/viscAhRemax  C      harmonic viscosity>sqrt(2*KE)*L/viscAhRemax
39  C  C
40  C     viscA4Remax is min value for grid point biharmonic Reynolds num  C     viscA4Remax is min value for grid point biharmonic Reynolds num
41  C      biharmonic viscosity>sqrt(2*KE)*L**3/16/viscA4Remax  C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4Remax
42  C  C
43  C     viscAhgridmax is CFL stability limiter for harmonic viscosity  C     viscAhgridmax is CFL stability limiter for harmonic viscosity
44  C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT  C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
# Line 58  C     viscC4LeithD=? Line 58  C     viscC4LeithD=?
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)
61  C     viscAhRemax>=1  C     viscAhRemax>=1, (<2 suppresses a computational mode)
62  C     viscA4Remax>=1  C     viscA4Remax>=1, (<2 suppresses a computational mode)
63  C     viscAhgridmax=1  C     viscAhgridmax=1
64  C     viscA4gridmax=1  C     viscA4gridmax=1
65  C     viscAhgrid<1  C     viscAhgrid<1
# Line 136  C     == Local variables == Line 136  C     == Local variables ==
136       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
137       &  .OR.(viscC2smag.NE.0.)       &  .OR.(viscC2smag.NE.0.)
138    
139        IF (harmonic) viscAhre_max=viscAhremax        IF ((harmonic).and.(viscAhremax.ne.0.)) THEN
140            viscAhre_max=sqrt(2. _d 0)/viscAhRemax
141          ELSE
142            viscAhre_max=0. _d 0
143          ENDIF
144    
145        biharmonic=        biharmonic=
146       &      (viscA4.NE.0.)       &      (viscA4.NE.0.)
# Line 147  C     == Local variables == Line 151  C     == Local variables ==
151       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
152       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
153    
154        IF (biharmonic) viscA4re_max=viscA4remax        IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN
155            viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax
156          ELSE
157            viscA4re_max=0. _d 0
158          ENDIF
159    
160        calcleith=        calcleith=
161       &      (viscC2leith.NE.0.)       &      (viscC2leith.NE.0.)
# Line 160  C     == Local variables == Line 168  C     == Local variables ==
168       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
169    
170        IF (deltaTmom.NE.0.) THEN        IF (deltaTmom.NE.0.) THEN
171         recip_dt=1d0/deltaTmom         recip_dt=1. _d 0/deltaTmom
172        ELSE        ELSE
173         recip_dt=0d0         recip_dt=0. _d 0
174        ENDIF        ENDIF
175    
176        IF (calcsmag) THEN        IF (calcsmag) THEN
177          smag2fac=(viscC2smag/pi)**2          smag2fac=(viscC2smag/pi)**2
178          smag4fac=0.125d0*(viscC4smag/pi)**2          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
179          ELSE
180            smag2fac=0. _d 0
181            smag4fac=0. _d 0
182        ENDIF        ENDIF
183    
184  C     - Viscosity  C     - Viscosity
# Line 177  C     - Viscosity Line 188  C     - Viscosity
188  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
189    
190  C These are (powers of) length scales  C These are (powers of) length scales
191           L2=2d0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))           L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
192           L3=(L2**1.5)           L3=(L2**1.5)
193           L4=(L2**2)           L4=(L2**2)
194           L5=(L2**2.5)           L5=(L2**2.5)
195    
196           L2rdt=0.25d0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
197    
198           L4rdt=recip_dt/( 6d0*(recip_DXF(I,J,bi,bj)**4           L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
199       &                       +recip_DYF(I,J,bi,bj)**4)       &                            +recip_DYF(I,J,bi,bj)**4)
200       &                   +8d0*((recip_DXF(I,J,bi,bj)       &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
201       &                        *recip_DYF(I,J,bi,bj))**2) )       &                             *recip_DYF(I,J,bi,bj))**2) )
202    
203  C Velocity Reynolds Scale  C Velocity Reynolds Scale
204           Uscl=sqrt(KE(i,j)*L2*0.5d0)/viscAhRe_max           Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
205           U4scl=0.125d0*L2*Uscl/viscA4Re_max           U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
206    
207           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
208  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
209            grdVrt=0.25d0*(            grdVrt=0.25 _d 0*(
210       &     ((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
211       &     +((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
212       &     +((vort3(i+1,j+1)-vort3(i,j+1))       &     +((vort3(i+1,j+1)-vort3(i,j+1))
# Line 205  C This is the vector magnitude of the vo Line 216  C This is the vector magnitude of the vo
216    
217  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
218  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
219            grdDiv=0.25d0*(            grdDiv=0.25 _d 0*(
220       &     ((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
221       &     +((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
222       &     +((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 213  C Using it in Leith serves to damp insta Line 224  C Using it in Leith serves to damp insta
224    
225            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
226       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
227            viscA4_DLth(i,j)=            viscA4_DLth(i,j)=0.125 _d 0*
228       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
229            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
230       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(viscC2leithD**2*grdDiv)*L3
231            viscA4_DLthd(i,j)=            viscA4_DLthd(i,j)=0.125 _d 0*
232       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(viscC4leithD**2*grdDiv)*L5
233           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
234  C but this approximation will work on cube  C but this approximation will work on cube
# Line 241  c (and differs by as much as 4X) Line 252  c (and differs by as much as 4X)
252  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
253            viscAh_Dlth(i,j)=            viscAh_Dlth(i,j)=
254       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3
255            viscA4_Dlth(i,j)=0.125d0*            viscA4_Dlth(i,j)=0.125 _d 0*
256       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5
257            viscAh_DlthD(i,j)=            viscAh_DlthD(i,j)=
258       &      ((viscC2leithD*grdDiv))*L3       &      ((viscC2leithD*grdDiv))*L3
259            viscA4_DlthD(i,j)=0.125d0*            viscA4_DlthD(i,j)=0.125 _d 0*
260       &      ((viscC4leithD*grdDiv))*L5       &      ((viscC4leithD*grdDiv))*L5
261           ELSE           ELSE
262            viscAh_Dlth(i,j)=0d0            viscAh_Dlth(i,j)=0. _d 0
263            viscA4_Dlth(i,j)=0d0            viscA4_Dlth(i,j)=0. _d 0
264            viscAh_DlthD(i,j)=0d0            viscAh_DlthD(i,j)=0. _d 0
265            viscA4_DlthD(i,j)=0d0            viscA4_DlthD(i,j)=0. _d 0
266           ENDIF           ENDIF
267    
268           IF (calcsmag) THEN           IF (calcsmag) THEN
269            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
270       &       *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
271       &       +0.25d0*(strain(i+1, j )**2+strain( i ,j+1)**2       &       +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2
272       &              +strain(i  , j )**2+strain(i+1,j+1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
273            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
274            viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)            viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
275           ELSE           ELSE
276            viscAh_DSmg(i,j)=0d0            viscAh_DSmg(i,j)=0. _d 0
277            viscA4_DSmg(i,j)=0d0            viscA4_DSmg(i,j)=0. _d 0
278           ENDIF           ENDIF
279    
280  C  Harmonic on Div.u points  C  Harmonic on Div.u points
# Line 284  C  BiHarmonic on Div.u points Line 295  C  BiHarmonic on Div.u points
295    
296  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
297  C These are (powers of) length scales  C These are (powers of) length scales
298           L2=2d0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))           L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
299           L3=(L2**1.5)           L3=(L2**1.5)
300           L4=(L2**2)           L4=(L2**2)
301           L5=(L2**2.5)           L5=(L2**2.5)
302    
303           L2rdt=0.25d0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
304           L4rdt=recip_dt/           L4rdt=recip_dt/
305       &     ( 6d0*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)       &     ( 6. _d 0*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)
306       &      +8d0*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))       &      +8. _d 0*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))
307    
308  C Velocity Reynolds Scale  C Velocity Reynolds Scale
309           Uscl=sqrt((KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))           Uscl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))
310       &         *L2*0.125d0)/viscAhRe_max       &         *L2)*viscAhRe_max
311           U4scl=0.125d0*L2*Uscl/viscA4Re_max           U4scl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1)))
312         &         *L3*viscA4Re_max
313    
314  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
315           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
316            grdVrt=0.25d0*(            grdVrt=0.25 _d 0*(
317       &     ((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
318       &     +((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
319       &     +((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
320       &     +((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)
321    
322  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
323            grdDiv=0.25d0*(            grdDiv=0.25 _d 0*(
324       &     ((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
325       &     +((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
326       &     +((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 316  C This is the vector magnitude of grad(d Line 328  C This is the vector magnitude of grad(d
328    
329            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
330       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
331            viscA4_ZLth(i,j)=            viscA4_ZLth(i,j)=0.125 _d 0*
332       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
333            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
334       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(viscC2leithD**2*grdDiv)*L3
335            viscA4_ZLthD(i,j)=            viscA4_ZLthD(i,j)=0.125 _d 0*
336       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(viscC4leithD**2*grdDiv)*L5
337    
338           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
# Line 337  C but this approximation will work on cu Line 349  C but this approximation will work on cu
349            grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
350       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))
351            grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
352       &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXG(i,j-1,bi,bj)))       &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj)))
353            grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
354       &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYG(i-1,j,bi,bj)))       &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj)))
355    
356            viscAh_ZLth(i,j)=(viscC2leith*grdVrt            viscAh_ZLth(i,j)=(viscC2leith*grdVrt
357       &                     +(viscC2leithD*grdDiv))*L3       &                     +(viscC2leithD*grdDiv))*L3
358            viscA4_ZLth(i,j)=(viscC4leith*grdVrt            viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt
359       &                     +(viscC4leithD*grdDiv))*L5       &                     +(viscC4leithD*grdDiv))*L5
360            viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3            viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3
361            viscA4_ZLthD(i,j)=((viscC4leithD*grdDiv))*L5            viscA4_ZLthD(i,j)=0.125 _d 0*((viscC4leithD*grdDiv))*L5
362           ELSE           ELSE
363            viscAh_ZLth(i,j)=0d0            viscAh_ZLth(i,j)=0. _d 0
364            viscA4_ZLth(i,j)=0d0            viscA4_ZLth(i,j)=0. _d 0
365            viscAh_ZLthD(i,j)=0d0            viscAh_ZLthD(i,j)=0. _d 0
366            viscA4_ZLthD(i,j)=0d0            viscA4_ZLthD(i,j)=0. _d 0
367           ENDIF           ENDIF
368    
369           IF (calcsmag) THEN           IF (calcsmag) THEN
370            viscAh_ZSmg(i,j)=L2            viscAh_ZSmg(i,j)=L2
371       &      *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
372       &        +0.25d0*(tension( i , j )**2+tension( i ,j-1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
373       &              +tension(i-1, j )**2+tension(i-1,j-1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
374            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
375            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
376           ENDIF           ENDIF
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22