/[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.37 - (hide annotations) (download)
Mon Apr 6 23:47:06 2009 UTC (15 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m
Changes since 1.36: +3 -2 lines
Overlaps had been forgotten in calculating ijk keys
(spotted by jmc using gfortran with check-bounds)

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

  ViewVC Help
Powered by ViewVC 1.1.22