/[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.9 by baylor, Wed Sep 21 23:17:34 2005 UTC revision 1.10 by jmc, Thu Sep 22 00:21:23 2005 UTC
# Line 137  C     == Local variables == Line 137  C     == Local variables ==
137       &  .OR.(viscC2smag.NE.0.)       &  .OR.(viscC2smag.NE.0.)
138    
139        IF ((harmonic).and.(viscAhremax.ne.0.)) THEN        IF ((harmonic).and.(viscAhremax.ne.0.)) THEN
140          viscAhre_max=sqrt(2d0)/viscAhRemax          viscAhre_max=sqrt(2. _d 0)/viscAhRemax
141        ELSE        ELSE
142          viscAhre_max=0d0          viscAhre_max=0. _d 0
143        ENDIF        ENDIF
144    
145        biharmonic=        biharmonic=
# Line 152  C     == Local variables == Line 152  C     == Local variables ==
152       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
153    
154        IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN        IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN
155          viscA4re_max=0.125d0*sqrt(2d0)/viscA4Remax          viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax
156        ELSE        ELSE
157          viscA4re_max=0d0          viscA4re_max=0. _d 0
158        ENDIF        ENDIF
159    
160        calcleith=        calcleith=
# Line 168  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        ELSE
180          smag2fac=0d0          smag2fac=0. _d 0
181          smag4fac=0d0          smag4fac=0. _d 0
182        ENDIF        ENDIF
183    
184  C     - Viscosity  C     - Viscosity
# Line 188  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)*viscAhRe_max           Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
# Line 206  C Velocity Reynolds Scale Line 206  C Velocity Reynolds Scale
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 216  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 224  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)=0.125d0*            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)=0.125d0*            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 252  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 295  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(0.25d0*(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)*viscAhRe_max       &         *L2)*viscAhRe_max
311           U4scl=sqrt(0.25d0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1)))           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       &         *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 328  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)=0.125d0*            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)=0.125d0*            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 349  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)=0.125d0*(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)=0.125d0*((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.9  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22