/[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.19 by baylor, Mon Oct 10 19:49:48 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/pi)**6*grad(Vort3)**2
21    C             +(viscC2leithD/pi)**6*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/pi)**6*grad(Vort3)**2
28    C                        +(viscC4leithD/pi)**6*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=1-3
55    C     viscC2LeithD=1-3
56    C     viscC4Leith=1-3
57    C     viscC4LeithD=1.5-3
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 leith2fac, leith4fac
95        _RL lenA2, lenAz2        _RL leithD2fac, leithD4fac
96        _RL Alin,Alth2,Alth4,grdVrt,grdDiv        _RL viscAhRe_max, viscA4Re_max
97        _RL vg2,vg4,vg4Min,vg4Max        _RL Alin,grdVrt,grdDiv, keZpt
98        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
99          _RL Uscl,U4scl
100          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104          _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105          _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106          _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107          _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108          _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109          _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110          _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112          _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115          _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116          _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117          _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118          _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119          _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122          LOGICAL calcLeith,calcSmag
123    
124        useVariableViscosity=        useVariableViscosity=
125       &      (viscAhGrid.NE.0.)       &      (viscAhGrid.NE.0.)
# Line 62  C     == Local variables == Line 140  C     == Local variables ==
140       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
141       &  .OR.(viscC2smag.NE.0.)       &  .OR.(viscC2smag.NE.0.)
142    
143          IF ((harmonic).and.(viscAhremax.ne.0.)) THEN
144            viscAhre_max=sqrt(2. _d 0)/viscAhRemax
145          ELSE
146            viscAhre_max=0. _d 0
147          ENDIF
148    
149        biharmonic=        biharmonic=
150       &      (viscA4.NE.0.)       &      (viscA4.NE.0.)
151       &  .OR.(viscA4D.NE.0.)       &  .OR.(viscA4D.NE.0.)
# Line 71  C     == Local variables == Line 155  C     == Local variables ==
155       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
156       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
157    
158          IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN
159            viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax
160          ELSE
161            viscA4re_max=0. _d 0
162          ENDIF
163    
164          calcleith=
165         &      (viscC2leith.NE.0.)
166         &  .OR.(viscC2leithD.NE.0.)
167         &  .OR.(viscC4leith.NE.0.)
168         &  .OR.(viscC4leithD.NE.0.)
169    
170          calcsmag=
171         &      (viscC2smag.NE.0.)
172         &  .OR.(viscC4smag.NE.0.)
173    
174        IF (deltaTmom.NE.0.) THEN        IF (deltaTmom.NE.0.) THEN
175         recip_dt=1./deltaTmom         recip_dt=1. _d 0/deltaTmom
176        ELSE        ELSE
177         recip_dt=0.         recip_dt=0. _d 0
178        ENDIF        ENDIF
179    
180        vg2=viscAhGrid*recip_dt        IF (calcsmag) THEN
181        vg2Min=viscAhGridMin*recip_dt          smag2fac=(viscC2smag/pi)**2
182        vg2Max=viscAhGridMax*recip_dt          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
183        vg4=viscA4Grid*recip_dt        ELSE
184        vg4Min=viscA4GridMin*recip_dt          smag2fac=0. _d 0
185        vg4Max=viscA4GridMax*recip_dt          smag4fac=0. _d 0
186          ENDIF
187    
188        smag2fac=(viscC2smag/pi)**2        IF (calcleith) THEN
189        smag4fac=0.125*(viscC4smag/pi)**2          IF (useFullLeith) THEN
190             leith2fac =(viscC2leith /pi)**6
191             leithD2fac=(viscC2leithD/pi)**6
192             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
193             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
194            ELSE
195             leith2fac =(viscC2leith /pi)**3
196             leithD2fac=(viscC2leithD/pi)**3
197             leith4fac =0.125 _d 0*(viscC4leith /pi)**3
198             leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
199            ENDIF
200          ELSE
201            leith2fac=0. _d 0
202            leith4fac=0. _d 0
203            leithD2fac=0. _d 0
204            leithD4fac=0. _d 0
205          ENDIF
206    
207  C     - Viscosity  C     - Viscosity
208        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
209    
210    C      horizontal gradient of horizontal divergence:
211           DO j=1-Oly,sNy+Oly
212            DO i=1-Olx,sNx+Olx
213              divDx(i,j) = 0.
214              divDy(i,j) = 0.
215            ENDDO
216           ENDDO
217           IF (calcleith) THEN
218    C-       gradient in x direction:
219    #ifndef ALLOW_AUTODIFF_TAMC
220             IF (useCubedSphereExchange) THEN
221    C        to compute d/dx(hDiv), fill corners with appropriate values:
222               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
223             ENDIF
224    #endif
225             DO j=2-Oly,sNy+Oly-1
226              DO i=2-Olx,sNx+Olx-1
227                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
228              ENDDO
229             ENDDO
230    
231    C-       gradient in y direction:
232    #ifndef ALLOW_AUTODIFF_TAMC
233             IF (useCubedSphereExchange) THEN
234    C        to compute d/dy(hDiv), fill corners with appropriate values:
235               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
236             ENDIF
237    #endif
238             DO j=2-Oly,sNy+Oly-1
239              DO i=2-Olx,sNx+Olx-1
240                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
241              ENDDO
242             ENDDO
243           ENDIF
244    
245         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
246          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
247  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
248    
249  C These are (powers of) length scales  C These are (powers of) length scales
250           L2=rA(i,j,bi,bj)           IF (useAreaViscLength) THEN
251              L2=rA(i,j,bi,bj)
252             ELSE
253              L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
254             ENDIF
255           L3=(L2**1.5)           L3=(L2**1.5)
256           L4=(L2**2)           L4=(L2**2)
257           L5=0.125*(L2**2.5)           L5=(L2**2.5)
258           IF (useAnisotropicViscAGridMax) THEN  
259            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
260       &                       +recip_DYF(I,J,bi,bj)**2) )  
261            L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4           IF (useAreaViscLength) THEN
262       &                       +recip_DYF(I,J,bi,bj)**4)            L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2
263       &                   +8.*((recip_DXF(I,J,bi,bj)           ELSE
264       &                        *recip_DYF(I,J,bi,bj))**2) )            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
265         &                            +recip_DYF(I,J,bi,bj)**4)
266         &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
267         &                             *recip_DYF(I,J,bi,bj))**2) )
268             ENDIF
269    
270    C Velocity Reynolds Scale
271             IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
272               Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
273             ELSE
274               Uscl=0.
275             ENDIF
276             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
277               U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
278             ELSE
279               U4scl=0.
280           ENDIF           ENDIF
281    
282           IF (useFullLeith) THEN           IF (useFullLeith.and.calcleith) THEN
283  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
284            grdVrt=0.25*(            grdVrt=0.25 _d 0*(
285       &     ((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
286       &     +((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
287       &     +((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))
288       &     +((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
289         &     +((vort3(i+1,j+1)-vort3(i+1,j))
290         &               *recip_DYG(i+1,j,bi,bj))**2)
291    
292  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
293  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
294             grdDiv=0.25*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
295       &      ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
296       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
297       &      +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
298       &      +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)  
299              viscAh_DLth(i,j)=
300             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
301       &          .NE. 0. ) THEN            viscA4_DLth(i,j)=
302              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
303             ELSE            viscAh_DLthd(i,j)=
304              Alth2=0. _d 0       &     sqrt(leithD2fac*grdDiv)*L3
305             ENDIF            viscA4_DLthd(i,j)=
306             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)       &     sqrt(leithD4fac*grdDiv)*L5
307       &          .NE. 0. ) THEN           ELSEIF (calcleith) THEN
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
308  C but this approximation will work on cube  C but this approximation will work on cube
309  c (and differs by as much as 4X)  c (and differs by as much as 4X)
310             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))
311             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
312       &      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)))
313             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
314       &      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)))
315             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
316       &      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)))
317    
318             grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
319             grdDiv=max(grdDiv,            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
320       &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i,j))   )
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj)))  
321    
322  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
323             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
324             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
325           ENDIF            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
326              viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
          IF (smag2fac.NE.0.) THEN  
           Asmag2=smag2fac*L2  
      &           *sqrt(tension(i,j)**2  
      &           +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2  
      &                 +strain(i-1, j )**2+strain( i ,j-1)**2))  
327           ELSE           ELSE
328            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
329              viscA4_Dlth(i,j)=0. _d 0
330              viscAh_DlthD(i,j)=0. _d 0
331              viscA4_DlthD(i,j)=0. _d 0
332           ENDIF           ENDIF
333    
334           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
335            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
336       &             *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
337       &             +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
338       &                   +strain(i-1, j )**2+strain( i ,j-1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
339              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
340              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
341           ELSE           ELSE
342            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
343              viscA4_DSmg(i,j)=0. _d 0
344           ENDIF           ENDIF
345    
346  C  Harmonic on Div.u points  C  Harmonic on Div.u points
347           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
348           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
349           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
350            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)
351            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
352           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))  
353    
354  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
355           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
356           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
357           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
358              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)
359              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
360           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))  
361    
362  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
363  C These are (powers of) length scales  C These are (powers of) length scales
364           L2=rAz(i,j,bi,bj)           IF (useAreaViscLength) THEN
365              L2=rAz(i,j,bi,bj)
366             ELSE
367              L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
368             ENDIF
369    
370           L3=(L2**1.5)           L3=(L2**1.5)
371           L4=(L2**2)           L4=(L2**2)
372           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  
373    
374  C This is the vector magnitude of the vorticity gradient squared           L2rdt=0.25 _d 0*recip_dt*L2
375           IF (useFullLeith) THEN           IF (useAreaViscLength) THEN
376             grdVrt=0.25*(            L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
377       &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2           ELSE
378       &      +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2            L4rdt=recip_dt/
379       &      +((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)
380       &      +((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))
381             ENDIF
382    
383  C This is the vector magnitude of grad(div.v) squared  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
384             grdDiv=0.25*(           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
385       &       ((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))
386       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                      +(KE(i-1,j)+KE(i,j-1)) )
387       &      +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2             IF ( keZpt.GT.0. ) THEN
388       &      +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)               Uscl = sqrt(keZpt*L2)*viscAhRe_max
389                 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  
390             ELSE             ELSE
391              Alth2=0. _d 0               Uscl =0.
392             ENDIF               U4scl=0.
            IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)  
      &          .NE. 0. ) THEN  
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
            ELSE  
             Alth4=0. _d 0  
393             ENDIF             ENDIF
394           ELSE           ELSE
395               Uscl =0.
396               U4scl=0.
397             ENDIF
398    
399    C This is the vector magnitude of the vorticity gradient squared
400             IF (useFullLeith.and.calcleith) THEN
401              grdVrt=0.25 _d 0*(
402         &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
403         &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
404         &     +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2
405         &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
406    
407    C This is the vector magnitude of grad(div.v) squared
408              grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
409         &                        + divDx(i,j)*divDx(i,j) )
410         &                     + (divDy(i-1,j)*divDy(i-1,j)
411         &                        + divDy(i,j)*divDy(i,j) )  )
412    
413              viscAh_ZLth(i,j)=
414         &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
415              viscA4_ZLth(i,j)=
416         &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
417              viscAh_ZLthD(i,j)=
418         &     sqrt(leithD2fac*grdDiv)*L3
419              viscA4_ZLthD(i,j)=
420         &     sqrt(leithD4fac*grdDiv)*L5
421    
422             ELSEIF (calcleith) THEN
423  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
424             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))
425             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
426       &       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)))
427             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
428       &       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)))
429             grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
430       &       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)))
431    
432             grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
433             grdDiv=max(grdDiv,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
434       &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
435             grdDiv=max(grdDiv,  
436       &      abs((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i-1,j,bi,bj)))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
437             grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
438       &      abs((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i,j-1,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
439              viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
 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  
   
          IF (smag2fac.NE.0.) THEN  
           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))  
440           ELSE           ELSE
441            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
442              viscA4_ZLth(i,j)=0. _d 0
443              viscAh_ZLthD(i,j)=0. _d 0
444              viscA4_ZLthD(i,j)=0. _d 0
445           ENDIF           ENDIF
446    
447           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
448            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
449       &             *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
450       &             +0.25*(tension( i , j )**2+tension( i ,j+1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
451       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
452           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
453            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
454           ENDIF           ENDIF
455    
456  C  Harmonic on Zeta points  C  Harmonic on Zeta points
457           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
458           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
459           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
460              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)
461              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
462           ELSE           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))
463           IF (vg2Max.GT.0.) THEN  
464              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
465              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
466           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
467           ENDIF           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
468           AlinMin=vg2Min*L2           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)
469           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
470             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))  
471          ENDDO          ENDDO
472         ENDDO         ENDDO
473        ELSE        ELSE
# Line 337  C  BiHarmonic on Zeta Points Line 487  C  BiHarmonic on Zeta Points
487         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
488         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
489         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
490    
491           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
492           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
493           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
494           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
495    
496           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
497           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
498           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
499           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
500    
501           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
502           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
503           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
504           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
505    
506           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
507         &   ,k,1,2,bi,bj,myThid)
508           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
509         &   ,k,1,2,bi,bj,myThid)
510           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
511         &   ,k,1,2,bi,bj,myThid)
512           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
513         &   ,k,1,2,bi,bj,myThid)
514    
515           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
516           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
517           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
518           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
519        ENDIF        ENDIF
520  #endif  #endif
521    
522        RETURN        RETURN
523        END        END
524    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22