KTH framework for Nek5000 toolboxes; testing version  0.0.1
postpro.f
Go to the documentation of this file.
1  subroutine load_fld(string)
2 
3  include 'SIZE'
4  include 'INPUT'
5  include 'RESTART'
6 
7  character string*(*)
8 
9  l = ltrunc(string,len(string))
10  if(l.gt.132) call exitti('invalid string length$',l)
11 
12  call blank (initc(1),132)
13  call chcopy (initc(1),string,l)
14  call setics
15 
16  return
17  end
18 c-----------------------------------------------------------------------
19  subroutine lambda2(l2)
20 c
21 c Generate Lambda-2 vortex of Jeong & Hussein, JFM '95
22 c
23  include 'SIZE'
24  include 'TOTAL'
25 
26  real l2(lx1,ly1,lz1,1)
27 
28  parameter(lxyz=lx1*ly1*lz1)
29 
30  real gije(lxyz,ldim,ldim)
31  real vv(ldim,ldim),ss(ldim,ldim),oo(ldim,ldim),w(ldim,ldim)
32  real lam(ldim)
33 
34  nxyz = lx1*ly1*lz1
35  n = nxyz*nelv
36 
37  do ie=1,nelv
38  ! Compute velocity gradient tensor
39  call comp_gije(gije,vx(1,1,1,ie),vy(1,1,1,ie),vz(1,1,1,ie),ie)
40 
41  do l=1,nxyz
42  ! decompose into symm. and antisymm. part
43  do j=1,ldim
44  do i=1,ldim
45  ss(i,j) = 0.5*(gije(l,i,j)+gije(l,j,i))
46  oo(i,j) = 0.5*(gije(l,i,j)-gije(l,j,i))
47  enddo
48  enddo
49 
50  call rzero(vv,ldim*ldim)
51  do j=1,ldim
52  do i=1,ldim
53  do k=1,ldim
54  vv(i,j) = vv(i,j) + ss(i,k)*ss(k,j) + oo(i,k)*oo(k,j)
55  enddo
56  enddo
57  enddo
58 
59 c Solve eigenvalue problemand sort
60 c eigenvalues in ascending order.
61  call find_lam3(lam,vv,w,ldim,ierr)
62 
63  l2(l,1,1,ie) = lam(2)
64  enddo
65  enddo
66 
67  ! smooth field
68  ifld = ifield
69  ifield = 1
70  wght = 0.5
71  ncut = 1
72  call filter_s0(l2,wght,ncut,'vortx')
73  ifield = ifld
74 
75  return
76  end
77 c-----------------------------------------------------------------------
78  subroutine find_lam3(lam,aa,w,ldim,ierr)
79  real aa(ldim,ldim),lam(ldim),w(ldim,ldim),lam2
80 c
81 c Use cubic eqn. to compute roots
82 c
83  common /ecmnr/ a,b,c,d,e,f,f2,ef,df,r0,r1
84  common /ecmni/ nr
85  common /ecmnl/ iffout,ifdefl
86  logical iffout,ifdefl
87 c
88 c
89  iffout = .false.
90  ierr = 0
91 c
92 c 2D case....
93 c
94 c
95  if (ldim.eq.2) then
96  a = aa(1,1)
97  b = aa(1,2)
98  c = aa(2,1)
99  d = aa(2,2)
100  aq = 1.
101  bq = -(a+d)
102  cq = a*d-c*b
103 c
104  call quadratic(x1,x2,aq,bq,cq,ierr)
105 c
106  lam(1) = min(x1,x2)
107  lam(2) = max(x1,x2)
108 c
109  return
110  endif
111 c
112 c Else ... 3D case....
113 c a d e
114 c Get symmetric 3x3 matrix d b f
115 c e f c
116 c
117  a = aa(1,1)
118  b = aa(2,2)
119  c = aa(3,3)
120  d = 0.5*(aa(1,2)+aa(2,1))
121  e = 0.5*(aa(1,3)+aa(3,1))
122  f = 0.5*(aa(2,3)+aa(3,2))
123  ef = e*f
124  df = d*f
125  f2 = f*f
126 c
127 c
128 c Use cubic eqn. to compute roots
129 c
130 c ax = a-x
131 c bx = b-x
132 c cx = c-x
133 c y = ax*(bx*cx-f2) - d*(d*cx-ef) + e*(df-e*bx)
134 c
135  a1 = -(a+b+c)
136  a2 = (a*b+b*c+a*c) - (d*d+e*e+f*f)
137  a3 = a*f*f + b*e*e + c*d*d - a*b*c - 2*d*e*f
138 c
139  call cubic (lam,a1,a2,a3,ierr)
140  call sort (lam,w,3)
141 c
142  return
143  end
144 c-----------------------------------------------------------------------
145  subroutine quadratic(x1,x2,a,b,c,ierr)
146 c
147 c Stable routine for computation of real roots of quadratic
148 c
149  ierr = 0
150  x1 = 0.
151  x2 = 0.
152 c
153  if (a.eq.0.) then
154  if (b.eq.0) then
155  if (c.ne.0) then
156 c write(6,10) x1,x2,a,b,c
157  ierr = 1
158  endif
159  return
160  endif
161  ierr = 2
162  x1 = -c/b
163 c write(6,11) x1,a,b,c
164  return
165  endif
166 c
167  d = b*b - 4.*a*c
168  if (d.lt.0) then
169  ierr = 1
170 c write(6,12) a,b,c,d
171  return
172  endif
173  if (d.gt.0) d = sqrt(d)
174 c
175  if (b.gt.0) then
176  x1 = -2.*c / ( d+b )
177  x2 = -( d+b ) / (2.*a)
178  else
179  x1 = ( d-b ) / (2.*a)
180  x2 = -2.*c / ( d-b )
181  endif
182 c
183  10 format('ERROR: Both a & b zero in routine quadratic NO ROOTS.'
184  $ ,1p5e12.4)
185  11 format('ERROR: a = 0 in routine quadratic, only one root.'
186  $ ,1p5e12.4)
187  12 format('ERROR: negative discriminate in routine quadratic.'
188  $ ,1p5e12.4)
189 c
190  return
191  end
192 c-----------------------------------------------------------------------
193  subroutine cubic(xo,ai1,ai2,ai3,ierr)
194  real xo(3),ai1,ai2,ai3
195  complex*16 x(3),a1,a2,a3,q,r,d,arg,t1,t2,t3,theta,sq,a13
196 c
197 c Compute real solutions to cubic root eqn. (Num. Rec. v. 1, p. 146)
198 c pff/Sang-Wook Lee Jan 19 , 2004
199 c
200 c Assumption is that all x's are *real*
201 c
202  real*8 twopi
203  save twopi
204  data twopi /6.283185307179586476925286766/
205 c
206  ierr = 0
207 c
208  zero = 0.
209  a1 = cmplx(ai1,zero)
210  a2 = cmplx(ai2,zero)
211  a3 = cmplx(ai3,zero)
212 c
213  q = (a1*a1 - 3*a2)/9.
214  if (q.eq.0) goto 999
215 c
216  r = (2*a1*a1*a1 - 9*a1*a2 + 27*a3)/54.
217 c
218  d = q*q*q - r*r
219 c
220 c if (d.lt.0) goto 999
221 c
222  arg = q*q*q
223  arg = sqrt(arg)
224  arg = r/arg
225 c
226  if (abs(arg).gt.1) goto 999
227  theta = acos(abs(arg))
228 c
229  t1 = theta / 3.
230  t2 = (theta + twopi) / 3.
231  t3 = (theta + 2.*twopi) / 3.
232 c
233  sq = -2.*sqrt(q)
234  a13 = a1/3.
235  x(1) = sq*cos(t1) - a13
236  x(2) = sq*cos(t2) - a13
237  x(3) = sq*cos(t3) - a13
238 c
239  xo(1) = real(x(1))
240  xo(2) = real(x(2))
241  xo(3) = real(x(3))
242 c
243  return
244 c
245  999 continue ! failed
246  ierr = 1
247  call rzero(x,3)
248 
249  return
250  end
251 c-----------------------------------------------------------------------
252  subroutine comp_gije(gije,u,v,w,e)
253 c
254 c du_i
255 c Compute the gradient tensor G_ij := ---- , for element e
256 c du_j
257 c
258  include 'SIZE'
259  include 'TOTAL'
260 
261  real gije(lx1*ly1*lz1,ldim,ldim)
262  real u (lx1*ly1*lz1)
263  real v (lx1*ly1*lz1)
264  real w (lx1*ly1*lz1)
265 
266  real ur (lx1*ly1*lz1)
267  real us (lx1*ly1*lz1)
268  real ut (lx1*ly1*lz1)
269 
270  integer e
271 
272  n = lx1-1 ! Polynomial degree
273  nxyz = lx1*ly1*lz1
274 
275  if (if3d) then ! 3D CASE
276 
277  do k=1,3
278  if (k.eq.1) call local_grad3(ur,us,ut,u,n,1,dxm1,dxtm1)
279  if (k.eq.2) call local_grad3(ur,us,ut,v,n,1,dxm1,dxtm1)
280  if (k.eq.3) call local_grad3(ur,us,ut,w,n,1,dxm1,dxtm1)
281 
282  do i=1,nxyz
283  dj = jacmi(i,e)
284 
285  ! d/dx
286  gije(i,k,1) = dj*(
287  $ ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e))
288  ! d/dy
289  gije(i,k,2) = dj*(
290  $ ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e))
291  ! d/dz
292  gije(i,k,3) = dj*(
293  $ ur(i)*rzm1(i,1,1,e)+us(i)*szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e))
294 
295  enddo
296  enddo
297 
298  elseif (ifaxis) then ! AXISYMMETRIC CASE
299  if(nid.eq.0) write(6,*)
300  & 'ABORT: comp_gije no axialsymmetric support for now'
301  call exitt
302  else ! 2D CASE
303 
304  do k=1,2
305  if (k.eq.1) call local_grad2(ur,us,u,n,1,dxm1,dxtm1)
306  if (k.eq.2) call local_grad2(ur,us,v,n,1,dxm1,dxtm1)
307  do i=1,nxyz
308  dj = jacmi(i,e)
309  ! d/dx
310  gije(i,k,1)=dj*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
311  ! d/dy
312  gije(i,k,2)=dj*(ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e))
313  enddo
314  enddo
315  endif
316 
317  return
318  end
319 c-----------------------------------------------------------------------
320  subroutine filter_s1(scalar,tf,nx,nel) ! filter scalar field
321 
322  include 'SIZE'
323 
324  parameter(lxyz=lx1*ly1*lz1)
325  real scalar(lxyz,1)
326  real fh(nx*nx),fht(nx*nx),tf(nx)
327 
328  real w1(lxyz,lelt)
329 
330 c Build 1D-filter based on the transfer function (tf)
331  call build_1d_filt(fh,fht,tf,nx,nio)
332 
333 c Filter scalar
334  call copy(w1,scalar,lxyz*nel)
335  do ie=1,nel
336  call tens3d1(scalar(1,ie),w1(1,ie),fh,fht,lx1,lx1) ! fh x fh x fh x scalar
337  enddo
338 
339  return
340  end
341 c-----------------------------------------------------------------------
342  subroutine filter_s0(scalar,wght,ncut,name5) ! filter scalar field
343 
344  include 'SIZE'
345  include 'TOTAL'
346 
347  real scalar(1)
348  character*5 name5
349 
350  parameter(l1=lx1*lx1)
351  real intdv(l1),intuv(l1),intdp(l1),intup(l1),intv(l1),intp(l1)
352  save intdv ,intuv ,intdp ,intup ,intv ,intp
353 
354  common /ctmp0/ intt
355  common /screv/ wk1,wk2
356  common /scrvh/ zgmv,wgtv,zgmp,wgtp,tmax(100)
357 
358  real intt (lx1,lx1)
359  real wk1 (lx1,lx1,lx1,lelt)
360  real wk2 (lx1,lx1,lx1)
361  real zgmv (lx1),wgtv(lx1),zgmp(lx1),wgtp(lx1)
362 
363 
364  integer icall
365  save icall
366  data icall /0/
367 
368  logical ifdmpflt
369 
370  imax = nid
371  imax = iglmax(imax,1)
372  jmax = iglmax(imax,1)
373 
374 c if (icall.eq.0) call build_new_filter(intv,zgm1,lx1,ncut,wght,nio)
375  call build_new_filter(intv,zgm1,lx1,ncut,wght,nio)
376 
377  icall = 1
378 
379  call filterq(scalar,intv,lx1,lz1,wk1,wk2,intt,if3d,fmax)
380  fmax = glmax(fmax,1)
381 
382  if (nio.eq.0) write(6,1) istep,fmax,name5
383  1 format(i8,' sfilt:',1pe12.4,a10)
384 
385  return
386  end
387 c-----------------------------------------------------------------------
388  subroutine tens3d1(v,u,f,ft,nv,nu) ! v = F x F x F x u
389 
390 c Note: this routine assumes that lx1=ly1=lz1
391 c
392  include 'SIZE'
393  include 'INPUT'
394 
395  parameter(lw=4*lx1*lx1*lz1)
396  common /ctensor/ w1(lw),w2(lw)
397 
398  real v(nv,nv,nv),u(nu,nu,nu)
399  real f(1),ft(1)
400 
401  if (nu*nu*nv.gt.lw) then
402  write(6,*) nid,nu,nv,lw,' ERROR in tens3d1. Increase lw.'
403  call exitt
404  endif
405 
406  if (if3d) then
407  nuv = nu*nv
408  nvv = nv*nv
409  call mxm(f,nv,u,nu,w1,nu*nu)
410  k=1
411  l=1
412  do iz=1,nu
413  call mxm(w1(k),nv,ft,nu,w2(l),nv)
414  k=k+nuv
415  l=l+nvv
416  enddo
417  call mxm(w2,nvv,ft,nu,v,nv)
418  else
419  call mxm(f ,nv,u,nu,w1,nu)
420  call mxm(w1,nv,ft,nu,v,nv)
421  endif
422  return
423  end
424 c-----------------------------------------------------------------------
425  subroutine build_1d_filt(fh,fht,trnsfr,nx,nio)
426 c
427 c This routing builds a 1D filter with transfer function diag()
428 c
429 c Here, nx = number of points
430 c
431  real fh(nx,nx),fht(nx,nx),trnsfr(nx)
432 c
433  parameter(lm=40)
434  parameter(lm2=lm*lm)
435  common /cfiltr/ phi(lm2),pht(lm2),diag(lm2),rmult(lm),lj(lm)
436  $ , zpts(lm)
437  real Lj
438 
439  common /cfilti/ indr(lm),indc(lm),ipiv(lm)
440 c
441  if (nx.gt.lm) then
442  write(6,*) 'ABORT in set_filt:',nx,lm
443  call exitt
444  endif
445 
446  call zwgll(zpts,rmult,nx)
447 
448  kj = 0
449  n = nx-1
450  do j=1,nx
451  z = zpts(j)
452  call legendre_poly(lj,z,n)
453  kj = kj+1
454  pht(kj) = lj(1)
455  kj = kj+1
456  pht(kj) = lj(2)
457  do k=3,nx
458  kj = kj+1
459  pht(kj) = lj(k)-lj(k-2)
460  enddo
461  enddo
462  call transpose (phi,nx,pht,nx)
463  call copy (pht,phi,nx*nx)
464  call gaujordf (pht,nx,nx,indr,indc,ipiv,ierr,rmult)
465 
466  call rzero(diag,nx*nx)
467  k=1
468  do i=1,nx
469  diag(k) = trnsfr(i)
470  k = k+(nx+1)
471  enddo
472 
473  call mxm (diag,nx,pht,nx,fh,nx) ! -1
474  call mxm (phi ,nx,fh,nx,pht,nx) ! V D V
475 
476  call copy (fh,pht,nx*nx)
477  call transpose (fht,nx,fh,nx)
478 
479  do k=1,nx*nx
480  pht(k) = 1.-diag(k)
481  enddo
482  np1 = nx+1
483  if (nio.eq.0) then
484  write(6,6) 'flt amp',(pht(k),k=1,nx*nx,np1)
485  write(6,6) 'flt trn',(diag(k),k=1,nx*nx,np1)
486  6 format(a8,16f7.4,6(/,8x,16f7.4))
487  endif
488 
489  return
490  end
491 c-----------------------------------------------------------------------
492  subroutine mag_tensor_e(mag,aije)
493 c
494 c Compute magnitude of tensor A_e for element e
495 c
496 c mag(A_e) = sqrt( 0.5 (A:A) )
497 c
498  include 'SIZE'
499  REAL mag (lx1*ly1*lz1)
500  REAL aije(lx1*ly1*lz1,ldim,ldim)
501 
502  nxyz = lx1*ly1*lz1
503 
504  call rzero(mag,nxyz)
505 
506  do 100 j=1,ldim
507  do 100 i=1,ldim
508  do 100 l=1,nxyz
509  mag(l) = mag(l) + 0.5*aije(l,i,j)*aije(l,i,j)
510  100 continue
511 
512  call vsqrt(mag,nxyz)
513 
514  return
515  end
516 c-----------------------------------------------------------------------
517  subroutine comp_sije(gije)
518 c
519 c Compute symmetric part of a tensor G_ij for element e
520 c
521  include 'SIZE'
522  include 'TOTAL'
523 
524  real gije(lx1*ly1*lz1,ldim,ldim)
525 
526  nxyz = lx1*ly1*lz1
527 
528  k = 1
529 
530  do j=1,ldim
531  do i=k,ldim
532  do l=1,nxyz
533  gije(l,i,j) = 0.5*(gije(l,i,j)+gije(l,j,i))
534  gije(l,j,i) = gije(l,i,j)
535  enddo
536  enddo
537  k = k + 1
538  enddo
539 
540  return
541  end
542 c-----------------------------------------------------------------------
543  subroutine map2reg(ur,n,u,nel)
544 c
545 c Map scalar field u() to regular n x n x n array ur
546 c
547  include 'SIZE'
548  real ur(1),u(lx1*ly1*lz1,1)
549 
550  integer e
551 
552  ldr = n**ldim
553 
554  k=1
555  do e=1,nel
556  if (ldim.eq.2) call map2reg_2di_e(ur(k),n,u(1,e),lx1)
557  if (ldim.eq.3) call map2reg_3di_e(ur(k),n,u(1,e),lx1)
558  k = k + ldr
559  enddo
560 
561  return
562  end
563 c-----------------------------------------------------------------------
564  subroutine map2reg_2di_e(uf,n,uc,m) ! Fine, uniform pt
565 
566  real uf(n,n),uc(m,m)
567 
568  parameter(l=50)
569  common /cmap2d/ j(l*l),jt(l*l),w(l*l),z(l)
570 
571  integer mo,no
572  save mo,no
573  data mo,no / 0,0 /
574 
575  if (m.gt.l) call exitti('map2reg_2di_e memory 1$',m)
576  if (n.gt.l) call exitti('map2reg_2di_e memory 2$',n)
577 
578  if (m.ne.mo .or. n.ne.no ) then
579 
580  call zwgll (z,w,m)
581  call zuni (w,n)
582 
583  call gen_int_gz(j,jt,w,n,z,m)
584 
585  endif
586 
587  call mxm(j,n,uc,m,w ,m)
588  call mxm(w,n,jt,m,uf,n)
589 
590  return
591  end
592 c-----------------------------------------------------------------------
593  subroutine map2reg_3di_e(uf,n,uc,m) ! Fine, uniform pt
594 
595  real uf(n,n,n),uc(m,m,m)
596 
597  parameter(l=16)
598  common /cmap3d/ j(l*l),jt(l*l),v(l*l*l),w(l*l*l),z(l)
599 
600  integer mo,no
601  save mo,no
602  data mo,no / 0,0 /
603 
604  if (m.gt.l) call exitti('map2reg_3di_e memory 1$',m)
605  if (n.gt.l) call exitti('map2reg_3di_e memory 2$',n)
606 
607  if (m.ne.mo .or. n.ne.no ) then
608 
609  call zwgll (z,w,m)
610  call zuni (w,n)
611 
612  call gen_int_gz(j,jt,w,n,z,m)
613 
614  endif
615 
616  mm = m*m
617  mn = m*n
618  nn = n*n
619 
620  call mxm(j,n,uc,m,v ,mm)
621  iv=1
622  iw=1
623  do k=1,m
624  call mxm(v(iv),n,jt,m,w(iw),n)
625  iv = iv+mn
626  iw = iw+nn
627  enddo
628  call mxm(w,nn,jt,m,uf,n)
629 
630  return
631  end
632 c-----------------------------------------------------------------------
633  subroutine gen_int_gz(j,jt,g,n,z,m)
634 
635 c Generate interpolater from m z points to n g points
636 
637 c j = interpolation matrix, mapping from z to g
638 c jt = transpose of interpolation matrix
639 c m = number of points on z grid
640 c n = number of points on g grid
641 
642  real j(n,m),jt(m,n),g(n),z(m)
643 
644  mpoly = m-1
645  do i=1,n
646  call fd_weights_full(g(i),z,mpoly,0,jt(1,i))
647  enddo
648 
649  call transpose(j,n,jt,m)
650 
651  return
652  end
653 c-----------------------------------------------------------------------
654  subroutine zuni(z,np)
655 c
656 c Generate equaly spaced np points on the interval [-1:1]
657 c
658  real z(1)
659 
660  dz = 2./(np-1)
661  z(1) = -1.
662  do i = 2,np-1
663  z(i) = z(i-1) + dz
664  enddo
665  z(np) = 1.
666 
667  return
668  end
669 c-----------------------------------------------------------------------
670  subroutine gen_re2(imid) ! Generate and output essential parts of .rea
671  ! And re2
672  ! Clobbers ccurve()
673  ! byte read is float size..
674  ! 4 wdsize
675  include 'SIZE'
676  include 'TOTAL'
677 
678  character*80 hdr
679  real*4 test
680  data test / 6.54321 /
681  integer ierr
682 
683 
684 c imid = 0 ! No midside node defs
685 c imid = 1 ! Midside defs where current curve sides don't exist
686 c imid = 2 ! All nontrivial midside node defs
687 
688  ierr = 0
689  if (nid.eq.0) then
690  call byte_open('newre2.re2' // char(0), ierr)
691  call blank(hdr,80)
692  if(wdsize.eq.8) then
693  write(hdr,112) nelgt,ldim,nelgv
694  else
695  write(hdr,111) nelgt,ldim,nelgv
696  endif
697  111 format('#v001',i9,i3,i9,' hdr')
698  112 format('#v002',i9,i3,i9,' hdr')
699  if(ierr.eq.0) call byte_write(hdr,20,ierr)
700  if(ierr.eq.0) call byte_write(test,1,ierr) !write endian discriminator
701  endif
702  call err_chk(ierr,'Error opening in gen_re2$')
703 
704  call gen_re2_xyz
705  call gen_re2_curve(imid) ! Clobbers ccurve()
706 
707  do ifld=1,nfield
708  call gen_re2_bc (ifld)
709  enddo
710 
711  if (nid.eq.0) call byte_close(ierr)
712  call err_chk(ierr,'Error closing in gen_re2$')
713 
714  return
715  end
716 c-----------------------------------------------------------------------
717  subroutine gen_re2_xyz
718  include 'SIZE'
719  include 'TOTAL'
720 
721  parameter(lv=2**ldim,lblock=1000)
722  common /scrns/ xyz(lv,ldim,lblock),wk(lv*ldim*lblock)
723  common /scruz/ igr(lblock)
724 
725  integer e,eb,eg,ierr,wdsiz2
726 
727  integer isym2pre(8) ! Symmetric-to-prenek vertex ordering
728  save isym2pre
729  data isym2pre / 1 , 2 , 4 , 3 , 5 , 6 , 8 , 7 /
730 
731  real*4 buf(50) ! nwds * 2 for double precision
732  real buf2(25) ! double precsn
733  equivalence(buf,buf2)
734 
735 
736  nxs = lx1-1
737  nys = ly1-1
738  nzs = lz1-1
739 
740  wdsiz2=4
741  if(wdsize.eq.8) wdsiz2=8
742  nblock = lv*ldim*lblock !memory size of data blocks
743 
744  ierr=0
745 
746  do eb=1,nelgt,lblock
747  nemax = min(eb+lblock-1,nelgt)
748  call rzero(xyz,nblock)
749  call izero(igr,lblock)
750  kb = 0
751  do eg=eb,nemax
752  mid = gllnid(eg)
753  e = gllel(eg)
754  kb = kb+1
755  l = 0
756  if (mid.eq.nid.and.if3d) then ! fill owning processor
757  igr(kb) = igroup(e)
758  do k=0,1
759  do j=0,1
760  do i=0,1
761  l=l+1
762  li=isym2pre(l)
763  xyz(li,1,kb) = xm1(1+i*nxs,1+j*nys,1+k*nzs,e)
764  xyz(li,2,kb) = ym1(1+i*nxs,1+j*nys,1+k*nzs,e)
765  xyz(li,3,kb) = zm1(1+i*nxs,1+j*nys,1+k*nzs,e)
766  enddo
767  enddo
768  enddo
769  elseif (mid.eq.nid) then ! 2D
770  igr(kb) = igroup(e)
771  do j=0,1
772  do i=0,1
773  l =l+1
774  li=isym2pre(l)
775  xyz(li,1,kb) = xm1(1+i*nxs,1+j*nys,1,e)
776  xyz(li,2,kb) = ym1(1+i*nxs,1+j*nys,1,e)
777  enddo
778  enddo
779  endif
780  enddo
781  call gop(xyz,wk,'+ ',nblock) ! Sum across all processors
782  call igop(igr,wk,'+ ',lblock) ! Sum across all processors
783 
784  if (nid.eq.0.and.ierr.eq.0) then
785  kb = 0
786  do eg=eb,nemax
787  kb = kb+1
788 
789  if(wdsiz2.eq.8) then
790  call rgrp=igr(kb)
791  call byte_write(rgrp,2,ierr)
792 
793  if(if3d) then
794  buf2(1) = xyz(1,1,kb)
795  buf2(2) = xyz(2,1,kb)
796  buf2(3) = xyz(3,1,kb)
797  buf2(4) = xyz(4,1,kb)
798 
799  buf2(5) = xyz(5,1,kb)
800  buf2(6) = xyz(6,1,kb)
801  buf2(7) = xyz(7,1,kb)
802  buf2(8) = xyz(8,1,kb)
803 
804  buf2(9) = xyz(1,2,kb)
805  buf2(10)= xyz(2,2,kb)
806  buf2(11)= xyz(3,2,kb)
807  buf2(12)= xyz(4,2,kb)
808 
809  buf2(13)= xyz(5,2,kb)
810  buf2(14)= xyz(6,2,kb)
811  buf2(15)= xyz(7,2,kb)
812  buf2(16)= xyz(8,2,kb)
813 
814  buf2(17)= xyz(1,3,kb)
815  buf2(18)= xyz(2,3,kb)
816  buf2(19)= xyz(3,3,kb)
817  buf2(20)= xyz(4,3,kb)
818 
819  buf2(21)= xyz(5,3,kb)
820  buf2(22)= xyz(6,3,kb)
821  buf2(23)= xyz(7,3,kb)
822  buf2(24)= xyz(8,3,kb)
823 
824  if(ierr.eq.0) call byte_write(buf,48,ierr)
825  else
826  buf2(1) = xyz(1,1,kb)
827  buf2(2) = xyz(2,1,kb)
828  buf2(3) = xyz(3,1,kb)
829  buf2(4) = xyz(4,1,kb)
830 
831  buf2(5) = xyz(1,2,kb)
832  buf2(6) = xyz(2,2,kb)
833  buf2(7) = xyz(3,2,kb)
834  buf2(8) = xyz(4,2,kb)
835 
836  if(ierr.eq.0) call byte_write(buf,16,ierr)
837  endif
838  else !!!! 4byte precision !!!!
839  call byte_write(igr(kb),1,ierr)
840  if (if3d) then
841 
842  buf(1) = xyz(1,1,kb)
843  buf(2) = xyz(2,1,kb)
844  buf(3) = xyz(3,1,kb)
845  buf(4) = xyz(4,1,kb)
846 
847  buf(5) = xyz(5,1,kb)
848  buf(6) = xyz(6,1,kb)
849  buf(7) = xyz(7,1,kb)
850  buf(8) = xyz(8,1,kb)
851 
852  buf(9) = xyz(1,2,kb)
853  buf(10) = xyz(2,2,kb)
854  buf(11) = xyz(3,2,kb)
855  buf(12) = xyz(4,2,kb)
856 
857  buf(13) = xyz(5,2,kb)
858  buf(14) = xyz(6,2,kb)
859  buf(15) = xyz(7,2,kb)
860  buf(16) = xyz(8,2,kb)
861 
862  buf(17) = xyz(1,3,kb)
863  buf(18) = xyz(2,3,kb)
864  buf(19) = xyz(3,3,kb)
865  buf(20) = xyz(4,3,kb)
866 
867  buf(21) = xyz(5,3,kb)
868  buf(22) = xyz(6,3,kb)
869  buf(23) = xyz(7,3,kb)
870  buf(24) = xyz(8,3,kb)
871 
872  if(ierr.eq.0) call byte_write(buf,24,ierr)
873 
874  else ! 2D
875  buf(1) = xyz(1,1,kb)
876  buf(2) = xyz(2,1,kb)
877  buf(3) = xyz(3,1,kb)
878  buf(4) = xyz(4,1,kb)
879 
880  buf(5) = xyz(1,2,kb)
881  buf(6) = xyz(2,2,kb)
882  buf(7) = xyz(3,2,kb)
883  buf(8) = xyz(4,2,kb)
884 
885  if(ierr.eq.0) call byte_write(buf,8,ierr)
886  endif
887  endif
888 
889  enddo
890  endif
891  enddo
892  call err_chk(ierr,'Error writing to newre2.re2 in gen_re2_xyz$')
893 
894  return
895  end
896 c-----------------------------------------------------------------------
897  subroutine gen_re2_curve(imid)
898 
899 c This routine is complex because we must first count number of
900 c nontrivial curved sides.
901 
902 c A two pass strategy is used: first count, then write
903 
904  include 'SIZE'
905  include 'TOTAL'
906 
907  integer e,eb,eg,wdsiz2
908  character*1 cc(4)
909 
910  real*4 buf(16) ! nwds * 2 for double precision
911  real buf2( 8) ! double precsn
912  equivalence(buf,buf2)
913 
914  parameter(lblock=500)
915  common /scrns/ vcurve(5,12,lblock),wk(5*12*lblock)
916  common /scruz/ icurve(12,lblock)
917 
918  wdsiz2=4
919  if(wdsize.eq.8) wdsiz2=8
920 
921  if (imid.gt.0) then
922 
923 c imid = 0 ! No midside node defs
924 c imid = 1 ! Midside defs where current curve sides don't exist
925 c imid = 2 ! All nontrivial midside node defs
926 
927  if (imid.eq.2) call blank(ccurve,12*lelt)
928 
929  do e=1,nelt
930  call gen_rea_midside_e(e)
931  enddo
932 
933  endif
934  nedge = 4 + 8*(ldim-2)
935 
936  ncurvn = 0
937  do e=1,nelt
938  do i=1,nedge
939  if (ccurve(i,e).ne.' ') ncurvn = ncurvn+1
940  enddo
941  enddo
942  ncurvn = iglsum(ncurvn,1)
943 
944  ierr=0
945 
946  if(nid.eq.0.and.wdsiz2.eq.8) then
947  rcurvn = ncurvn
948  call byte_write(rcurvn,2,ierr)
949  elseif(nid.eq.0) then
950  call byte_write(ncurvn,1,ierr)
951  endif
952 
953  do eb=1,nelgt,lblock
954 
955  nemax = min(eb+lblock-1,nelgt)
956  call izero(icurve,12*lblock)
957  call rzero(vcurve,60*lblock)
958 
959  kb = 0
960  do eg=eb,nemax
961  mid = gllnid(eg)
962  e = gllel(eg)
963  kb = kb+1
964  if (mid.eq.nid) then ! fill owning processor
965  do i=1,nedge
966  icurve(i,kb) = 0
967  if (ccurve(i,e).eq.'C') icurve(i,kb) = 1
968  if (ccurve(i,e).eq.'s') icurve(i,kb) = 2
969  if (ccurve(i,e).eq.'m') icurve(i,kb) = 3
970  call copy(vcurve(1,i,kb),curve(1,i,e),5)
971  enddo
972  endif
973  enddo
974  call igop(icurve,wk,'+ ',12*lblock) ! Sum across all processors
975  call gop(vcurve,wk,'+ ',60*lblock) ! Sum across all processors
976 
977  if (nid.eq.0) then
978  kb = 0
979  do eg=eb,nemax
980  kb = kb+1
981 
982  do i=1,nedge
983  ii = icurve(i,kb) ! equivalenced to s4
984  if (ii.ne.0) then
985  if (ii.eq.1) cc(1)='C'
986  if (ii.eq.2) cc(1)='s'
987  if (ii.eq.3) cc(1)='m'
988 
989  if(wdsiz2.eq.8) then
990 
991  buf2(1) = eg
992  buf2(2) = i
993  call copy (buf2(3),vcurve(1,i,kb),5)!real*8 write
994  call blank(buf2(8),8)
995  call chcopy(buf2(8),cc,4)
996  iz = 16
997  else
998  call icopy(buf(1),eg,1)
999  call icopy(buf(2), i,1)
1000  call copyx4(buf(3),vcurve(1,i,kb),5) !real*4 write
1001  call blank(buf(8),4)
1002  call chcopy(buf(8),cc,4)
1003  iz = 8
1004  endif
1005 
1006  if(ierr.eq.0) call byte_write(buf,iz,ierr)
1007  endif
1008  enddo
1009  enddo
1010  endif
1011 
1012  enddo
1013  call err_chk(ierr,'Error writing to newre2.re2 in gen_re2_curve$')
1014 
1015  return
1016  end
1017 c-----------------------------------------------------------------------
1018  subroutine gen_re2_bc (ifld)
1019 
1020  include 'SIZE'
1021  include 'TOTAL'
1022 
1023  integer e,eb,eg,wdsiz2
1024 
1025  parameter(lblock=500)
1026  common /scrns/ vbc(5,6,lblock),wk(5*6*lblock)
1027  common /scruz/ ibc(6,lblock)
1028 
1029  character*1 s4(4)
1030  character*3 s3
1031  integer i4
1032  equivalence(i4,s4)
1033  equivalence(s3,s4)
1034 
1035  character*1 chtemp
1036  save chtemp
1037  data chtemp /' '/ ! For mesh bcs
1038 
1039  real*4 buf(16) ! nwds * 2 for double precision
1040  real buf2( 8) ! double precsn
1041  equivalence(buf,buf2)
1042 
1043 
1044  nface = 2*ldim
1045  ierr = 0
1046  nbc = 0
1047  rbc = 0
1048 
1049  wdsiz2=4
1050  if(wdsize.eq.8) wdsiz2=8
1051 
1052  nlg = nelg(ifld)
1053 
1054  if (ifld.eq.1.and..not.ifflow) then ! NO B.C.'s for this field
1055  if(nid.eq.0.and.wdsiz2.eq.4) call byte_write(nbc,1,ierr)
1056  if(nid.eq.0.and.wdsiz2.eq.8) call byte_write(rbc,2,ierr)
1057  call err_chk(ierr,'Error writing to newre2.re2 in gen_re2_bc$')
1058  return
1059  endif
1060 
1061  do ii = 1,nelt
1062  do jj = 1,nface
1063  if(cbc(jj,ii,ifld).ne.'E ')nbc=nbc+1
1064  enddo
1065  enddo
1066  call igop(nbc,wk,'+ ', 1 ) ! Sum across all processors
1067  if(nid.eq.0.and.wdsiz2.eq.8) then
1068  rbc = nbc
1069  call byte_write(rbc,2,ierr)
1070  elseif(nid.eq.0) then
1071  call byte_write(nbc,1,ierr)
1072  endif
1073 
1074  do eb=1,nlg,lblock
1075  nemax = min(eb+lblock-1,nlg)
1076  call izero(ibc, 6*lblock)
1077  call rzero(vbc,30*lblock)
1078  kb = 0
1079  do eg=eb,nemax
1080  mid = gllnid(eg)
1081  e = gllel(eg)
1082  kb = kb+1
1083  if (mid.eq.nid) then ! fill owning processor
1084  do i=1,nface
1085  i4 = 0
1086  call chcopy(s4,cbc(i,e,ifld),3)
1087  ibc(i,kb) = i4
1088  call copy(vbc(1,i,kb),bc(1,i,e,ifld),5)
1089  enddo
1090  endif
1091  enddo
1092  call igop(ibc,wk,'+ ', 6*lblock) ! Sum across all processors
1093  call gop(vbc,wk,'+ ',30*lblock) ! Sum across all processors
1094 
1095  if (nid.eq.0) then
1096  kb = 0
1097  do eg=eb,nemax
1098  kb = kb+1
1099 
1100  do i=1,nface
1101  i4 = ibc(i,kb) ! equivalenced to s4
1102  if (s3.ne.'E ') then
1103 
1104  if(wdsiz2.eq.8) then
1105  buf2(1)=eg
1106  buf2(2)=i
1107  call copy (buf2(3),vbc(1,i,eg),5)
1108  call blank (buf2(8),8)
1109  call chcopy (buf2(8),s3,3)
1110  if(nlg.ge.1000000) then
1111  call icopy(i_vbc,vbc(1,i,eg),1)
1112  buf2(3)=i_vbc
1113  endif
1114  iz=16
1115  else
1116  call icopy (buf(1),eg,1)
1117  call icopy (buf(2),i,1)
1118  call copyx4 (buf(3),vbc(1,i,eg),5)
1119  call blank (buf(8),4)
1120  if(nlg.ge.1000000)call icopy(buf(3),vbc(1,i,eg),1)
1121  call chcopy (buf(8),s3,3)
1122  iz=8
1123  endif
1124 
1125  if(ierr.eq.0) call byte_write (buf,iz,ierr)
1126 
1127  endif
1128  enddo
1129  enddo
1130  endif
1131  enddo
1132  call err_chk(ierr,'Error writing to newre2.re2 in gen_re2_bc$')
1133 
1134  return
1135  end
1136 c-----------------------------------------------------------------------
1137  subroutine gen_rea(imid) ! Generate and output essential parts of .rea
1138  ! Clobbers ccurve()
1139  include 'SIZE'
1140  include 'TOTAL'
1141 
1142 c imid = 0 ! No midside node defs
1143 c imid = 1 ! Midside defs where current curve sides don't exist
1144 c imid = 2 ! All nontrivial midside node defs
1145 
1146  if (nid.eq.0) open(unit=10,file='newrea.out',status='unknown') ! clobbers existing file
1147 
1148  call gen_rea_xyz
1149 
1150  call gen_rea_curve(imid) ! Clobbers ccurve()
1151 
1152  if (nid.eq.0) write(10,*)' ***** BOUNDARY CONDITIONS *****'
1153  do ifld=1,nfield
1154  call gen_rea_bc (ifld)
1155  enddo
1156 
1157  if (nid.eq.0) close(10)
1158 
1159  return
1160  end
1161 c-----------------------------------------------------------------------
1162  subroutine gen_rea_xyz
1163  include 'SIZE'
1164  include 'TOTAL'
1165 
1166  parameter(lv=2**ldim,lblock=1000)
1167  common /scrns/ xyz(lv,ldim,lblock),wk(lv*ldim*lblock)
1168  common /scruz/ igr(lblock)
1169 
1170  integer e,eb,eg
1171  character*1 letapt
1172 
1173  integer isym2pre(8) ! Symmetric-to-prenek vertex ordering
1174  save isym2pre
1175  data isym2pre / 1 , 2 , 4 , 3 , 5 , 6 , 8 , 7 /
1176 
1177  letapt = 'a'
1178  numapt = 1
1179 
1180  nxs = lx1-1
1181  nys = ly1-1
1182  nzs = lz1-1
1183  nblock = lv*ldim*lblock
1184 
1185  if (nid.eq.0)
1186  $ write(10,'(i12,i3,i12,'' NEL,ldim,NELV'')') nelgt,ldim,nelgv
1187 
1188  do eb=1,nelgt,lblock
1189  nemax = min(eb+lblock-1,nelgt)
1190  call rzero(xyz,nblock)
1191  call izero(igr,lblock)
1192  kb = 0
1193  do eg=eb,nemax
1194  mid = gllnid(eg)
1195  e = gllel(eg)
1196  kb = kb+1
1197  l = 0
1198  if (mid.eq.nid.and.if3d) then ! fill owning processor
1199  igr(kb) = igroup(e)
1200  do k=0,1
1201  do j=0,1
1202  do i=0,1
1203  l=l+1
1204  li=isym2pre(l)
1205  xyz(li,1,kb) = xm1(1+i*nxs,1+j*nys,1+k*nzs,e)
1206  xyz(li,2,kb) = ym1(1+i*nxs,1+j*nys,1+k*nzs,e)
1207  xyz(li,3,kb) = zm1(1+i*nxs,1+j*nys,1+k*nzs,e)
1208  enddo
1209  enddo
1210  enddo
1211  elseif (mid.eq.nid) then ! 2D
1212  igr(kb) = igroup(e)
1213  do j=0,1
1214  do i=0,1
1215  l =l+1
1216  li=isym2pre(l)
1217  xyz(li,1,kb) = xm1(1+i*nxs,1+j*nys,1,e)
1218  xyz(li,2,kb) = ym1(1+i*nxs,1+j*nys,1,e)
1219  enddo
1220  enddo
1221  endif
1222  enddo
1223  call gop(xyz,wk,'+ ',nblock) ! Sum across all processors
1224  call igop(igr,wk,'+ ',lblock) ! Sum across all processors
1225 
1226  if (nid.eq.0) then
1227  kb = 0
1228  do eg=eb,nemax
1229  kb = kb+1
1230 
1231  write(10,'(a15,i12,a2,i5,a1,a10,i6)')
1232  $ ' ELEMENT ',eg,' [',numapt,letapt,'] GROUP',igr(kb)
1233 
1234  if (if3d) then
1235 
1236  write(10,'(4g15.7)')(xyz(ic,1,kb),ic=1,4)
1237  write(10,'(4g15.7)')(xyz(ic,2,kb),ic=1,4)
1238  write(10,'(4g15.7)')(xyz(ic,3,kb),ic=1,4)
1239 
1240  write(10,'(4g15.7)')(xyz(ic,1,kb),ic=5,8)
1241  write(10,'(4g15.7)')(xyz(ic,2,kb),ic=5,8)
1242  write(10,'(4g15.7)')(xyz(ic,3,kb),ic=5,8)
1243 
1244  else ! 2D
1245 
1246  write(10,'(4g15.7)')(xyz(ic,1,kb),ic=1,4)
1247  write(10,'(4g15.7)')(xyz(ic,2,kb),ic=1,4)
1248 
1249  endif
1250 
1251  enddo
1252  endif
1253  enddo
1254 
1255  return
1256  end
1257 c-----------------------------------------------------------------------
1258  subroutine gen_rea_curve(imid)
1259 
1260 c This routine is complex because we must first count number of
1261 c nontrivial curved sides.
1262 
1263 c A two pass strategy is used: first count, then write
1264 
1265  include 'SIZE'
1266  include 'TOTAL'
1267 
1268  integer e,eb,eg
1269  character*1 cc
1270 
1271  parameter(lblock=500)
1272  common /scrns/ vcurve(5,12,lblock),wk(5*12*lblock)
1273  common /scruz/ icurve(12,lblock)
1274 
1275  if (imid.gt.0) then
1276 
1277 c imid = 0 ! No midside node defs
1278 c imid = 1 ! Midside defs where current curve sides don't exist
1279 c imid = 2 ! All nontrivial midside node defs
1280 
1281  if (imid.eq.2) call blank(ccurve,12*lelt)
1282 
1283  do e=1,nelt
1284  call gen_rea_midside_e(e)
1285  enddo
1286 
1287  endif
1288 
1289  nedge = 4 + 8*(ldim-2)
1290 
1291  ncurvn = 0
1292  do e=1,nelt
1293  do i=1,nedge
1294  if (ccurve(i,e).ne.' ') ncurvn = ncurvn+1
1295  enddo
1296  enddo
1297  ncurvn = iglsum(ncurvn,1)
1298 
1299  if (nid.eq.0) then
1300  WRITE(10,*)' ***** CURVED SIDE DATA *****'
1301  WRITE(10,'(I10,A20,A33)') ncurvn,' Curved sides follow',
1302  $ ' IEDGE,IEL,CURVE(I),I=1,5, CCURVE'
1303  endif
1304 
1305  do eb=1,nelgt,lblock
1306 
1307  nemax = min(eb+lblock-1,nelgt)
1308  call izero(icurve,12*lblock)
1309  call rzero(vcurve,60*lblock)
1310 
1311  kb = 0
1312  do eg=eb,nemax
1313  mid = gllnid(eg)
1314  e = gllel(eg)
1315  kb = kb+1
1316  if (mid.eq.nid) then ! fill owning processor
1317  do i=1,nedge
1318  icurve(i,kb) = 0
1319  if (ccurve(i,e).eq.'C') icurve(i,kb) = 1
1320  if (ccurve(i,e).eq.'s') icurve(i,kb) = 2
1321  if (ccurve(i,e).eq.'m') icurve(i,kb) = 3
1322  call copy(vcurve(1,i,kb),curve(1,i,e),5)
1323  enddo
1324  endif
1325  enddo
1326  call igop(icurve,wk,'+ ',12*lblock) ! Sum across all processors
1327  call gop(vcurve,wk,'+ ',60*lblock) ! Sum across all processors
1328 
1329  if (nid.eq.0) then
1330  kb = 0
1331  do eg=eb,nemax
1332  kb = kb+1
1333 
1334  do i=1,nedge
1335  ii = icurve(i,kb) ! equivalenced to s4
1336  if (ii.ne.0) then
1337  if (ii.eq.1) cc='C'
1338  if (ii.eq.2) cc='s'
1339  if (ii.eq.3) cc='m'
1340  if (nelgt.lt.1000) then
1341  write(10,'(i3,i3,5g14.6,1x,a1)') i,eg,
1342  $ (vcurve(k,i,kb),k=1,5),cc
1343  elseif (nelgt.lt.1000000) then
1344  write(10,'(i2,i6,5g14.6,1x,a1)') i,eg,
1345  $ (vcurve(k,i,kb),k=1,5),cc
1346  else
1347  write(10,'(i2,i12,5g14.6,1x,a1)') i,eg,
1348  $ (vcurve(k,i,kb),k=1,5),cc
1349  endif
1350  endif
1351  enddo
1352  enddo
1353  endif
1354 
1355  enddo
1356 
1357  return
1358  end
1359 c-----------------------------------------------------------------------
1360  subroutine gen_rea_bc (ifld)
1361 
1362  include 'SIZE'
1363  include 'TOTAL'
1364 
1365  integer e,eb,eg
1366 
1367  parameter(lblock=500)
1368  common /scrns/ vbc(5,6,lblock),wk(5*6*lblock)
1369  common /scruz/ ibc(6,lblock)
1370 
1371  character*1 s4(4)
1372  character*3 s3
1373  integer i4
1374  equivalence(i4,s4)
1375  equivalence(s3,s4)
1376 
1377  character*1 chtemp
1378  save chtemp
1379  data chtemp /' '/ ! For mesh bcs
1380 
1381  nface = 2*ldim
1382 
1383  nlg = nelg(ifld)
1384 
1385  if (ifld.eq.1.and..not.ifflow) then ! NO B.C.'s for this field
1386  if (nid.eq.0) write(10,*)
1387  $ ' ***** NO FLUID BOUNDARY CONDITIONS *****'
1388  return
1389  elseif (ifld.eq.1.and.nid.eq.0) then ! NO B.C.'s for this field
1390  write(10,*) ' ***** FLUID BOUNDARY CONDITIONS *****'
1391  elseif (ifld.ge.2.and.nid.eq.0) then ! NO B.C.'s for this field
1392  write(10,*) ' ***** THERMAL BOUNDARY CONDITIONS *****'
1393  endif
1394 
1395  do eb=1,nlg,lblock
1396  nemax = min(eb+lblock-1,nlg)
1397  call izero(ibc, 6*lblock)
1398  call rzero(vbc,30*lblock)
1399  kb = 0
1400  do eg=eb,nemax
1401  mid = gllnid(eg)
1402  e = gllel(eg)
1403  kb = kb+1
1404  if (mid.eq.nid) then ! fill owning processor
1405  do i=1,nface
1406  i4 = 0
1407  call chcopy(s4,cbc(i,e,ifld),3)
1408  ibc(i,kb) = i4
1409  call copy(vbc(1,i,kb),bc(1,i,e,ifld),5)
1410  enddo
1411  endif
1412  enddo
1413  call igop(ibc,wk,'+ ', 6*lblock) ! Sum across all processors
1414  call gop(vbc,wk,'+ ',30*lblock) ! Sum across all processors
1415 
1416  if (nid.eq.0) then
1417  kb = 0
1418  do eg=eb,nemax
1419  kb = kb+1
1420 
1421  do i=1,nface
1422  i4 = ibc(i,kb) ! equivalenced to s4
1423 
1424 c chtemp=' '
1425 c if (ifld.eq.1 .or. (ifld.eq.2 .and. .not. ifflow))
1426 c $ chtemp = cbc(i,kb,0)
1427 
1428  if (nlg.lt.1000) then
1429  write(10,'(a1,a3,2i3,5g14.6)')
1430  $ chtemp,s3,eg,i,(vbc(ii,i,kb),ii=1,5)
1431  elseif (nlg.lt.100000) then
1432  write(10,'(a1,a3,i5,i1,5g14.6)')
1433  $ chtemp,s3,eg,i,(vbc(ii,i,kb),ii=1,5)
1434  elseif (nlg.lt.1000000) then
1435  write(10,'(a1,a3,i6,5g14.6)')
1436  $ chtemp,s3,eg,(vbc(ii,i,kb),ii=1,5)
1437  else
1438  write(10,'(a1,a3,i12,5g18.11)')
1439  $ chtemp,s3,eg,(vbc(ii,i,kb),ii=1,5)
1440  endif
1441  enddo
1442  enddo
1443  endif
1444  enddo
1445 
1446  return
1447  end
1448 c-----------------------------------------------------------------------
1449  subroutine gen_rea_midside_e(e)
1450 
1451  include 'SIZE'
1452  include 'TOTAL'
1453 
1454  common /scrns/ x3(27),y3(27),z3(27),xyz(3,3)
1455  character*1 ccrve(12)
1456  integer e,edge
1457 
1458  integer e3(3,12)
1459  save e3
1460  data e3 / 1, 2, 3, 3, 6, 9, 9, 8, 7, 7, 4, 1
1461  $ , 19,20,21, 21,24,27, 27,26,25, 25,22,19
1462  $ , 1,10,19, 3,12,21, 9,18,27, 7,16,25 /
1463 
1464  real len
1465 
1466  call chcopy(ccrve,ccurve(1,e),12)
1467 
1468  call map2reg(x3,3,xm1(1,1,1,e),1) ! Map to 3x3x3 array
1469  call map2reg(y3,3,ym1(1,1,1,e),1)
1470  if (if3d) call map2reg(z3,3,zm1(1,1,1,e),1)
1471 
1472 c Take care of spherical curved face defn
1473  if (ccurve(5,e).eq.'s') then
1474  call chcopy(ccrve(1),'ssss',4) ! face 5
1475  call chcopy(ccrve(5),' ',1) ! face 5
1476  endif
1477  if (ccurve(6,e).eq.'s') then
1478  call chcopy(ccrve(5),'ssss',4) ! face 6
1479  endif
1480 
1481  tol = 1.e-4
1482  tol2 = tol**2
1483  nedge = 4 + 8*(ldim-2)
1484 
1485  do i=1,nedge
1486  if (ccrve(i).eq.' ') then
1487  do j=1,3
1488  xyz(1,j)=x3(e3(j,i))
1489  xyz(2,j)=y3(e3(j,i))
1490  xyz(3,j)=z3(e3(j,i))
1491  enddo
1492  len = 0.
1493  h = 0.
1494  do j=1,ldim
1495  xmid = .5*(xyz(j,1)+xyz(j,3))
1496  h = h + (xyz(j,2)-xmid)**2
1497  len = len + (xyz(j,3)-xyz(j,1))**2
1498  enddo
1499  if (h.gt.tol2*len) ccurve(i,e) = 'm'
1500  if (h.gt.tol2*len) call copy(curve(1,i,e),xyz(1,2),ldim)
1501  endif
1502  enddo
1503 
1504  return
1505  end
1506 c-----------------------------------------------------------------------
1507  subroutine hpts
1508 c
1509 c evaluate velocity, temperature, pressure and ps-scalars
1510 c for list of points and dump results
1511 c note: read/write on rank0 only
1512 c
1513 c ASSUMING LHIS IS MAX NUMBER OF POINTS TO READ IN ON ONE PROCESSOR
1514 
1515  include 'SIZE'
1516  include 'TOTAL'
1517 
1518  parameter(nfldm=ldim+ldimt+1)
1519 
1520  real pts, fieldout, dist, rst
1521  common /c_hptsr/ pts(ldim,lhis)
1522  $ , fieldout(nfldm,lhis)
1523  $ , dist(lhis)
1524  $ , rst(lhis*ldim)
1525 
1526  common /nekmpi/ nidd,npp,nekcomm,nekgroup,nekreal
1527 
1528  integer rcode, elid, proc
1529  common /c_hptsi/ rcode(lhis),elid(lhis),proc(lhis)
1530 
1531  common /scrcg/ pm1(lx1,ly1,lz1,lelv) ! mapped pressure
1532  common /outtmp/ wrk(lx1*ly1*lz1*lelt,nfldm)
1533 
1534 
1535  logical iffind
1536 
1537  integer icalld,npoints,npts
1538  save icalld,npoints,npts
1539  data icalld /0/
1540  data npoints /0/
1541 
1542  save inth_hpts
1543 
1544  nxyz = lx1*ly1*lz1
1545  ntot = nxyz*nelt
1546  nbuff = lhis ! point to be read in on 1 proc.
1547 
1548  toldist = 5e-6
1549 
1550  if(nio.eq.0) write(6,*) 'dump history points'
1551 
1552  if(icalld.eq.0) then
1553  npts = lhis ! number of points per proc
1554  call hpts_in(pts,npts,npoints)
1555 
1556  tol = 5e-13
1557  n = lx1*ly1*lz1*lelt
1558  npt_max = 128
1559  nxf = 2*lx1 ! fine mesh for bb-test
1560  nyf = 2*ly1
1561  nzf = 2*lz1
1562  bb_t = 0.01 ! relative size to expand bounding boxes by
1563  call fgslib_findpts_setup(inth_hpts,nekcomm,np,ldim,
1564  & xm1,ym1,zm1,lx1,ly1,lz1,
1565  & nelt,nxf,nyf,nzf,bb_t,n,n,
1566  & npt_max,tol)
1567  endif
1568 
1569 
1570  call prepost_map(0) ! maps axisymm and pressure
1571 
1572  ! pack working array
1573  nflds = 0
1574  if(ifvo) then
1575  call copy(wrk(1,1),vx,ntot)
1576  call copy(wrk(1,2),vy,ntot)
1577  if(if3d) call copy(wrk(1,3),vz,ntot)
1578  nflds = ldim
1579  endif
1580  if(ifpo) then
1581  nflds = nflds + 1
1582  call copy(wrk(1,nflds),pm1,ntot)
1583  endif
1584  if(ifto) then
1585  nflds = nflds + 1
1586  call copy(wrk(1,nflds),t,ntot)
1587  endif
1588  do i = 1,ldimt
1589  if(ifpsco(i)) then
1590  nflds = nflds + 1
1591  call copy(wrk(1,nflds),t(1,1,1,1,i+1),ntot)
1592  endif
1593  enddo
1594 
1595  ! interpolate
1596  if(icalld.eq.0) then
1597  call fgslib_findpts(inth_hpts,rcode,1,
1598  & proc,1,
1599  & elid,1,
1600  & rst,ldim,
1601  & dist,1,
1602  & pts(1,1),ldim,
1603  & pts(2,1),ldim,
1604  & pts(3,1),ldim,npts)
1605 
1606  nfail = 0
1607  do i=1,npts
1608  ! check return code
1609  if(rcode(i).eq.1) then
1610  if(sqrt(dist(i)).gt.toldist) then
1611  nfail = nfail + 1
1612  IF (nfail.LE.5) WRITE(6,'(a,1p4e15.7)')
1613  & ' WARNING: point on boundary or outside the mesh xy[z]d^2:'
1614  & ,(pts(k,i),k=1,ldim),dist(i)
1615  endif
1616  elseif(rcode(i).eq.2) then
1617  nfail = nfail + 1
1618  if (nfail.le.5) write(6,'(a,1p3e15.7)')
1619  & ' WARNING: point not within mesh xy[z]: !',
1620  & (pts(k,i),k=1,ldim)
1621  endif
1622  enddo
1623  icalld = 1
1624  endif
1625 
1626 
1627  ! evaluate input field at given points
1628  do ifld = 1,nflds
1629  call fgslib_findpts_eval(inth_hpts,fieldout(ifld,1),nfldm,
1630  & rcode,1,
1631  & proc,1,
1632  & elid,1,
1633  & rst,ldim,npts,
1634  & wrk(1,ifld))
1635  enddo
1636  ! write interpolation results to hpts.out
1637  call hpts_out(fieldout,nflds,nfldm,npoints,nbuff)
1638 
1639  call prepost_map(1) ! maps back axisymm arrays
1640 
1641  return
1642  end
1643 c-----------------------------------------------------------------------
1644  subroutine buffer_in(buffer,npp,npoints,nbuf)
1645 
1646  include 'SIZE'
1647  include 'INPUT'
1648  include 'PARALLEL'
1649 
1650  real buffer(ldim,nbuf)
1651 
1652  ierr = 0
1653  if(nid.eq.0) then
1654  write(6,*) 'reading history points'
1655  open(50,file=hisfle,status='old',err=100)
1656  read(50,*,err=100) npoints
1657  goto 101
1658  100 ierr = 1
1659  101 continue
1660  endif
1661  ierr=iglsum(ierr,1)
1662  if(ierr.gt.0) then
1663  if(nio.eq.0)
1664  & write(6,*) 'Cannot open history file in subroutine hpts()'
1665  call exitt
1666  endif
1667 
1668  call bcast(npoints,isize)
1669  if(npoints.gt.(lhis-1)*np) then
1670  if(nid.eq.0) write(6,*) 'ABORT: Increase lhis in SIZE!'
1671  call exitt
1672  endif
1673  if(nid.eq.0) write(6,*) 'found ', npoints, ' points'
1674 
1675 
1676  npass = npoints/nbuf +1 !number of passes to cover all pts
1677  n0 = mod(npoints,nbuf)!remainder
1678  if(n0.eq.0) then
1679  npass = npass-1
1680  n0 = nbuf
1681  endif
1682 
1683  len = wdsize*ldim*nbuf
1684  if (nid.gt.0.and.nid.lt.npass) msg_id=irecv(nid,buffer,len)
1685  call nekgsync
1686 
1687  npp=0
1688  if(nid.eq.0) then
1689  i1 = nbuf
1690  do ipass = 1,npass
1691  if(ipass.eq.npass) i1 = n0
1692  do i = 1,i1
1693  read(50,*) (buffer(j,i),j=1,ldim)
1694  enddo
1695  if(ipass.lt.npass)call csend(ipass,buffer,len,ipass,0)
1696  enddo
1697  npp = n0
1698  elseif (nid.lt.npass) then !processors receiving data
1699  call msgwait(msg_id)
1700  npp=nbuf
1701  endif
1702 
1703  return
1704  end
1705 c-----------------------------------------------------------------------
1706  subroutine hpts_in(pts,npts,npoints)
1707 c npts=local count; npoints=total count
1708 
1709  include 'SIZE'
1710  include 'PARALLEL'
1711 
1712  parameter(lt2=2*lx1*ly1*lz1*lelt)
1713  common /scrns/ xyz(ldim,lt2)
1714  common /scruz/ mid(lt2) ! Target proc id
1715  real pts(ldim,npts)
1716 
1717  if (lt2.gt.npts) then
1718 
1719  call buffer_in(xyz,npp,npoints,lt2)
1720  if(npoints.gt.np*npts) then
1721  if(nid.eq.0)write(6,*)'ABORT in hpts(): npoints > NP*lhis!!'
1722  if(nid.eq.0)write(6,*)'Change SIZE: ',np,npts,npoints
1723  call exitt
1724  endif
1725 
1726  npmax = (npoints/npts)
1727  if(mod(npoints,npts).eq.0) npmax=npmax+1
1728 
1729  if(nid.gt.0.and.npp.gt.0) then
1730  npts_b = lt2*(nid-1) ! # pts offset(w/o 0)
1731  nprc_b = npts_b/npts ! # proc offset(w/o 0)
1732 
1733  istart = mod(npts_b,npts) ! istart-->npts pts left
1734  ip = nprc_b + 1 ! PID offset
1735  icount = istart ! point offset
1736  elseif(nid.eq.0) then
1737  npts0 = mod1(npoints,lt2) ! Node 0 pts
1738  npts_b = npoints - npts0 ! # pts before Node 0
1739  nprc_b = npts_b/npts
1740 
1741  istart = mod(npts_b,npts)
1742  ip = nprc_b + 1
1743  icount = istart
1744  endif
1745 
1746  do i =1,npp
1747  icount = icount + 1
1748  if(ip.gt.npmax) ip = 0
1749  mid(i) = ip
1750  if (icount.eq.npts) then
1751  ip = ip+1
1752  icount = 0
1753  endif
1754  enddo
1755 
1756  call fgslib_crystal_tuple_transfer
1757  & (cr_h,npp,lt2,mid,1,pts,0,xyz,ldim,1)
1758 
1759  call copy(pts,xyz,ldim*npp)
1760  else
1761  call buffer_in(pts,npp,npoints,npts)
1762  endif
1763  npts = npp
1764 
1765 
1766  return
1767  end
1768 c-----------------------------------------------------------------------
1769  subroutine hpts_out(fieldout,nflds,nfldm,npoints,nbuff)
1770 
1771  include 'SIZE'
1772  include 'TOTAL'
1773 
1774  real buf(nfldm,nbuff),fieldout(nfldm,nbuff)
1775 
1776  len = wdsize*nfldm*nbuff
1777 
1778 
1779  npass = npoints/nbuff + 1
1780  il = mod(npoints,nbuff)
1781  if(il.eq.0) then
1782  il = nbuff
1783  npass = npass-1
1784  endif
1785 
1786  do ipass = 1,npass
1787 
1788  call nekgsync
1789 
1790  if(ipass.lt.npass) then
1791  if(nid.eq.0) then
1792  call crecv(ipass,buf,len)
1793  do ip = 1,nbuff
1794  write(50,'(1p20E15.7)') time,
1795  & (buf(i,ip), i=1,nflds)
1796  enddo
1797  elseif(nid.eq.ipass) then
1798  call csend(ipass,fieldout,len,0,nid)
1799  endif
1800 
1801  else !ipass.eq.npass
1802 
1803  if(nid.eq.0) then
1804  do ip = 1,il
1805  write(50,'(1p20E15.7)') time,
1806  & (fieldout(i,ip), i=1,nflds)
1807  enddo
1808  endif
1809 
1810  endif
1811  enddo
1812 
1813  return
1814  end
1815 c-----------------------------------------------------------------------
1816  subroutine gen_rea_full(imid) ! Generate and output essential parts of .rea
1817  ! Clobbers ccurve()
1818  include 'SIZE'
1819  include 'TOTAL'
1820 
1821 c imid = 0 ! No midside node defs
1822 c imid = 1 ! Midside defs where current curve sides don't exist
1823 c imid = 2 ! All nontrivial midside node defs
1824 
1825  if (nid.eq.0) open(unit=10,file='newrea.rea',status='unknown') ! clobbers existing file
1826 
1827  call gen_rea_top
1828 
1829  call gen_rea_xyz
1830 
1831  call gen_rea_curve(imid) ! Clobbers ccurve()
1832 
1833  if (nid.eq.0) write(10,*)' ***** BOUNDARY CONDITIONS *****'
1834  do ifld=1,nfield
1835  call gen_rea_bc (ifld)
1836  enddo
1837 
1838  call gen_rea_bottom
1839 
1840  if (nid.eq.0) close(10)
1841 
1842  return
1843  end
1844 c-----------------------------------------------------------------------
1845  subroutine gen_rea_top
1846 
1847  include 'SIZE'
1848  include 'INPUT'
1849  include 'PARALLEL'
1850  include 'CTIMER'
1851 
1852  logical ifbswap,ifre2,parfound
1853  character*132 string
1854  character*72 string2
1855  integer idum(3*numsts+3)
1856  integer paramval
1857 
1858  ierr = 0
1859  call flush_io
1860 
1861  if(nid.eq.0) then
1862  write(6,'(A,A)') ' Reading ', reafle
1863  open (unit=9,file=reafle,status='old', iostat=ierr)
1864  endif
1865 
1866  call bcast(ierr,isize)
1867  if (ierr .gt. 0) call exitti('Cannot open .rea file!$',1)
1868 
1869 
1870  IF(nid.EQ.0) THEN
1871  READ(9,'(a)') string2
1872  write(10,'(a)') string2
1873  READ(9,'(a)') string2
1874  write(10,*) string2
1875  READ(9,'(a)') string2
1876  write(10,*) string2
1877  READ(9,'(a)') string2
1878  write(10,*) string2
1879 
1880  READ(string2,*) paramval
1881 
1882 c WRITE PARAMETERS
1883  DO 20 i=1,paramval
1884  READ(9,'(a)') string2
1885  write(10,*) string2
1886  20 CONTINUE
1887 
1888 c LINES OF PASSIVE SCALARS
1889  read(9,'(a)') string2
1890  write(10,*) string2
1891 
1892  read(string2,*) paramval
1893  if (paramval.gt.0) then
1894  do i=1,paramval
1895  READ(9,'(a)') string2
1896  write(10,*) string2
1897  enddo
1898  endif
1899 
1900 c LINES OF LOGICAL SWITCHES
1901  read(9,'(a)') string2
1902  write(10,*) string2
1903 
1904  read(string2,*) paramval
1905  if (paramval.gt.0) then
1906  do i=1,paramval
1907  READ(9,'(a)') string2
1908  write(10,*) string2
1909  enddo
1910  endif
1911 
1912 c LAST TWO LINES BEFORE ELEMENT DATA BEGINS
1913  do i=1,2
1914  READ(9,'(a)') string2
1915  write(10,*) string2
1916  enddo
1917 
1918 
1919  ENDIF
1920 
1921 
1922  return
1923  end
1924 c-----------------------------------------------------------------------
1925  subroutine gen_rea_bottom
1926 
1927  include 'SIZE'
1928  include 'INPUT'
1929  include 'PARALLEL'
1930  include 'CTIMER'
1931 
1932  logical ifbswap,ifre2,parfound
1933  character*132 string
1934  character*72 string2
1935  integer idum(3*numsts+3)
1936  integer paramval,i,j
1937  integer n1,n2,n3
1938 
1939 c IGNORE ELEMENT DATA
1940 c IGNORE XY DATA
1941 
1942  IF(nid.EQ.0) THEN
1943  READ(9,'(a)') string2
1944  READ(string2,*) n1,n2,n3
1945  if (n1.lt.0) goto 1001
1946  do i=1,nelgt
1947  READ(9,'(a)') string2
1948  do j=1,2+(ldim-2)*4
1949  READ(9,'(a)') string2
1950  enddo
1951  enddo
1952 c CURVE SIDE DATA
1953  READ(9,'(a)') string2
1954  READ(9,*) paramval
1955  if (paramval.gt.0) then
1956  do i=1,paramval
1957  READ(9,'(a)') string2
1958  enddo
1959  endif
1960 c BOUNDARY CONDITIONS
1961  READ(9,'(a)') string2
1962 c FLUID
1963  READ(9,'(a)') string2
1964  if (ifflow) then
1965  do i=1,nelgv*2*ldim
1966  READ(9,'(a)') string2
1967  enddo
1968  endif
1969 c Thermal
1970  READ(9,'(a)') string2
1971  if (ifheat) then
1972  do i=1,nelgt*2*ldim
1973  READ(9,'(a)') string2
1974  enddo
1975  else
1976  write(10,*) string2
1977  endif
1978 
1979  1001 continue
1980 
1981 c PRESOLVE
1982  READ(9,'(a)') string2
1983  write(10,*) string2
1984 c INITIAL CONDITIONS
1985  READ(9,'(a)') string2
1986  write(10,*) string2
1987  read(string2,*) paramval
1988  if (paramval.gt.0) then
1989  do i=1,paramval
1990  READ(9,'(a)') string2
1991  write(10,*) string2
1992  enddo
1993  endif
1994 c DRIVE FORCE
1995  READ(9,'(a)') string2
1996  write(10,*) string2
1997  READ(9,'(a)') string2
1998  write(10,*) string2
1999  read(string2,*) paramval
2000 
2001  if (paramval.gt.0) then
2002  do i=1,paramval
2003  READ(9,'(a)') string2
2004  write(10,*) string2
2005  enddo
2006  endif
2007 
2008 c VARIABLE PROPERTY DATA
2009  read(9,'(a)') string2
2010  write(10,*) string2
2011  read(9,'(a)') string2
2012  write(10,*) string2
2013  read(string2,*) paramval
2014  if (paramval.gt.0) then
2015  do i=1,paramval
2016  READ(9,'(a)') string2
2017  write(10,*) string2
2018  enddo
2019  endif
2020 
2021 c HISTORY AND INTEGRAL DATA
2022  read(9,'(a)') string2
2023  write(10,*) string2
2024  read(9,'(a)') string2
2025  write(10,*) string2
2026 
2027  read(string2,*) paramval
2028  if (paramval.gt.0) then
2029  do i=1,paramval
2030  READ(9,'(a)') string2
2031  write(10,*) string2
2032  enddo
2033  endif
2034 
2035 c OUPUT FIELD SPECIFICATION
2036  read(9,'(a)') string2
2037  write(10,*) string2
2038 
2039  read(9,'(a)') string2
2040  write(10,*) string2
2041  read(string2,*) paramval
2042  if (paramval.gt.0) then
2043  do i=1,paramval
2044  READ(9,'(a)') string2
2045  write(10,*) string2
2046  enddo
2047  endif
2048 
2049 c LAST FOUR LINES
2050  do i=1,5
2051  READ(9,'(a)') string2
2052  write(10,*) string2
2053  enddo
2054 
2055 
2056  ENDIF
2057 
2058  if (nid.eq.0) close(9)
2059 
2060  return
2061  end
2062 c-----------------------------------------------------------------------
byte_open
#define byte_open
Definition: byte.c:35
exitt
void exitt()
Definition: comm_mpi.f:604
byte_write
#define byte_write
Definition: byte.c:39
byte_close
#define byte_close
Definition: byte.c:36
crecv
subroutine crecv(mtype, buf, lenm)
Definition: comm_mpi.f:313
igop
subroutine igop(x, w, op, n)
Definition: comm_mpi.f:247
msgwait
subroutine msgwait(imsg)
Definition: comm_mpi.f:489
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
irecv
function irecv(msgtag, x, len)
Definition: comm_mpi.f:471
err_chk
subroutine err_chk(ierr, string)
Definition: comm_mpi.f:558
gllel
integer function gllel(ieg)
Definition: dprocmap.f:183
gllnid
integer function gllnid(ieg)
Definition: dprocmap.f:161
fd_weights_full
subroutine fd_weights_full(xx, x, n, m, c)
Definition: fast3d.f:1293
setics
subroutine setics
Definition: ic.f:3
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
iglsum
function iglsum(a, n)
Definition: math.f:926
copy
subroutine copy(a, b, n)
Definition: math.f:260
blank
subroutine blank(A, N)
Definition: math.f:19
iglmax
function iglmax(a, n)
Definition: math.f:913
vsqrt
subroutine vsqrt(A, N)
Definition: math.f:41
izero
subroutine izero(a, n)
Definition: math.f:215
ltrunc
function ltrunc(string, l)
Definition: math.f:494
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
sort
subroutine sort(a, ind, n)
Definition: math.f:1325
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
build_new_filter
subroutine build_new_filter(intv, zpts, nx, kut, wght, nio)
Definition: navier5.f:801
gaujordf
subroutine gaujordf(a, m, n, indr, indc, ipiv, ierr, rmult)
Definition: navier5.f:668
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
legendre_poly
subroutine legendre_poly(L, x, N)
Definition: navier5.f:785
filterq
subroutine filterq(v, f, nx, nz, w1, w2, ft, if3d, dmax)
Definition: navier5.f:167
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
hpts
subroutine hpts
Definition: postpro.f:1508
quadratic
subroutine quadratic(x1, x2, a, b, c, ierr)
Definition: postpro.f:146
gen_rea_top
subroutine gen_rea_top
Definition: postpro.f:1846
hpts_in
subroutine hpts_in(pts, npts, npoints)
Definition: postpro.f:1707
comp_gije
subroutine comp_gije(gije, u, v, w, e)
Definition: postpro.f:253
tens3d1
subroutine tens3d1(v, u, f, ft, nv, nu)
Definition: postpro.f:389
buffer_in
subroutine buffer_in(buffer, npp, npoints, nbuf)
Definition: postpro.f:1645
build_1d_filt
subroutine build_1d_filt(fh, fht, trnsfr, nx, nio)
Definition: postpro.f:426
gen_int_gz
subroutine gen_int_gz(j, jt, g, n, z, m)
Definition: postpro.f:634
gen_rea_full
subroutine gen_rea_full(imid)
Definition: postpro.f:1817
lambda2
subroutine lambda2(l2)
Definition: postpro.f:20
gen_rea_bottom
subroutine gen_rea_bottom
Definition: postpro.f:1926
cubic
subroutine cubic(xo, ai1, ai2, ai3, ierr)
Definition: postpro.f:194
gen_rea_bc
subroutine gen_rea_bc(ifld)
Definition: postpro.f:1361
gen_rea_curve
subroutine gen_rea_curve(imid)
Definition: postpro.f:1259
filter_s1
subroutine filter_s1(scalar, tf, nx, nel)
Definition: postpro.f:321
gen_rea_midside_e
subroutine gen_rea_midside_e(e)
Definition: postpro.f:1450
gen_re2_xyz
subroutine gen_re2_xyz
Definition: postpro.f:718
gen_re2
subroutine gen_re2(imid)
Definition: postpro.f:671
find_lam3
subroutine find_lam3(lam, aa, w, ldim, ierr)
Definition: postpro.f:79
gen_re2_bc
subroutine gen_re2_bc(ifld)
Definition: postpro.f:1019
map2reg_3di_e
subroutine map2reg_3di_e(uf, n, uc, m)
Definition: postpro.f:594
hpts_out
subroutine hpts_out(fieldout, nflds, nfldm, npoints, nbuff)
Definition: postpro.f:1770
map2reg_2di_e
subroutine map2reg_2di_e(uf, n, uc, m)
Definition: postpro.f:565
filter_s0
subroutine filter_s0(scalar, wght, ncut, name5)
Definition: postpro.f:343
comp_sije
subroutine comp_sije(gije)
Definition: postpro.f:518
load_fld
subroutine load_fld(string)
Definition: postpro.f:2
gen_rea_xyz
subroutine gen_rea_xyz
Definition: postpro.f:1163
gen_re2_curve
subroutine gen_re2_curve(imid)
Definition: postpro.f:898
map2reg
subroutine map2reg(ur, n, u, nel)
Definition: postpro.f:544
gen_rea
subroutine gen_rea(imid)
Definition: postpro.f:1138
zuni
subroutine zuni(z, np)
Definition: postpro.f:655
mag_tensor_e
subroutine mag_tensor_e(mag, aije)
Definition: postpro.f:493
prepost_map
subroutine prepost_map(isave)
Definition: prepost.f:136
copyx4
subroutine copyx4(a, b, n)
Definition: prepost.f:560
zwgll
subroutine zwgll(Z, W, NP)
Definition: speclib.f:108
flush_io
subroutine flush_io
Definition: subs1.f:1561