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