/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_visc.F
ViewVC logotype

Annotation of /MITgcm/pkg/mom_common/mom_calc_visc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.46 - (hide annotations) (download)
Thu Aug 1 20:13:38 2013 UTC (10 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n
Changes since 1.45: +31 -33 lines
- always set horiz. viscosity arrays to background value before calling
  MOM_CALC_VISC (in MOM_FLUXFORM & MOM_VECINV) and call S/R MOM_CALC_VISC
  only when using variable horiz. viscosity (useVariableVisc=T);

1 jmc 1.46 C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_visc.F,v 1.45 2013/07/28 21:04:25 jmc Exp $
2 jmc 1.14 C $Name: $
3 baylor 1.1
4     #include "MOM_COMMON_OPTIONS.h"
5    
6 jmc 1.45 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7     CBOP
8     C !ROUTINE: MOM_CALC_VISC
9 baylor 1.5
10 jmc 1.45 C !INTERFACE:
11 baylor 1.1 SUBROUTINE MOM_CALC_VISC(
12     I bi,bj,k,
13     O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
14 jmc 1.12 I hDiv,vort3,tension,strain,KE,hFacZ,
15 baylor 1.1 I myThid)
16    
17 jmc 1.45 C !DESCRIPTION:
18 baylor 1.5 C Calculate horizontal viscosities (L is typical grid width)
19     C harmonic viscosity=
20     C viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
21     C +0.25*L**2*viscAhGrid/deltaT
22 baylor 1.17 C +sqrt((viscC2leith/pi)**6*grad(Vort3)**2
23     C +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3
24 baylor 1.5 C +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
25     C
26     C biharmonic viscosity=
27     C viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
28     C +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
29 baylor 1.17 C +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2
30     C +(viscC4leithD/pi)**6*grad(hDiv)**2)
31 baylor 1.5 C +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
32     C
33     C Note that often 0.125*L**2 is the scale between harmonic and
34     C biharmonic (see Griffies and Hallberg (2000))
35     C This allows the same value of the coefficient to be used
36     C for roughly similar results with biharmonic and harmonic
37     C
38     C LIMITERS -- limit min and max values of viscosities
39 jmc 1.32 C viscAhReMax is min value for grid point harmonic Reynolds num
40     C harmonic viscosity>sqrt(2*KE)*L/viscAhReMax
41 baylor 1.5 C
42 jmc 1.32 C viscA4ReMax is min value for grid point biharmonic Reynolds num
43     C biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4ReMax
44 baylor 1.5 C
45     C viscAhgridmax is CFL stability limiter for harmonic viscosity
46     C harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
47     C
48     C viscA4gridmax is CFL stability limiter for biharmonic viscosity
49     C biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)
50     C
51     C viscAhgridmin and viscA4gridmin are lower limits for viscosity:
52 cnh 1.25 C harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT
53     C biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)
54    
55 baylor 1.5 C RECOMMENDED VALUES
56 baylor 1.18 C viscC2Leith=1-3
57     C viscC2LeithD=1-3
58     C viscC4Leith=1-3
59     C viscC4LeithD=1.5-3
60 jmc 1.32 C viscC2smag=2.2-4 (Griffies and Hallberg,2000)
61 baylor 1.5 C 0.2-0.9 (Smagorinsky,1993)
62 jmc 1.32 C viscC4smag=2.2-4 (Griffies and Hallberg,2000)
63     C viscAhReMax>=1, (<2 suppresses a computational mode)
64     C viscA4ReMax>=1, (<2 suppresses a computational mode)
65 baylor 1.5 C viscAhgridmax=1
66     C viscA4gridmax=1
67     C viscAhgrid<1
68     C viscA4grid<1
69     C viscAhgridmin<<1
70     C viscA4gridmin<<1
71 baylor 1.1
72 jmc 1.45 C !USES:
73     IMPLICIT NONE
74    
75 baylor 1.1 C == Global variables ==
76     #include "SIZE.h"
77     #include "GRID.h"
78     #include "EEPARAMS.h"
79     #include "PARAMS.h"
80 jmc 1.45 #include "MOM_VISC.h"
81 heimbach 1.33 #ifdef ALLOW_AUTODIFF_TAMC
82     #include "tamc.h"
83     #include "tamc_keys.h"
84     #endif /* ALLOW_AUTODIFF_TAMC */
85 baylor 1.1
86 jmc 1.45 C !INPUT/OUTPUT PARAMETERS:
87     C myThid :: my thread Id number
88 baylor 1.1 INTEGER bi,bj,k
89     _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90     _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91     _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92     _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93     _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94     _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95     _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96     _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97     _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     INTEGER myThid
100 jmc 1.45 CEOP
101 baylor 1.1
102 jmc 1.45 C !LOCAL VARIABLES:
103 jmc 1.44 INTEGER i,j
104 jmc 1.29 #ifdef ALLOW_NONHYDROSTATIC
105 jmc 1.45 _RL shiftAh, shiftA4
106 jmc 1.29 #endif
107 jmc 1.42 #ifdef ALLOW_AUTODIFF_TAMC
108 heimbach 1.34 INTEGER lockey_1, lockey_2
109 jmc 1.42 #endif
110 baylor 1.5 _RL smag2fac, smag4fac
111 baylor 1.17 _RL leith2fac, leith4fac
112     _RL leithD2fac, leithD4fac
113 baylor 1.6 _RL viscAhRe_max, viscA4Re_max
114 jmc 1.15 _RL Alin,grdVrt,grdDiv, keZpt
115 jmc 1.45 _RL L2, L3, L5, L2rdt, L4rdt, recip_dt
116 baylor 1.5 _RL Uscl,U4scl
117 jmc 1.16 _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118     _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119 jmc 1.20 _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120     _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 baylor 1.5 _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122     _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123     _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124     _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125     _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126     _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127     _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128     _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129     _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130     _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131     _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132     _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133     _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134     _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135     _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136     _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137     _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138     _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139     _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140     _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141 jmc 1.32 LOGICAL calcLeith, calcSmag
142 baylor 1.1
143 heimbach 1.33 #ifdef ALLOW_AUTODIFF_TAMC
144     act1 = bi - myBxLo(myThid)
145     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
146     act2 = bj - myByLo(myThid)
147     max2 = myByHi(myThid) - myByLo(myThid) + 1
148     act3 = myThid - 1
149     max3 = nTx*nTy
150     act4 = ikey_dynamics - 1
151     ikey = (act1 + 1) + act2*max1
152     & + act3*max1*max2
153     & + act4*max1*max2*max3
154 heimbach 1.34 lockey_1 = (ikey-1)*Nr + k
155 heimbach 1.33 #endif /* ALLOW_AUTODIFF_TAMC */
156    
157 jmc 1.32 C-- Set flags which are used in this S/R and elsewhere :
158 jmc 1.45 C useVariableVisc, useHarmonicVisc and useBiharmonicVisc
159     C are now set early on (in S/R SET_PARAMS)
160 baylor 1.1
161 jmc 1.46 c IF ( useVariableVisc ) THEN
162 jmc 1.45 C---- variable viscosity :
163 baylor 1.1
164 jmc 1.45 recip_dt = 1. _d 0
165     IF ( deltaTmom.NE.0. ) recip_dt = 1. _d 0/deltaTmom
166 jmc 1.32
167 jmc 1.45 IF ( useHarmonicVisc .AND. viscAhReMax.NE.0. ) THEN
168 jmc 1.32 viscAhRe_max=SQRT(2. _d 0)/viscAhReMax
169     ELSE
170     viscAhRe_max=0. _d 0
171     ENDIF
172    
173 jmc 1.45 IF ( useBiharmonicVisc .AND. viscA4ReMax.NE.0. ) THEN
174 jmc 1.32 viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax
175     ELSE
176     viscA4Re_max=0. _d 0
177     ENDIF
178 baylor 1.5
179 jmc 1.32 calcLeith=
180 baylor 1.5 & (viscC2leith.NE.0.)
181     & .OR.(viscC2leithD.NE.0.)
182     & .OR.(viscC4leith.NE.0.)
183     & .OR.(viscC4leithD.NE.0.)
184    
185 jmc 1.32 calcSmag=
186 baylor 1.5 & (viscC2smag.NE.0.)
187     & .OR.(viscC4smag.NE.0.)
188    
189 jmc 1.32 IF (calcSmag) THEN
190 baylor 1.5 smag2fac=(viscC2smag/pi)**2
191 jmc 1.10 smag4fac=0.125 _d 0*(viscC4smag/pi)**2
192 jmc 1.32 ELSE
193 jmc 1.10 smag2fac=0. _d 0
194     smag4fac=0. _d 0
195 jmc 1.32 ENDIF
196 baylor 1.1
197 jmc 1.32 IF (calcLeith) THEN
198 baylor 1.17 IF (useFullLeith) THEN
199 baylor 1.19 leith2fac =(viscC2leith /pi)**6
200 baylor 1.17 leithD2fac=(viscC2leithD/pi)**6
201 baylor 1.19 leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
202 baylor 1.17 leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
203     ELSE
204 baylor 1.19 leith2fac =(viscC2leith /pi)**3
205 baylor 1.17 leithD2fac=(viscC2leithD/pi)**3
206 baylor 1.19 leith4fac =0.125 _d 0*(viscC4leith /pi)**3
207     leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
208 baylor 1.17 ENDIF
209 jmc 1.32 ELSE
210 baylor 1.17 leith2fac=0. _d 0
211     leith4fac=0. _d 0
212     leithD2fac=0. _d 0
213     leithD4fac=0. _d 0
214 jmc 1.32 ENDIF
215 baylor 1.17
216 jmc 1.45 DO j=1-OLy,sNy+OLy
217     DO i=1-OLx,sNx+OLx
218 jmc 1.46 C- viscosity arrays have been initialised everywhere before calling this S/R
219     c viscAh_D(i,j) = viscAhD
220     c viscAh_Z(i,j) = viscAhZ
221     c viscA4_D(i,j) = viscA4D
222     c viscA4_Z(i,j) = viscA4Z
223 jmc 1.45
224 heimbach 1.21 visca4_zsmg(i,j) = 0. _d 0
225     viscah_zsmg(i,j) = 0. _d 0
226 jmc 1.45
227 heimbach 1.21 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 jmc 1.45
232 heimbach 1.21 viscAh_DSmg(i,j) = 0. _d 0
233     viscA4_DSmg(i,j) = 0. _d 0
234 jmc 1.45
235 heimbach 1.21 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 jmc 1.32 ENDDO
241 jmc 1.16
242 jmc 1.46 C- Initialise to zero gradient of vorticity & divergence:
243 jmc 1.45 DO j=1-OLy,sNy+OLy
244     DO i=1-OLx,sNx+OLx
245 jmc 1.16 divDx(i,j) = 0.
246     divDy(i,j) = 0.
247 jmc 1.20 vrtDx(i,j) = 0.
248     vrtDy(i,j) = 0.
249 jmc 1.16 ENDDO
250     ENDDO
251 jmc 1.20
252 jmc 1.46 IF ( calcLeith ) THEN
253     C-- horizontal gradient of horizontal divergence:
254 jmc 1.16 C- gradient in x direction:
255     IF (useCubedSphereExchange) THEN
256     C to compute d/dx(hDiv), fill corners with appropriate values:
257 jmc 1.36 CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
258 jmc 1.27 & hDiv, bi,bj, myThid )
259 jmc 1.16 ENDIF
260 jmc 1.45 DO j=2-OLy,sNy+OLy-1
261     DO i=2-OLx,sNx+OLx-1
262 jmc 1.44 divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_dxC(i,j,bi,bj)
263 jmc 1.16 ENDDO
264     ENDDO
265    
266     C- gradient in y direction:
267     IF (useCubedSphereExchange) THEN
268     C to compute d/dy(hDiv), fill corners with appropriate values:
269 jmc 1.36 CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
270 jmc 1.27 & hDiv, bi,bj, myThid )
271 jmc 1.16 ENDIF
272 jmc 1.45 DO j=2-OLy,sNy+OLy-1
273     DO i=2-OLx,sNx+OLx-1
274 jmc 1.44 divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_dyC(i,j,bi,bj)
275 jmc 1.16 ENDDO
276     ENDDO
277 jmc 1.20
278 jmc 1.46 C-- horizontal gradient of vertical vorticity:
279 jmc 1.20 C- gradient in x direction:
280 jmc 1.45 DO j=2-OLy,sNy+OLy
281     DO i=2-OLx,sNx+OLx-1
282 jmc 1.20 vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
283 jmc 1.44 & *recip_dxG(i,j,bi,bj)
284 jmc 1.20 & *maskS(i,j,k,bi,bj)
285 jmc 1.44 #ifdef ALLOW_OBCS
286     & *maskInS(i,j,bi,bj)
287     #endif
288 jmc 1.20 ENDDO
289     ENDDO
290     C- gradient in y direction:
291 jmc 1.45 DO j=2-OLy,sNy+OLy-1
292     DO i=2-OLx,sNx+OLx
293 jmc 1.20 vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
294 jmc 1.44 & *recip_dyG(i,j,bi,bj)
295 jmc 1.20 & *maskW(i,j,k,bi,bj)
296 jmc 1.44 #ifdef ALLOW_OBCS
297     & *maskInW(i,j,bi,bj)
298     #endif
299 jmc 1.20 ENDDO
300     ENDDO
301    
302 jmc 1.46 C-- end if calcLeith
303 jmc 1.16 ENDIF
304    
305 jmc 1.45 DO j=2-OLy,sNy+OLy-1
306     DO i=2-OLx,sNx+OLx-1
307 baylor 1.1 CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
308 baylor 1.5
309 heimbach 1.34 #ifdef ALLOW_AUTODIFF_TAMC
310     # ifndef AUTODIFF_DISABLE_LEITH
311 dfer 1.38 lockey_2 = i+olx + (sNx+2*olx)*(j+oly-1)
312 heimbach 1.37 & + (sNx+2*olx)*(sNy+2*oly)*(lockey_1-1)
313 jmc 1.44 CADJ STORE viscA4_ZSmg(i,j)
314 heimbach 1.34 CADJ & = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
315 jmc 1.44 CADJ STORE viscAh_ZSmg(i,j)
316 heimbach 1.34 CADJ & = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
317     # endif
318     #endif /* ALLOW_AUTODIFF_TAMC */
319    
320 jmc 1.32 C These are (powers of) length scales
321 dfer 1.40 L2 = L2_D(i,j,bi,bj)
322     L2rdt = 0.25 _d 0*recip_dt*L2
323     L3 = L3_D(i,j,bi,bj)
324     L4rdt = L4rdt_D(i,j,bi,bj)
325     L5 = (L2*L3)
326 baylor 1.5
327 mlosch 1.41 #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
328 baylor 1.5 C Velocity Reynolds Scale
329 jmc 1.15 IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
330 jmc 1.32 Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
331 jmc 1.15 ELSE
332     Uscl=0.
333     ENDIF
334     IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
335 jmc 1.32 U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
336 jmc 1.15 ELSE
337     U4scl=0.
338     ENDIF
339 mlosch 1.41 #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
340 baylor 1.5
341 heimbach 1.33 #ifndef AUTODIFF_DISABLE_LEITH
342 jmc 1.32 IF (useFullLeith.AND.calcLeith) THEN
343 baylor 1.1 C This is the vector magnitude of the vorticity gradient squared
344 jmc 1.20 grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
345     & + vrtDx(i,j)*vrtDx(i,j) )
346     & + (vrtDy(i+1,j)*vrtDy(i+1,j)
347     & + vrtDy(i,j)*vrtDy(i,j) ) )
348 baylor 1.1
349     C This is the vector magnitude of grad (div.v) squared
350     C Using it in Leith serves to damp instabilities in w.
351 jmc 1.16 grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
352     & + divDx(i,j)*divDx(i,j) )
353     & + (divDy(i,j+1)*divDy(i,j+1)
354     & + divDy(i,j)*divDy(i,j) ) )
355 baylor 1.5
356     viscAh_DLth(i,j)=
357 jmc 1.32 & SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
358 baylor 1.17 viscA4_DLth(i,j)=
359 jmc 1.32 & SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
360 baylor 1.5 viscAh_DLthd(i,j)=
361 jmc 1.32 & SQRT(leithD2fac*grdDiv)*L3
362 baylor 1.17 viscA4_DLthd(i,j)=
363 jmc 1.32 & SQRT(leithD4fac*grdDiv)*L5
364     ELSEIF (calcLeith) THEN
365 jmc 1.45 C but this approximation will work on cube (and differs by as much as 4X)
366 jmc 1.32 grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )
367     grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
368     grdVrt=MAX( grdVrt, ABS(vrtDy(i,j)) )
369 baylor 1.5
370 jmc 1.45 C This approximation is good to the same order as above...
371 jmc 1.32 grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
372     grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
373     grdDiv=MAX( grdDiv, ABS(divDy(i,j)) )
374 baylor 1.1
375 baylor 1.17 viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
376     viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
377     viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
378     viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
379 baylor 1.1 ELSE
380 jmc 1.10 viscAh_Dlth(i,j)=0. _d 0
381     viscA4_Dlth(i,j)=0. _d 0
382     viscAh_DlthD(i,j)=0. _d 0
383     viscA4_DlthD(i,j)=0. _d 0
384 baylor 1.1 ENDIF
385    
386 jmc 1.32 IF (calcSmag) THEN
387 baylor 1.5 viscAh_DSmg(i,j)=L2
388 jmc 1.32 & *SQRT(tension(i,j)**2
389 jmc 1.10 & +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2
390     & +strain(i , j )**2+strain(i+1,j+1)**2))
391 baylor 1.5 viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
392     viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
393 baylor 1.1 ELSE
394 jmc 1.10 viscAh_DSmg(i,j)=0. _d 0
395     viscA4_DSmg(i,j)=0. _d 0
396 baylor 1.1 ENDIF
397 heimbach 1.33 #endif /* AUTODIFF_DISABLE_LEITH */
398 baylor 1.1
399     C Harmonic on Div.u points
400 baylor 1.5 Alin=viscAhD+viscAhGrid*L2rdt
401     & +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
402 gforget 1.43 #ifdef ALLOW_3D_VISCAH
403     & +viscAhDfld(i,j,k,bi,bj)
404     #endif
405 jmc 1.32 viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
406     viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
407     viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
408     viscAh_D(i,j)=MIN(viscAh_DMax(i,j),viscAh_D(i,j))
409 baylor 1.1
410     C BiHarmonic on Div.u points
411 baylor 1.5 Alin=viscA4D+viscA4Grid*L4rdt
412     & +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
413 gforget 1.43 #ifdef ALLOW_3D_VISCA4
414     & +viscA4Dfld(i,j,k,bi,bj)
415 jmc 1.44 #endif
416 jmc 1.32 viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
417     viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
418     viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
419     viscA4_D(i,j)=MIN(viscA4_DMax(i,j),viscA4_D(i,j))
420 baylor 1.1
421     CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
422 jmc 1.32 C These are (powers of) length scales
423 dfer 1.40 L2 = L2_Z(i,j,bi,bj)
424     L2rdt = 0.25 _d 0*recip_dt*L2
425     L3 = L3_Z(i,j,bi,bj)
426     L4rdt = L4rdt_Z(i,j,bi,bj)
427     L5 = (L2*L3)
428 baylor 1.5
429 mlosch 1.41 #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
430 jmc 1.15 C Velocity Reynolds Scale (Pb here at CS-grid corners !)
431     IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
432     keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
433     & +(KE(i-1,j)+KE(i,j-1)) )
434     IF ( keZpt.GT.0. ) THEN
435 jmc 1.32 Uscl = SQRT(keZpt*L2)*viscAhRe_max
436     U4scl= SQRT(keZpt)*L3*viscA4Re_max
437 jmc 1.15 ELSE
438     Uscl =0.
439     U4scl=0.
440     ENDIF
441     ELSE
442     Uscl =0.
443     U4scl=0.
444     ENDIF
445 mlosch 1.41 #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
446 baylor 1.1
447 heimbach 1.33 #ifndef AUTODIFF_DISABLE_LEITH
448 baylor 1.1 C This is the vector magnitude of the vorticity gradient squared
449 jmc 1.32 IF (useFullLeith.AND.calcLeith) THEN
450 jmc 1.20 grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
451     & + vrtDx(i,j)*vrtDx(i,j) )
452     & + (vrtDy(i,j-1)*vrtDy(i,j-1)
453     & + vrtDy(i,j)*vrtDy(i,j) ) )
454 baylor 1.1
455     C This is the vector magnitude of grad(div.v) squared
456 jmc 1.16 grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
457     & + divDx(i,j)*divDx(i,j) )
458     & + (divDy(i-1,j)*divDy(i-1,j)
459     & + divDy(i,j)*divDy(i,j) ) )
460 baylor 1.5
461     viscAh_ZLth(i,j)=
462 jmc 1.32 & SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
463 baylor 1.17 viscA4_ZLth(i,j)=
464 jmc 1.32 & SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
465 baylor 1.5 viscAh_ZLthD(i,j)=
466 jmc 1.32 & SQRT(leithD2fac*grdDiv)*L3
467 baylor 1.17 viscA4_ZLthD(i,j)=
468 jmc 1.32 & SQRT(leithD4fac*grdDiv)*L5
469 baylor 1.5
470 jmc 1.32 ELSEIF (calcLeith) THEN
471 baylor 1.1 C but this approximation will work on cube (and differs by 4X)
472 jmc 1.32 grdVrt=MAX( ABS(vrtDx(i-1,j)), ABS(vrtDx(i,j)) )
473     grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
474     grdVrt=MAX( grdVrt, ABS(vrtDy(i,j)) )
475    
476     grdDiv=MAX( ABS(divDx(i,j)), ABS(divDx(i,j-1)) )
477     grdDiv=MAX( grdDiv, ABS(divDy(i,j)) )
478     grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) )
479 baylor 1.5
480 baylor 1.17 viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
481     viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
482     viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
483     viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
484 baylor 1.1 ELSE
485 jmc 1.10 viscAh_ZLth(i,j)=0. _d 0
486     viscA4_ZLth(i,j)=0. _d 0
487     viscAh_ZLthD(i,j)=0. _d 0
488     viscA4_ZLthD(i,j)=0. _d 0
489 baylor 1.1 ENDIF
490    
491 jmc 1.32 IF (calcSmag) THEN
492 baylor 1.5 viscAh_ZSmg(i,j)=L2
493 jmc 1.32 & *SQRT(strain(i,j)**2
494 jmc 1.10 & +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
495     & +tension(i-1, j )**2+tension(i-1,j-1)**2))
496 baylor 1.5 viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
497     viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
498 baylor 1.1 ENDIF
499 heimbach 1.33 #endif /* AUTODIFF_DISABLE_LEITH */
500 baylor 1.1
501     C Harmonic on Zeta points
502 baylor 1.5 Alin=viscAhZ+viscAhGrid*L2rdt
503     & +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
504 jmc 1.44 #ifdef ALLOW_3D_VISCAH
505 gforget 1.43 & +viscAhZfld(i,j,k,bi,bj)
506 jmc 1.44 #endif
507 jmc 1.32 viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
508     viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
509     viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
510     viscAh_Z(i,j)=MIN(viscAh_ZMax(i,j),viscAh_Z(i,j))
511 baylor 1.5
512     C BiHarmonic on Zeta points
513     Alin=viscA4Z+viscA4Grid*L4rdt
514     & +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
515 gforget 1.43 #ifdef ALLOW_3D_VISCA4
516     & +viscA4Zfld(i,j,k,bi,bj)
517 jmc 1.44 #endif
518 jmc 1.32 viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
519     viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
520     viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
521     viscA4_Z(i,j)=MIN(viscA4_ZMax(i,j),viscA4_Z(i,j))
522 baylor 1.1 ENDDO
523     ENDDO
524 jmc 1.32
525 jmc 1.45 #ifdef ALLOW_NONHYDROSTATIC
526     IF ( nonHydrostatic ) THEN
527 jmc 1.46 C-- Pass Viscosities to calc_gw (if constant, not necessary)
528 jmc 1.45
529     IF ( k.LT.Nr ) THEN
530     C Prepare for next level (next call)
531     DO j=1-OLy,sNy+OLy
532     DO i=1-OLx,sNx+OLx
533     viscAh_W(i,j,k+1,bi,bj) = halfRL*viscAh_D(i,j)
534     viscA4_W(i,j,k+1,bi,bj) = halfRL*viscA4_D(i,j)
535     ENDDO
536     ENDDO
537     ENDIF
538    
539     shiftAh = viscAhW - viscAhD
540     shiftA4 = viscA4W - viscA4D
541     IF ( k.EQ.1 ) THEN
542     C These values dont get used
543     DO j=1-OLy,sNy+OLy
544     DO i=1-OLx,sNx+OLx
545     viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_D(i,j)
546     viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_D(i,j)
547     ENDDO
548     ENDDO
549     ELSE
550     C Note that previous call of this function has already added half.
551     DO j=1-OLy,sNy+OLy
552     DO i=1-OLx,sNx+OLx
553     viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_W(i,j,k,bi,bj)
554     & + halfRL*viscAh_D(i,j)
555     viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_W(i,j,k,bi,bj)
556     & + halfRL*viscA4_D(i,j)
557     ENDDO
558     ENDDO
559     ENDIF
560    
561     ENDIF
562     #endif /* ALLOW_NONHYDROSTATIC */
563    
564 jmc 1.46 c ELSE
565 jmc 1.45 C---- use constant viscosity (useVariableVisc=F):
566 jmc 1.46 c DO j=1-OLy,sNy+OLy
567     c DO i=1-OLx,sNx+OLx
568     c viscAh_D(i,j) = viscAhD
569     c viscAh_Z(i,j) = viscAhZ
570     c viscA4_D(i,j) = viscA4D
571     c viscA4_Z(i,j) = viscA4Z
572     c ENDDO
573     c ENDDO
574 jmc 1.32 C---- variable/constant viscosity : end if/else block
575 jmc 1.46 c ENDIF
576 baylor 1.1
577     #ifdef ALLOW_DIAGNOSTICS
578     IF (useDiagnostics) THEN
579     CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid)
580     CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
581     CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
582     CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
583 baylor 1.5
584     CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
585     CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
586     CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
587     CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
588    
589     CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
590     CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
591     CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
592     CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
593    
594     CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
595     CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
596     CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
597     CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
598    
599 jmc 1.46 CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD',
600     & k,1,2,bi,bj,myThid)
601     CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD',
602     & k,1,2,bi,bj,myThid)
603     CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD',
604     & k,1,2,bi,bj,myThid)
605     CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD',
606     & k,1,2,bi,bj,myThid)
607 baylor 1.5
608     CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
609     CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
610     CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
611     CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
612 baylor 1.1 ENDIF
613     #endif
614    
615     RETURN
616     END

  ViewVC Help
Powered by ViewVC 1.1.22