/[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.11 by baylor, Thu Sep 22 14:00:59 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,Alinmin,grdVrt,grdDiv
       _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.)
165        vg4Max=viscA4GridMax*recip_dt  
166          calcsmag=
167         &      (viscC2smag.NE.0.)
168         &  .OR.(viscC4smag.NE.0.)
169    
170        smag2fac=(viscC2smag/pi)**2        IF (deltaTmom.NE.0.) THEN
171        smag4fac=0.125*(viscC4smag/pi)**2         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             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
213             U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
214    
215             IF (useFullLeith.and.calcleith) THEN
216  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
217            grdVrt=0.25*(            grdVrt=0.25 _d 0*(
218       &     ((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
219       &     +((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
220       &     +((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))
221       &     +((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
222         &     +((vort3(i+1,j+1)-vort3(i+1,j))
223         &               *recip_DYG(i+1,j,bi,bj))**2)
224    
225  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
226  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
227             grdDiv=0.25*(            grdDiv=0.25 _d 0*(
228       &      ((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
229       &      +((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
230       &      +((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
231       &      +((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)
232    
233             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_DLth(i,j)=
234       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
235              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_DLth(i,j)=0.125 _d 0*
236             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
237              Alth2=0. _d 0            viscAh_DLthd(i,j)=
238             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
239             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_DLthd(i,j)=0.125 _d 0*
240       &          .NE. 0. ) THEN       &     sqrt(viscC4leithD**2*grdDiv)*L5
241              Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5           ELSEIF (calcleith) THEN
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
242  C but this approximation will work on cube  C but this approximation will work on cube
243  c (and differs by as much as 4X)  c (and differs by as much as 4X)
244             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))
245             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
246       &      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)))
247             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
248       &      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)))
249             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
250       &      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)))
251    
252             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))
253             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
254       &      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)))
255             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
256       &      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)))
257             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
258       &      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)))
259    
260  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
261             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            viscAh_Dlth(i,j)=
262             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3
263           ENDIF            viscA4_Dlth(i,j)=0.125 _d 0*
264         &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5
265           IF (smag2fac.NE.0.) THEN            viscAh_DlthD(i,j)=
266            Asmag2=smag2fac*L2       &      ((viscC2leithD*grdDiv))*L3
267       &           *sqrt(tension(i,j)**2            viscA4_DlthD(i,j)=0.125 _d 0*
268       &           +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))  
269           ELSE           ELSE
270            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
271              viscA4_Dlth(i,j)=0. _d 0
272              viscAh_DlthD(i,j)=0. _d 0
273              viscA4_DlthD(i,j)=0. _d 0
274           ENDIF           ENDIF
275    
276           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
277            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
278       &             *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
279       &             +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
280       &                   +strain(i-1, j )**2+strain( i ,j-1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
281              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
282              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
283           ELSE           ELSE
284            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
285              viscA4_DSmg(i,j)=0. _d 0
286           ENDIF           ENDIF
287    
288  C  Harmonic on Div.u points  C  Harmonic on Div.u points
289           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
290           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
291           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
292            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)
293            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
294           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))  
295    
296  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
297           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
298           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
299           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
300              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)
301              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
302           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))  
303    
304  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
305  C These are (powers of) length scales  C These are (powers of) length scales
306           L2=rAz(i,j,bi,bj)           IF (useAreaViscLength) THEN
307              L2=RaZ(i,j,bi,bj)
308             ELSE
309              L2=2._d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
310             ENDIF
311    
312           L3=(L2**1.5)           L3=(L2**1.5)
313           L4=(L2**2)           L4=(L2**2)
314           L5=0.125*(L2**2.5)           L5=(L2**2.5)
315           IF (useAnisotropicViscAGridMax) THEN  
316           L2rdt=recip_dt/( 2.*(recip_DXV(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
317       &                       +recip_DYU(I,J,bi,bj)**2) )           IF (useAreaViscLength) THEN
318           L4rdt=recip_dt/( 6.*(recip_DXV(I,J,bi,bj)**4            L4rdt=0.125 _d 0*recip_dt*RaZ(i,j,bi,bj)**2
319       &                       +recip_DYU(I,J,bi,bj)**4)           ELSE
320       &                   +8.*((recip_DXV(I,J,bi,bj)            L4rdt=recip_dt/
321       &                        *recip_DYU(I,J,bi,bj))**2) )       &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
322         &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
323           ENDIF           ENDIF
324    
325    C Velocity Reynolds Scale
326             Uscl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))
327         &         *L2)*viscAhRe_max
328             U4scl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1)))
329         &         *L3*viscA4Re_max
330    
331  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
332           IF (useFullLeith) THEN           IF (useFullLeith.and.calcleith) THEN
333             grdVrt=0.25*(            grdVrt=0.25 _d 0*(
334       &      ((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
335       &      +((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
336       &      +((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
337       &      +((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)
338    
339  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
340             grdDiv=0.25*(            grdDiv=0.25 _d 0*(
341       &       ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2
342       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2
343       &      +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2       &     +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2
344       &      +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)       &     +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)
345    
346             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_ZLth(i,j)=
347       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
348              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_ZLth(i,j)=0.125 _d 0*
349             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
350              Alth2=0. _d 0            viscAh_ZLthD(i,j)=
351             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
352             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_ZLthD(i,j)=0.125 _d 0*
353       &          .NE. 0. ) THEN       &     sqrt(viscC4leithD**2*grdDiv)*L5
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
 C but this approximation will work on cube (and differs by 4X)  
            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))  
            grdVrt=max(grdVrt,  
      &       abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))  
            grdVrt=max(grdVrt,  
      &       abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))  
            grdVrt=max(grdVrt,  
      &       abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))  
   
            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i-1,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i,j-1,bi,bj)))  
   
 C This if statement is just to prevent bitwise changes when leithd=0  
            Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3  
            Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
          ENDIF  
354    
355           IF (smag2fac.NE.0.) THEN           ELSEIF (calcleith) THEN
356            Asmag2=smag2fac*L2  C but this approximation will work on cube (and differs by 4X)
357       &           *sqrt(strain(i,j)**2            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
358       &              +0.25*(tension( i , j )**2+tension( i ,j+1)**2            grdVrt=max(grdVrt,
359       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
360              grdVrt=max(grdVrt,
361         &     abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
362              grdVrt=max(grdVrt,
363         &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
364    
365              grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))
366              grdDiv=max(grdDiv,
367         &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))
368              grdDiv=max(grdDiv,
369         &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj)))
370              grdDiv=max(grdDiv,
371         &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj)))
372    
373              viscAh_ZLth(i,j)=(viscC2leith*grdVrt
374         &                     +(viscC2leithD*grdDiv))*L3
375              viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt
376         &                     +(viscC4leithD*grdDiv))*L5
377              viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3
378              viscA4_ZLthD(i,j)=0.125 _d 0*((viscC4leithD*grdDiv))*L5
379           ELSE           ELSE
380            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
381              viscA4_ZLth(i,j)=0. _d 0
382              viscAh_ZLthD(i,j)=0. _d 0
383              viscA4_ZLthD(i,j)=0. _d 0
384           ENDIF           ENDIF
385    
386           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
387            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
388       &             *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
389       &             +0.25*(tension( i , j )**2+tension( i ,j+1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
390       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
391           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
392            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
393           ENDIF           ENDIF
394    
395  C  Harmonic on Zeta points  C  Harmonic on Zeta points
396           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
397           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
398           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
399              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)
400              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
401           ELSE           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))
402           IF (vg2Max.GT.0.) THEN  
403              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
404              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
405           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
406           ENDIF           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
407           AlinMin=vg2Min*L2           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)
408           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
409             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))  
410          ENDDO          ENDDO
411         ENDDO         ENDDO
412        ELSE        ELSE
# Line 337  C  BiHarmonic on Zeta Points Line 426  C  BiHarmonic on Zeta Points
426         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
427         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
428         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
429    
430           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
431           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
432           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
433           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
434    
435           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
436           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
437           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
438           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
439    
440           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
441           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
442           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
443           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
444    
445           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
446         &   ,k,1,2,bi,bj,myThid)
447           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
448         &   ,k,1,2,bi,bj,myThid)
449           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
450         &   ,k,1,2,bi,bj,myThid)
451           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
452         &   ,k,1,2,bi,bj,myThid)
453    
454           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
455           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
456           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
457           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
458        ENDIF        ENDIF
459  #endif  #endif
460    
461        RETURN        RETURN
462        END        END
463    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22