/[MITgcm]/MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F
ViewVC logotype

Contents of /MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F

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


Revision 1.21 - (show annotations) (download)
Thu Mar 24 14:32:46 2005 UTC (19 years, 2 months ago) by baylor
Branch: MAIN
Changes since 1.20: +33 -5 lines
Add approximate grdDiv formula (useFullLeith=.FALSE. for cubed sphere) capability to modified Leith scheme.

1 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F,v 1.20 2005/03/23 18:21:06 adcroft Exp $
2 C $Name: $
3
4 #include "MOM_VECINV_OPTIONS.h"
5
6 SUBROUTINE MOM_VI_HDISSIP(
7 I bi,bj,k,
8 I hDiv,vort3,hFacZ,dStar,zStar,
9 O uDissip,vDissip,
10 I myThid)
11
12 cph(
13 cph The following line was commented in the argument list
14 cph TAMC cannot digest commented lines within continuing lines
15 c I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
16 cph)
17
18 IMPLICIT NONE
19 C
20 C Calculate horizontal dissipation terms
21 C [del^2 - del^4] (u,v)
22 C
23
24 C == Global variables ==
25 #include "SIZE.h"
26 #include "GRID.h"
27 #include "EEPARAMS.h"
28 #include "PARAMS.h"
29
30 C == Routine arguments ==
31 INTEGER bi,bj,k
32 _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
35 _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36 _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
37 _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38 _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39 INTEGER myThid
40
41 C == Local variables ==
42 _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43 _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44 _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45 _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46 INTEGER I,J
47 _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4
48 _RL Alin,AlinMin,AlinMax,Alth2,Alth4,grdVrt,grdDiv
49 _RL vg2,vg2Min,vg2Max,vg4,vg4Min,vg4Max
50 _RL recip_dt,L2,L3,L4,L5
51 LOGICAL useVariableViscosity
52
53 useVariableViscosity=
54 & (viscAhGrid.NE.0.)
55 & .OR.(viscA4Grid.NE.0.)
56 & .OR.(viscC2leith.NE.0.)
57 & .OR.(viscC2leithD.NE.0.)
58 & .OR.(viscC4leith.NE.0.)
59 & .OR.(viscC4leithD.NE.0.)
60
61 IF (deltaTmom.NE.0.) THEN
62 recip_dt=1./deltaTmom
63 ELSE
64 recip_dt=0.
65 ENDIF
66
67 vg2=viscAhGrid*recip_dt
68 vg2Min=viscAhGridMin*recip_dt
69 vg2Max=viscAhGridMax*recip_dt
70 vg4=viscA4Grid*recip_dt
71 vg4Min=viscA4GridMin*recip_dt
72 vg4Max=viscA4GridMax*recip_dt
73
74 C - Viscosity
75 IF (useVariableViscosity) THEN
76 DO j=2-Oly,sNy+Oly-1
77 DO i=2-Olx,sNx+Olx-1
78 C These are (powers of) length scales used in the Leith viscosity calculation
79 L2=rA(i,j,bi,bj)
80 L3=(L2**1.5)
81 L4=(L2**2)
82 L5=0.125*(L2**2.5)
83
84
85 IF (useFullLeith) THEN
86 C This is the vector magnitude of the vorticity gradient squared
87 grdVrt=0.25*(
88 & ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
89 & +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
90 & +((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj))**2
91 & +((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)
92
93 C This is the vector magnitude of grad (div.v) squared
94 C Using it in Leith serves to damp instabilities in w.
95 grdDiv=0.25*(
96 & ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2
97 & +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2
98 & +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2
99 & +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)
100
101 Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
102 Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
103 ELSE
104 C but this approximation will work on cube
105 c (and differs by as much as 4X)
106 grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
107 grdVrt=max(grdVrt,
108 & abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
109 grdVrt=max(grdVrt,
110 & abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))
111 grdVrt=max(grdVrt,
112 & abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))
113
114 grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))
115 grdDiv=max(grdDiv,
116 & abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))
117 grdDiv=max(grdDiv,
118 & abs((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj)))
119 grdDiv=max(grdDiv,
120 & abs((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj)))
121
122 C This if statement is just to prevent bitwise changes when leithd=0
123 IF ((viscC2leithd.eq.0).and.(viscc4leithd.eq.0)) THEN
124 Alth2=viscC2leith*grdVrt*L3
125 Alth4=viscC4leith*grdVrt*L5
126 ELSE
127 Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
128 Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
129 ENDIF
130 ENDIF
131
132 C Harmonic Modified Leith on Div.u points
133 Alin=viscAhD+vg2*L2+Alth2
134 viscAh_D(i,j)=min(viscAhMax,Alin)
135 IF (vg2Max.GT.0.) THEN
136 AlinMax=vg2Max*L2
137 viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))
138 ENDIF
139 AlinMin=vg2Min*L2
140 viscAh_D(i,j)=max(AlinMin,viscAh_D(i,j))
141
142 C BiHarmonic Modified Leith on Div.u points
143 Alin=viscA4D+vg4*L4+Alth4
144 viscA4_D(i,j)=min(viscA4Max,Alin)
145 IF (vg4Max.GT.0.) THEN
146 AlinMax=vg4Max*L4
147 viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))
148 ENDIF
149 AlinMin=vg4Min*L4
150 viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))
151
152 C These are (powers of) length scales used in the Leith viscosity calculation
153 L2=rAz(i,j,bi,bj)
154 L3=(L2**1.5)
155 L4=(L2**2)
156 L5=0.125*(L2**2.5)
157
158 C This is the vector magnitude of the vorticity gradient squared
159 IF (useFullLeith) THEN
160 grdVrt=0.25*(
161 & ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
162 & +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
163 & +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2
164 & +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
165
166 C This is the vector magnitude of grad(div.v) squared
167 grdDiv=0.25*(
168 & ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2
169 & +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2
170 & +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2
171 & +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)
172
173 Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
174 Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
175 ELSE
176 C but this approximation will work on cube (and differs by as much as 4X)
177 grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
178 grdVrt=max(grdVrt,
179 & abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
180 grdVrt=max(grdVrt,
181 & abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
182 grdVrt=max(grdVrt,
183 & abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
184
185 grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))
186 grdDiv=max(grdDiv,
187 & abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))
188 grdDiv=max(grdDiv,
189 & abs((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i-1,j,bi,bj)))
190 grdDiv=max(grdDiv,
191 & abs((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i,j-1,bi,bj)))
192
193 C This if statement is just to prevent bitwise changes when leithd=0
194 IF ((viscC2leithd.eq.0).and.(viscc4leithd.eq.0)) THEN
195 Alth2=viscC2leith*grdVrt*L3
196 Alth4=viscC4leith*grdVrt*L5
197 ELSE
198 Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
199 Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5
200 ENDIF
201 ENDIF
202
203 C Harmonic Modified Leith on Zeta points
204 Alin=viscAhZ+vg2*L2+Alth2
205 viscAh_Z(i,j)=min(viscAhMax,Alin)
206 IF (vg2Max.GT.0.) THEN
207 AlinMax=vg2Max*L2
208 viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))
209 ENDIF
210 AlinMin=vg2Min*L2
211 viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))
212
213 C BiHarmonic Modified Leith on Zeta Points
214 Alin=viscA4Z+vg4*L4+Alth4
215 viscA4_Z(i,j)=min(viscA4Max,Alin)
216 IF (vg4Max.GT.0.) THEN
217 AlinMax=vg4Max*L4
218 viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))
219 ENDIF
220 AlinMin=vg4Min*L4
221 viscA4_Z(i,j)=max(AlinMin,viscA4_Z(i,j))
222 ENDDO
223 ENDDO
224 ELSE
225 DO j=1-Oly,sNy+Oly
226 DO i=1-Olx,sNx+Olx
227 viscAh_D(i,j)=viscAhD
228 viscAh_Z(i,j)=viscAhZ
229 viscA4_D(i,j)=viscA4D
230 viscA4_Z(i,j)=viscA4Z
231 ENDDO
232 ENDDO
233 ENDIF
234
235 C - Laplacian terms
236 IF ( viscC2leith.NE.0. .OR. viscAhGrid.NE.0.
237 & .OR. viscAhD.NE.0. .OR. viscAhZ.NE.0. ) THEN
238 DO j=2-Oly,sNy+Oly-1
239 DO i=2-Olx,sNx+Olx-1
240
241 Dim=hDiv( i ,j-1)
242 Dij=hDiv( i , j )
243 Dmj=hDiv(i-1, j )
244 Zip=hFacZ( i ,j+1)*vort3( i ,j+1)
245 Zij=hFacZ( i , j )*vort3( i , j )
246 Zpj=hFacZ(i+1, j )*vort3(i+1, j )
247
248 C This bit scales the harmonic dissipation operator to be proportional
249 C to the grid-cell area over the time-step. viscAh is then non-dimensional
250 C and should be less than 1/8, for example viscAh=0.01
251 IF (useVariableViscosity) THEN
252 Dij=Dij*viscAh_D(i,j)
253 Dim=Dim*viscAh_D(i,j-1)
254 Dmj=Dmj*viscAh_D(i-1,j)
255 Zij=Zij*viscAh_Z(i,j)
256 Zip=Zip*viscAh_Z(i,j+1)
257 Zpj=Zpj*viscAh_Z(i+1,j)
258 uD2 = (
259 & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
260 & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
261 vD2 = (
262 & recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
263 & *cosFacV(j,bi,bj)
264 & +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
265 ELSE
266 uD2 = viscAhD*
267 & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
268 & - viscAhZ*recip_hFacW(i,j,k,bi,bj)*
269 & ( Zip-Zij )*recip_DYG(i,j,bi,bj)
270 vD2 = viscAhZ*recip_hFacS(i,j,k,bi,bj)*
271 & cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
272 & + viscAhD* ( Dij-Dim )*recip_DYC(i,j,bi,bj)
273 ENDIF
274
275 uDissip(i,j) = uD2
276 vDissip(i,j) = vD2
277
278 ENDDO
279 ENDDO
280 ELSE
281 DO j=2-Oly,sNy+Oly-1
282 DO i=2-Olx,sNx+Olx-1
283 uDissip(i,j) = 0.
284 vDissip(i,j) = 0.
285 ENDDO
286 ENDDO
287 ENDIF
288
289 C - Bi-harmonic terms
290 IF ( viscC4leith.NE.0. .OR. viscA4Grid.NE.0.
291 & .OR. viscA4D.NE.0. .OR. viscA4Z.NE.0. ) THEN
292 DO j=2-Oly,sNy+Oly-1
293 DO i=2-Olx,sNx+Olx-1
294
295 #ifdef MOM_VI_ORIGINAL_VISCA4
296 Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)
297 Dij=dyF( i , j ,bi,bj)*dStar( i , j )
298 Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )
299
300 Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)
301 Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )
302 Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )
303 #else
304 Dim=dStar( i ,j-1)
305 Dij=dStar( i , j )
306 Dmj=dStar(i-1, j )
307
308 Zip=hFacZ( i ,j+1)*zStar( i ,j+1)
309 Zij=hFacZ( i , j )*zStar( i , j )
310 Zpj=hFacZ(i+1, j )*zStar(i+1, j )
311 #endif
312
313 C This bit scales the harmonic dissipation operator to be proportional
314 C to the grid-cell area over the time-step. viscAh is then non-dimensional
315 C and should be less than 1/8, for example viscAh=0.01
316 IF (useVariableViscosity) THEN
317 Dij=Dij*viscA4_D(i,j)
318 Dim=Dim*viscA4_D(i,j-1)
319 Dmj=Dmj*viscA4_D(i-1,j)
320 Zij=Zij*viscA4_Z(i,j)
321 Zip=Zip*viscA4_Z(i,j+1)
322 Zpj=Zpj*viscA4_Z(i+1,j)
323
324 #ifdef MOM_VI_ORIGINAL_VISCA4
325 uD4 = recip_rAw(i,j,bi,bj)*(
326 & ( (Dij-Dmj)*cosFacU(j,bi,bj) )
327 & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij ) )
328 vD4 = recip_rAs(i,j,bi,bj)*(
329 & recip_hFacS(i,j,k,bi,bj)*( (Zpj-Zij)*cosFacV(j,bi,bj) )
330 & + ( Dij-Dim ) )
331 ELSE
332 uD4 = recip_rAw(i,j,bi,bj)*(
333 & viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )
334 & -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )
335 vD4 = recip_rAs(i,j,bi,bj)*(
336 & recip_hFacS(i,j,k,bi,bj)*viscA4*( (Zpj-Zij)*cosFacV(j,bi,bj) )
337 & + viscA4*( Dij-Dim ) )
338 #else /* MOM_VI_ORIGINAL_VISCA4 */
339 uD4 = (
340 & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
341 & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
342 vD4 = (
343 & recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
344 & *cosFacV(j,bi,bj)
345 & +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
346 ELSE
347 uD4 = viscA4D*
348 & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
349 & - viscA4Z*recip_hFacW(i,j,k,bi,bj)*
350 & ( Zip-Zij )*recip_DYG(i,j,bi,bj)
351 vD4 = viscA4Z*recip_hFacS(i,j,k,bi,bj)*
352 & cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
353 & + viscA4D* ( Dij-Dim )*recip_DYC(i,j,bi,bj)
354 #endif /* MOM_VI_ORIGINAL_VISCA4 */
355 ENDIF
356
357 uDissip(i,j) = uDissip(i,j) - uD4
358 vDissip(i,j) = vDissip(i,j) - vD4
359
360 ENDDO
361 ENDDO
362 ENDIF
363
364 #ifdef ALLOW_DIAGNOSTICS
365 IF (useDiagnostics) THEN
366 CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAH ',k,1,2,bi,bj,myThid)
367 CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4 ',k,1,2,bi,bj,myThid)
368 ENDIF
369 #endif
370
371 RETURN
372 END

  ViewVC Help
Powered by ViewVC 1.1.22