/[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.1 by baylor, Fri Sep 16 19:33:59 2005 UTC revision 1.15 by jmc, Mon Oct 3 21:43:03 2005 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "MOM_COMMON_OPTIONS.h"  #include "MOM_COMMON_OPTIONS.h"
5    
6    
7        SUBROUTINE MOM_CALC_VISC(        SUBROUTINE MOM_CALC_VISC(
8       I        bi,bj,k,       I        bi,bj,k,
9       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
10       O        harmonic,biharmonic,useVariableViscosity,       O        harmonic,biharmonic,useVariableViscosity,
11       I        hDiv,vort3,tension,strain,KE,       I        hDiv,vort3,tension,strain,KE,hFacZ,
12       I        myThid)       I        myThid)
13    
14        IMPLICIT NONE        IMPLICIT NONE
15    C
16    C     Calculate horizontal viscosities (L is typical grid width)
17    C     harmonic viscosity=
18    C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
19    C       +0.25*L**2*viscAhGrid/deltaT
20    C       +sqrt(viscC2leith**2*grad(Vort3)**2
21    C             +viscC2leithD**2*grad(hDiv)**2)*L**3
22    C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
23    C
24    C     biharmonic viscosity=
25    C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
26    C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
27    C       +0.125*L**5*sqrt(viscC4leith**2*grad(Vort3)**2
28    C                        +viscC4leithD**2*grad(hDiv)**2)
29    C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
30    C
31    C     Note that often 0.125*L**2 is the scale between harmonic and
32    C     biharmonic (see Griffies and Hallberg (2000))
33    C     This allows the same value of the coefficient to be used
34    C     for roughly similar results with biharmonic and harmonic
35    C
36    C     LIMITERS -- limit min and max values of viscosities
37    C     viscAhRemax is min value for grid point harmonic Reynolds num
38    C      harmonic viscosity>sqrt(2*KE)*L/viscAhRemax
39    C
40    C     viscA4Remax is min value for grid point biharmonic Reynolds num
41    C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4Remax
42    C
43    C     viscAhgridmax is CFL stability limiter for harmonic viscosity
44    C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
45    C
46    C     viscA4gridmax is CFL stability limiter for biharmonic viscosity
47    C      biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)
48    C
49    C     viscAhgridmin and viscA4gridmin are lower limits for viscosity:
50    C      harmonic viscosity>0.25*viscAhgridmax*L**2/deltaT
51    C      biharmonic viscosity>viscA4gridmax*L**4/32/deltaT (approx)
52    C
53    C     RECOMMENDED VALUES
54    C     viscC2Leith=?
55    C     viscC2LeithD=?
56    C     viscC4Leith=?
57    C     viscC4LeithD=?
58    C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
59    C               0.2-0.9 (Smagorinsky,1993)
60    C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
61    C     viscAhRemax>=1, (<2 suppresses a computational mode)
62    C     viscA4Remax>=1, (<2 suppresses a computational mode)
63    C     viscAhgridmax=1
64    C     viscA4gridmax=1
65    C     viscAhgrid<1
66    C     viscA4grid<1
67    C     viscAhgridmin<<1
68    C     viscA4gridmin<<1
69    
70  C     == Global variables ==  C     == Global variables ==
71  #include "SIZE.h"  #include "SIZE.h"
# Line 35  C     == Routine arguments == Line 90  C     == Routine arguments ==
90    
91  C     == Local variables ==  C     == Local variables ==
92        INTEGER I,J        INTEGER I,J
93        _RL ASmag2, ASmag4, smag2fac, smag4fac        _RL smag2fac, smag4fac
94        _RL vg2Min, vg2Max, AlinMax, AlinMin        _RL viscAhRe_max, viscA4Re_max
95        _RL lenA2, lenAz2        _RL Alin,grdVrt,grdDiv, keZpt
       _RL Alin,Alth2,Alth4,grdVrt,grdDiv  
       _RL vg2,vg4,vg4Min,vg4Max  
96        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
97          _RL Uscl,U4scl
98          _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99          _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100          _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104          _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105          _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106          _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107          _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108          _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109          _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110          _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112          _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115          _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116          _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117          _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118          LOGICAL calcLeith,calcSmag
119    
120        useVariableViscosity=        useVariableViscosity=
121       &      (viscAhGrid.NE.0.)       &      (viscAhGrid.NE.0.)
# Line 62  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).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.)
147       &  .OR.(viscA4D.NE.0.)       &  .OR.(viscA4D.NE.0.)
# Line 71  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 (deltaTmom.NE.0.) THEN        IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN
155         recip_dt=1./deltaTmom          viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax
156        ELSE        ELSE
157         recip_dt=0.          viscA4re_max=0. _d 0
158        ENDIF        ENDIF
159    
160        vg2=viscAhGrid*recip_dt        calcleith=
161        vg2Min=viscAhGridMin*recip_dt       &      (viscC2leith.NE.0.)
162        vg2Max=viscAhGridMax*recip_dt       &  .OR.(viscC2leithD.NE.0.)
163        vg4=viscA4Grid*recip_dt       &  .OR.(viscC4leith.NE.0.)
164        vg4Min=viscA4GridMin*recip_dt       &  .OR.(viscC4leithD.NE.0.)
       vg4Max=viscA4GridMax*recip_dt  
165    
166        smag2fac=(viscC2smag/pi)**2        calcsmag=
167        smag4fac=0.125*(viscC4smag/pi)**2       &      (viscC2smag.NE.0.)
168         &  .OR.(viscC4smag.NE.0.)
169    
170          IF (deltaTmom.NE.0.) THEN
171           recip_dt=1. _d 0/deltaTmom
172          ELSE
173           recip_dt=0. _d 0
174          ENDIF
175    
176          IF (calcsmag) THEN
177            smag2fac=(viscC2smag/pi)**2
178            smag4fac=0.125 _d 0*(viscC4smag/pi)**2
179          ELSE
180            smag2fac=0. _d 0
181            smag4fac=0. _d 0
182          ENDIF
183    
184  C     - Viscosity  C     - Viscosity
185        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
186         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
187          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
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=rA(i,j,bi,bj)           IF (useAreaViscLength) THEN
192              L2=rA(i,j,bi,bj)
193             ELSE
194              L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
195             ENDIF
196           L3=(L2**1.5)           L3=(L2**1.5)
197           L4=(L2**2)           L4=(L2**2)
198           L5=0.125*(L2**2.5)           L5=(L2**2.5)
199           IF (useAnisotropicViscAGridMax) THEN  
200            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
201       &                       +recip_DYF(I,J,bi,bj)**2) )  
202            L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4           IF (useAreaViscLength) THEN
203       &                       +recip_DYF(I,J,bi,bj)**4)            L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2
204       &                   +8.*((recip_DXF(I,J,bi,bj)           ELSE
205       &                        *recip_DYF(I,J,bi,bj))**2) )            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
206         &                            +recip_DYF(I,J,bi,bj)**4)
207         &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
208         &                             *recip_DYF(I,J,bi,bj))**2) )
209           ENDIF           ENDIF
210    
211           IF (useFullLeith) THEN  C Velocity Reynolds Scale
212             IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
213               Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
214             ELSE
215               Uscl=0.
216             ENDIF
217             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
218               U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
219             ELSE
220               U4scl=0.
221             ENDIF
222    
223             IF (useFullLeith.and.calcleith) THEN
224  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
225            grdVrt=0.25*(            grdVrt=0.25 _d 0*(
226       &     ((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
227       &     +((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
228       &     +((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))
229       &     +((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
230         &     +((vort3(i+1,j+1)-vort3(i+1,j))
231         &               *recip_DYG(i+1,j,bi,bj))**2)
232    
233  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
234  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
235             grdDiv=0.25*(            grdDiv=0.25 _d 0*(
236       &      ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2
237       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2
238       &      +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2       &     +((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2
239       &      +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2)
240    
241             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_DLth(i,j)=
242       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
243              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_DLth(i,j)=0.125 _d 0*
244             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
245              Alth2=0. _d 0            viscAh_DLthd(i,j)=
246             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
247             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_DLthd(i,j)=0.125 _d 0*
248       &          .NE. 0. ) THEN       &     sqrt(viscC4leithD**2*grdDiv)*L5
249              Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5           ELSEIF (calcleith) THEN
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
250  C but this approximation will work on cube  C but this approximation will work on cube
251  c (and differs by as much as 4X)  c (and differs by as much as 4X)
252             grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
253             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
254       &      abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
255             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
256       &      abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))       &     abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))
257             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
258       &      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)))
259    
260             grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))
261             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
262       &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))       &     abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj)))
263             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
264       &      abs((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj)))       &     abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)))
265             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
266       &      abs((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj)))       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))
267    
268  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
269             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            viscAh_Dlth(i,j)=
270             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3
271           ENDIF            viscA4_Dlth(i,j)=0.125 _d 0*
272         &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5
273           IF (smag2fac.NE.0.) THEN            viscAh_DlthD(i,j)=
274            Asmag2=smag2fac*L2       &      ((viscC2leithD*grdDiv))*L3
275       &           *sqrt(tension(i,j)**2            viscA4_DlthD(i,j)=0.125 _d 0*
276       &           +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2       &      ((viscC4leithD*grdDiv))*L5
      &                 +strain(i-1, j )**2+strain( i ,j-1)**2))  
277           ELSE           ELSE
278            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
279              viscA4_Dlth(i,j)=0. _d 0
280              viscAh_DlthD(i,j)=0. _d 0
281              viscA4_DlthD(i,j)=0. _d 0
282           ENDIF           ENDIF
283    
284           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
285            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
286       &             *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
287       &             +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2       &       +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2
288       &                   +strain(i-1, j )**2+strain( i ,j-1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
289              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
290              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
291           ELSE           ELSE
292            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
293              viscA4_DSmg(i,j)=0. _d 0
294           ENDIF           ENDIF
295    
296  C  Harmonic on Div.u points  C  Harmonic on Div.u points
297           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
298           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
299           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
300            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)
301            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
302           ELSE           viscAh_D(i,j)=min(viscAh_DMax(i,j),viscAh_D(i,j))
           IF (vg2Max.GT.0.) THEN  
            AlinMax=vg2Max*L2  
            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))  
           ENDIF  
          ENDIF  
          AlinMin=vg2Min*L2  
          viscAh_D(i,j)=max(AlinMin,viscAh_D(i,j))  
303    
304  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
305           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
306           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
307           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
308              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)
309              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
310           ELSE           viscA4_D(i,j)=min(viscA4_DMax(i,j),viscA4_D(i,j))
          IF (vg4Max.GT.0.) THEN  
             AlinMax=vg4Max*L4  
             viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))  
          ENDIF  
          ENDIF  
          AlinMin=vg4Min*L4  
          viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))  
311    
312  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
313  C These are (powers of) length scales  C These are (powers of) length scales
314           L2=rAz(i,j,bi,bj)           IF (useAreaViscLength) THEN
315              L2=rAz(i,j,bi,bj)
316             ELSE
317              L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
318             ENDIF
319    
320           L3=(L2**1.5)           L3=(L2**1.5)
321           L4=(L2**2)           L4=(L2**2)
322           L5=0.125*(L2**2.5)           L5=(L2**2.5)
          IF (useAnisotropicViscAGridMax) THEN  
          L2rdt=recip_dt/( 2.*(recip_DXV(I,J,bi,bj)**2  
      &                       +recip_DYU(I,J,bi,bj)**2) )  
          L4rdt=recip_dt/( 6.*(recip_DXV(I,J,bi,bj)**4  
      &                       +recip_DYU(I,J,bi,bj)**4)  
      &                   +8.*((recip_DXV(I,J,bi,bj)  
      &                        *recip_DYU(I,J,bi,bj))**2) )  
          ENDIF  
323    
324  C This is the vector magnitude of the vorticity gradient squared           L2rdt=0.25 _d 0*recip_dt*L2
325           IF (useFullLeith) THEN           IF (useAreaViscLength) THEN
326             grdVrt=0.25*(            L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
327       &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2           ELSE
328       &      +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2            L4rdt=recip_dt/
329       &      +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2       &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
330       &      +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)       &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
331             ENDIF
332    
333  C This is the vector magnitude of grad(div.v) squared  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
334             grdDiv=0.25*(           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
335       &       ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
336       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                      +(KE(i-1,j)+KE(i,j-1)) )
337       &      +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2             IF ( keZpt.GT.0. ) THEN
338       &      +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)               Uscl = sqrt(keZpt*L2)*viscAhRe_max
339                 U4scl= sqrt(keZpt)*L3*viscA4Re_max
            IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)  
      &          .NE. 0. ) THEN  
             Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3  
            ELSE  
             Alth2=0. _d 0  
            ENDIF  
            IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)  
      &          .NE. 0. ) THEN  
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
340             ELSE             ELSE
341              Alth4=0. _d 0               Uscl =0.
342                 U4scl=0.
343             ENDIF             ENDIF
344           ELSE           ELSE
345               Uscl =0.
346               U4scl=0.
347             ENDIF
348    
349    C This is the vector magnitude of the vorticity gradient squared
350             IF (useFullLeith.and.calcleith) THEN
351              grdVrt=0.25 _d 0*(
352         &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
353         &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
354         &     +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2
355         &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
356    
357    C This is the vector magnitude of grad(div.v) squared
358              grdDiv=0.25 _d 0*(
359         &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2
360         &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2
361         &     +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2
362         &     +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)
363    
364              viscAh_ZLth(i,j)=
365         &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
366              viscA4_ZLth(i,j)=0.125 _d 0*
367         &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
368              viscAh_ZLthD(i,j)=
369         &     sqrt(viscC2leithD**2*grdDiv)*L3
370              viscA4_ZLthD(i,j)=0.125 _d 0*
371         &     sqrt(viscC4leithD**2*grdDiv)*L5
372    
373             ELSEIF (calcleith) THEN
374  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
375             grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
376             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
377       &       abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
378             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
379       &       abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))       &     abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
380             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
381       &       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)))
382    
383             grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))
384             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
385       &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))
386             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
387       &      abs((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i-1,j,bi,bj)))       &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj)))
388             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
389       &      abs((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i,j-1,bi,bj)))       &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj)))
390    
391  C This if statement is just to prevent bitwise changes when leithd=0            viscAh_ZLth(i,j)=(viscC2leith*grdVrt
392             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3       &                     +(viscC2leithD*grdDiv))*L3
393             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5            viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt
394           ENDIF       &                     +(viscC4leithD*grdDiv))*L5
395              viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3
396           IF (smag2fac.NE.0.) THEN            viscA4_ZLthD(i,j)=0.125 _d 0*((viscC4leithD*grdDiv))*L5
           Asmag2=smag2fac*L2  
      &           *sqrt(strain(i,j)**2  
      &              +0.25*(tension( i , j )**2+tension( i ,j+1)**2  
      &                +tension(i+1, j )**2+tension(i+1,j+1)**2))  
397           ELSE           ELSE
398            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
399              viscA4_ZLth(i,j)=0. _d 0
400              viscAh_ZLthD(i,j)=0. _d 0
401              viscA4_ZLthD(i,j)=0. _d 0
402           ENDIF           ENDIF
403    
404           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
405            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
406       &             *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
407       &             +0.25*(tension( i , j )**2+tension( i ,j+1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
408       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
409           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
410            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
411           ENDIF           ENDIF
412    
413  C  Harmonic on Zeta points  C  Harmonic on Zeta points
414           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
415           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
416           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
417              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)
418              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
419           ELSE           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))
420           IF (vg2Max.GT.0.) THEN  
421              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
422              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
423           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
424           ENDIF           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
425           AlinMin=vg2Min*L2           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)
426           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
427             viscA4_Z(i,j)=min(viscA4_ZMax(i,j),viscA4_Z(i,j))
 C  BiHarmonic on Zeta Points  
          Alin=viscA4Z+vg4*L4+Alth4+Asmag4  
          viscA4_Z(i,j)=min(viscA4Max,Alin)  
          IF (useAnisotropicViscAGridMax) THEN  
             AlinMax=viscA4GridMax*L4rdt  
             viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))  
          ELSE  
          IF (vg4Max.GT.0.) THEN  
             AlinMax=vg4Max*L4  
             viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))  
          ENDIF  
          ENDIF  
          AlinMin=vg4Min*L4  
          viscA4_Z(i,j)=max(AlinMin,viscA4_Z(i,j))  
428          ENDDO          ENDDO
429         ENDDO         ENDDO
430        ELSE        ELSE
# Line 337  C  BiHarmonic on Zeta Points Line 444  C  BiHarmonic on Zeta Points
444         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
445         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
446         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
447    
448           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
449           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
450           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
451           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
452    
453           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
454           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
455           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
456           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
457    
458           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
459           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
460           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
461           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
462    
463           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
464         &   ,k,1,2,bi,bj,myThid)
465           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
466         &   ,k,1,2,bi,bj,myThid)
467           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
468         &   ,k,1,2,bi,bj,myThid)
469           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
470         &   ,k,1,2,bi,bj,myThid)
471    
472           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
473           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
474           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
475           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
476        ENDIF        ENDIF
477  #endif  #endif
478    
479        RETURN        RETURN
480        END        END
481    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22