/[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.6 by baylor, Tue Sep 20 21:01:30 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/2/viscAhRemax
39    C
40    C     viscA4Remax is min value for grid point biharmonic Reynolds num
41    C      biharmonic viscosity>sqrt(2*KE)*L**3/16/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
62    C     viscA4Remax>=1
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) viscAhre_max=viscAhremax
140    
141        biharmonic=        biharmonic=
142       &      (viscA4.NE.0.)       &      (viscA4.NE.0.)
143       &  .OR.(viscA4D.NE.0.)       &  .OR.(viscA4D.NE.0.)
# Line 71  C     == Local variables == Line 147  C     == Local variables ==
147       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
148       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
149    
150          IF (biharmonic) viscA4re_max=viscA4remax
151    
152          calcleith=
153         &      (viscC2leith.NE.0.)
154         &  .OR.(viscC2leithD.NE.0.)
155         &  .OR.(viscC4leith.NE.0.)
156         &  .OR.(viscC4leithD.NE.0.)
157    
158          calcsmag=
159         &      (viscC2smag.NE.0.)
160         &  .OR.(viscC4smag.NE.0.)
161    
162        IF (deltaTmom.NE.0.) THEN        IF (deltaTmom.NE.0.) THEN
163         recip_dt=1./deltaTmom         recip_dt=1./deltaTmom
164        ELSE        ELSE
165         recip_dt=0.         recip_dt=0.
166        ENDIF        ENDIF
167    
168        vg2=viscAhGrid*recip_dt        IF (calcsmag) THEN
169        vg2Min=viscAhGridMin*recip_dt          smag2fac=(viscC2smag/pi)**2
170        vg2Max=viscAhGridMax*recip_dt          smag4fac=0.125*(viscC4smag/pi)**2
171        vg4=viscA4Grid*recip_dt        ENDIF
       vg4Min=viscA4GridMin*recip_dt  
       vg4Max=viscA4GridMax*recip_dt  
   
       smag2fac=(viscC2smag/pi)**2  
       smag4fac=0.125*(viscC4smag/pi)**2  
172    
173  C     - Viscosity  C     - Viscosity
174        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
175         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
176          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
177  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
178    
179  C These are (powers of) length scales  C These are (powers of) length scales
180           L2=rA(i,j,bi,bj)           L2=2./((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
181           L3=(L2**1.5)           L3=(L2**1.5)
182           L4=(L2**2)           L4=(L2**2)
183           L5=0.125*(L2**2.5)           L5=(L2**2.5)
184           IF (useAnisotropicViscAGridMax) THEN  
185            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2           L2rdt=0.25*recip_dt*L2
186       &                       +recip_DYF(I,J,bi,bj)**2) )  
187            L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4           L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4
188       &                       +recip_DYF(I,J,bi,bj)**4)       &                       +recip_DYF(I,J,bi,bj)**4)
189       &                   +8.*((recip_DXF(I,J,bi,bj)       &                   +8.*((recip_DXF(I,J,bi,bj)
190       &                        *recip_DYF(I,J,bi,bj))**2) )       &                        *recip_DYF(I,J,bi,bj))**2) )
          ENDIF  
191    
192           IF (useFullLeith) THEN  C Velocity Reynolds Scale
193             Uscl=sqrt(KE(i,j)*L2*0.5)/viscAhRe_max
194             U4scl=0.125*L2*Uscl/viscA4Re_max
195    
196             IF (useFullLeith.and.calcleith) THEN
197  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
198            grdVrt=0.25*(            grdVrt=0.25*(
199       &     ((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
# Line 116  C This is the vector magnitude of the vo Line 203  C This is the vector magnitude of the vo
203    
204  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
205  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
206             grdDiv=0.25*(            grdDiv=0.25*(
207       &      ((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
208       &      +((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
209       &      +((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
210       &      +((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)
211    
212             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_DLth(i,j)=
213       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
214              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_DLth(i,j)=
215             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
216              Alth2=0. _d 0            viscAh_DLthd(i,j)=
217             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
218             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_DLthd(i,j)=
219       &          .NE. 0. ) THEN       &     sqrt(viscC4leithD**2*grdDiv)*L5
220              Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5           ELSEIF (calcleith) THEN
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
221  C but this approximation will work on cube  C but this approximation will work on cube
222  c (and differs by as much as 4X)  c (and differs by as much as 4X)
223             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))
224             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
225       &      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)))
226             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
227       &      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)))
228             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
229       &      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)))
230    
231             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))
232             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
233       &      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)))
234             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
235       &      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)))
236             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
237       &      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)))
238    
239  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
240             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            viscAh_Dlth(i,j)=
241             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3
242           ENDIF            viscA4_Dlth(i,j)=0.125*
243         &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5
244           IF (smag2fac.NE.0.) THEN            viscAh_DlthD(i,j)=
245            Asmag2=smag2fac*L2       &      ((viscC2leithD*grdDiv))*L3
246       &           *sqrt(tension(i,j)**2            viscA4_DlthD(i,j)=0.125*
247       &           +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))  
248           ELSE           ELSE
249            Asmag2=0d0            viscAh_Dlth(i,j)=0d0
250              viscA4_Dlth(i,j)=0d0
251              viscAh_DlthD(i,j)=0d0
252              viscA4_DlthD(i,j)=0d0
253           ENDIF           ENDIF
254    
255           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
256            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
257       &             *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
258       &             +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2       &       +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2
259       &                   +strain(i-1, j )**2+strain( i ,j-1)**2))       &              +strain(i  , j )**2+strain(i+1,j+1)**2))
260              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
261              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
262           ELSE           ELSE
263            Asmag4=0d0            viscAh_DSmg(i,j)=0d0
264              viscA4_DSmg(i,j)=0d0
265           ENDIF           ENDIF
266    
267  C  Harmonic on Div.u points  C  Harmonic on Div.u points
268           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
269           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
270           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
271            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)
272            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
273           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))  
274    
275  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
276           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
277           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
278           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
279              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)
280              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
281           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))  
282    
283  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
284  C These are (powers of) length scales  C These are (powers of) length scales
285           L2=rAz(i,j,bi,bj)           L2=2./((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
286           L3=(L2**1.5)           L3=(L2**1.5)
287           L4=(L2**2)           L4=(L2**2)
288           L5=0.125*(L2**2.5)           L5=(L2**2.5)
289           IF (useAnisotropicViscAGridMax) THEN  
290           L2rdt=recip_dt/( 2.*(recip_DXV(I,J,bi,bj)**2           L2rdt=0.25*recip_dt*L2
291       &                       +recip_DYU(I,J,bi,bj)**2) )           L4rdt=recip_dt/
292           L4rdt=recip_dt/( 6.*(recip_DXV(I,J,bi,bj)**4       &     ( 6.*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)
293       &                       +recip_DYU(I,J,bi,bj)**4)       &      +8.*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))
294       &                   +8.*((recip_DXV(I,J,bi,bj)  
295       &                        *recip_DYU(I,J,bi,bj))**2) )  C Velocity Reynolds Scale
296           ENDIF           Uscl=sqrt((KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))*L2*0.125)/
297         &         viscAhRe_max
298             U4scl=0.125*L2*Uscl/viscA4Re_max
299    
300  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
301           IF (useFullLeith) THEN           IF (useFullLeith.and.calcleith) THEN
302             grdVrt=0.25*(            grdVrt=0.25*(
303       &      ((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
304       &      +((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
305       &      +((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
306       &      +((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)
307    
308  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
309             grdDiv=0.25*(            grdDiv=0.25*(
310       &       ((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
311       &      +((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
312       &      +((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
313       &      +((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)
314    
315             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_ZLth(i,j)=
316       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
317              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_ZLth(i,j)=
318             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
319              Alth2=0. _d 0            viscAh_ZLthD(i,j)=
320             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
321             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_ZLthD(i,j)=
322       &          .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  
323    
324           IF (smag2fac.NE.0.) THEN           ELSEIF (calcleith) THEN
325            Asmag2=smag2fac*L2  C but this approximation will work on cube (and differs by 4X)
326       &           *sqrt(strain(i,j)**2            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
327       &              +0.25*(tension( i , j )**2+tension( i ,j+1)**2            grdVrt=max(grdVrt,
328       &                +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)))
329              grdVrt=max(grdVrt,
330         &     abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
331              grdVrt=max(grdVrt,
332         &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
333    
334              grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))
335              grdDiv=max(grdDiv,
336         &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))
337              grdDiv=max(grdDiv,
338         &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXG(i,j-1,bi,bj)))
339              grdDiv=max(grdDiv,
340         &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYG(i-1,j,bi,bj)))
341    
342              viscAh_ZLth(i,j)=(viscC2leith*grdVrt
343         &                     +(viscC2leithD*grdDiv))*L3
344              viscA4_ZLth(i,j)=(viscC4leith*grdVrt
345         &                     +(viscC4leithD*grdDiv))*L5
346              viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3
347              viscA4_ZLthD(i,j)=((viscC4leithD*grdDiv))*L5
348           ELSE           ELSE
349            Asmag2=0d0            viscAh_ZLth(i,j)=0d0
350              viscA4_ZLth(i,j)=0d0
351              viscAh_ZLthD(i,j)=0d0
352              viscA4_ZLthD(i,j)=0d0
353           ENDIF           ENDIF
354    
355           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
356            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
357       &             *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
358       &             +0.25*(tension( i , j )**2+tension( i ,j+1)**2       &        +0.25*(tension( i , j )**2+tension( i ,j-1)**2
359       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &              +tension(i-1, j )**2+tension(i-1,j-1)**2))
360           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
361            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
362           ENDIF           ENDIF
363    
364  C  Harmonic on Zeta points  C  Harmonic on Zeta points
365           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
366           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
367           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
368              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)
369              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
370           ELSE           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))
371           IF (vg2Max.GT.0.) THEN  
372              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
373              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
374           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
375           ENDIF           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
376           AlinMin=vg2Min*L2           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)
377           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
378             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))  
379          ENDDO          ENDDO
380         ENDDO         ENDDO
381        ELSE        ELSE
# Line 337  C  BiHarmonic on Zeta Points Line 395  C  BiHarmonic on Zeta Points
395         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
396         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
397         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
398    
399           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
400           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
401           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
402           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
403    
404           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
405           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
406           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
407           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
408    
409           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
410           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
411           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
412           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
413    
414           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD',k,1,2,bi,bj,myThid)
415           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD',k,1,2,bi,bj,myThid)
416           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD',k,1,2,bi,bj,myThid)
417           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD',k,1,2,bi,bj,myThid)
418    
419           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
420           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
421           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
422           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
423        ENDIF        ENDIF
424  #endif  #endif
425    
426        RETURN        RETURN
427        END        END
428    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22