Context Navigation

navier5.f

main

Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c

Property mode set to 100755

File size: 91.0 KB

Line
1	c-----------------------------------------------------------------------
2	subroutine q_filter(wght)
3	c
4	c filter vx,vy,vz, and p by simple interpolation
5	c
6	include 'SIZE'
7	include 'TOTAL'
8	c
9	c
10	c These are the dimensions that we interpolate onto for v and p:
11	parameter(lxv=lx1-1)
12	parameter(lxp=lx2-1)
13	c
14	real intdv(lx1,lx1)
15	real intuv(lx1,lx1)
16	real intdp(lx1,lx1)
17	real intup(lx1,lx1)
18	real intv(lx1,lx1)
19	real intp(lx1,lx1)
20	c
21	save intdv
22	save intuv
23	save intdp
24	save intup
25	save intv
26	save intp
27
28	common /ctmp0/ intw,intt
29	$ , wk1,wk2
30	$ , zgmv,wgtv,zgmp,wgtp,tmax(100),omax(103)
31
32	real intw(lx1,lx1)
33	real intt(lx1,lx1)
34	real wk1 (lx1,lx1,lx1,lelt)
35	real wk2 (lx1,lx1,lx1)
36	real zgmv(lx1),wgtv(lx1),zgmp(lx1),wgtp(lx1)
37	c
38	c outpost arrays
39	parameter (lt=lx1ly1lz1*lelv)
40	common /scruz/ w1(lt),w2(lt),w3(lt),wt(lt)
41
42	character*18 sfmt
43
44	integer icalld
45	save icalld
46	data icalld /0/
47
48	logical if_fltv
49
50	ncut = param(101)+1
51
52	if(wght.le.0) return
53	if(ifaxis) call exitti('Filtering not supported w/ IFAXIS!$',1)
54	if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'apply q_filter ',
55	$ ncut, wght
56
57	imax = nid
58	imax = iglmax(imax,1)
59	jmax = iglmax(imax,1)
60	if (icalld.eq.0) then
61	icalld = 1
62	call build_new_filter(intv,zgm1,lx1,ncut,wght,nio)
63	elseif (icalld.lt.0) then ! old (std.) filter
64	icalld = 1
65	call zwgll(zgmv,wgtv,lxv)
66	call igllm(intuv,intw,zgmv,zgm1,lxv,lx1,lxv,lx1)
67	call igllm(intdv,intw,zgm1,zgmv,lx1,lxv,lx1,lxv)
68	c
69	call zwgl (zgmp,wgtp,lxp)
70	call iglm (intup,intw,zgmp,zgm2,lxp,lx2,lxp,lx2)
71	call iglm (intdp,intw,zgm2,zgmp,lx2,lxp,lx2,lxp)
72	c
73	c Multiply up and down interpolation into single operator
74	c
75	call mxm(intup,lx2,intdp,lxp,intp,lx2)
76	call mxm(intuv,lx1,intdv,lxv,intv,lx1)
77	c
78	c Weight the filter to make it a smooth (as opposed to truncated)
79	c decay in wave space
80
81	w0 = 1.-wght
82	call ident(intup,lx2)
83	call add2sxy(intp,wght,intup,w0,lx2*lx2)
84
85	call ident (intuv,lx1)
86	call add2sxy (intv ,wght,intuv,w0,lx1*lx1)
87
88	endif
89
90	ifldt = ifield
91	ifield = 1
92
93	if_fltv = .false.
94	if ( ifflow .and. .not. ifmhd ) if_fltv = .true.
95	if ( ifmhd ) if_fltv = .true.
96	if ( .not.ifbase ) if_fltv = .false. ! base-flow frozen
97	if ( .not. iffilter(1) ) if_fltv = .false.
98
99	if ( if_fltv ) then
100	call filterq(vx,intv,lx1,lz1,wk1,wk2,intt,if3d,umax)
101	call filterq(vy,intv,lx1,lz1,wk1,wk2,intt,if3d,vmax)
102	if (if3d)
103	$ call filterq(vz,intv,lx1,lz1,wk1,wk2,intt,if3d,wmax)
104
105	if (ifsplit.and..not.iflomach)
106	$ call filterq(pr,intv,lx1,lz1,wk1,wk2,intt,if3d,pmax)
107
108	if (ifpert) then
109	do j=1,npert
110	call filterq(vxp(1,j),intv,lx1,lz1,wk1,wk2,intt,if3d,umax)
111	call filterq(vyp(1,j),intv,lx1,lz1,wk1,wk2,intt,if3d,vmax)
112	if (if3d)
113	$ call filterq(vzp(1,j),intv,lx1,lz1,wk1,wk2,intt,if3d,wmax)
114	enddo
115	endif
116	endif
117
118	if (ifmhd.and.ifield.eq.ifldmhd) then
119	call filterq(bx,intv,lx1,lz1,wk1,wk2,intt,if3d,umax)
120	call filterq(by,intv,lx1,lz1,wk1,wk2,intt,if3d,vmax)
121	if (if3d)
122	$ call filterq(bz,intv,lx1,lz1,wk1,wk2,intt,if3d,wmax)
123	endif
124
125	mmax = 0
126	if (ifflow) then
127	c pmax = glmax(pmax,1)
128	c omax(1) = glmax(umax,1)
129	c omax(2) = glmax(vmax,1)
130	c omax(3) = glmax(wmax,1)
131	mmax = ldim
132	endif
133
134	nfldt = 1+npscal
135	do ifld=1,nfldt
136	ifield = ifld + 1
137	if (.not. iffilter(ifield)) goto 10
138	call filterq(t(1,1,1,1,ifld),intv,
139	$ lx1,lz1,wk1,wk2,intt,if3d,tmax(ifld))
140	if (ifpert) then
141	do j=1,npert
142	call filterq(tp(1,j,1),intv,lx1,lz1,wk1,wk2,intt,if3d,
143	$ wmax)
144	enddo
145	endif
146	10 mmax = mmax+1
147	c omax(mmax) = glmax(tmax(ifld),1)
148	enddo
149
150	c if (nio.eq.0) then
151	c if (if3d) then
152	c write(6,1) istep,ifield,umax,vmax,wmax
153	c else
154	c write(6,1) istep,ifield,umax,vmax
155	c endif
156	c 1 format(4x,i7,i3,' qfilt:',1p3e12.4)
157	c if(ifheat)
158	c & write(6,'(1p50e12.4)') (tmax(k),k=1,nfldt)
159	c endif
160
161	ifield = ifldt ! RESTORE ifield
162
163	return
164	end
165	c-----------------------------------------------------------------------
166	subroutine filterq(v,f,nx,nz,w1,w2,ft,if3d,dmax)
167	c
168	include 'SIZE'
169	include 'TSTEP'
170
171	real v(nxnxnz,nelt),w1(1),w2(1)
172	logical if3d
173	c
174	real f(nx,nx),ft(nx,nx)
175	c
176	integer e
177	c
178	call transpose(ft,nx,f,nx)
179	c
180	nxyz=nxnxnz
181	dmax = 0.
182
183	if (nio.eq.0 .and. loglevel.gt.2) write(6,*) 'call filterq',ifield
184	nel = nelfld(ifield)
185
186	if (if3d) then
187	do e=1,nel
188	c Filter
189	call copy(w2,v(1,e),nxyz)
190	call mxm(f,nx,w2,nx,w1,nx*nx)
191	i=1
192	j=1
193	do k=1,nx
194	call mxm(w1(i),nx,ft,nx,w2(j),nx)
195	i = i+nx*nx
196	j = j+nx*nx
197	enddo
198	call mxm (w2,nx*nx,ft,nx,w1,nx)
199	call sub3(w2,v(1,e),w1,nxyz)
200	call copy(v(1,e),w1,nxyz)
201	smax = vlamax(w2,nxyz)
202	dmax = max(dmax,abs(smax))
203	enddo
204	c
205	else
206	do e=1,nel
207	c Filter
208	call copy(w1,v(1,e),nxyz)
209	call mxm(f ,nx,w1,nx,w2,nx)
210	call mxm(w2,nx,ft,nx,w1,nx)
211	c
212	call sub3(w2,v(1,e),w1,nxyz)
213	call copy(v(1,e),w1,nxyz)
214	smax = vlamax(w2,nxyz)
215	dmax = max(dmax,abs(smax))
216	enddo
217	endif
218	c
219	return
220	end
221	c-----------------------------------------------------------------------
222	subroutine outmatx(a,m,n,io,name)
223	real a(m*n)
224	character*4 name
225	c
226	open(unit=io,file=name)
227	do i=1,m*n
228	write(io,1) a(i)
229	enddo
230	1 format(1p1e22.13)
231	close(unit=io)
232	c
233	return
234	end
235	c-----------------------------------------------------------------------
236	subroutine mappr(pm1,pm2,pa,pb)
237	c
238	INCLUDE 'SIZE'
239	INCLUDE 'TOTAL'
240	real pm1(lx1,ly1,lz1,lelv),pm2(lx2,ly2,lz2,lelv)
241	$ ,pa (lx1,ly2,lz2) ,pb (lx1,ly1,lz2)
242	c
243	C Map the pressure onto the velocity mesh
244	C
245	NGLOB1 = lx1ly1lz1*NELV
246	NYZ2 = ly2*lz2
247	NXY1 = lx1*ly1
248	NXYZ = lx1ly1lz1
249	C
250	IF (IFSPLIT) THEN
251	CALL COPY(PM1,PM2,NGLOB1)
252	ELSE
253	DO 1000 IEL=1,NELV
254	CALL MXM (IXM21,lx1,PM2(1,1,1,IEL),lx2,pa (1,1,1),NYZ2)
255	DO 100 IZ=1,lz2
256	CALL MXM (PA(1,1,IZ),lx1,IYTM21,ly2,PB(1,1,IZ),ly1)
257	100 CONTINUE
258	CALL MXM (PB(1,1,1),NXY1,IZTM21,lz2,PM1(1,1,1,IEL),lz1)
259	1000 CONTINUE
260
261	C Average the pressure on elemental boundaries
262	IFIELD=1
263	CALL DSSUM (PM1,lx1,ly1,lz1)
264	CALL COL2 (PM1,VMULT,NGLOB1)
265
266	ENDIF
267	C
268	C
269	return
270	end
271	c
272	c-----------------------------------------------------------------------
273	function facint_a(a,area,f,e)
274	c Integrate areal array a() on face f of element e. 27 June, 2012 pff
275
276	c f is in the preprocessor notation
277
278	include 'SIZE'
279	include 'TOPOL'
280	real a(lx1,lz1,6,lelt),area(lx1,lz1,6,lelt)
281
282	integer e,f
283
284	sum=0.0
285	do i=1,lx1*lz1
286	sum = sum + a(i,1,f,e)*area(i,1,f,e)
287	enddo
288
289	facint_a = sum
290
291	return
292	end
293	c-----------------------------------------------------------------------
294	function facint_v(a,area,f,e)
295	c Integrate volumetric array a() on face f of element e
296
297	c f is in the preprocessor notation
298	c fd is the dssum notation.
299	c 27 June, 2012 PFF
300
301	include 'SIZE'
302	include 'TOPOL'
303	real a(lx1,ly1,lz1,lelt),area(lx1,lz1,6,lelt)
304
305	integer e,f,fd
306
307	call dsset(lx1,ly1,lz1) ! set counters
308	fd = eface1(f)
309	js1 = skpdat(1,fd)
310	jf1 = skpdat(2,fd)
311	jskip1 = skpdat(3,fd)
312	js2 = skpdat(4,fd)
313	jf2 = skpdat(5,fd)
314	jskip2 = skpdat(6,fd)
315
316	sum=0.0
317	i = 0
318	do 100 j2=js2,jf2,jskip2
319	do 100 j1=js1,jf1,jskip1
320	i = i+1
321	sum = sum + a(j1,j2,1,e)*area(i,1,f,e)
322	100 continue
323
324	facint_v = sum
325
326	return
327	end
328	c-----------------------------------------------------------------------
329	function facint(a,b,area,ifc,ie)
330	c
331	C
332	C Take the dot product of A and B on the surface IFACE1 of element IE.
333	C
334	C IFACE1 is in the preprocessor notation
335	C IFACE is the dssum notation.
336	C 5 Jan 1989 15:12:22 PFF
337	C
338	INCLUDE 'SIZE'
339	INCLUDE 'TOPOL'
340	DIMENSION A (LX1,LY1,LZ1,lelv)
341	$ ,B (lx1,lz1,6,lelv)
342	$ ,area (lx1,lz1,6,lelv)
343	C
344	C Set up counters
345	C
346	CALL DSSET(lx1,ly1,lz1)
347	IFACE = EFACE1(IFC)
348	JS1 = SKPDAT(1,IFACE)
349	JF1 = SKPDAT(2,IFACE)
350	JSKIP1 = SKPDAT(3,IFACE)
351	JS2 = SKPDAT(4,IFACE)
352	JF2 = SKPDAT(5,IFACE)
353	JSKIP2 = SKPDAT(6,IFACE)
354	C
355	SUM=0.0
356	I = 0
357	DO 100 J2=JS2,JF2,JSKIP2
358	DO 100 J1=JS1,JF1,JSKIP1
359	I = I+1
360	SUM = SUM + A(J1,J2,1,ie)B(I,1,ifc,ie)area(I,1,ifc,ie)
361	c SUM = SUM + A(J1,J2,1,ie)B(J1,J2,1,ie)area(I,1,ifc,ie)
362	100 CONTINUE
363	C
364	facint = SUM
365	c
366	return
367	end
368	c-----------------------------------------------------------------------
369	function facint2(a,b,c,area,ifc,ie)
370	include 'SIZE'
371	include 'TOPOL'
372	dimension a (lx1,ly1,lz1,lelv)
373	$ , b (lx1,lz1,6,lelv)
374	$ , c (lx1,lz1,6,lelv)
375	$ , area (lx1,lz1,6,lelv)
376	call dsset(lx1,ly1,lz1)
377	iface = eface1(ifc)
378	js1 = skpdat(1,iface)
379	jf1 = skpdat(2,iface)
380	jskip1 = skpdat(3,iface)
381	js2 = skpdat(4,iface)
382	jf2 = skpdat(5,iface)
383	jskip2 = skpdat(6,iface)
384	sum=0.0
385	i=0
386	do j2=js2,jf2,jskip2
387	do j1=js1,jf1,jskip1
388	i=i+1
389	sum=sum+a(j1,j2,1,ie)b(i,1,ifc,ie)c(i,1,ifc,ie)
390	$ *area(i,1,ifc,ie)
391	end do
392	end do
393	facint2=sum
394	return
395	end
396	c-----------------------------------------------------------------------
397	subroutine out_csrmats(acsr,ia,ja,n,name9)
398	real acsr(1)
399	integer ia(1),ja(1)
400	c
401	character*9 name9
402	character*9 s(16)
403	c
404	nnz = ia(n+1)-ia(1)
405	c
406	write(6,1) name9,n,nnz
407	1 format(/,'CSR Mat:',a9,3x,'n =',i5,3x,'nnz =',i5,/)
408	c
409	n16 = min(n,16)
410	n29 = min(n,29)
411	do i=1,n29
412	call blank(s,9*16)
413	n1 = ia(i)
414	n2 = ia(i+1)-1
415	do jj=n1,n2
416	j = ja (jj)
417	a = acsr(jj)
418	if (a.ne.0..and.j.le.n16) write(s(j),9) a
419	enddo
420	write(6,16) (s(k),k=1,n16)
421	enddo
422	9 format(f9.4)
423	16 format(16a9)
424	c
425	return
426	end
427	c-----------------------------------------------------------------------
428	subroutine local_grad3(ur,us,ut,u,N,e,D,Dt)
429	c Output: ur,us,ut Input:u,N,e,D,Dt
430	real ur(0:N,0:N,0:N),us(0:N,0:N,0:N),ut(0:N,0:N,0:N)
431	real u (0:N,0:N,0:N,1)
432	real D (0:N,0:N),Dt(0:N,0:N)
433	integer e
434	c
435	m1 = N+1
436	m2 = m1*m1
437	c
438	call mxm(D ,m1,u(0,0,0,e),m1,ur,m2)
439	do k=0,N
440	call mxm(u(0,0,k,e),m1,Dt,m1,us(0,0,k),m1)
441	enddo
442	call mxm(u(0,0,0,e),m2,Dt,m1,ut,m1)
443	c
444	return
445	end
446	c-----------------------------------------------------------------------
447	subroutine local_grad2(ur,us,u,N,e,D,Dt)
448	c Output: ur,us Input:u,N,e,D,Dt
449	real ur(0:N,0:N),us(0:N,0:N)
450	real u (0:N,0:N,1)
451	real D (0:N,0:N),Dt(0:N,0:N)
452	integer e
453	c
454	m1 = N+1
455	c
456	call mxm(D ,m1,u(0,0,e),m1,ur,m1)
457	call mxm(u(0,0,e),m1,Dt,m1,us,m1)
458	c
459	return
460	end
461	c-----------------------------------------------------------------------
462	subroutine gradm1(ux,uy,uz,u)
463	c
464	c Compute gradient of T -- mesh 1 to mesh 1 (vel. to vel.)
465	c
466	include 'SIZE'
467	include 'DXYZ'
468	include 'GEOM'
469	include 'INPUT'
470	include 'TSTEP'
471	c
472	parameter (lxyz=lx1ly1lz1)
473	real ux(lxyz,1),uy(lxyz,1),uz(lxyz,1),u(lxyz,1)
474
475	common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
476
477	integer e
478
479	nxyz = lx1ly1lz1
480	ntot = nxyz*nelt
481
482	N = lx1-1
483	do e=1,nelt
484	if (if3d) then
485	call local_grad3(ur,us,ut,u,N,e,dxm1,dxtm1)
486	do i=1,lxyz
487	ux(i,e) = jacmi(i,e)(ur(i)rxm1(i,1,1,e)
488	$ + us(i)*sxm1(i,1,1,e)
489	$ + ut(i)*txm1(i,1,1,e) )
490	uy(i,e) = jacmi(i,e)(ur(i)rym1(i,1,1,e)
491	$ + us(i)*sym1(i,1,1,e)
492	$ + ut(i)*tym1(i,1,1,e) )
493	uz(i,e) = jacmi(i,e)(ur(i)rzm1(i,1,1,e)
494	$ + us(i)*szm1(i,1,1,e)
495	$ + ut(i)*tzm1(i,1,1,e) )
496	enddo
497	else
498	if (ifaxis) call setaxdy (ifrzer(e))
499	call local_grad2(ur,us,u,N,e,dxm1,dytm1)
500	do i=1,lxyz
501	ux(i,e) =jacmi(i,e)(ur(i)rxm1(i,1,1,e)
502	$ + us(i)*sxm1(i,1,1,e) )
503	uy(i,e) =jacmi(i,e)(ur(i)rym1(i,1,1,e)
504	$ + us(i)*sym1(i,1,1,e) )
505	enddo
506	endif
507	enddo
508	c
509	return
510	end
511	c-----------------------------------------------------------------------
512	subroutine comp_vort3(vort,work1,work2,u,v,w)
513
514	include 'SIZE'
515	include 'TOTAL'
516
517	parameter(lt=lx1ly1lz1*lelv)
518	real vort(lt,3),work1(1),work2(1),u(1),v(1),w(1)
519
520	parameter(lx=lx1ly1lz1)
521	real ur(lx),us(lx),ut(lx)
522	$ ,vr(lx),vs(lx),vt(lx)
523	$ ,wr(lx),ws(lx),wt(lx)
524	common /ctmp0/ ur,us,ut,vr,vs,vt,wr,ws,wt
525
526	integer e
527	real jacmil
528
529	ntot = lx1ly1lz1*nelt
530	nxyz = lx1ly1lz1
531	nx = lx1 - 1 ! Polynomial degree
532
533	if (ldim.eq.3) then
534	k=0
535	do e=1,nelt
536	call local_grad3(ur,us,ut,u,nx,e,dxm1,dxtm1)
537	call local_grad3(vr,vs,vt,v,nx,e,dxm1,dxtm1)
538	call local_grad3(wr,ws,wt,w,nx,e,dxm1,dxtm1)
539	do i=1,lx
540	jacmil = jacmi(i,e)
541	c vux=ur(i)rxm1(i,1,1,e)+us(i)sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e)
542	vuy=ur(i)rym1(i,1,1,e)+us(i)sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e)
543	vuz=ur(i)rzm1(i,1,1,e)+us(i)szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e)
544	vvx=vr(i)rxm1(i,1,1,e)+vs(i)sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e)
545	c vvy=vr(i)rym1(i,1,1,e)+vs(i)sym1(i,1,1,e)+vt(i)*tym1(i,1,1,e)
546	vvz=vr(i)rzm1(i,1,1,e)+vs(i)szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e)
547	vwx=wr(i)rxm1(i,1,1,e)+ws(i)sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e)
548	vwy=wr(i)rym1(i,1,1,e)+ws(i)sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e)
549	c vwz=wr(i)rzm1(i,1,1,e)+ws(i)szm1(i,1,1,e)+wt(i)*tzm1(i,1,1,e)
550
551	k = k+1
552	vort(k,1) = (vwy-vvz)*jacmil
553	vort(k,2) = (vuz-vwx)*jacmil
554	vort(k,3) = (vvx-vuy)*jacmil
555	c write(6,*) i,jacmil,vuy,vvx,k,e,' vort'
556	enddo
557	enddo
558
559	else ! 2D
560
561	k=0
562	do e=1,nelt
563	call local_grad2(ur,us,u,nx,e,dxm1,dxtm1)
564	call local_grad2(vr,vs,v,nx,e,dxm1,dxtm1)
565	do i=1,lx
566	c vux=ur(i)rxm1(i,1,1,e)+us(i)sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e)
567	vuy=ur(i)rym1(i,1,1,e)+us(i)sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e)
568	vvx=vr(i)rxm1(i,1,1,e)+vs(i)sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e)
569	c vvy=vr(i)rym1(i,1,1,e)+vs(i)sym1(i,1,1,e)+vt(i)*tym1(i,1,1,e)
570
571	k = k+1
572	vort(k,1) = (vvx-vuy)*jacmi(i,e)
573	enddo
574	enddo
575	endif
576	c
577	c Avg at bndry
578	c
579	ifielt = ifield
580	ifield = 1
581	if (if3d) then
582	do idim=1,ldim
583	call col2 (vort(1,idim),bm1,ntot)
584	call dssum (vort(1,idim),lx1,ly1,lz1)
585	call col2 (vort(1,idim),binvm1,ntot)
586	enddo
587	else
588	call col2 (vort,bm1,ntot)
589	call dssum (vort,lx1,ly1,lz1)
590	call col2 (vort,binvm1,ntot)
591	endif
592	ifield = ifielt
593	c
594	return
595	end
596	c-----------------------------------------------------------------------
597	subroutine surface_int(sint,sarea,a,e,f)
598	C
599	include 'SIZE'
600	include 'GEOM'
601	include 'PARALLEL'
602	include 'TOPOL'
603	real a(lx1,ly1,lz1,1)
604
605	integer e,f
606
607	call dsset(lx1,ly1,lz1)
608
609	iface = eface1(f)
610	js1 = skpdat(1,iface)
611	jf1 = skpdat(2,iface)
612	jskip1 = skpdat(3,iface)
613	js2 = skpdat(4,iface)
614	jf2 = skpdat(5,iface)
615	jskip2 = skpdat(6,iface)
616
617	sarea = 0.
618	sint = 0.
619	i = 0
620
621	do 100 j2=js2,jf2,jskip2
622	do 100 j1=js1,jf1,jskip1
623	i = i+1
624	sarea = sarea+area(i,1,f,e)
625	sint = sint +area(i,1,f,e)*a(j1,j2,1,e)
626	100 continue
627
628	return
629	end
630	c-----------------------------------------------------------------------
631	subroutine surface_flux(dq,qx,qy,qz,e,f,w)
632
633	include 'SIZE'
634	include 'GEOM'
635	include 'INPUT'
636	include 'PARALLEL'
637	include 'TOPOL'
638	parameter (l=lx1ly1lz1)
639
640	real qx(l,1),qy(l,1),qz(l,1),w(lx1,ly1,lz1)
641	integer e,f
642
643	call faccl3 (w,qx(1,e),unx(1,1,f,e),f)
644	call faddcl3 (w,qy(1,e),uny(1,1,f,e),f)
645	if (if3d) call faddcl3 (w,qz(1,e),unz(1,1,f,e),f)
646
647	call dsset(lx1,ly1,lz1)
648	iface = eface1(f)
649	js1 = skpdat(1,iface)
650	jf1 = skpdat(2,iface)
651	jskip1 = skpdat(3,iface)
652	js2 = skpdat(4,iface)
653	jf2 = skpdat(5,iface)
654	jskip2 = skpdat(6,iface)
655
656	dq = 0
657	i = 0
658	do 100 j2=js2,jf2,jskip2
659	do 100 j1=js1,jf1,jskip1
660	i = i+1
661	dq = dq + area(i,1,f,e)*w(j1,j2,1)
662	100 continue
663
664	return
665	end
666	c-----------------------------------------------------------------------
667	subroutine gaujordf(a,m,n,indr,indc,ipiv,ierr,rmult)
668	C
669	C Gauss-Jordan matrix inversion with full pivoting
670	c
671	c Num. Rec. p. 30, 2nd Ed., Fortran
672	c
673	C
674	C a is an m x n matrix
675	C rmult is a work array of dimension m
676	C
677	c
678	real a(m,n),rmult(m)
679	integer indr(m),indc(n),ipiv(n)
680
681	c call outmat(a,m,n,'ab4',n)
682	c do i=1,m
683	c write(6,1) (a(i,j),j=1,n)
684	c enddo
685	c 1 format('mat: ',1p6e12.4)
686
687	ierr = 0
688	eps = 1.e-9
689	call izero(ipiv,m)
690
691	do k=1,m
692	amx=0.
693	do i=1,m ! Pivot search
694	if (ipiv(i).ne.1) then
695	do j=1,m
696	if (ipiv(j).eq.0) then
697	if (abs(a(i,j)).ge.amx) then
698	amx = abs(a(i,j))
699	ir = i
700	jc = j
701	endif
702	elseif (ipiv(j).gt.1) then
703	ierr = -ipiv(j)
704	return
705	endif
706	enddo
707	endif
708	enddo
709	ipiv(jc) = ipiv(jc) + 1
710	c
711	c Swap rows
712	if (ir.ne.jc) then
713	do j=1,n
714	tmp = a(ir,j)
715	a(ir,j) = a(jc,j)
716	a(jc,j) = tmp
717	enddo
718	endif
719	indr(k)=ir
720	indc(k)=jc
721	c write(6 ,*) k,' Piv:',jc,a(jc,jc)
722	c write(28,*) k,' Piv:',jc,a(jc,jc)
723	if (abs(a(jc,jc)).lt.eps) then
724	write(6 ,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
725	write(28,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
726	ierr = jc
727	call exitt
728	return
729	endif
730	piv = 1./a(jc,jc)
731	a(jc,jc)=1.
732	do j=1,n
733	a(jc,j) = a(jc,j)*piv
734	enddo
735	c
736	do j=1,n
737	work = a(jc,j)
738	a(jc,j) = a(1 ,j)
739	a(1 ,j) = work
740	enddo
741	do i=2,m
742	rmult(i) = a(i,jc)
743	a(i,jc) = 0.
744	enddo
745	c
746	do j=1,n
747	do i=2,m
748	a(i,j) = a(i,j) - rmult(i)*a(1,j)
749	enddo
750	enddo
751	c
752	do j=1,n
753	work = a(jc,j)
754	a(jc,j) = a(1 ,j)
755	a(1 ,j) = work
756	enddo
757	c
758	c do i=1,m
759	c if (i.ne.jc) then
760	c rmult = a(i,jc)
761	c a(i,jc) = 0.
762	c do j=1,n
763	c a(i,j) = a(i,j) - rmult*a(jc,j)
764	c enddo
765	c endif
766	c enddo
767	c
768	enddo
769	c
770	c Unscramble matrix
771	do j=m,1,-1
772	if (indr(j).ne.indc(j)) then
773	do i=1,m
774	tmp=a(i,indr(j))
775	a(i,indr(j))=a(i,indc(j))
776	a(i,indc(j))=tmp
777	enddo
778	endif
779	enddo
780	c
781	return
782	end
783	c-----------------------------------------------------------------------
784	subroutine legendre_poly(L,x,N)
785	c
786	c Evaluate Legendre polynomials of degrees 0-N at point x
787	c
788	real L(0:N)
789	c
790	L(0) = 1.
791	L(1) = x
792	c
793	do j=2,N
794	L(j) = ( (2j-1) x * L(j-1) - (j-1) * L(j-2) ) / j
795	enddo
796	c
797	return
798	end
799	c-----------------------------------------------------------------------
800	subroutine build_new_filter(intv,zpts,nx,kut,wght,nio)
801	c
802	c This routing builds a 1D filter with a transfer function that
803	c looks like:
804	c
805	c
806	c ^
807	c d_k \|
808	c \| \|
809	c 1 \|__________ _v_
810	c \| -_
811	c \| \ wght
812	c \| \ ___
813	c \| \| ^
814	c 0 \|-------------\|---\|>
815	c
816	c 0 c N k-->
817	c
818	c Where c := N-kut is the point below which d_k = 1.
819	c
820	c
821	c
822	c Here, nx = number of points
823	c
824	real intv(nx,nx),zpts(nx)
825	c
826	parameter (lm=40)
827	parameter (lm2=lm*lm)
828	real phi(lm2),pht(lm2),diag(lm2),rmult(lm),Lj(lm)
829	integer indr(lm),indc(lm),ipiv(lm)
830	c
831	if (nx.gt.lm) then
832	write(6,*) 'ABORT in build_new_filter:',nx,lm
833	call exitt
834	endif
835	c
836	kj = 0
837	n = nx-1
838	do j=1,nx
839	z = zpts(j)
840	call legendre_poly(Lj,z,n)
841	kj = kj+1
842	pht(kj) = Lj(1)
843	kj = kj+1
844	pht(kj) = Lj(2)
845	do k=3,nx
846	kj = kj+1
847	pht(kj) = Lj(k)-Lj(k-2)
848	enddo
849	enddo
850	call transpose (phi,nx,pht,nx)
851	call copy (pht,phi,nx*nx)
852	call gaujordf (pht,nx,nx,indr,indc,ipiv,ierr,rmult)
853	c
854	c Set up transfer function
855	c
856	call ident (diag,nx)
857	c
858	k0 = nx-kut
859	do k=k0+1,nx
860	kk = k+nx*(k-1)
861	amp = wght(k-k0)(k-k0)/(kut*kut) ! quadratic growth
862	diag(kk) = 1.-amp
863	enddo
864	c
865	call mxm (diag,nx,pht,nx,intv,nx) ! -1
866	call mxm (phi ,nx,intv,nx,pht,nx) ! V D V
867	call copy (intv,pht,nx*nx)
868	c
869	do k=1,nx*nx
870	pht(k) = 1.-diag(k)
871	enddo
872	np1 = nx+1
873	if (nio.eq.0) then
874	write(6,6) 'filt amp',(pht (k),k=1,nx*nx,np1)
875	write(6,6) 'filt trn',(diag(k),k=1,nx*nx,np1)
876	6 format(a8,16f7.4,6(/,8x,16f7.4))
877	endif
878	c
879	return
880	end
881	c-----------------------------------------------------------------------
882	subroutine avg_all
883	c
884	c This routine computes running averages E(X),E(X^2),E(X*Y)
885	c and outputs to avg.fld, rms.fld, and rm2.fld for all
886	c fields.
887	c
888	c E denotes the expected value operator and X,Y two
889	c real valued random variables.
890	c
891	c variances and covariances can be computed in a post-processing step:
892	c
893	c var(X) := E(X^X) - E(X)*E(X)
894	c cov(X,Y) := E(XY) - E(X)E(Y)
895	c
896	c Note: The E-operator is linear, in the sense that the expected
897	c value is given by E(X) = 1/N * sum[ E(X)_i ], where E(X)_i
898	c is the expected value of the sub-ensemble i (i=1...N).
899	c
900	include 'SIZE'
901	include 'TOTAL'
902	include 'AVG'
903
904	logical ifverbose
905	integer icalld
906	save icalld
907	data icalld /0/
908
909	if (ax1.ne.lx1 .or. ay1.ne.ly1 .or. az1.ne.lz1) then
910	if(nid.eq.0) write(6,*)
911	$ 'ABORT: wrong size of ax1,ay1,az1 in avg_all(), check SIZE!'
912	call exitt
913	endif
914	if (ax2.ne.lx2 .or. ay2.ne.ay2 .or. az2.ne.lz2) then
915	if(nid.eq.0) write(6,*)
916	$ 'ABORT: wrong size of ax2,ay2,az2 in avg_all(), check SIZE!'
917	call exitt
918	endif
919
920	ntot = lx1ly1lz1*nelv
921	ntott = lx1ly1lz1*nelt
922	nto2 = lx2ly2lz2*nelv
923
924	! initialization
925	if (icalld.eq.0) then
926	icalld = icalld + 1
927	atime = 0.
928	timel = time
929
930	call rzero(uavg,ntot)
931	call rzero(vavg,ntot)
932	call rzero(wavg,ntot)
933	call rzero(pavg,nto2)
934	do i = 1,ldimt
935	call rzero(tavg(1,1,1,1,i),ntott)
936	enddo
937
938	call rzero(urms,ntot)
939	call rzero(vrms,ntot)
940	call rzero(wrms,ntot)
941	call rzero(prms,nto2)
942	do i = 1,ldimt
943	call rzero(trms(1,1,1,1,i),ntott)
944	enddo
945
946	call rzero(vwms,ntot)
947	call rzero(wums,ntot)
948	call rzero(uvms,ntot)
949	endif
950
951	dtime = time - timel
952	atime = atime + dtime
953
954	! dump freq
955	iastep = param(68)
956	if (iastep.eq.0) iastep=param(15) ! same as iostep
957	if (iastep.eq.0) iastep=500
958
959	ifverbose=.false.
960	if (istep.le.10) ifverbose=.true.
961	if (mod(istep,iastep).eq.0) ifverbose=.true.
962
963	if (atime.ne.0..and.dtime.ne.0.) then
964	if(nio.eq.0) write(6,*) 'Compute statistics ...'
965	beta = dtime/atime
966	alpha = 1.-beta
967	! compute averages E(X)
968	call avg1 (uavg,vx,alpha,beta,ntot ,'um ',ifverbose)
969	call avg1 (vavg,vy,alpha,beta,ntot ,'vm ',ifverbose)
970	call avg1 (wavg,vz,alpha,beta,ntot ,'wm ',ifverbose)
971	call avg1 (pavg,pr,alpha,beta,nto2 ,'prm ',ifverbose)
972	call avg1 (tavg,t ,alpha,beta,ntott,'tm ',ifverbose)
973	do i = 2,ldimt
974	call avg1 (tavg(1,1,1,1,i),t(1,1,1,1,i),alpha,beta,
975	& ntott,'psav',ifverbose)
976	enddo
977
978	! compute averages E(X^2)
979	call avg2 (urms,vx,alpha,beta,ntot ,'ums ',ifverbose)
980	call avg2 (vrms,vy,alpha,beta,ntot ,'vms ',ifverbose)
981	call avg2 (wrms,vz,alpha,beta,ntot ,'wms ',ifverbose)
982	call avg2 (prms,pr,alpha,beta,nto2 ,'prms',ifverbose)
983	call avg2 (trms,t ,alpha,beta,ntott,'tms ',ifverbose)
984	do i = 2,ldimt
985	call avg2 (trms(1,1,1,1,i),t(1,1,1,1,i),alpha,beta,
986	& ntott,'psms',ifverbose)
987	enddo
988
989	! compute averages E(X*Y) (for now just for the velocities)
990	call avg3 (uvms,vx,vy,alpha,beta,ntot,'uvm ',ifverbose)
991	call avg3 (vwms,vy,vz,alpha,beta,ntot,'vwm ',ifverbose)
992	call avg3 (wums,vz,vx,alpha,beta,ntot,'wum ',ifverbose)
993	endif
994	c
995	c-----------------------------------------------------------------------
996	if ( (mod(istep,iastep).eq.0.and.istep.gt.1) .or.lastep.eq.1) then
997
998	time_temp = time
999	time = atime ! Output the duration of this avg
1000	dtmp = param(63)
1001	param(63) = 1 ! Enforce 64-bit output
1002
1003	call outpost2(uavg,vavg,wavg,pavg,tavg,ldimt,'avg')
1004	call outpost2(urms,vrms,wrms,prms,trms,ldimt,'rms')
1005	call outpost (uvms,vwms,wums,prms,trms, 'rm2')
1006
1007	param(63) = dtmp
1008	atime = 0.
1009	time = time_temp ! Restore clock
1010
1011	endif
1012	c
1013	timel = time
1014	c
1015	return
1016	end
1017	c-----------------------------------------------------------------------
1018	subroutine avg1(avg,f,alpha,beta,n,name,ifverbose)
1019	include 'SIZE'
1020	include 'TSTEP'
1021	c
1022	real avg(n),f(n)
1023	character*4 name
1024	logical ifverbose
1025	c
1026	do k=1,n
1027	avg(k) = alphaavg(k) + betaf(k)
1028	enddo
1029	c
1030	if (ifverbose) then
1031	avgmax = glmax(avg,n)
1032	avgmin = glmin(avg,n)
1033	if (nio.eq.0) write(6,1) istep,time,avgmin,avgmax
1034	$ ,alpha,beta,name
1035	1 format(i9,1p5e13.5,1x,a4,' av1mnx')
1036	endif
1037	c
1038	return
1039	end
1040	c-----------------------------------------------------------------------
1041	subroutine avg2(avg,f,alpha,beta,n,name,ifverbose)
1042	include 'SIZE'
1043	include 'TSTEP'
1044	c
1045	real avg(n),f(n)
1046	character*4 name
1047	logical ifverbose
1048	c
1049	do k=1,n
1050	avg(k) = alphaavg(k) + betaf(k)*f(k)
1051	enddo
1052	c
1053	if (ifverbose) then
1054	avgmax = glmax(avg,n)
1055	avgmin = glmin(avg,n)
1056	if (nio.eq.0) write(6,1) istep,time,avgmin,avgmax
1057	$ ,alpha,beta,name
1058	1 format(i9,1p5e13.5,1x,a4,' av2mnx')
1059	endif
1060	c
1061	return
1062	end
1063	c-----------------------------------------------------------------------
1064	subroutine avg3(avg,f,g,alpha,beta,n,name,ifverbose)
1065	include 'SIZE'
1066	include 'TSTEP'
1067	c
1068	real avg(n),f(n),g(n)
1069	character*4 name
1070	logical ifverbose
1071	c
1072	do k=1,n
1073	avg(k) = alphaavg(k) + betaf(k)*g(k)
1074	enddo
1075	c
1076	if (ifverbose) then
1077	avgmax = glmax(avg,n)
1078	avgmin = glmin(avg,n)
1079	if (nio.eq.0) write(6,1) istep,time,avgmin,avgmax
1080	$ ,alpha,beta,name
1081	1 format(i9,1p5e13.5,1x,a4,' av3mnx')
1082	endif
1083	c
1084	return
1085	end
1086	c-----------------------------------------------------------------------
1087	subroutine build_legend_transform(Lj,Ljt,zpts,nx)
1088	c
1089	real Lj(nxnx),Ljt(nxnx),zpts(nx)
1090	c
1091	parameter (lm=90)
1092	integer indr(lm),indc(lm),ipiv(lm)
1093	real rmult(lm)
1094	c
1095	if (nx.gt.lm) then
1096	write(6,*) 'ABORT in build_legend_transform:',nx,lm
1097	call exitt
1098	endif
1099	c
1100	j = 1
1101	n = nx-1
1102	do i=1,nx
1103	z = zpts(i)
1104	call legendre_poly(Lj(j),z,n) ! Return Lk(z), k=0,...,n
1105	j = j+nx
1106	enddo
1107	call transpose1(Lj,nx)
1108	c call outmat(Lj,nx,nx,'Lj ',n)
1109	c call exitt
1110	call gaujordf (Lj,nx,nx,indr,indc,ipiv,ierr,rmult)
1111	call transpose (Ljt,nx,Lj,nx)
1112	c
1113	return
1114	end
1115	c-----------------------------------------------------------------------
1116	subroutine local_err_est(err,u,nx,Lj,Ljt,uh,w,if3d)
1117	c
1118	c Local error estimates for u_e
1119	c
1120	include 'SIZE'
1121	real err(5,2),u(1),uh(nx,nx,nx),w(1),Lj(1),Ljt(1)
1122	logical if3d
1123	c
1124	call rzero(err,10)
1125	c
1126	nxyz = nx**ldim
1127	utot = vlsc2(u,u,nxyz)
1128	if (utot.eq.0) return
1129	c
1130	call tensr3(uh,nx,u,nx,Lj,Ljt,Ljt,w) ! Go to Legendre space
1131	c
1132	c
1133	c Get energy in low modes
1134	c
1135	m = nx-2
1136	c
1137	if (if3d) then
1138	amp2_l = 0.
1139	do k=1,m
1140	do j=1,m
1141	do i=1,m
1142	amp2_l = amp2_l + uh(i,j,k)**2
1143	enddo
1144	enddo
1145	enddo
1146	c
1147	c Energy in each spatial direction
1148	c
1149	amp2_t = 0
1150	do k=m+1,nx
1151	do j=1,m
1152	do i=1,m
1153	amp2_t = amp2_t + uh(i,j,k)**2
1154	enddo
1155	enddo
1156	enddo
1157	c
1158	amp2_s = 0
1159	do k=1,m
1160	do j=m+1,nx
1161	do i=1,m
1162	amp2_s = amp2_s + uh(i,j,k)**2
1163	enddo
1164	enddo
1165	enddo
1166	c
1167	amp2_r = 0
1168	do k=1,m
1169	do j=1,m
1170	do i=m+1,nx
1171	amp2_r = amp2_r + uh(i,j,k)**2
1172	enddo
1173	enddo
1174	enddo
1175	c
1176	amp2_h = 0
1177	do k=m+1,nx
1178	do j=m+1,nx
1179	do i=m+1,nx
1180	amp2_h = amp2_h + uh(i,j,k)**2
1181	enddo
1182	enddo
1183	enddo
1184	c
1185	etot = amp2_l + amp2_r + amp2_s + amp2_t + amp2_h
1186	c
1187	relr = amp2_r / (amp2_r + amp2_l)
1188	rels = amp2_s / (amp2_s + amp2_l)
1189	relt = amp2_t / (amp2_t + amp2_l)
1190	relh = (amp2_r + amp2_s + amp2_t + amp2_h) / etot
1191	c
1192	else
1193	k = 1
1194	amp2_l = 0.
1195	do j=1,m
1196	do i=1,m
1197	amp2_l = amp2_l + uh(i,j,k)**2
1198	enddo
1199	enddo
1200	if (amp2_l.eq.0) return
1201	c
1202	c Energy in each spatial direction
1203	c
1204	amp2_t = 0
1205	c
1206	amp2_s = 0
1207	do j=m+1,nx
1208	do i=1,m
1209	amp2_s = amp2_s + uh(i,j,k)**2
1210	enddo
1211	enddo
1212	c
1213	amp2_r = 0
1214	do j=1,m
1215	do i=m+1,nx
1216	amp2_r = amp2_r + uh(i,j,k)**2
1217	enddo
1218	enddo
1219	c
1220	amp2_h = 0
1221	do j=m+1,nx
1222	do i=m+1,nx
1223	amp2_h = amp2_h + uh(i,j,k)**2
1224	enddo
1225	enddo
1226	c
1227	etot = amp2_l + amp2_r + amp2_s + amp2_h
1228	c
1229	relr = amp2_r / (amp2_r + amp2_l)
1230	rels = amp2_s / (amp2_s + amp2_l)
1231	relt = 0
1232	relh = (amp2_r + amp2_s + amp2_h) / etot
1233	c
1234	endif
1235	c
1236	err (1,1) = sqrt(relr)
1237	err (2,1) = sqrt(rels)
1238	if (if3d) err (3,1) = sqrt(relt)
1239	err (4,1) = sqrt(relh)
1240	err (5,1) = sqrt(etot)
1241	c
1242	err (1,2) = sqrt(amp2_r)
1243	err (2,2) = sqrt(amp2_s)
1244	if (if3d) err (3,2) = sqrt(amp2_t)
1245	err (4,2) = sqrt(amp2_h)
1246	err (5,2) = sqrt(utot)
1247	c
1248	return
1249	end
1250	c-----------------------------------------------------------------------
1251	subroutine get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
1252	integer ex,ey,ez,eg
1253	c
1254	nelxy = nelx*nely
1255	c
1256	ez = 1 + (eg-1)/nelxy
1257	ey = mod1 (eg,nelxy)
1258	ey = 1 + (ey-1)/nelx
1259	ex = mod1 (eg,nelx)
1260	c
1261	return
1262	end
1263	c-----------------------------------------------------------------------
1264	subroutine dump_header2d(excode,nx,ny,nlx,nly,ierr)
1265
1266	include 'SIZE'
1267	include 'TOTAL'
1268
1269	character*2 excode(15)
1270
1271	real*4 test_pattern
1272
1273	character*1 fhdfle1(132)
1274	character*132 fhdfle
1275	equivalence (fhdfle,fhdfle1)
1276
1277	jstep = istep
1278	if (jstep.gt.10000) jstep = jstep / 10
1279	if (jstep.gt.10000) jstep = jstep / 10
1280	if (jstep.gt.10000) jstep = jstep / 10
1281	if (jstep.gt.10000) jstep = jstep / 10
1282	if (jstep.gt.10000) jstep = jstep / 10
1283
1284	nlxy = nlx*nly
1285	nzz = 1
1286
1287	c write(6,'(4i4,1PE14.7,i5,1x,15a2,1x,a12)')
1288	c $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
1289	c $ 'NELT,NX,NY,N'
1290	c
1291	p66 = 0.
1292	ierr= 0
1293	IF (p66.EQ.1.0) THEN
1294	C unformatted i/o
1295	WRITE(24)
1296	$ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15)
1297	ELSEIF (p66.EQ.3.0) THEN
1298	C formatted i/o to header file
1299	WRITE(27,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
1300	$ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
1301	$ 'NELT,NX,NY,N'
1302	ELSEIF (p66.eq.4.0) THEN
1303	C formatted i/o to header file
1304	WRITE(fhdfle,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
1305	$ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
1306	$ ' 4 NELT,NX,NY,N'
1307	call byte_write(fhdfle,20,ierr)
1308	ELSEIF (p66.eq.5.0) THEN
1309	C formatted i/o to header file
1310	WRITE(fhdfle,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
1311	$ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
1312	$ ' 8 NELT,NX,NY,N'
1313	call byte_write(fhdfle,20,ierr)
1314	ELSE
1315	C formatted i/o
1316	WRITE(24,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
1317	$ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
1318	$ 'NELT,NX,NY,N'
1319	ENDIF
1320	C cdrror is a dummy cerror value for now.
1321	CDRROR=0.0
1322	IF (p66.EQ.1.0) THEN
1323	C unformatted i/o
1324	WRITE(24)(CDRROR,I=1,nlxy)
1325	ELSEIF (p66.eq.3. .or. p66.eq.4.0) then
1326	C write byte-ordering test pattern to byte file...
1327	test_pattern = 6.54321
1328	call byte_write(test_pattern,1,ierr)
1329	ELSEIF (p66.eq.5.) then
1330	test_pattern8 = 6.54321
1331	call byte_write(test_pattern8,2,ierr)
1332	ELSE
1333	C formatted i/o
1334	WRITE(24,'(6G11.4)')(CDRROR,I=1,nlxy)
1335	ENDIF
1336
1337	return
1338	end
1339	c-----------------------------------------------------------------------
1340	subroutine outfld2d_p(u,v,w,nx,ny,nlx,nly,name,ifld,jid,npido,ir)
1341
1342	include 'SIZE'
1343	include 'TOTAL'
1344
1345	integer icalld
1346	save icalld
1347	data icalld /0/
1348
1349	real u(nxnynlx*nly)
1350	real v(nxnynlx*nly)
1351	real w(nxnynlx*nly)
1352	character*4 name
1353
1354	character*2 excode(15)
1355	character*12 fm
1356	character*20 outfile
1357
1358	character*1 slash,dot
1359	save slash,dot
1360	data slash,dot / '/' , '.' /
1361
1362	icalld = icalld+1
1363
1364	call blank(excode,30)
1365	excode(4) = 'U '
1366	excode(5) = ' '
1367	excode(6) = 'T '
1368	nthings = 3
1369	ir = 0 !error code for dump_header2d
1370
1371	call blank(outfile,20)
1372	if (npido.lt.100) then
1373	if (ifld.lt.100) then
1374	write(outfile,22) jid,slash,name,ifld
1375	22 format('B',i2.2,a1,a4,'.fld',i2.2)
1376	elseif (ifld.lt.1000) then
1377	write(outfile,23) jid,slash,name,ifld
1378	23 format('B',i2.2,a1,a4,'.fld',i3)
1379	elseif (ifld.lt.10000) then
1380	write(outfile,24) jid,slash,name,ifld
1381	24 format('B',i2.2,a1,a4,'.fld',i4)
1382	elseif (ifld.lt.100000) then
1383	write(outfile,25) jid,slash,name,ifld
1384	25 format('B',i2.2,a1,a4,'.fld',i5)
1385	elseif (ifld.lt.1000000) then
1386	write(outfile,26) jid,slash,name,ifld
1387	26 format('B',i2.2,a1,a4,'.fld',i6)
1388	endif
1389	else
1390	if (ifld.lt.100) then
1391	write(outfile,32) jid,slash,name,ifld
1392	32 format('B',i3.3,a1,a4,'.fld',i2.2)
1393	elseif (ifld.lt.1000) then
1394	write(outfile,33) jid,slash,name,ifld
1395	33 format('B',i3.3,a1,a4,'.fld',i3)
1396	elseif (ifld.lt.10000) then
1397	write(outfile,34) jid,slash,name,ifld
1398	34 format('B',i3.3,a1,a4,'.fld',i4)
1399	elseif (ifld.lt.100000) then
1400	write(outfile,35) jid,slash,name,ifld
1401	35 format('B',i3.3,a1,a4,'.fld',i5)
1402	elseif (ifld.lt.1000000) then
1403	write(outfile,36) jid,slash,name,ifld
1404	36 format('B',i3.3,a1,a4,'.fld',i6)
1405	endif
1406	endif
1407
1408	if (icalld.le.4) write(6,*) nid,outfile,' OPEN',nlx,nly
1409	open(unit=24,file=outfile,status='unknown')
1410	call dump_header2d(excode,nx,ny,nlx,nly,ir)
1411
1412	n = nxnynlx*nly
1413	write(fm,10) nthings
1414	write(24,fm) (u(i),v(i),w(i),i=1,n)
1415	10 format('(1p',i1,'e14.6)')
1416	close(24)
1417
1418	return
1419	end
1420	c-----------------------------------------------------------------------
1421	subroutine outfld2d(u,v,w,nx,ny,nlx,nly,name,ifld)
1422
1423	include 'SIZE'
1424	include 'TOTAL'
1425
1426	real u(nxnynlx*nly)
1427	real v(nxnynlx*nly)
1428	real w(nxnynlx*nly)
1429	character*3 name
1430
1431	character*2 excode(15)
1432	character*12 fm
1433	character*20 outfile
1434
1435	c if (istep.le.10) write(6,*) nid,' in out2d:',iz
1436
1437	call blank(excode,30)
1438	c
1439	c excode(1) = 'X '
1440	c excode(2) = 'Y '
1441	c excode(3) = 'U '
1442	c excode(4) = 'V '
1443	c excode(5) = 'P '
1444	c excode(6) = 'T '
1445	c
1446	excode(4) = 'U '
1447	excode(5) = ' '
1448	excode(6) = 'T '
1449	nthings = 3
1450
1451	ierr = 0
1452	if (nid.eq.0) then
1453	call blank(outfile,20)
1454	if (ifld.lt.100) then
1455	write(outfile,2) name,ifld
1456	2 format(a3,'2d.fld',i2.2)
1457	elseif (ifld.lt.1000) then
1458	write(outfile,3) name,ifld
1459	3 format(a3,'2d.fld',i3)
1460	elseif (ifld.lt.10000) then
1461	write(outfile,4) name,ifld
1462	4 format(a3,'2d.fld',i4)
1463	elseif (ifld.lt.100000) then
1464	write(outfile,5) name,ifld
1465	5 format(a3,'2d.fld',i5)
1466	elseif (ifld.lt.1000000) then
1467	write(outfile,6) name,ifld
1468	6 format(a3,'2d.fld',i6)
1469	endif
1470	open(unit=24,file=outfile,status='unknown')
1471	call dump_header2d(excode,nx,ny,nlx,nly,ierr)
1472
1473	n = nxnynlx*nly
1474	write(fm,10) nthings
1475	c write(6,*) fm
1476	c call exitt
1477	write(24,fm) (u(i),v(i),w(i),i=1,n)
1478	10 format('(1p',i1,'e14.6)')
1479	c 10 format('''(1p',i1,'e15.7)''')
1480	c 10 format(1p7e15.7)
1481	c
1482	close(24)
1483	endif
1484	call err_chk(ierr,'Error using byte_write,outfld2d. Abort.$')
1485
1486	return
1487	end
1488	c-----------------------------------------------------------------------
1489	subroutine drgtrq(dgtq,xm0,ym0,zm0,sij,pm1,visc,f,e)
1490	c
1491	INCLUDE 'SIZE'
1492	INCLUDE 'GEOM'
1493	INCLUDE 'INPUT'
1494	INCLUDE 'TOPOL'
1495	INCLUDE 'TSTEP'
1496	c
1497	real dgtq(3,4)
1498	real xm0 (lx1,ly1,lz1,lelt)
1499	real ym0 (lx1,ly1,lz1,lelt)
1500	real zm0 (lx1,ly1,lz1,lelt)
1501	real sij (lx1,ly1,lz1,3*ldim-3,lelv)
1502	real pm1 (lx1,ly1,lz1,lelv)
1503	real visc(lx1,ly1,lz1,lelv)
1504	c
1505	real dg(3,2)
1506	c
1507	integer f,e,pf
1508	real n1,n2,n3
1509	c
1510	call dsset(lx1,ly1,lz1) ! set up counters
1511	pf = eface1(f) ! convert from preproc. notation
1512	js1 = skpdat(1,pf)
1513	jf1 = skpdat(2,pf)
1514	jskip1 = skpdat(3,pf)
1515	js2 = skpdat(4,pf)
1516	jf2 = skpdat(5,pf)
1517	jskip2 = skpdat(6,pf)
1518	C
1519	call rzero(dgtq,12)
1520	c
1521	if (if3d.or.ifaxis) then
1522	i = 0
1523	a = 0
1524	do j2=js2,jf2,jskip2
1525	do j1=js1,jf1,jskip1
1526	i = i+1
1527	n1 = unx(i,1,f,e)*area(i,1,f,e)
1528	n2 = uny(i,1,f,e)*area(i,1,f,e)
1529	n3 = unz(i,1,f,e)*area(i,1,f,e)
1530	a = a + area(i,1,f,e)
1531	c
1532	v = visc(j1,j2,1,e)
1533	c
1534	s11 = sij(j1,j2,1,1,e)
1535	s21 = sij(j1,j2,1,4,e)
1536	s31 = sij(j1,j2,1,6,e)
1537	c
1538	s12 = sij(j1,j2,1,4,e)
1539	s22 = sij(j1,j2,1,2,e)
1540	s32 = sij(j1,j2,1,5,e)
1541	c
1542	s13 = sij(j1,j2,1,6,e)
1543	s23 = sij(j1,j2,1,5,e)
1544	s33 = sij(j1,j2,1,3,e)
1545	c
1546	dg(1,1) = pm1(j1,j2,1,e)*n1 ! pressure drag
1547	dg(2,1) = pm1(j1,j2,1,e)*n2
1548	dg(3,1) = pm1(j1,j2,1,e)*n3
1549	c
1550	dg(1,2) = -v(s11n1 + s12n2 + s13n3) ! viscous drag
1551	dg(2,2) = -v(s21n1 + s22n2 + s23n3)
1552	dg(3,2) = -v(s31n1 + s32n2 + s33n3)
1553	c
1554	r1 = xm0(j1,j2,1,e)
1555	r2 = ym0(j1,j2,1,e)
1556	r3 = zm0(j1,j2,1,e)
1557	c
1558	do l=1,2
1559	do k=1,3
1560	dgtq(k,l) = dgtq(k,l) + dg(k,l)
1561	enddo
1562	enddo
1563	c
1564	dgtq(1,3) = dgtq(1,3) + (r2dg(3,1)-r3dg(2,1)) ! pressure
1565	dgtq(2,3) = dgtq(2,3) + (r3dg(1,1)-r1dg(3,1)) ! torque
1566	dgtq(3,3) = dgtq(3,3) + (r1dg(2,1)-r2dg(1,1))
1567	c
1568	dgtq(1,4) = dgtq(1,4) + (r2dg(3,2)-r3dg(2,2)) ! viscous
1569	dgtq(2,4) = dgtq(2,4) + (r3dg(1,2)-r1dg(3,2)) ! torque
1570	dgtq(3,4) = dgtq(3,4) + (r1dg(2,2)-r2dg(1,2))
1571	enddo
1572	enddo
1573
1574	else ! 2D
1575
1576	i = 0
1577	a = 0
1578	do j2=js2,jf2,jskip2
1579	do j1=js1,jf1,jskip1
1580	i = i+1
1581	n1 = unx(i,1,f,e)*area(i,1,f,e)
1582	n2 = uny(i,1,f,e)*area(i,1,f,e)
1583	a = a + area(i,1,f,e)
1584	v = visc(j1,j2,1,e)
1585
1586	s11 = sij(j1,j2,1,1,e)
1587	s12 = sij(j1,j2,1,3,e)
1588	s21 = sij(j1,j2,1,3,e)
1589	s22 = sij(j1,j2,1,2,e)
1590
1591	dg(1,1) = pm1(j1,j2,1,e)*n1 ! pressure drag
1592	dg(2,1) = pm1(j1,j2,1,e)*n2
1593	dg(3,1) = 0
1594
1595	dg(1,2) = -v(s11n1 + s12*n2) ! viscous drag
1596	dg(2,2) = -v(s21n1 + s22*n2)
1597	dg(3,2) = 0.
1598
1599	r1 = xm0(j1,j2,1,e)
1600	r2 = ym0(j1,j2,1,e)
1601	r3 = 0.
1602
1603	do l=1,2
1604	do k=1,3
1605	dgtq(k,l) = dgtq(k,l) + dg(k,l)
1606	enddo
1607	enddo
1608
1609	dgtq(1,3) = 0! dgtq(1,3) + (r2dg(3,1)-r3dg(2,1)) ! pressure
1610	dgtq(2,3) = 0! dgtq(2,3) + (r3dg(1,1)-r1dg(3,1)) ! torque
1611	dgtq(3,3) = dgtq(3,3) + (r1dg(2,1)-r2dg(1,1))
1612
1613	dgtq(1,4) = 0! dgtq(1,4) + (r2dg(3,2)-r3dg(2,2)) ! viscous
1614	dgtq(2,4) = 0! dgtq(2,4) + (r3dg(1,2)-r1dg(3,2)) ! torque
1615	dgtq(3,4) = dgtq(3,4) + (r1dg(2,2)-r2dg(1,2))
1616	enddo
1617	enddo
1618	endif
1619
1620	return
1621	end
1622	c-----------------------------------------------------------------------
1623	subroutine torque_calc(scale,x0,ifdout,iftout)
1624	c
1625	c Compute torque about point x0
1626	c
1627	c Scale is a user-supplied multiplier so that results may be
1628	c scaled to any convenient non-dimensionalization.
1629	c
1630	c
1631	INCLUDE 'SIZE'
1632	INCLUDE 'TOTAL'
1633
1634	common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
1635	$ , scale_vf(3)
1636
1637	c
1638	real x0(3),w1(0:maxobj)
1639	logical ifdout,iftout
1640	c
1641	common /scrns/ sij (lx1ly1lz16lelv)
1642	common /scrcg/ pm1 (lx1,ly1,lz1,lelv)
1643	common /scrsf/ xm0(lx1,ly1,lz1,lelt)
1644	$, ym0(lx1,ly1,lz1,lelt)
1645	$, zm0(lx1,ly1,lz1,lelt)
1646	c
1647	parameter (lr=lx1ly1lz1)
1648	common /scruz/ ur(lr),us(lr),ut(lr)
1649	$ , vr(lr),vs(lr),vt(lr)
1650	$ , wr(lr),ws(lr),wt(lr)
1651	c
1652	n = lx1ly1lz1*nelv
1653	c
1654	call mappr(pm1,pr,xm0,ym0) ! map pressure onto Mesh 1
1655	c
1656	c Add mean_pressure_gradient.X to p:
1657
1658	if (param(55).ne.0) then
1659	dpdx_mean = -scale_vf(1)
1660	dpdy_mean = -scale_vf(2)
1661	dpdz_mean = -scale_vf(3)
1662	endif
1663
1664	call add2s2(pm1,xm1,dpdx_mean,n) ! Doesn't work if object is cut by
1665	call add2s2(pm1,ym1,dpdy_mean,n) ! periodicboundary. In this case,
1666	call add2s2(pm1,zm1,dpdz_mean,n) ! set ._mean=0 and compensate in
1667	c
1668	c Compute sij
1669	c
1670	nij = 3
1671	if (if3d.or.ifaxis) nij=6
1672	call comp_sij(sij,nij,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)
1673	c
1674	c
1675	c Fill up viscous array w/ default
1676	c
1677	if (istep.lt.1) call cfill(vdiff,param(2),n)
1678	c
1679	call cadd2(xm0,xm1,-x0(1),n)
1680	call cadd2(ym0,ym1,-x0(2),n)
1681	call cadd2(zm0,zm1,-x0(3),n)
1682	c
1683	x1min=glmin(xm0(1,1,1,1),n)
1684	x2min=glmin(ym0(1,1,1,1),n)
1685	x3min=glmin(zm0(1,1,1,1),n)
1686	c
1687	x1max=glmax(xm0(1,1,1,1),n)
1688	x2max=glmax(ym0(1,1,1,1),n)
1689	x3max=glmax(zm0(1,1,1,1),n)
1690	c
1691	do i=0,maxobj
1692	dragpx(i) = 0 ! BIG CODE :}
1693	dragvx(i) = 0
1694	dragx (i) = 0
1695	dragpy(i) = 0
1696	dragvy(i) = 0
1697	dragy (i) = 0
1698	dragpz(i) = 0
1699	dragvz(i) = 0
1700	dragz (i) = 0
1701	torqpx(i) = 0
1702	torqvx(i) = 0
1703	torqx (i) = 0
1704	torqpy(i) = 0
1705	torqvy(i) = 0
1706	torqy (i) = 0
1707	torqpz(i) = 0
1708	torqvz(i) = 0
1709	torqz (i) = 0
1710	enddo
1711	c
1712	c
1713	ifield = 1
1714	do iobj = 1,nobj
1715	memtot = nmember(iobj)
1716	do mem = 1,memtot
1717	ieg = object(iobj,mem,1)
1718	ifc = object(iobj,mem,2)
1719	if (gllnid(ieg).eq.nid) then ! this processor has a contribution
1720	ie = gllel(ieg)
1721	call drgtrq(dgtq,xm0,ym0,zm0,sij,pm1,vdiff,ifc,ie)
1722
1723	call cmult(dgtq,scale,12)
1724
1725	dragpx(iobj) = dragpx(iobj) + dgtq(1,1) ! pressure
1726	dragpy(iobj) = dragpy(iobj) + dgtq(2,1)
1727	dragpz(iobj) = dragpz(iobj) + dgtq(3,1)
1728
1729	dragvx(iobj) = dragvx(iobj) + dgtq(1,2) ! viscous
1730	dragvy(iobj) = dragvy(iobj) + dgtq(2,2)
1731	dragvz(iobj) = dragvz(iobj) + dgtq(3,2)
1732
1733	torqpx(iobj) = torqpx(iobj) + dgtq(1,3) ! pressure
1734	torqpy(iobj) = torqpy(iobj) + dgtq(2,3)
1735	torqpz(iobj) = torqpz(iobj) + dgtq(3,3)
1736
1737	torqvx(iobj) = torqvx(iobj) + dgtq(1,4) ! viscous
1738	torqvy(iobj) = torqvy(iobj) + dgtq(2,4)
1739	torqvz(iobj) = torqvz(iobj) + dgtq(3,4)
1740	endif
1741	enddo
1742	enddo
1743	c
1744	c Sum contributions from all processors
1745	c
1746	call gop(dragpx,w1,'+ ',maxobj+1)
1747	call gop(dragpy,w1,'+ ',maxobj+1)
1748	call gop(dragpz,w1,'+ ',maxobj+1)
1749	call gop(dragvx,w1,'+ ',maxobj+1)
1750	call gop(dragvy,w1,'+ ',maxobj+1)
1751	call gop(dragvz,w1,'+ ',maxobj+1)
1752	c
1753	call gop(torqpx,w1,'+ ',maxobj+1)
1754	call gop(torqpy,w1,'+ ',maxobj+1)
1755	call gop(torqpz,w1,'+ ',maxobj+1)
1756	call gop(torqvx,w1,'+ ',maxobj+1)
1757	call gop(torqvy,w1,'+ ',maxobj+1)
1758	call gop(torqvz,w1,'+ ',maxobj+1)
1759
1760	do i=1,nobj
1761	dragx(i) = dragpx(i) + dragvx(i)
1762	dragy(i) = dragpy(i) + dragvy(i)
1763	dragz(i) = dragpz(i) + dragvz(i)
1764
1765	torqx(i) = torqpx(i) + torqvx(i)
1766	torqy(i) = torqpy(i) + torqvy(i)
1767	torqz(i) = torqpz(i) + torqvz(i)
1768
1769	dragpx(0) = dragpx (0) + dragpx (i)
1770	dragvx(0) = dragvx (0) + dragvx (i)
1771	dragx (0) = dragx (0) + dragx (i)
1772
1773	dragpy(0) = dragpy (0) + dragpy (i)
1774	dragvy(0) = dragvy (0) + dragvy (i)
1775	dragy (0) = dragy (0) + dragy (i)
1776
1777	dragpz(0) = dragpz (0) + dragpz (i)
1778	dragvz(0) = dragvz (0) + dragvz (i)
1779	dragz (0) = dragz (0) + dragz (i)
1780
1781	torqpx(0) = torqpx (0) + torqpx (i)
1782	torqvx(0) = torqvx (0) + torqvx (i)
1783	torqx (0) = torqx (0) + torqx (i)
1784
1785	torqpy(0) = torqpy (0) + torqpy (i)
1786	torqvy(0) = torqvy (0) + torqvy (i)
1787	torqy (0) = torqy (0) + torqy (i)
1788
1789	torqpz(0) = torqpz (0) + torqpz (i)
1790	torqvz(0) = torqvz (0) + torqvz (i)
1791	torqz (0) = torqz (0) + torqz (i)
1792	enddo
1793
1794	do i=1,nobj
1795	if (nio.eq.0) then
1796	if (if3d.or.ifaxis) then
1797	if (ifdout) then
1798	write(6,6) istep,time,dragx(i),dragpx(i),dragvx(i),i,'dragx'
1799	write(6,6) istep,time,dragy(i),dragpy(i),dragvy(i),i,'dragy'
1800	write(6,6) istep,time,dragz(i),dragpz(i),dragvz(i),i,'dragz'
1801	endif
1802	if (iftout) then
1803	write(6,6) istep,time,torqx(i),torqpx(i),torqvx(i),i,'torqx'
1804	write(6,6) istep,time,torqy(i),torqpy(i),torqvy(i),i,'torqy'
1805	write(6,6) istep,time,torqz(i),torqpz(i),torqvz(i),i,'torqz'
1806	endif
1807	else
1808	if (ifdout) then
1809	write(6,6) istep,time,dragx(i),dragpx(i),dragvx(i),i,'dragx'
1810	write(6,6) istep,time,dragy(i),dragpy(i),dragvy(i),i,'dragy'
1811	endif
1812	if (iftout) then
1813	write(6,6) istep,time,torqz(i),torqpz(i),torqvz(i),i,'torqz'
1814	endif
1815	endif
1816	endif
1817	6 format(i8,1p4e19.11,2x,i5,a5)
1818	enddo
1819
1820	return
1821	end
1822	c-----------------------------------------------------------------------
1823	subroutine comp_sij(sij,nij,u,v,w,ur,us,ut,vr,vs,vt,wr,ws,wt)
1824	c du_i du_j
1825	c Compute the stress tensor S_ij := ---- + ----
1826	c du_j du_i
1827	c
1828	include 'SIZE'
1829	include 'TOTAL'
1830	c
1831	integer e
1832	c
1833	real sij(lx1ly1lz1,nij,lelv)
1834	real u (lx1ly1lz1,lelv)
1835	real v (lx1ly1lz1,lelv)
1836	real w (lx1ly1lz1,lelv)
1837	real ur (1) , us (1) , ut (1)
1838	$ , vr (1) , vs (1) , vt (1)
1839	$ , wr (1) , ws (1) , wt (1)
1840
1841	real j ! Inverse Jacobian
1842
1843	n = lx1-1 ! Polynomial degree
1844	nxyz = lx1ly1lz1
1845
1846	if (if3d) then ! 3D CASE
1847	do e=1,nelv
1848	call local_grad3(ur,us,ut,u,N,e,dxm1,dxtm1)
1849	call local_grad3(vr,vs,vt,v,N,e,dxm1,dxtm1)
1850	call local_grad3(wr,ws,wt,w,N,e,dxm1,dxtm1)
1851
1852	do i=1,nxyz
1853
1854	j = jacmi(i,e)
1855
1856	sij(i,1,e) = j* ! du/dx + du/dx
1857	$ 2(ur(i)rxm1(i,1,1,e)+us(i)sxm1(i,1,1,e)+ut(i)txm1(i,1,1,e))
1858
1859	sij(i,2,e) = j* ! dv/dy + dv/dy
1860	$ 2(vr(i)rym1(i,1,1,e)+vs(i)sym1(i,1,1,e)+vt(i)tym1(i,1,1,e))
1861
1862	sij(i,3,e) = j* ! dw/dz + dw/dz
1863	$ 2(wr(i)rzm1(i,1,1,e)+ws(i)szm1(i,1,1,e)+wt(i)tzm1(i,1,1,e))
1864
1865	sij(i,4,e) = j* ! du/dy + dv/dx
1866	$ (ur(i)rym1(i,1,1,e)+us(i)sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e) +
1867	$ vr(i)rxm1(i,1,1,e)+vs(i)sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e) )
1868
1869	sij(i,5,e) = j* ! dv/dz + dw/dy
1870	$ (wr(i)rym1(i,1,1,e)+ws(i)sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e) +
1871	$ vr(i)rzm1(i,1,1,e)+vs(i)szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e) )
1872
1873	sij(i,6,e) = j* ! du/dz + dw/dx
1874	$ (ur(i)rzm1(i,1,1,e)+us(i)szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e) +
1875	$ wr(i)rxm1(i,1,1,e)+ws(i)sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e) )
1876
1877	enddo
1878	enddo
1879
1880	elseif (ifaxis) then ! AXISYMMETRIC CASE
1881
1882	c
1883	c Notation: ( 2 x Acheson, p. 353)
1884	c Cylindrical
1885	c Nek5k Coordinates
1886	c
1887	c x z
1888	c y r
1889	c z theta
1890	c
1891
1892	do e=1,nelv
1893	call setaxdy ( ifrzer(e) ) ! change dytm1 if on-axis
1894	call local_grad2(ur,us,u,N,e,dxm1,dytm1)
1895	call local_grad2(vr,vs,v,N,e,dxm1,dytm1)
1896	call local_grad2(wr,ws,w,N,e,dxm1,dytm1)
1897
1898	do i=1,nxyz
1899	j = jacmi(i,e)
1900	r = ym1(i,1,1,e) ! Cyl. Coord:
1901
1902	sij(i,1,e) = j* ! du/dx + du/dx ! e_zz
1903	$ 2(ur(i)rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
1904
1905	sij(i,2,e) = j* ! dv/dy + dv/dy ! e_rr
1906	$ 2(vr(i)rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
1907
1908	if (r.gt.0) then ! e_@@
1909	sij(i,3,e) = v(i,e)/r ! v / r
1910	else
1911	sij(i,3,e) = j* ! L'Hopital's rule: e_@@ = dv/dr
1912	$ 2(vr(i)rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
1913	endif
1914
1915	sij(i,4,e) = j* ! du/dy + dv/dx ! e_zr
1916	$ ( ur(i)rym1(i,1,1,e)+us(i)sym1(i,1,1,e) +
1917	$ vr(i)rxm1(i,1,1,e)+vs(i)sxm1(i,1,1,e) )
1918
1919	if (r.gt.0) then ! e_r@
1920	sij(i,5,e) = j* ! dw/dy
1921	$ ( wr(i)rym1(i,1,1,e)+ws(i)sym1(i,1,1,e) )
1922	$ - w(i,e) / r
1923	else
1924	sij(i,5,e) = 0
1925	endif
1926
1927	sij(i,6,e) = j* ! dw/dx ! e_@z
1928	$ ( wr(i)rxm1(i,1,1,e)+ws(i)sxm1(i,1,1,e) )
1929
1930	enddo
1931	enddo
1932
1933	else ! 2D CASE
1934
1935	do e=1,nelv
1936	call local_grad2(ur,us,u,N,e,dxm1,dxtm1)
1937	call local_grad2(vr,vs,v,N,e,dxm1,dxtm1)
1938
1939	do i=1,nxyz
1940	j = jacmi(i,e)
1941
1942	sij(i,1,e) = j* ! du/dx + du/dx
1943	$ 2(ur(i)rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
1944
1945	sij(i,2,e) = j* ! dv/dy + dv/dy
1946	$ 2(vr(i)rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
1947
1948	sij(i,3,e) = j* ! du/dy + dv/dx
1949	$ (ur(i)rym1(i,1,1,e)+us(i)sym1(i,1,1,e) +
1950	$ vr(i)rxm1(i,1,1,e)+vs(i)sxm1(i,1,1,e) )
1951
1952	enddo
1953	enddo
1954	endif
1955	return
1956	end
1957	c-----------------------------------------------------------------------
1958	subroutine auto_averager(fname127) ! simple average of files
1959
1960	c This routine reads files specificed of file.list and averages
1961	c them with uniform weight
1962	c
1963	c Note that it relies on scravg and scruz common blocks. pff 12/7/14
1964	c
1965
1966	include 'SIZE'
1967	include 'TOTAL'
1968
1969	character*127 fname127
1970	character*1 f1(127)
1971
1972	parameter (lt=lx1ly1lz1*lelt)
1973	common /scravg/ ua(lt),va(lt),wa(lt),pa(lt)
1974	common /scrsf/ ta(lt,ldimt)
1975
1976	character*1 s1(127)
1977	equivalence (s1,initc) ! equivalence to initial condition
1978
1979	if (nid.eq.0) then
1980	ib=indx1(fname127,' ',1)-1
1981	call chcopy(f1,fname127,ib)
1982	write(6,2) (f1(k),k=1,ib)
1983	2 format('Open file: ',127a1)
1984	endif
1985
1986	ierr = 0
1987	if (nid.eq.0) open(77,file=fname127,status='old',err=199)
1988	ierr = iglmax(ierr,1)
1989	if (ierr.gt.0) goto 199
1990	n = lx1ly1lz1*nelt
1991	n2= lx2ly2lz2*nelt
1992
1993	call rzero (ua,n)
1994	call rzero (va,n)
1995	call rzero (wa,n)
1996	call rzero (pa,n2)
1997	do k=1,npscal+1
1998	call rzero (ta(1,k),n)
1999	enddo
2000
2001	icount = 0
2002	do ipass=1,9999999
2003
2004	call blank(initc,127)
2005	initc(1) = 'done '
2006	if (nid.eq.0) read(77,127,end=998) initc(1)
2007	998 call bcast(initc,127)
2008	127 format(a127)
2009
2010	iblank = indx1(initc,' ',1)-1
2011	if (nio.eq.0) write(6,1) ipass,(s1(k),k=1,iblank)
2012	1 format(i8,'Reading: ',127a1)
2013
2014	if (indx1(initc,'done ',5).eq.0) then ! We're not done
2015
2016	nfiles = 1
2017	call restart(nfiles) ! Note -- time is reset.
2018
2019	call opadd2 (ua,va,wa,vx,vy,vz)
2020	call add2 (pa,pr,n2)
2021	do k=1,npscal+1
2022	call add2(ta(1,k),t(1,1,1,1,k),n)
2023	enddo
2024	icount = icount+1
2025
2026	else
2027	goto 999
2028	endif
2029
2030	enddo
2031
2032	999 continue ! clean up averages
2033	if (nid.eq.0) close(77)
2034
2035	scale = 1./icount
2036	call cmult2(vx,ua,scale,n)
2037	call cmult2(vy,va,scale,n)
2038	call cmult2(vz,wa,scale,n)
2039	call cmult2(pr,pa,scale,n2)
2040	do k=1,npscal+1
2041	call cmult2(t(1,1,1,1,k),ta(1,k),scale,n)
2042	enddo
2043	return
2044
2045	199 continue ! exception handle for file not found
2046	ierr = 1
2047	if (nid.eq.0) ierr = iglmax(ierr,1)
2048	call exitti('Auto averager did not find list file.$',ierr)
2049
2050	return
2051	end
2052	c-----------------------------------------------------------------------
2053	subroutine outmesh
2054	include 'SIZE'
2055	include 'TOTAL'
2056	integer e,eg
2057
2058	common /cmesh/ xt(2**ldim,ldim)
2059
2060	len = wdsizeldim(2**ldim)
2061
2062	if (nid.eq.0) open(unit=29,file='rea.new')
2063
2064	do eg=1,nelgt
2065	call nekgsync() ! belt
2066	jnid = gllnid(eg)
2067	e = gllel (eg)
2068	mtype = e
2069	if (jnid.eq.0 .and. nid.eq.0) then
2070	call get_el(xt,xm1(1,1,1,e),ym1(1,1,1,e),zm1(1,1,1,e))
2071	call out_el(xt,eg)
2072	elseif (nid.eq.0) then
2073	call crecv2(mtype,xt,len,jnid)
2074	call out_el(xt,eg)
2075	elseif (jnid.eq.nid) then
2076	call get_el(xt,xm1(1,1,1,e),ym1(1,1,1,e),zm1(1,1,1,e))
2077	call csend(mtype,xt,len,0,0)
2078	endif
2079	call nekgsync() ! suspenders
2080	enddo
2081
2082	if (nid.eq.0) close(29)
2083	call nekgsync()
2084	call exitt
2085
2086	return
2087	end
2088	c-----------------------------------------------------------------------
2089	subroutine out_el(xt,e)
2090	include 'SIZE'
2091	include 'TOTAL'
2092
2093	real xt(2**ldim,ldim)
2094	integer e
2095
2096	integer ed(8)
2097	save ed
2098	data ed / 1,2,4,3 , 5,6,8,7 /
2099
2100	write(29,1) e
2101	write(29,2) ((xt(ed(k),j),k=1,4),j=1,ldim)
2102	write(29,2) ((xt(ed(k),j),k=5,8),j=1,ldim)
2103
2104	1 format(12x,'ELEMENT',i6,' [ 1 ] GROUP 0')
2105	2 format(1p4e18.10)
2106
2107	return
2108	end
2109	c-----------------------------------------------------------------------
2110	subroutine get_el(xt,x,y,z)
2111	include 'SIZE'
2112	include 'TOTAL'
2113
2114	real xt(2**ldim,ldim)
2115	real x(lx1,ly1,lz1),y(lx1,ly1,lz1),z(lx1,ly1,lz1)
2116
2117	l = 0
2118	do k=1,lz1,max(lz1-1,1)
2119	do j=1,ly1,ly1-1
2120	do i=1,lx1,lx1-1
2121	l = l+1
2122	xt(l,1) = x(i,j,k)
2123	xt(l,2) = y(i,j,k)
2124	xt(l,3) = z(i,j,k)
2125	enddo
2126	enddo
2127	enddo
2128
2129	return
2130	end
2131	c-----------------------------------------------------------------------
2132	subroutine shear_calc_max(strsmx,scale,x0,ifdout,iftout)
2133	c
2134	c Compute maximum shear stress on objects
2135	c
2136	c Scale is a user-supplied multiplier so that results may be
2137	c scaled to any convenient non-dimensionalization.
2138	c
2139	c
2140	INCLUDE 'SIZE'
2141	INCLUDE 'TOTAL'
2142
2143	real strsmx(maxobj),x0(3),w1(0:maxobj)
2144	logical ifdout,iftout
2145
2146	common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
2147	$ , scale_vf(3)
2148
2149
2150	common /scrns/ sij (lx1ly1lz16lelv)
2151	common /scrcg/ pm1 (lx1,ly1,lz1,lelv)
2152	common /scrsf/ xm0(lx1,ly1,lz1,lelt)
2153	$, ym0(lx1,ly1,lz1,lelt)
2154	$, zm0(lx1,ly1,lz1,lelt)
2155
2156	parameter (lr=lx1ly1lz1)
2157	common /scruz/ ur(lr),us(lr),ut(lr)
2158	$ , vr(lr),vs(lr),vt(lr)
2159	$ , wr(lr),ws(lr),wt(lr)
2160
2161
2162	n = lx1ly1lz1*nelv
2163
2164	call mappr(pm1,pr,xm0,ym0) ! map pressure onto Mesh 1
2165
2166	c Add mean_pressure_gradient.X to p:
2167
2168	if (param(55).ne.0) then
2169	dpdx_mean = -scale_vf(1)
2170	dpdy_mean = -scale_vf(2)
2171	dpdz_mean = -scale_vf(3)
2172	endif
2173
2174	call add2s2(pm1,xm1,dpdx_mean,n) ! Doesn't work if object is cut by
2175	call add2s2(pm1,ym1,dpdy_mean,n) ! periodicboundary. In this case,
2176	call add2s2(pm1,zm1,dpdz_mean,n) ! set ._mean=0 and compensate in
2177	c
2178	c Compute sij
2179	c
2180	nij = 3
2181	if (if3d.or.ifaxis) nij=6
2182	call comp_sij(sij,nij,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)
2183	c
2184	c
2185	c Fill up viscous array w/ default
2186	c
2187	if (istep.lt.1) call cfill(vdiff,param(2),n)
2188	c
2189	call cadd2(xm0,xm1,-x0(1),n)
2190	call cadd2(ym0,ym1,-x0(2),n)
2191	call cadd2(zm0,zm1,-x0(3),n)
2192	c
2193	x1min=glmin(xm0(1,1,1,1),n)
2194	x2min=glmin(ym0(1,1,1,1),n)
2195	x3min=glmin(zm0(1,1,1,1),n)
2196	c
2197	x1max=glmax(xm0(1,1,1,1),n)
2198	x2max=glmax(ym0(1,1,1,1),n)
2199	x3max=glmax(zm0(1,1,1,1),n)
2200	c
2201	call rzero(strsmx,maxobj)
2202
2203	ifield = 1
2204
2205	strsmx(ii) = 0.
2206	do iobj = 1,nobj
2207	memtot = nmember(iobj)
2208	do mem = 1,memtot
2209	ieg = object(iobj,mem,1)
2210	ifc = object(iobj,mem,2)
2211	if (gllnid(ieg).eq.nid) then ! this processor has a contribution
2212	ie = gllel(ieg)
2213	call get_strsmax(strsmxl,xm0,ym0,zm0,sij,pm1,vdiff,ifc,ie)
2214	call cmult(strsmxl,scale,1)
2215	strsmx(ii)=max(strsmx(ii),strsmxl)
2216	endif
2217	enddo
2218	enddo
2219	c
2220	c Max contributions over all processors
2221	c
2222	call gop(strsmx,w1,'M ',maxobj)
2223
2224	return
2225	end
2226	c-----------------------------------------------------------------------
2227	subroutine get_strsmax(strsmax,xm0,ym0,zm0,sij,pm1,visc,f,e)
2228	c
2229	INCLUDE 'SIZE'
2230	INCLUDE 'GEOM'
2231	INCLUDE 'INPUT'
2232	INCLUDE 'TOPOL'
2233	INCLUDE 'TSTEP'
2234	c
2235	real dgtq(3,4)
2236	real xm0 (lx1,ly1,lz1,lelt)
2237	real ym0 (lx1,ly1,lz1,lelt)
2238	real zm0 (lx1,ly1,lz1,lelt)
2239	real sij (lx1,ly1,lz1,3*ldim-3,lelv)
2240	real pm1 (lx1,ly1,lz1,lelv)
2241	real visc(lx1,ly1,lz1,lelv)
2242
2243	integer f,e,pf
2244	real n1,n2,n3
2245
2246	call dsset(lx1,ly1,lz1) ! set up counters
2247	pf = eface1(f) ! convert from preproc. notation
2248	js1 = skpdat(1,pf)
2249	jf1 = skpdat(2,pf)
2250	jskip1 = skpdat(3,pf)
2251	js2 = skpdat(4,pf)
2252	jf2 = skpdat(5,pf)
2253	jskip2 = skpdat(6,pf)
2254
2255	if (if3d.or.ifaxis) then
2256	i = 0
2257	strsmax = 0
2258	do j2=js2,jf2,jskip2
2259	do j1=js1,jf1,jskip1
2260	i = i+1
2261	n1 = unx(i,1,f,e)
2262	n2 = uny(i,1,f,e)
2263	n3 = unz(i,1,f,e)
2264	c
2265	v = visc(j1,j2,1,e)
2266	c
2267	s11 = sij(j1,j2,1,1,e)
2268	s21 = sij(j1,j2,1,4,e)
2269	s31 = sij(j1,j2,1,6,e)
2270	c
2271	s12 = sij(j1,j2,1,4,e)
2272	s22 = sij(j1,j2,1,2,e)
2273	s32 = sij(j1,j2,1,5,e)
2274
2275	s13 = sij(j1,j2,1,6,e)
2276	s23 = sij(j1,j2,1,5,e)
2277	s33 = sij(j1,j2,1,3,e)
2278
2279	stress1 = -v(s11n1 + s12n2 + s13n3) ! viscous drag
2280	stress2 = -v(s21n1 + s22n2 + s23n3)
2281	stress3 = -v(s31n1 + s32n2 + s33n3)
2282
2283	strsnrm = stress1stress1+stress2stress2+stress3*stress3
2284	strsmax = max(strsmax,strsnrm)
2285
2286	enddo
2287	enddo
2288
2289	else ! 2D
2290
2291	i = 0
2292	strsmax = 0
2293	do j2=js2,jf2,jskip2
2294	do j1=js1,jf1,jskip1
2295	i = i+1
2296	n1 = unx(i,1,f,e)*area(i,1,f,e)
2297	n2 = uny(i,1,f,e)*area(i,1,f,e)
2298	v = visc(j1,j2,1,e)
2299
2300	s11 = sij(j1,j2,1,1,e)
2301	s12 = sij(j1,j2,1,3,e)
2302	s21 = sij(j1,j2,1,3,e)
2303	s22 = sij(j1,j2,1,2,e)
2304
2305	stress1 = -v(s11n1 + s12*n2) ! viscous drag
2306	stress2 = -v(s21n1 + s22*n2)
2307
2308	strsnrm = stress1stress1+stress2stress2
2309	strsmax = max(strsmax,strsnrm)
2310
2311	enddo
2312	enddo
2313
2314	endif
2315
2316	if (strsmax.gt.0) strsmax = sqrt(strsmax)
2317
2318	return
2319	end
2320	c-----------------------------------------------------------------------
2321	subroutine fix_geom ! fix up geometry irregularities
2322
2323	include 'SIZE'
2324	include 'TOTAL'
2325
2326	parameter (lt = lx1ly1lz1)
2327	common /scrns/ xb(lt,lelt),yb(lt,lelt),zb(lt,lelt)
2328	common /scruz/ tmsk(lt,lelt),tmlt(lt,lelt),w1(lt),w2(lt)
2329
2330	integer e,f
2331	character*3 cb
2332
2333	n = lx1ly1lz1*nelt
2334	nxyz = lx1ly1lz1
2335	nfaces = 2*ldim
2336	ifield = 1 ! velocity field
2337	if (nelgv.ne.nelgt .or. .not.ifflow) ifield = 2 ! temperature field
2338
2339	call rone (tmlt,n)
2340	call dssum (tmlt,lx1,ly1,lz1) ! denominator
2341
2342	call rone (tmsk,n)
2343	do e=1,nelt ! fill mask where bc is periodic
2344	do f=1,nfaces ! so we don't translate periodic bcs (z only)
2345	cb =cbc(f,e,ifield)
2346	if (cb.eq.'P ') call facev (tmsk,e,f,0.0,lx1,ly1,lz1)
2347	enddo
2348	enddo
2349	call dsop(tmsk,'* ',lx1,ly1,lz1)
2350	call dsop(tmsk,'* ',lx1,ly1,lz1)
2351	call dsop(tmsk,'* ',lx1,ly1,lz1)
2352
2353	do kpass = 1,ldim+1 ! This doesn't work for 2D, yet.
2354	! Extra pass is just to test convergence
2355
2356	c call opcopy (xb,yb,zb,xm1,ym1,zm1) ! Must use WHOLE field,
2357	c call opdssum(xb,yb,zb) ! not just fluid domain.
2358	call copy (xb,xm1,n)
2359	call copy (yb,ym1,n)
2360	call copy (zb,zm1,n)
2361	call dssum (xb,lx1,ly1,lz1)
2362	call dssum (yb,lx1,ly1,lz1)
2363	call dssum (zb,lx1,ly1,lz1)
2364
2365	xm = 0.
2366	ym = 0.
2367	zm = 0.
2368
2369	do e=1,nelt
2370	do i=1,nxyz ! compute averages of geometry
2371	s = 1./tmlt(i,e)
2372	xb(i,e) = s*xb(i,e)
2373	yb(i,e) = s*yb(i,e)
2374	zb(i,e) = s*zb(i,e)
2375
2376	xb(i,e) = xb(i,e) - xm1(i,1,1,e) ! local displacements
2377	yb(i,e) = yb(i,e) - ym1(i,1,1,e)
2378	zb(i,e) = zb(i,e) - zm1(i,1,1,e)
2379	xb(i,e) = xb(i,e)*tmsk(i,e)
2380	yb(i,e) = yb(i,e)*tmsk(i,e)
2381	zb(i,e) = zb(i,e)*tmsk(i,e)
2382
2383	xm = max(xm,abs(xb(i,e)))
2384	ym = max(ym,abs(yb(i,e)))
2385	zm = max(zm,abs(zb(i,e)))
2386	enddo
2387
2388	if (kpass.le.ldim) then
2389	call gh_face_extend(xb(1,e),zgm1,lx1,kpass,w1,w2)
2390	call gh_face_extend(yb(1,e),zgm1,lx1,kpass,w1,w2)
2391	call gh_face_extend(zb(1,e),zgm1,lx1,kpass,w1,w2)
2392	endif
2393
2394	enddo
2395
2396	if (kpass.le.ldim) then
2397	call add2(xm1,xb,n)
2398	call add2(ym1,yb,n)
2399	call add2(zm1,zb,n)
2400	endif
2401
2402	xx = glamax(xb,n)
2403	yx = glamax(yb,n)
2404	zx = glamax(zb,n)
2405
2406	xm = glmax(xm,1)
2407	ym = glmax(ym,1)
2408	zm = glmax(zm,1)
2409
2410	if (nio.eq.0) write(6,1) xm,ym,zm,xx,yx,zx,kpass
2411	1 format(1p6e12.4,' xyz repair',i2)
2412
2413	enddo
2414
2415	param(59) = 1. ! ifdef = .true.
2416	call geom_reset(1) ! reset metrics, etc.
2417
2418	return
2419	end
2420	c-----------------------------------------------------------------------
2421	subroutine gh_face_extend(x,zg,n,gh_type,e,v)
2422	include 'SIZE'
2423
2424	real x(1),zg(1),e(1),v(1)
2425	integer gh_type
2426
2427	if (ldim.eq.2) then
2428	call gh_face_extend_2d(x,zg,n,gh_type,e,v)
2429	else
2430	call gh_face_extend_3d(x,zg,n,gh_type,e,v)
2431	endif
2432
2433	return
2434	end
2435	c-----------------------------------------------------------------------
2436	subroutine gh_face_extend_2d(x,zg,n,gh_type,e,v)
2437	c
2438	c Extend 2D faces into interior via gordon hall
2439	c
2440	c gh_type: 1 - vertex only
2441	c 2 - vertex and faces
2442	c
2443	c
2444	real x(n,n)
2445	real zg(n)
2446	real e(n,n)
2447	real v(n,n)
2448	integer gh_type
2449	c
2450	c Build vertex interpolant
2451	c
2452	ntot=n*n
2453	call rzero(v,ntot)
2454	do jj=1,n,n-1
2455	do ii=1,n,n-1
2456	do j=1,n
2457	do i=1,n
2458	si = 0.5((n-ii)(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
2459	sj = 0.5((n-jj)(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
2460	v(i,j) = v(i,j) + sisjx(ii,jj)
2461	enddo
2462	enddo
2463	enddo
2464	enddo
2465	if (gh_type.eq.1) then
2466	call copy(x,v,ntot)
2467	return
2468	endif
2469
2470
2471	c Extend 4 edges
2472	call rzero(e,ntot)
2473	c
2474	c x-edges
2475	c
2476	do jj=1,n,n-1
2477	do j=1,n
2478	do i=1,n
2479	hj = 0.5((n-jj)(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
2480	e(i,j) = e(i,j) + hj*(x(i,jj)-v(i,jj))
2481	enddo
2482	enddo
2483	enddo
2484	c
2485	c y-edges
2486	c
2487	do ii=1,n,n-1
2488	do j=1,n
2489	do i=1,n
2490	hi = 0.5((n-ii)(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
2491	e(i,j) = e(i,j) + hi*(x(ii,j)-v(ii,j))
2492	enddo
2493	enddo
2494	enddo
2495
2496	call add3(x,e,v,ntot)
2497
2498	return
2499	end
2500	c-----------------------------------------------------------------------
2501	subroutine gh_face_extend_3d(x,zg,n,gh_type,e,v)
2502	c
2503	c Extend faces into interior via gordon hall
2504	c
2505	c gh_type: 1 - vertex only
2506	c 2 - vertex and edges
2507	c 3 - vertex, edges, and faces
2508	c
2509	c
2510	real x(n,n,n)
2511	real zg(n)
2512	real e(n,n,n)
2513	real v(n,n,n)
2514	integer gh_type
2515	c
2516	c Build vertex interpolant
2517	c
2518	ntot=nnn
2519	call rzero(v,ntot)
2520	do kk=1,n,n-1
2521	do jj=1,n,n-1
2522	do ii=1,n,n-1
2523	do k=1,n
2524	do j=1,n
2525	do i=1,n
2526	si = 0.5((n-ii)(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
2527	sj = 0.5((n-jj)(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
2528	sk = 0.5((n-kk)(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
2529	v(i,j,k) = v(i,j,k) + sisjsk*x(ii,jj,kk)
2530	enddo
2531	enddo
2532	enddo
2533	enddo
2534	enddo
2535	enddo
2536	if (gh_type.eq.1) then
2537	call copy(x,v,ntot)
2538	return
2539	endif
2540	c
2541	c
2542	c Extend 12 edges
2543	call rzero(e,ntot)
2544	c
2545	c x-edges
2546	c
2547	do kk=1,n,n-1
2548	do jj=1,n,n-1
2549	do k=1,n
2550	do j=1,n
2551	do i=1,n
2552	hj = 0.5((n-jj)(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
2553	hk = 0.5((n-kk)(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
2554	e(i,j,k) = e(i,j,k) + hjhk(x(i,jj,kk)-v(i,jj,kk))
2555	enddo
2556	enddo
2557	enddo
2558	enddo
2559	enddo
2560	c
2561	c y-edges
2562	c
2563	do kk=1,n,n-1
2564	do ii=1,n,n-1
2565	do k=1,n
2566	do j=1,n
2567	do i=1,n
2568	hi = 0.5((n-ii)(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
2569	hk = 0.5((n-kk)(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
2570	e(i,j,k) = e(i,j,k) + hihk(x(ii,j,kk)-v(ii,j,kk))
2571	enddo
2572	enddo
2573	enddo
2574	enddo
2575	enddo
2576	c
2577	c z-edges
2578	c
2579	do jj=1,n,n-1
2580	do ii=1,n,n-1
2581	do k=1,n
2582	do j=1,n
2583	do i=1,n
2584	hi = 0.5((n-ii)(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
2585	hj = 0.5((n-jj)(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
2586	e(i,j,k) = e(i,j,k) + hihj(x(ii,jj,k)-v(ii,jj,k))
2587	enddo
2588	enddo
2589	enddo
2590	enddo
2591	enddo
2592	c
2593	call add2(e,v,ntot)
2594	c
2595	if (gh_type.eq.2) then
2596	call copy(x,e,ntot)
2597	return
2598	endif
2599	c
2600	c Extend faces
2601	c
2602	call rzero(v,ntot)
2603	c
2604	c x-edges
2605	c
2606	do ii=1,n,n-1
2607	do k=1,n
2608	do j=1,n
2609	do i=1,n
2610	hi = 0.5((n-ii)(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
2611	v(i,j,k) = v(i,j,k) + hi*(x(ii,j,k)-e(ii,j,k))
2612	enddo
2613	enddo
2614	enddo
2615	enddo
2616	c
2617	c y-edges
2618	c
2619	do jj=1,n,n-1
2620	do k=1,n
2621	do j=1,n
2622	do i=1,n
2623	hj = 0.5((n-jj)(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
2624	v(i,j,k) = v(i,j,k) + hj*(x(i,jj,k)-e(i,jj,k))
2625	enddo
2626	enddo
2627	enddo
2628	enddo
2629	c
2630	c z-edges
2631	c
2632	do kk=1,n,n-1
2633	do k=1,n
2634	do j=1,n
2635	do i=1,n
2636	hk = 0.5((n-kk)(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
2637	v(i,j,k) = v(i,j,k) + hk*(x(i,j,kk)-e(i,j,kk))
2638	enddo
2639	enddo
2640	enddo
2641	enddo
2642	c
2643	call add2(v,e,ntot)
2644	call copy(x,v,ntot)
2645
2646	return
2647	end
2648	c-----------------------------------------------------------------------
2649	function ran1(idum)
2650	c
2651	integer idum,ia,im,iq,ir,ntab,ndiv
2652	real ran1,am,eps,rnmx
2653	c
2654	parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836)
2655	parameter (ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
2656	c
2657	c Numerical Rec. in Fortran, 2nd eD. P. 271
2658	c
2659	integer j,k
2660	integer iv(ntab),iy
2661	save iv,iy
2662	data iv,iy /ntab*0,0/
2663	c
2664	if (idum.le.0.or.iy.eq.0) then
2665	idum=max(-idum,1)
2666	do j=ntab+8,1,-1
2667	k = idum/iq
2668	idum = ia(idum-kiq)-ir*k
2669	if(idum.lt.0) idum = idum+im
2670	if (j.le.ntab) iv(j) = idum
2671	enddo
2672	iy = iv(1)
2673	endif
2674	k = idum/iq
2675	idum = ia(idum-kiq)-ir*k
2676	if(idum.lt.0) idum = idum+im
2677	j = 1+iy/ndiv
2678	iy = iv(j)
2679	iv(j) = idum
2680	ran1 = min(am*iy,rnmx)
2681	c ran1 = cos(ran1*1.e8)
2682
2683	return
2684	end
2685	c-----------------------------------------------------------------------
2686	subroutine rand_fld_h1(x)
2687
2688	include 'SIZE'
2689	real x(1)
2690
2691	n=lx1ly1lz1*nelt
2692	id = n
2693	do i=1,n
2694	x(i) = ran1(id)
2695	enddo
2696	call dsavg(x)
2697
2698	return
2699	end
2700	c-----------------------------------------------------------------------
2701	subroutine rescale_x (x,x0,x1)
2702	include 'SIZE'
2703	real x(1)
2704
2705	n = lx1ly1lz1*nelt
2706	xmin = glmin(x,n)
2707	xmax = glmax(x,n)
2708
2709	if (xmax.le.xmin) return
2710
2711	scale = (x1-x0)/(xmax-xmin)
2712	do i=1,n
2713	x(i) = x0 + scale*(x(i)-xmin)
2714	enddo
2715
2716	return
2717	end
2718	c-----------------------------------------------------------------------
2719	subroutine build_filter(f,diag,nx)
2720	include 'SIZE'
2721
2722	real f(nx,nx),diag(nx),zpts(nx)
2723
2724	parameter (lm=4*lx1) ! Totally arbitrary
2725	parameter (lm2=lm*lm)
2726
2727	common /cfilt1/ phi,pht,ft,rmult,Lj,gpts,indr,indc,ipiv
2728	real phi(lm2),pht(lm2),ft(lm2),rmult(lm),Lj(lm),gpts(lm)
2729	integer indr(lm),indc(lm),ipiv(lm)
2730
2731	integer nxl
2732	save nxl
2733	data nxl / -9 /
2734
2735	if (nx.gt.lm) call exitti('ABORT in build_filter:$',nx)
2736
2737	if (nx.ne.nxl) then
2738
2739	nxl = nx
2740
2741	call zwgll (gpts,f,nx) ! Get nx GLL points
2742
2743	kj = 0
2744	n = nx-1
2745	do j=1,nx
2746	z = gpts(j)
2747	call legendre_poly(Lj,z,n)
2748	kj = kj+1
2749	pht(kj) = Lj(1)
2750	kj = kj+1
2751	pht(kj) = Lj(2)
2752	do k=3,nx
2753	kj = kj+1
2754	pht(kj) = Lj(k)-Lj(k-2)
2755	enddo
2756	enddo
2757
2758	call transpose (phi,nx,pht,nx)
2759	call copy (pht,phi,nx*nx)
2760	call gaujordf (pht,nx,nx,indr,indc,ipiv,ierr,rmult)
2761
2762	endif ! End of save section
2763
2764	ij=0
2765	do j=1,nx
2766	do i=1,nx
2767	ij = ij+1
2768	ft(ij) = diag(i)*pht(ij)
2769	enddo
2770	enddo
2771	! -1
2772	call mxm (phi,nx,ft,nx,f,nx) ! V D V
2773
2774	return
2775	end
2776	c-----------------------------------------------------------------------
2777	subroutine g_filter(u,diag,ifld)
2778	c
2779	c Generalized filter: F(u) with F = J^T D J, where D=diag(diag)
2780	c
2781	include 'SIZE'
2782	include 'TOTAL'
2783
2784	real u(1),diag(1)
2785
2786	parameter (lxx=lx1lx1,lxyz=lx1ly1*lz1)
2787	common /ctmp0/ f(lxx),wk1(lxyz),wk2(lxyz),wk3(lxyz)
2788
2789	ifldt = ifield
2790	ifield = ifld
2791
2792	call build_filter(f,diag,lx1)
2793	call filterq(u,f,lx1,lz1,wk1,wk2,wk3,if3d,umax)
2794
2795	ifield = ifldt
2796
2797	return
2798	end
2799	c-----------------------------------------------------------------------
2800	subroutine cut_off_filter(u,mx,ifld) ! mx=max saved mode
2801	c
2802	c Generalized filter: F(u) with F = J^T D J, where D=diag(diag)
2803	c
2804	include 'SIZE'
2805	include 'TOTAL'
2806
2807	real u(1)
2808
2809	parameter (lxx=lx1lx1,lxyz=lx1ly1*lz1)
2810	common /ctmp0/ f(lxx),wk1(lxyz),wk2(lxyz),wk3(lxyz),diag(lx1)
2811
2812	ifldt = ifield
2813	ifield = ifld
2814
2815	call rone(diag,lx1)
2816	do i=mx+1,lx1
2817	diag(i)=0.
2818	enddo
2819
2820	call build_filter(f,diag,lx1)
2821	call filterq(u,f,lx1,lz1,wk1,wk2,wk3,if3d,umax)
2822
2823	ifield = ifldt
2824
2825	return
2826	end
2827	c-----------------------------------------------------------------------
2828	subroutine filter_d2(v,nx,nz,wgt,ifd4)
2829
2830	include 'SIZE'
2831	include 'TOTAL'
2832
2833	parameter (lt=lx1ly1lz1)
2834	real v(lt,nelt)
2835	logical ifd4
2836
2837	common /ctmp1/ w(lt,lelt),ur(lt),us(lt),ut(lt),w1(2*lt)
2838
2839	integer e
2840
2841	n = lx1ly1lz1*nelt
2842	nn = nx-1
2843	nel = nelfld(ifield)
2844
2845	bmax = glamax(v,n)
2846
2847	if (if3d) then
2848	do e=1,nel
2849	call local_grad3(ur,us,ut,v(1,e),nn,1,dxm1,dxtm1)
2850	do i=1,lt
2851	ur(i) = ur(i)*w3m1(i,1,1)
2852	us(i) = us(i)*w3m1(i,1,1)
2853	ut(i) = ut(i)*w3m1(i,1,1)
2854	enddo
2855	call local_grad3_t(w(1,e),ur,us,ut,nn,1,dxm1,dxtm1,w1)
2856	enddo
2857	call dsavg(w) !NOTE STILL NEED BC TREATMENT !
2858
2859	if (ifd4) then
2860	wght = 20./(lx1**4)
2861	do e=1,nel
2862	do i=1,lt
2863	w(i,e) = wght*w(i,e)/w3m1(i,1,1)
2864	enddo
2865	call local_grad3(ur,us,ut,w(1,e),nn,1,dxm1,dxtm1)
2866	do i=1,lt
2867	ur(i) = ur(i)*w3m1(i,1,1)
2868	us(i) = us(i)*w3m1(i,1,1)
2869	ut(i) = ut(i)*w3m1(i,1,1)
2870	enddo
2871	call local_grad3_t(w(1,e),ur,us,ut,nn,1,dxm1,dxtm1,w1)
2872	enddo
2873	call dsavg(w) !NOTE STILL NEED BC TREATMENT !
2874	endif
2875
2876	wght = wgt/(lx1**4)
2877	do e=1,nel
2878	do i=1,lt
2879	v(i,e) = v(i,e) - wght*w(i,e)/w3m1(i,1,1)
2880	enddo
2881	enddo
2882
2883	else ! 2D
2884
2885	do e=1,nel
2886	call local_grad2(ur,us,v(1,e),nn,1,dxm1,dxtm1)
2887	do i=1,lt
2888	ur(i) = ur(i)*w3m1(i,1,1)
2889	us(i) = us(i)*w3m1(i,1,1)
2890	enddo
2891	call local_grad2_t(w(1,e),ur,us,nn,1,dxm1,dxtm1,w1)
2892	enddo
2893	call dsavg(w) !NOTE STILL NEED BC TREATMENT !
2894
2895	if (ifd4) then
2896	wght = 200./(lx1**4)
2897	do e=1,nel
2898	do i=1,lt
2899	w(i,e) = wght*w(i,e)/w3m1(i,1,1)
2900	enddo
2901	call local_grad2(ur,us,w(1,e),nn,1,dxm1,dxtm1)
2902	do i=1,lt
2903	ur(i) = ur(i)*w3m1(i,1,1)
2904	us(i) = us(i)*w3m1(i,1,1)
2905	enddo
2906	call local_grad2_t(w(1,e),ur,us,nn,1,dxm1,dxtm1,w1)
2907	enddo
2908	call dsavg(w) !NOTE STILL NEED BC TREATMENT !
2909	endif
2910
2911	wght = wgt/(lx1**4)
2912	do e=1,nel
2913	do i=1,lt
2914	v(i,e) = v(i,e) - wght*w(i,e)/w3m1(i,1,1)
2915	enddo
2916	enddo
2917
2918	endif
2919
2920	vmax = glamax(v,n)
2921	if (nio.eq.0) write(6,1) istep,time,vmax,bmax,' filter max'
2922	1 format(i9,1p3e12.4,a11)
2923
2924	return
2925	end
2926	c-------------------------------------------------------------------------
2927	function dist3d(a,b,c,x,y,z)
2928
2929	d = (a-x)2 + (b-y)2 + (c-z)**2
2930
2931	dist3d = 0.
2932	if (d.gt.0) dist3d = sqrt(d)
2933
2934	return
2935	end
2936	c-----------------------------------------------------------------------
2937	function dist2d(a,b,x,y)
2938
2939	d = (a-x)2 + (b-y)2
2940
2941	dist2d = 0.
2942	if (d.gt.0) dist2d = sqrt(d)
2943
2944	return
2945	end
2946	c-----------------------------------------------------------------------
2947	subroutine domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
2948
2949	include 'SIZE'
2950	include 'TOTAL'
2951
2952	n = lx1ly1lz1*nelt
2953
2954	xmin = glmin(xm1,n)
2955	xmax = glmax(xm1,n)
2956	ymin = glmin(ym1,n)
2957	ymax = glmax(ym1,n)
2958	if (if3d) then
2959	zmin = glmin(zm1,n)
2960	zmax = glmax(zm1,n)
2961	else
2962	zmin = 0.
2963	zmax = 0.
2964	endif
2965
2966	return
2967	end
2968	c-----------------------------------------------------------------------
2969	subroutine cheap_dist(d,ifld,b)
2970
2971	c Finds a pseudo-distance function.
2972
2973	c INPUT: ifld - field type for which distance function is to be found.
2974	c ifld = 1 for velocity
2975	c ifld = 2 for temperature, etc.
2976
2977	c OUTPUT: d = "path" distance to nearest wall
2978
2979	c This approach has a significant advantage that it works for
2980	c periodict boundary conditions, whereas most other approaches
2981	c will not.
2982
2983	include 'SIZE'
2984	include 'GEOM' ! Coordinates
2985	include 'INPUT' ! cbc()
2986	include 'TSTEP' ! nelfld
2987	include 'PARALLEL' ! gather-scatter handle for field "ifld"
2988
2989	real d(lx1,ly1,lz1,lelt)
2990
2991	character*3 b ! Boundary condition of interest
2992
2993	integer e,eg,f
2994
2995	nel = nelfld(ifld)
2996	n = lx1ly1lz1*nel
2997
2998	call domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
2999
3000	xmn = min(xmin,ymin)
3001	xmx = max(xmax,ymax)
3002	if (if3d) xmn = min(xmn ,zmin)
3003	if (if3d) xmx = max(xmx ,zmax)
3004
3005	big = 10*(xmx-xmn)
3006	call cfill(d,big,n)
3007
3008
3009	nface = 2*ldim
3010	do e=1,nel ! Set d=0 on walls
3011	do f=1,nface
3012	if (cbc(f,e,ifld).eq.b) call facev(d,e,f,0.,lx1,ly1,lz1)
3013	enddo
3014	enddo
3015
3016	do ipass=1,10000
3017	dmax = 0
3018	nchange = 0
3019	do e=1,nel
3020	do k=1,lz1
3021	do j=1,ly1
3022	do i=1,lx1
3023	i0=max( 1,i-1)
3024	j0=max( 1,j-1)
3025	k0=max( 1,k-1)
3026	i1=min(lx1,i+1)
3027	j1=min(ly1,j+1)
3028	k1=min(lz1,k+1)
3029	do kk=k0,k1
3030	do jj=j0,j1
3031	do ii=i0,i1
3032
3033	if (if3d) then
3034	dtmp = d(ii,jj,kk,e) + dist3d(
3035	$ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e),zm1(ii,jj,kk,e)
3036	$ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e),zm1(i ,j ,k ,e))
3037	else
3038	dtmp = d(ii,jj,kk,e) + dist2d(
3039	$ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e)
3040	$ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e))
3041	endif
3042
3043	if (dtmp.lt.d(i,j,k,e)) then
3044	d(i,j,k,e) = dtmp
3045	nchange = nchange+1
3046	dmax = max(dmax,d(i,j,k,e))
3047	endif
3048	enddo
3049	enddo
3050	enddo
3051
3052	enddo
3053	enddo
3054	enddo
3055	enddo
3056	call fgslib_gs_op(gsh_fld(ifld),d,1,3,0) ! min over all elements
3057	nchange = iglsum(nchange,1)
3058	dmax = glmax(dmax,1)
3059	if (nio.eq.0.and.loglevel.gt.2) write(6,1) ipass,nchange,dmax,b
3060	1 format(i9,i12,1pe12.4,' max distance b: ',a3)
3061	if (nchange.eq.0) goto 1000
3062	enddo
3063	1000 return
3064	end
3065	c-----------------------------------------------------------------------
3066	subroutine distf(d,ifld,b,dmin,emin,xn,yn,zn)
3067
3068	c Generate a distance function to boundary with bc "b".
3069	c This approach does not yet work with periodic boundary conditions.
3070
3071	c INPUT: ifld - field type for which distance function is to be found.
3072	c ifld = 1 for velocity
3073	c ifld = 2 for temperature, etc.
3074
3075	c OUTPUT: d = distance to nearest boundary with boundary condition "b"
3076
3077	c Work arrays: dmin,emin,xn,yn,zn
3078
3079	include 'SIZE'
3080	include 'GEOM' ! Coordinates
3081	include 'INPUT' ! cbc()
3082	include 'TSTEP' ! nelfld
3083	include 'PARALLEL' ! gather-scatter handle for field "ifld"
3084
3085	real d(lx1,ly1,lz1,lelt)
3086	character*3 b
3087
3088	real dmin(lx1,ly1,lz1,lelt),emin(lx1,ly1,lz1,lelt)
3089	real xn(lx1,ly1,lz1,lelt),yn(lx1,ly1,lz1,lelt)
3090	real zn(lx1,ly1,lz1,lelt)
3091
3092
3093	integer e,eg,f
3094
3095	nel = nelfld(ifld)
3096	n = lx1ly1lz1*nel
3097
3098	call domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
3099
3100	xmn = min(xmin,ymin)
3101	xmx = max(xmax,ymax)
3102	if (if3d) xmn = min(xmn ,zmin)
3103	if (if3d) xmx = max(xmx ,zmax)
3104
3105	big = 10*(xmx-xmn)
3106	call cfill (d,big,n)
3107
3108	call opcopy(xn,yn,zn,xm1,ym1,zm1)
3109
3110	nface = 2*ldim
3111	do e=1,nel ! Set d=0 on walls
3112	do f=1,nface
3113	if (cbc(f,e,1).eq.b) call facev(d,e,f,0.,lx1,ly1,lz1)
3114	enddo
3115	enddo
3116
3117	nxyz = lx1ly1lz1
3118
3119	do ipass=1,10000
3120	dmax = 0
3121	nchange = 0
3122	do e=1,nel
3123	do k=1,lz1
3124	do j=1,ly1
3125	do i=1,lx1
3126	i0=max( 1,i-1)
3127	j0=max( 1,j-1)
3128	k0=max( 1,k-1)
3129	i1=min(lx1,i+1)
3130	j1=min(ly1,j+1)
3131	k1=min(lz1,k+1)
3132	do kk=k0,k1
3133	do jj=j0,j1
3134	do ii=i0,i1
3135
3136	dself = d(i,j,k,e)
3137	dneigh = d(ii,jj,kk,e)
3138	if (dneigh.lt.dself) then ! check neighbor's nearest point
3139	d2 = (xm1(i,j,k,e)-xn(ii,jj,kk,e))**2
3140	$ + (ym1(i,j,k,e)-yn(ii,jj,kk,e))**2
3141	if (if3d) d2 = d2 + (zm1(i,j,k,e)-zn(ii,jj,kk,e))**2
3142	if (d2.gt.0) d2 = sqrt(d2)
3143	if (d2.lt.dself) then
3144	nchange = nchange+1
3145	d (i,j,k,e) = d2
3146	xn(i,j,k,e) = xn(ii,jj,kk,e)
3147	yn(i,j,k,e) = yn(ii,jj,kk,e)
3148	zn(i,j,k,e) = zn(ii,jj,kk,e)
3149	dmax = max(dmax,d(i,j,k,e))
3150	endif
3151	endif
3152	enddo
3153	enddo
3154	enddo
3155
3156	enddo
3157	enddo
3158	enddo
3159
3160	re = lglel(e)
3161	call cfill(emin(1,1,1,e),re,nxyz)
3162	call copy (dmin(1,1,1,e),d(1,1,1,e),nxyz)
3163
3164	enddo
3165	nchange = iglsum(nchange,1)
3166
3167	call fgslib_gs_op(gsh_fld(ifld),dmin,1,3,0) ! min over all elements
3168
3169
3170	nchange2=0
3171	do e=1,nel
3172	do i=1,nxyz
3173	if (dmin(i,1,1,e).ne.d(i,1,1,e)) then
3174	nchange2 = nchange2+1
3175	emin(i,1,1,e) = 0 ! Flag
3176	endif
3177	enddo
3178	enddo
3179	call copy(d,dmin,n) ! Ensure updated distance
3180	nchange2 = iglsum(nchange2,1)
3181	nchange = nchange + nchange2
3182	call fgslib_gs_op(gsh_fld(ifld),emin,1,4,0) ! max over all elements
3183
3184	do e=1,nel ! Propagate nearest wall points
3185	do i=1,nxyz
3186	eg = emin(i,1,1,e)
3187	if (eg.ne.lglel(e)) then
3188	xn(i,1,1,e) = 0
3189	yn(i,1,1,e) = 0
3190	zn(i,1,1,e) = 0
3191	endif
3192	enddo
3193	enddo
3194	call fgslib_gs_op(gsh_fld(ifld),xn,1,1,0) ! Sum over all elements to
3195	call fgslib_gs_op(gsh_fld(ifld),yn,1,1,0) ! convey nearest point
3196	call fgslib_gs_op(gsh_fld(ifld),zn,1,1,0) ! to shared neighbor.
3197
3198	dmax = glmax(dmax,1)
3199	if (nio.eq.0) write(6,1) ipass,nchange,dmax
3200	1 format(i9,i12,1pe12.4,' max wall distance 2')
3201	if (nchange.eq.0) goto 1000
3202	enddo
3203	1000 continue
3204
3205	c wgt = 0.3
3206	c call filter_d2(d,lx1,lz1,wgt,.true.)
3207
3208	return
3209	end
3210	c-----------------------------------------------------------------------
3211	subroutine turb_outflow(d,m1,rq,uin)
3212
3213	c . Set div U > 0 in elements with 'O ' bc.
3214	c
3215	c . rq is nominally the ratio of Qout/Qin and is typically 1.5
3216	c
3217	c . uin is normally zero, unless your flow is zero everywhere
3218	c
3219	c . d and m1 are work arrays of size (lx1,ly1,lz1,lelt), assumed persistant
3220	c
3221	c This routine may or may not work with multiple outlets --- it has
3222	c not been tested for this case.
3223	c
3224	c
3225	c TYPICAL USAGE -- ADD THESE LINES TO userchk() in your .usr file:
3226	c (uncommented)
3227	c
3228	c common /myoutflow/ d(lx1,ly1,lz1,lelt),m1(lx1ly1lz1,lelt)
3229	c real m1
3230	c rq = 2.
3231	c uin = 0.
3232	c call turb_outflow(d,m1,rq,uin)
3233	c
3234
3235	include 'SIZE'
3236	include 'TOTAL'
3237
3238	real d(lx2,ly2,lz2,lelt),m1(lx1ly1lz1,lelt)
3239
3240	parameter (lw = 3lx1ly1*lz1)
3241	common /ctmp1/ w(lw)
3242
3243	integer icalld,noutf,e,f
3244	save icalld,noutf
3245	data icalld,noutf /0,0/
3246
3247	real ddmax,cso
3248	save ddmax,cso
3249	logical ifout
3250
3251	character*3 b
3252
3253	n = lx1ly1lz1*nelv
3254	n2 = lx2ly2lz2*nelv
3255	nxyz = lx1ly1lz1
3256	nxyz2 = lx2ly2lz2
3257
3258	if (icalld.eq.0) then
3259	icalld = 1
3260
3261	b = 'O '
3262	call cheap_dist(m1,1,b)
3263
3264	call rzero (d,n2)
3265
3266	ddmax = 0.
3267	noutf = 0
3268
3269	do e=1,nelv
3270	ifout = .false.
3271	do f=1,2*ldim
3272	if (cbc(f,e,1).eq.b) ifout = .true. ! outflow
3273	if (cbc(f,e,1).eq.b) noutf = noutf+1
3274	enddo
3275	if (ifout) then
3276	if (lx2.lt.lx1) then ! Map distance function to Gauss
3277	call maph1_to_l2(d(1,1,1,e),lx2,m1(1,e),lx1,if3d,w,lw)
3278	else
3279	call copy(d(1,1,1,e),m1(1,e),nxyz)
3280	endif
3281	dmax = vlmax(m1(1,e),nxyz)
3282	ddmax = max(ddmax,dmax)
3283	call rzero(m1(1,e),nxyz) ! mask points at outflow
3284	else
3285	call rone (m1(1,e),nxyz)
3286	endif
3287	enddo
3288
3289	ddmax = glamax(ddmax,1)
3290
3291	do e=1,nelv
3292	ifout = .false.
3293	do f=1,2*ldim
3294	if (cbc(f,e,1).eq.b) ifout = .true. ! outflow
3295	enddo
3296	if (ifout) then
3297	do i=1,nxyz2
3298	d(i,1,1,e) = (ddmax - d(i,1,1,e))/ddmax
3299	enddo
3300	endif
3301	enddo
3302	noutf = iglsum(noutf,1)
3303	endif
3304
3305	if (noutf.eq.0) return
3306
3307	if (uin.ne.0) then ! Use user-supplied characteristic velocity
3308	ubar = uin
3309	else
3310	ubar = glsc3(vx,bm1,m1,n) ! Masked average
3311	vbar = glsc3(vy,bm1,m1,n)
3312	wbar = glsc3(vz,bm1,m1,n)
3313	volu = glsc2(bm1,m1,n)
3314	ubar = abs(ubar)+abs(vbar)
3315	if (if3d) ubar = abs(ubar)+abs(wbar)
3316	ubar = ubar/volu
3317	endif
3318
3319	cs = 3(rq-1.)(ubar/ddmax)
3320	if (.not.ifsplit) cs = cs*bd(1)/dt
3321
3322	do i=1,n2
3323	usrdiv(i,1,1,1) = cs(d(i,1,1,1)*2)
3324	enddo
3325
3326	return
3327	end
3328	c-----------------------------------------------------------------------
3329	subroutine add_temp(f2tbc,nbc,npass)
3330
3331	c add multiple passive scalar fields (npass new ones)
3332
3333	include 'SIZE'
3334	include 'TOTAL'
3335
3336	character*3 f2tbc(2,nbc)
3337
3338	do i=1,npass
3339	call add_temp_1(f2tbc,nbc)
3340	enddo
3341
3342	igeom = 2
3343	call setup_topo ! Set gs handles and multiplicity
3344
3345	return
3346	end
3347	c-----------------------------------------------------------------------
3348	subroutine add_temp_1(f2tbc,nbc)
3349
3350	c
3351	c TYPICAL USAGE: Add the below to usrdat().
3352	c
3353	c parameter (lbc=10) ! Maximum number of bc types
3354	c character*3 f2tbc(2,lbc)
3355	c
3356	c f2tbc(1,1) = 'W ' ! Any 'W ' bc is swapped to ft2bc(2,1)
3357	c f2tbc(2,1) = 'I '
3358	c
3359	c f2tbc(1,2) = 'v ' ! Any 'v ' bc is swapped to ft2bc(2,2)
3360	c f2tbc(2,2) = 't '
3361	c
3362	c nbc = 2 ! Number of boundary condition pairings (e.g., W-->t)
3363	c do i=1,ldimt-1
3364	c call add_temp(f2tbc,nbc)
3365	c enddo
3366
3367	include 'SIZE'
3368	include 'TOTAL'
3369	character*3 f2tbc(2,nbc)
3370
3371	integer e,f
3372
3373	nfld=nfield+1
3374
3375	if (nio.eq.0) write(6,*) 'add temp: ',nfld,nfield,istep
3376
3377	nelfld(nfld) = nelfld(nfield)
3378	nel = nelfld(nfield)
3379	call copy (bc(1,1,1,nfld),bc(1,1,1,nfield),30*nel)
3380	call chcopy(cbc(1,1,nfld),cbc(1,1,nfield),36nel)
3381
3382	do k=1,3
3383	cpfld(nfld,k)=cpfld(nfield,k)
3384	call copy (cpgrp(-5,nfld,k),cpgrp(-5,nfield,k),16)
3385	enddo
3386	call icopy(matype(-5,nfld),matype(-5,nfield),16)
3387
3388	param(7) = param(1) ! rhoCP = rho
3389	param(8) = param(2) ! conduct = dyn. visc
3390
3391	ifheat = .true.
3392	ifadvc(nfld) = .true.
3393	iftmsh(nfld) = .true.
3394	ifvarp(nfld) = ifvarp(nfield)
3395	ifdeal(nfld) = ifdeal(nfield)
3396	ifprojfld(nfld) = ldimt_proj.ge.(nfld-1).and.ifprojfld(nfield)
3397
3398	if (nfld.eq.2) ifto = .true.
3399	if (nfld.gt.2) ifpsco(nfld-2) = .true.
3400	if (nfld.gt.2) npscal = npscal+1
3401
3402	ifldmhd = npscal + 3
3403
3404	nfield = nfld
3405
3406	nface = 2*ldim
3407	do k=1,nbc ! BC conversion
3408	do e=1,nelfld(nfld)
3409	do f=1,nface
3410	if (cbc(f,e,nfld).eq.f2tbc(1,k)) cbc(f,e,nfld)=f2tbc(2,k)
3411	enddo
3412	enddo
3413	enddo
3414
3415
3416	return
3417	end
3418	c-----------------------------------------------------------------------
3419	subroutine planar_avg(ua,u,hndl)
3420	c
3421	c Examples:
3422	c
3423	c ! average field in z and then in x
3424	c idir = 3 ! z
3425	c call gtpp_gs_setup(igs_z,nelx*nely,1,nelz,idir)
3426	c call planar_avg(uavg_z,u,igs_z)
3427	c idir = 1 ! x
3428	c call gtpp_gs_setup(igs_x,nelx,nely,nelz,idir)
3429	c call planar_avg(uavg_xz,uavg,igs_z)
3430	c
3431	c Note, mesh has to be extruded in idir (tensor product)
3432	c
3433	include 'SIZE'
3434	include 'TOTAL'
3435
3436	real ua(*)
3437	real u (*)
3438
3439	common /scrcg/ wrk(lx1ly1lz1*lelv)
3440
3441	n = nx1ny1nz1*nelv
3442
3443	call copy(wrk,bm1,n) ! Set the averaging weights
3444	call fgslib_gs_op(hndl,wrk,1,1,0) ! Sum weights
3445	call invcol1(wrk,n)
3446
3447	do i=1,n
3448	ua(i) = bm1(i,1,1,1)u(i)wrk(i)
3449	enddo
3450
3451	call fgslib_gs_op(hndl,ua,1,1,0) ! Sum weighted values
3452
3453	return
3454	end

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CIVL/examples/fortran/nek5000/core/navier5.f

Download in other formats: