/[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.3 by baylor, Tue Sep 20 13:09:10 2005 UTC revision 1.10 by jmc, Thu Sep 22 00:21:23 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          IF (deltaTmom.NE.0.) THEN
171           recip_dt=1. _d 0/deltaTmom
172          ELSE
173           recip_dt=0. _d 0
174          ENDIF
175    
176        smag2fac=(viscC2smag/pi)**2        IF (calcsmag) THEN
177        smag4fac=0.125*(viscC4smag/pi)**2          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)           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=0.125*(L2**2.5)           L5=(L2**2.5)
195           IF (useAnisotropicViscAGridMax) THEN  
196            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
197       &                       +recip_DYF(I,J,bi,bj)**2) )  
198            L4rdt=recip_dt/( 6.*(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       &                   +8.*((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) )
          ENDIF  
202    
203           IF (useFullLeith) THEN  C Velocity Reynolds Scale
204             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
205             U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
206    
207             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.25*(            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))*recip_DXG(i,j+1,bi,bj))**2       &     +((vort3(i+1,j+1)-vort3(i,j+1))
213       &     +((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
214         &     +((vort3(i+1,j+1)-vort3(i+1,j))
215         &               *recip_DYG(i+1,j,bi,bj))**2)
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.25*(            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
223       &      +((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)
224    
225             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_DLth(i,j)=
226       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
227              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_DLth(i,j)=0.125 _d 0*
228             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
229              Alth2=0. _d 0            viscAh_DLthd(i,j)=
230             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
231             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_DLthd(i,j)=0.125 _d 0*
232       &          .NE. 0. ) THEN       &     sqrt(viscC4leithD**2*grdDiv)*L5
233              Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5           ELSEIF (calcleith) THEN
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
234  C but this approximation will work on cube  C but this approximation will work on cube
235  c (and differs by as much as 4X)  c (and differs by as much as 4X)
236             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))
237             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
238       &      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)))
239             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
240       &      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)))
241             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
242       &      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)))
243    
244             grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))
245             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
246       &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj)))       &     abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj)))
247             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
248       &      abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)))       &     abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)))
249             grdDiv=max(grdDiv,            grdDiv=max(grdDiv,
250       &      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)))
251    
252  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
253             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            viscAh_Dlth(i,j)=
254             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3
255           ENDIF            viscA4_Dlth(i,j)=0.125 _d 0*
256         &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5
257           IF (smag2fac.NE.0.) THEN            viscAh_DlthD(i,j)=
258            Asmag2=smag2fac*L2       &      ((viscC2leithD*grdDiv))*L3
259       &           *sqrt(tension(i,j)**2            viscA4_DlthD(i,j)=0.125 _d 0*
260       &           +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2       &      ((viscC4leithD*grdDiv))*L5
      &                 +strain(i  , j )**2+strain(i+1,j+1)**2))  
261           ELSE           ELSE
262            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
263              viscA4_Dlth(i,j)=0. _d 0
264              viscAh_DlthD(i,j)=0. _d 0
265              viscA4_DlthD(i,j)=0. _d 0
266           ENDIF           ENDIF
267    
268           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
269            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
270       &             *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
271       &           +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
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)
274              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
275           ELSE           ELSE
276            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
277              viscA4_DSmg(i,j)=0. _d 0
278           ENDIF           ENDIF
279    
280  C  Harmonic on Div.u points  C  Harmonic on Div.u points
281           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
282           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
283           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
284            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)
285            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
286           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))  
287    
288  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
289           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
290           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
291           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
292              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)
293              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
294           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))  
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=rAz(i,j,bi,bj)           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=0.125*(L2**2.5)           L5=(L2**2.5)
302           IF (useAnisotropicViscAGridMax) THEN  
303           L2rdt=recip_dt/( 2.*(recip_DXV(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
304       &                       +recip_DYU(I,J,bi,bj)**2) )           L4rdt=recip_dt/
305           L4rdt=recip_dt/( 6.*(recip_DXV(I,J,bi,bj)**4       &     ( 6. _d 0*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)
306       &                       +recip_DYU(I,J,bi,bj)**4)       &      +8. _d 0*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))
307       &                   +8.*((recip_DXV(I,J,bi,bj)  
308       &                        *recip_DYU(I,J,bi,bj))**2) )  C Velocity Reynolds Scale
309           ENDIF           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
311             U4scl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1)))
312         &         *L3*viscA4Re_max
313    
314  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
315           IF (useFullLeith) THEN           IF (useFullLeith.and.calcleith) THEN
316             grdVrt=0.25*(            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.25*(            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
327       &      +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)       &     +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)
328    
329             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)            viscAh_ZLth(i,j)=
330       &          .NE. 0. ) THEN       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
331              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3            viscA4_ZLth(i,j)=0.125 _d 0*
332             ELSE       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
333              Alth2=0. _d 0            viscAh_ZLthD(i,j)=
334             ENDIF       &     sqrt(viscC2leithD**2*grdDiv)*L3
335             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)            viscA4_ZLthD(i,j)=0.125 _d 0*
336       &          .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,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXG(i,j-1,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYG(i-1,j,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  
337    
338           IF (smag2fac.NE.0.) THEN           ELSEIF (calcleith) THEN
339            Asmag2=smag2fac*L2  C but this approximation will work on cube (and differs by 4X)
340       &           *sqrt(strain(i,j)**2            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
341       &              +0.25*(tension( i , j )**2+tension( i ,j-1)**2            grdVrt=max(grdVrt,
342       &                    +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)))
343              grdVrt=max(grdVrt,
344         &     abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
345              grdVrt=max(grdVrt,
346         &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
347    
348              grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))
349              grdDiv=max(grdDiv,
350         &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))
351              grdDiv=max(grdDiv,
352         &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj)))
353              grdDiv=max(grdDiv,
354         &     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
357         &                     +(viscC2leithD*grdDiv))*L3
358              viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt
359         &                     +(viscC4leithD*grdDiv))*L5
360              viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3
361              viscA4_ZLthD(i,j)=0.125 _d 0*((viscC4leithD*grdDiv))*L5
362           ELSE           ELSE
363            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
364              viscA4_ZLth(i,j)=0. _d 0
365              viscAh_ZLthD(i,j)=0. _d 0
366              viscA4_ZLthD(i,j)=0. _d 0
367           ENDIF           ENDIF
368    
369           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
370            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
371       &             *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
372       &              +0.25*(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           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
375            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
376           ENDIF           ENDIF
377    
378  C  Harmonic on Zeta points  C  Harmonic on Zeta points
379           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
380           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
381           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
382              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)
383              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
384           ELSE           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))
385           IF (vg2Max.GT.0.) THEN  
386              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
387              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
388           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
389           ENDIF           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
390           AlinMin=vg2Min*L2           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)
391           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
392             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))  
393          ENDDO          ENDDO
394         ENDDO         ENDDO
395        ELSE        ELSE
# Line 338  C  BiHarmonic on Zeta Points Line 409  C  BiHarmonic on Zeta Points
409         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
410         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
411         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
412    
413           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
414           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
415           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
416           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
417    
418           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
419           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
420           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
421           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
422    
423           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
424           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
425           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
426           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
427    
428           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
429         &   ,k,1,2,bi,bj,myThid)
430           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
431         &   ,k,1,2,bi,bj,myThid)
432           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
433         &   ,k,1,2,bi,bj,myThid)
434           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
435         &   ,k,1,2,bi,bj,myThid)
436    
437           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
438           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
439           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
440           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
441        ENDIF        ENDIF
442  #endif  #endif
443    
444        RETURN        RETURN
445        END        END
446    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22