/[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.24 by mlosch, Thu Nov 23 09:49:36 2006 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"
72  #include "GRID.h"  #include "GRID.h"
73  #include "EEPARAMS.h"  #include "EEPARAMS.h"
74  #include "PARAMS.h"  #include "PARAMS.h"
75    #ifdef ALLOW_NONHYDROSTATIC
76    #include "NH_VARS.h"
77    #endif
78    
79  C     == Routine arguments ==  C     == Routine arguments ==
80        INTEGER bi,bj,k        INTEGER bi,bj,k
# Line 35  C     == Routine arguments == Line 93  C     == Routine arguments ==
93    
94  C     == Local variables ==  C     == Local variables ==
95        INTEGER I,J        INTEGER I,J
96        _RL ASmag2, ASmag4, smag2fac, smag4fac        INTEGER kp1
97        _RL vg2Min, vg2Max, AlinMax, AlinMin        _RL smag2fac, smag4fac
98        _RL lenA2, lenAz2        _RL leith2fac, leith4fac
99        _RL Alin,Alth2,Alth4,grdVrt,grdDiv        _RL leithD2fac, leithD4fac
100        _RL vg2,vg4,vg4Min,vg4Max        _RL viscAhRe_max, viscA4Re_max
101          _RL Alin,grdVrt,grdDiv, keZpt
102        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
103          _RL Uscl,U4scl
104          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108          _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109          _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110          _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112          _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115          _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116          _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117          _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118          _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119          _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122          _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123          _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124          _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127          _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128          LOGICAL calcLeith,calcSmag
129    
130        useVariableViscosity=        useVariableViscosity=
131       &      (viscAhGrid.NE.0.)       &      (viscAhGrid.NE.0.)
# Line 62  C     == Local variables == Line 146  C     == Local variables ==
146       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
147       &  .OR.(viscC2smag.NE.0.)       &  .OR.(viscC2smag.NE.0.)
148    
149          IF ((harmonic).and.(viscAhremax.ne.0.)) THEN
150            viscAhre_max=sqrt(2. _d 0)/viscAhRemax
151          ELSE
152            viscAhre_max=0. _d 0
153          ENDIF
154    
155        biharmonic=        biharmonic=
156       &      (viscA4.NE.0.)       &      (viscA4.NE.0.)
157       &  .OR.(viscA4D.NE.0.)       &  .OR.(viscA4D.NE.0.)
# Line 71  C     == Local variables == Line 161  C     == Local variables ==
161       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
162       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
163    
164          IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN
165            viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax
166          ELSE
167            viscA4re_max=0. _d 0
168          ENDIF
169    
170          calcleith=
171         &      (viscC2leith.NE.0.)
172         &  .OR.(viscC2leithD.NE.0.)
173         &  .OR.(viscC4leith.NE.0.)
174         &  .OR.(viscC4leithD.NE.0.)
175    
176          calcsmag=
177         &      (viscC2smag.NE.0.)
178         &  .OR.(viscC4smag.NE.0.)
179    
180        IF (deltaTmom.NE.0.) THEN        IF (deltaTmom.NE.0.) THEN
181         recip_dt=1./deltaTmom         recip_dt=1. _d 0/deltaTmom
182        ELSE        ELSE
183         recip_dt=0.         recip_dt=0. _d 0
184        ENDIF        ENDIF
185    
186        vg2=viscAhGrid*recip_dt        IF (calcsmag) THEN
187        vg2Min=viscAhGridMin*recip_dt          smag2fac=(viscC2smag/pi)**2
188        vg2Max=viscAhGridMax*recip_dt          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
189        vg4=viscA4Grid*recip_dt        ELSE
190        vg4Min=viscA4GridMin*recip_dt          smag2fac=0. _d 0
191        vg4Max=viscA4GridMax*recip_dt          smag4fac=0. _d 0
192          ENDIF
193    
194        smag2fac=(viscC2smag/pi)**2        IF (calcleith) THEN
195        smag4fac=0.125*(viscC4smag/pi)**2          IF (useFullLeith) THEN
196             leith2fac =(viscC2leith /pi)**6
197             leithD2fac=(viscC2leithD/pi)**6
198             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
199             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
200            ELSE
201             leith2fac =(viscC2leith /pi)**3
202             leithD2fac=(viscC2leithD/pi)**3
203             leith4fac =0.125 _d 0*(viscC4leith /pi)**3
204             leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
205            ENDIF
206          ELSE
207            leith2fac=0. _d 0
208            leith4fac=0. _d 0
209            leithD2fac=0. _d 0
210            leithD4fac=0. _d 0
211          ENDIF
212    
213    #ifdef ALLOW_AUTODIFF_TAMC
214          IF ( calcLeith .OR. calcSmag ) THEN
215           STOP 'calcLeith or calcSmag not implemented for ADJOINT'
216          ENDIF
217    #endif
218          DO j=1-Oly,sNy+Oly
219            DO i=1-Olx,sNx+Olx
220              viscAh_D(i,j)=viscAhD
221              viscAh_Z(i,j)=viscAhZ
222              viscA4_D(i,j)=viscA4D
223              viscA4_Z(i,j)=viscA4Z
224    c
225              visca4_zsmg(i,j) = 0. _d 0
226              viscah_zsmg(i,j) = 0. _d 0
227    c
228              viscAh_Dlth(i,j) = 0. _d 0
229              viscA4_Dlth(i,j) = 0. _d 0
230              viscAh_DlthD(i,j)= 0. _d 0
231              viscA4_DlthD(i,j)= 0. _d 0
232    c
233              viscAh_DSmg(i,j) = 0. _d 0
234              viscA4_DSmg(i,j) = 0. _d 0
235    c
236              viscAh_ZLth(i,j) = 0. _d 0
237              viscA4_ZLth(i,j) = 0. _d 0
238              viscAh_ZLthD(i,j)= 0. _d 0
239              viscA4_ZLthD(i,j)= 0. _d 0
240            ENDDO
241          ENDDO
242    
243  C     - Viscosity  C     - Viscosity
244        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
245    
246    C-     Initialise to zero gradient of vorticity & divergence:
247           DO j=1-Oly,sNy+Oly
248            DO i=1-Olx,sNx+Olx
249              divDx(i,j) = 0.
250              divDy(i,j) = 0.
251              vrtDx(i,j) = 0.
252              vrtDy(i,j) = 0.
253            ENDDO
254           ENDDO
255    
256           IF (calcleith) THEN
257    C      horizontal gradient of horizontal divergence:
258    
259    C-       gradient in x direction:
260    #ifndef ALLOW_AUTODIFF_TAMC
261             IF (useCubedSphereExchange) THEN
262    C        to compute d/dx(hDiv), fill corners with appropriate values:
263               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
264             ENDIF
265    #endif
266             DO j=2-Oly,sNy+Oly-1
267              DO i=2-Olx,sNx+Olx-1
268                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
269              ENDDO
270             ENDDO
271    
272    C-       gradient in y direction:
273    #ifndef ALLOW_AUTODIFF_TAMC
274             IF (useCubedSphereExchange) THEN
275    C        to compute d/dy(hDiv), fill corners with appropriate values:
276               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
277             ENDIF
278    #endif
279             DO j=2-Oly,sNy+Oly-1
280              DO i=2-Olx,sNx+Olx-1
281                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
282              ENDDO
283             ENDDO
284    
285    C      horizontal gradient of vertical vorticity:
286    C-       gradient in x direction:
287             DO j=2-Oly,sNy+Oly
288              DO i=2-Olx,sNx+Olx-1
289                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
290         &                  *recip_DXG(i,j,bi,bj)
291         &                  *maskS(i,j,k,bi,bj)
292              ENDDO
293             ENDDO
294    C-       gradient in y direction:
295             DO j=2-Oly,sNy+Oly-1
296              DO i=2-Olx,sNx+Olx
297                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
298         &                  *recip_DYG(i,j,bi,bj)
299         &                  *maskW(i,j,k,bi,bj)
300              ENDDO
301             ENDDO
302    
303           ENDIF
304    
305         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
306          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
307  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
308    
309  C These are (powers of) length scales  C These are (powers of) length scales
310           L2=rA(i,j,bi,bj)           IF (useAreaViscLength) THEN
311              L2=rA(i,j,bi,bj)
312             ELSE
313              L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
314             ENDIF
315           L3=(L2**1.5)           L3=(L2**1.5)
316           L4=(L2**2)           L4=(L2**2)
317           L5=0.125*(L2**2.5)           L5=(L2**2.5)
318           IF (useAnisotropicViscAGridMax) THEN  
319            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
320       &                       +recip_DYF(I,J,bi,bj)**2) )  
321            L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4           IF (useAreaViscLength) THEN
322       &                       +recip_DYF(I,J,bi,bj)**4)            L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2
323       &                   +8.*((recip_DXF(I,J,bi,bj)           ELSE
324       &                        *recip_DYF(I,J,bi,bj))**2) )            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
325         &                            +recip_DYF(I,J,bi,bj)**4)
326         &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
327         &                             *recip_DYF(I,J,bi,bj))**2) )
328           ENDIF           ENDIF
329    
330           IF (useFullLeith) THEN  C Velocity Reynolds Scale
331             IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
332               Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
333             ELSE
334               Uscl=0.
335             ENDIF
336             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
337               U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
338             ELSE
339               U4scl=0.
340             ENDIF
341    
342    #ifndef ALLOW_AUTODIFF_TAMC
343             IF (useFullLeith.and.calcleith) THEN
344  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
345            grdVrt=0.25*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
346       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
347       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
348       &     +((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj))**2       &                        + vrtDy(i,j)*vrtDy(i,j) )  )
      &     +((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)  
349    
350  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
351  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
352             grdDiv=0.25*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
353       &      ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
354       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
355       &      +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
356       &      +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)  
357              viscAh_DLth(i,j)=
358             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
359       &          .NE. 0. ) THEN            viscA4_DLth(i,j)=
360              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
361             ELSE            viscAh_DLthd(i,j)=
362              Alth2=0. _d 0       &     sqrt(leithD2fac*grdDiv)*L3
363             ENDIF            viscA4_DLthd(i,j)=
364             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)       &     sqrt(leithD4fac*grdDiv)*L5
365       &          .NE. 0. ) THEN           ELSEIF (calcleith) THEN
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
366  C but this approximation will work on cube  C but this approximation will work on cube
367  c (and differs by as much as 4X)  c (and differs by as much as 4X)
368             grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=max( abs(vrtDx(i,j+1)), abs(vrtDx(i,j)) )
369             grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
370       &      abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
            grdVrt=max(grdVrt,  
      &      abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))  
            grdVrt=max(grdVrt,  
      &      abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,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)-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)))  
371    
372  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
373             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
374             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
375           ENDIF            grdDiv=max( grdDiv, abs(divDy(i,j))   )
376    
377           IF (smag2fac.NE.0.) THEN            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
378            Asmag2=smag2fac*L2            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
379       &           *sqrt(tension(i,j)**2            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
380       &           +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
      &                 +strain(i-1, j )**2+strain( i ,j-1)**2))  
381           ELSE           ELSE
382            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
383              viscA4_Dlth(i,j)=0. _d 0
384              viscAh_DlthD(i,j)=0. _d 0
385              viscA4_DlthD(i,j)=0. _d 0
386           ENDIF           ENDIF
387    
388           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
389            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
390       &             *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
391       &             +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
392       &                   +strain(i-1, j )**2+strain( i ,j-1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
393              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
394              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
395           ELSE           ELSE
396            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
397              viscA4_DSmg(i,j)=0. _d 0
398           ENDIF           ENDIF
399    #endif /* ALLOW_AUTODIFF_TAMC */
400    
401  C  Harmonic on Div.u points  C  Harmonic on Div.u points
402           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
403           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
404           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
405            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)
406            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
407           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))  
408    
409  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
410           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
411           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
412           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
413              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)
414              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
415           ELSE           viscA4_D(i,j)=min(viscA4_DMax(i,j),viscA4_D(i,j))
416           IF (vg4Max.GT.0.) THEN  
417              AlinMax=vg4Max*L4  #ifdef ALLOW_NONHYDROSTATIC
418              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))  C     /* Pass Viscosities to calc_gw, if constant, not necessary */
419           ENDIF  
420           ENDIF           kp1  = MIN(k+1,Nr)
421           AlinMin=vg4Min*L4  
422           viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))           if (k .eq. 1) then
423              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
424              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
425    
426              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j) /* These values dont get used */
427              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
428             else
429    C Note that previous call of this function has already added half.
430              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
431              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
432    
433              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
434              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
435             endif
436    #endif /* ALLOW_NONHYDROSTATIC */
437    
438  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
439  C These are (powers of) length scales  C These are (powers of) length scales
440           L2=rAz(i,j,bi,bj)           IF (useAreaViscLength) THEN
441              L2=rAz(i,j,bi,bj)
442             ELSE
443              L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
444             ENDIF
445    
446           L3=(L2**1.5)           L3=(L2**1.5)
447           L4=(L2**2)           L4=(L2**2)
448           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  
449    
450  C This is the vector magnitude of the vorticity gradient squared           L2rdt=0.25 _d 0*recip_dt*L2
451           IF (useFullLeith) THEN           IF (useAreaViscLength) THEN
452             grdVrt=0.25*(            L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
453       &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2           ELSE
454       &      +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2            L4rdt=recip_dt/
455       &      +((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)
456       &      +((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))
457             ENDIF
458    
459  C This is the vector magnitude of grad(div.v) squared  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
460             grdDiv=0.25*(           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
461       &       ((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))
462       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                      +(KE(i-1,j)+KE(i,j-1)) )
463       &      +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2             IF ( keZpt.GT.0. ) THEN
464       &      +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)               Uscl = sqrt(keZpt*L2)*viscAhRe_max
465                 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  
466             ELSE             ELSE
467              Alth2=0. _d 0               Uscl =0.
468             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  
469             ENDIF             ENDIF
470           ELSE           ELSE
471               Uscl =0.
472               U4scl=0.
473             ENDIF
474    
475    #ifndef ALLOW_AUTODIFF_TAMC
476    C This is the vector magnitude of the vorticity gradient squared
477             IF (useFullLeith.and.calcleith) THEN
478              grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
479         &                        + vrtDx(i,j)*vrtDx(i,j) )
480         &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
481         &                        + vrtDy(i,j)*vrtDy(i,j) )  )
482    
483    C This is the vector magnitude of grad(div.v) squared
484              grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
485         &                        + divDx(i,j)*divDx(i,j) )
486         &                     + (divDy(i-1,j)*divDy(i-1,j)
487         &                        + divDy(i,j)*divDy(i,j) )  )
488    
489              viscAh_ZLth(i,j)=
490         &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
491              viscA4_ZLth(i,j)=
492         &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
493              viscAh_ZLthD(i,j)=
494         &     sqrt(leithD2fac*grdDiv)*L3
495              viscA4_ZLthD(i,j)=
496         &     sqrt(leithD4fac*grdDiv)*L5
497    
498             ELSEIF (calcleith) THEN
499  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
500             grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=max( abs(vrtDx(i-1,j)), abs(vrtDx(i,j)) )
501             grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
502       &       abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
503             grdVrt=max(grdVrt,  
504       &       abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
505             grdVrt=max(grdVrt,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
506       &       abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
507    
508             grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
509             grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
510       &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
511             grdDiv=max(grdDiv,            viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
      &      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  
   
          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))  
512           ELSE           ELSE
513            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
514              viscA4_ZLth(i,j)=0. _d 0
515              viscAh_ZLthD(i,j)=0. _d 0
516              viscA4_ZLthD(i,j)=0. _d 0
517           ENDIF           ENDIF
518    
519           IF (smag4fac.NE.0.) THEN           IF (calcsmag) THEN
520            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
521       &             *sqrt(strain(i,j)**2       &      *sqrt(strain(i,j)**2
522       &             +0.25*(tension( i , j )**2+tension( i ,j+1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
523       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
524           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
525            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
526           ENDIF           ENDIF
527    #endif /* ALLOW_AUTODIFF_TAMC */
528    
529  C  Harmonic on Zeta points  C  Harmonic on Zeta points
530           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
531           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
532           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)
533              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)
534              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)
535           ELSE           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))
536           IF (vg2Max.GT.0.) THEN  
537              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
538              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
539           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
540           ENDIF           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)
541           AlinMin=vg2Min*L2           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)
542           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
543             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))  
544          ENDDO          ENDDO
545         ENDDO         ENDDO
546        ELSE        ELSE
# Line 337  C  BiHarmonic on Zeta Points Line 560  C  BiHarmonic on Zeta Points
560         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
561         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
562         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
563    #ifdef ALLOW_NONHYDROSTATIC
564           CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',k,1,2,bi,bj,myThid)
565           CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',k,1,2,bi,bj,myThid)
566    #endif
567    
568           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
569           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
570           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
571           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
572    
573           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
574           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
575           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
576           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
577    
578           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
579           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
580           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
581           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
582    
583           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
584         &   ,k,1,2,bi,bj,myThid)
585           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
586         &   ,k,1,2,bi,bj,myThid)
587           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
588         &   ,k,1,2,bi,bj,myThid)
589           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
590         &   ,k,1,2,bi,bj,myThid)
591    
592           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
593           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
594           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
595           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
596        ENDIF        ENDIF
597  #endif  #endif
598    
599        RETURN        RETURN
600        END        END
601    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22