KTH framework for Nek5000 toolboxes; testing version  0.0.1
navier5.f
Go to the documentation of this file.
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=lx1*ly1*lz1*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(nx*nx*nz,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=nx*nx*nz
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 = lx1*ly1*lz1*nelv
246  nyz2 = ly2*lz2
247  nxy1 = lx1*ly1
248  nxyz = lx1*ly1*lz1
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=lx1*ly1*lz1)
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 = lx1*ly1*lz1
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=lx1*ly1*lz1*lelv)
518  real vort(lt,3),work1(1),work2(1),u(1),v(1),w(1)
519 
520  parameter(lx=lx1*ly1*lz1)
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 = lx1*ly1*lz1*nelt
530  nxyz = lx1*ly1*lz1
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=lx1*ly1*lz1)
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) = ( (2*j-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(X*Y) - 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 = lx1*ly1*lz1*nelv
921  ntott = lx1*ly1*lz1*nelt
922  nto2 = lx2*ly2*lz2*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) = alpha*avg(k) + beta*f(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) = alpha*avg(k) + beta*f(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) = alpha*avg(k) + beta*f(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(nx*nx),Ljt(nx*nx),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(nx*ny*nlx*nly)
1350  real v(nx*ny*nlx*nly)
1351  real w(nx*ny*nlx*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 = nx*ny*nlx*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(nx*ny*nlx*nly)
1427  real v(nx*ny*nlx*nly)
1428  real w(nx*ny*nlx*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 = nx*ny*nlx*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*(s11*n1 + s12*n2 + s13*n3) ! viscous drag
1551  dg(2,2) = -v*(s21*n1 + s22*n2 + s23*n3)
1552  dg(3,2) = -v*(s31*n1 + s32*n2 + s33*n3)
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) + (r2*dg(3,1)-r3*dg(2,1)) ! pressure
1565  dgtq(2,3) = dgtq(2,3) + (r3*dg(1,1)-r1*dg(3,1)) ! torque
1566  dgtq(3,3) = dgtq(3,3) + (r1*dg(2,1)-r2*dg(1,1))
1567 c
1568  dgtq(1,4) = dgtq(1,4) + (r2*dg(3,2)-r3*dg(2,2)) ! viscous
1569  dgtq(2,4) = dgtq(2,4) + (r3*dg(1,2)-r1*dg(3,2)) ! torque
1570  dgtq(3,4) = dgtq(3,4) + (r1*dg(2,2)-r2*dg(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*(s11*n1 + s12*n2) ! viscous drag
1596  dg(2,2) = -v*(s21*n1 + 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) + (r2*dg(3,1)-r3*dg(2,1)) ! pressure
1610  dgtq(2,3) = 0! dgtq(2,3) + (r3*dg(1,1)-r1*dg(3,1)) ! torque
1611  dgtq(3,3) = dgtq(3,3) + (r1*dg(2,1)-r2*dg(1,1))
1612 
1613  dgtq(1,4) = 0! dgtq(1,4) + (r2*dg(3,2)-r3*dg(2,2)) ! viscous
1614  dgtq(2,4) = 0! dgtq(2,4) + (r3*dg(1,2)-r1*dg(3,2)) ! torque
1615  dgtq(3,4) = dgtq(3,4) + (r1*dg(2,2)-r2*dg(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(lx1*ly1*lz1*6*lelv)
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=lx1*ly1*lz1)
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 = lx1*ly1*lz1*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(lx1*ly1*lz1,nij,lelv)
1834  real u (lx1*ly1*lz1,lelv)
1835  real v (lx1*ly1*lz1,lelv)
1836  real w (lx1*ly1*lz1,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 = lx1*ly1*lz1
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=lx1*ly1*lz1*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 = lx1*ly1*lz1*nelt
1991  n2= lx2*ly2*lz2*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 = wdsize*ldim*(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(lx1*ly1*lz1*6*lelv)
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=lx1*ly1*lz1)
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 = lx1*ly1*lz1*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*(s11*n1 + s12*n2 + s13*n3) ! viscous drag
2280  stress2 = -v*(s21*n1 + s22*n2 + s23*n3)
2281  stress3 = -v*(s31*n1 + s32*n2 + s33*n3)
2282 
2283  strsnrm = stress1*stress1+stress2*stress2+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*(s11*n1 + s12*n2) ! viscous drag
2306  stress2 = -v*(s21*n1 + s22*n2)
2307 
2308  strsnrm = stress1*stress1+stress2*stress2
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 = lx1*ly1*lz1)
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 = lx1*ly1*lz1*nelt
2334  nxyz = lx1*ly1*lz1
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) + si*sj*x(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=n*n*n
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) + si*sj*sk*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) + hj*hk*(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) + hi*hk*(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) + hi*hj*(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-k*iq)-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-k*iq)-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=lx1*ly1*lz1*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 = lx1*ly1*lz1*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=lx1*lx1,lxyz=lx1*ly1*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=lx1*lx1,lxyz=lx1*ly1*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=lx1*ly1*lz1)
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 = lx1*ly1*lz1*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 = lx1*ly1*lz1*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 = lx1*ly1*lz1*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 = lx1*ly1*lz1*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 = lx1*ly1*lz1
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(lx1*ly1*lz1,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(lx1*ly1*lz1,lelt)
3239 
3240  parameter(lw = 3*lx1*ly1*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 = lx1*ly1*lz1*nelv
3254  n2 = lx2*ly2*lz2*nelv
3255  nxyz = lx1*ly1*lz1
3256  nxyz2 = lx2*ly2*lz2
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),3*6*nel)
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(lx1*ly1*lz1*lelv)
3440 
3441  n = nx1*ny1*nz1*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
exitt
void exitt()
Definition: comm_mpi.f:604
byte_write
#define byte_write
Definition: byte.c:39
ident
subroutine ident(a, n)
Definition: calcz.f:147
crecv2
subroutine crecv2(mtype, buf, lenm, jnid)
Definition: comm_mpi.f:333
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
nekgsync
subroutine nekgsync()
Definition: comm_mpi.f:502
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
csend
subroutine csend(mtype, buf, len, jnid, jpid)
Definition: comm_mpi.f:303
bcast
subroutine bcast(buf, len)
Definition: comm_mpi.f:431
err_chk
subroutine err_chk(ierr, string)
Definition: comm_mpi.f:558
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
setup_topo
subroutine setup_topo
Definition: connect1.f:3
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
gllel
integer function gllel(ieg)
Definition: dprocmap.f:183
gllnid
integer function gllnid(ieg)
Definition: dprocmap.f:161
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
tensr3
subroutine tensr3(v, nv, u, nu, A, Bt, Ct, w)
Definition: fasts.f:126
dsavg
subroutine dsavg(u)
Definition: ic.f:1878
restart
subroutine restart(nfiles)
Definition: ic.f:416
geom_reset
subroutine geom_reset(icall)
Definition: ic.f:1802
indx1
integer function indx1(S1, S2, L2)
Definition: ic.f:1357
glmin
function glmin(a, n)
Definition: math.f:973
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
mod1
function mod1(i, n)
Definition: math.f:509
icopy
subroutine icopy(a, b, n)
Definition: math.f:289
glsc2
function glsc2(x, y, n)
Definition: math.f:794
iglsum
function iglsum(a, n)
Definition: math.f:926
glsc3
function glsc3(a, b, mult, n)
Definition: math.f:776
vlamax
real function vlamax(vec, n)
Definition: math.f:406
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
add2sxy
subroutine add2sxy(x, a, y, b, n)
Definition: math.f:1574
blank
subroutine blank(A, N)
Definition: math.f:19
add3
subroutine add3(a, b, c, n)
Definition: math.f:644
iglmax
function iglmax(a, n)
Definition: math.f:913
izero
subroutine izero(a, n)
Definition: math.f:215
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
cadd2
subroutine cadd2(a, b, const, n)
Definition: math.f:346
transpose1
subroutine transpose1(a, n)
Definition: math.f:1605
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
glamax
real function glamax(a, n)
Definition: math.f:874
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rone
subroutine rone(a, n)
Definition: math.f:230
invcol1
subroutine invcol1(a, n)
Definition: math.f:62
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
local_grad3_t
subroutine local_grad3_t(u, ur, us, ut, N, e, D, Dt, w)
Definition: navier1.f:4669
opcopy
subroutine opcopy(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2474
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
local_grad2_t
subroutine local_grad2_t(u, ur, us, N, e, D, Dt, w)
Definition: navier1.f:4694
facint_v
function facint_v(a, area, f, e)
Definition: navier5.f:295
build_new_filter
subroutine build_new_filter(intv, zpts, nx, kut, wght, nio)
Definition: navier5.f:801
rand_fld_h1
subroutine rand_fld_h1(x)
Definition: navier5.f:2687
cheap_dist
subroutine cheap_dist(d, ifld, b)
Definition: navier5.f:2970
comp_sij
subroutine comp_sij(sij, nij, u, v, w, ur, us, ut, vr, vs, vt, wr, ws, wt)
Definition: navier5.f:1824
build_legend_transform
subroutine build_legend_transform(Lj, Ljt, zpts, nx)
Definition: navier5.f:1088
drgtrq
subroutine drgtrq(dgtq, xm0, ym0, zm0, sij, pm1, visc, f, e)
Definition: navier5.f:1490
avg1
subroutine avg1(avg, f, alpha, beta, n, name, ifverbose)
Definition: navier5.f:1019
domain_size
subroutine domain_size(xmin, xmax, ymin, ymax, zmin, zmax)
Definition: navier5.f:2948
outmatx
subroutine outmatx(a, m, n, io, name)
Definition: navier5.f:223
facint2
function facint2(a, b, c, area, ifc, ie)
Definition: navier5.f:370
turb_outflow
subroutine turb_outflow(d, m1, rq, uin)
Definition: navier5.f:3212
add_temp
subroutine add_temp(f2tbc, nbc, npass)
Definition: navier5.f:3330
gaujordf
subroutine gaujordf(a, m, n, indr, indc, ipiv, ierr, rmult)
Definition: navier5.f:668
planar_avg
subroutine planar_avg(ua, u, hndl)
Definition: navier5.f:3420
shear_calc_max
subroutine shear_calc_max(strsmx, scale, x0, ifdout, iftout)
Definition: navier5.f:2133
facint
function facint(a, b, area, ifc, ie)
Definition: navier5.f:330
q_filter
subroutine q_filter(wght)
Definition: navier5.f:3
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
mappr
subroutine mappr(pm1, pm2, pa, pb)
Definition: navier5.f:237
out_el
subroutine out_el(xt, e)
Definition: navier5.f:2090
dump_header2d
subroutine dump_header2d(excode, nx, ny, nlx, nly, ierr)
Definition: navier5.f:1265
rescale_x
subroutine rescale_x(x, x0, x1)
Definition: navier5.f:2702
cut_off_filter
subroutine cut_off_filter(u, mx, ifld)
Definition: navier5.f:2801
avg_all
subroutine avg_all
Definition: navier5.f:883
g_filter
subroutine g_filter(u, diag, ifld)
Definition: navier5.f:2778
auto_averager
subroutine auto_averager(fname127)
Definition: navier5.f:1959
gh_face_extend_2d
subroutine gh_face_extend_2d(x, zg, n, gh_type, e, v)
Definition: navier5.f:2437
out_csrmats
subroutine out_csrmats(acsr, ia, ja, n, name9)
Definition: navier5.f:398
filter_d2
subroutine filter_d2(v, nx, nz, wgt, ifd4)
Definition: navier5.f:2829
get_el
subroutine get_el(xt, x, y, z)
Definition: navier5.f:2111
legendre_poly
subroutine legendre_poly(L, x, N)
Definition: navier5.f:785
dist2d
function dist2d(a, b, x, y)
Definition: navier5.f:2938
get_exyz
subroutine get_exyz(ex, ey, ez, eg, nelx, nely, nelz)
Definition: navier5.f:1252
gh_face_extend_3d
subroutine gh_face_extend_3d(x, zg, n, gh_type, e, v)
Definition: navier5.f:2502
outfld2d_p
subroutine outfld2d_p(u, v, w, nx, ny, nlx, nly, name, ifld, jid, npido, ir)
Definition: navier5.f:1341
avg2
subroutine avg2(avg, f, alpha, beta, n, name, ifverbose)
Definition: navier5.f:1042
facint_a
function facint_a(a, area, f, e)
Definition: navier5.f:274
fix_geom
subroutine fix_geom
Definition: navier5.f:2322
local_err_est
subroutine local_err_est(err, u, nx, Lj, Ljt, uh, w, if3d)
Definition: navier5.f:1117
filterq
subroutine filterq(v, f, nx, nz, w1, w2, ft, if3d, dmax)
Definition: navier5.f:167
outfld2d
subroutine outfld2d(u, v, w, nx, ny, nlx, nly, name, ifld)
Definition: navier5.f:1422
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
ran1
real function ran1(idum)
Definition: navier5.f:2650
dist3d
function dist3d(a, b, c, x, y, z)
Definition: navier5.f:2928
build_filter
subroutine build_filter(f, diag, nx)
Definition: navier5.f:2720
outmesh
subroutine outmesh
Definition: navier5.f:2054
torque_calc
subroutine torque_calc(scale, x0, ifdout, iftout)
Definition: navier5.f:1624
distf
subroutine distf(d, ifld, b, dmin, emin, xn, yn, zn)
Definition: navier5.f:3067
comp_vort3
subroutine comp_vort3(vort, work1, work2, u, v, w)
Definition: navier5.f:513
avg3
subroutine avg3(avg, f, g, alpha, beta, n, name, ifverbose)
Definition: navier5.f:1065
surface_int
subroutine surface_int(sint, sarea, a, e, f)
Definition: navier5.f:598
gradm1
subroutine gradm1(ux, uy, uz, u)
Definition: navier5.f:463
gh_face_extend
subroutine gh_face_extend(x, zg, n, gh_type, e, v)
Definition: navier5.f:2422
surface_flux
subroutine surface_flux(dq, qx, qy, qz, e, f, w)
Definition: navier5.f:632
add_temp_1
subroutine add_temp_1(f2tbc, nbc)
Definition: navier5.f:3349
get_strsmax
subroutine get_strsmax(strsmax, xm0, ym0, zm0, sij, pm1, visc, f, e)
Definition: navier5.f:2228
maph1_to_l2
subroutine maph1_to_l2(a, na, b, nb, if3d, w, ldw)
Definition: navier8.f:1127
outpost
subroutine outpost(v1, v2, v3, vp, vt, name3)
Definition: prepost.f:1378
outpost2
subroutine outpost2(v1, v2, v3, vp, vt, nfldt, name3)
Definition: prepost.f:1394
zwgl
subroutine zwgl(Z, W, NP)
Definition: speclib.f:91
zwgll
subroutine zwgll(Z, W, NP)
Definition: speclib.f:108
iglm
subroutine iglm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2)
Definition: speclib.f:1077
igllm
subroutine igllm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2)
Definition: speclib.f:1102
faddcl3
subroutine faddcl3(a, b, c, iface1)
Definition: subs1.f:978
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342
faccl3
subroutine faccl3(a, b, c, iface1)
Definition: subs1.f:944
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341