KTH framework for Nek5000 toolboxes; testing version  0.0.1
navier1.f
Go to the documentation of this file.
1  subroutine plan1 (igeom)
2 C-------------------------------------------------------------------------
3 C
4 C Compute pressure and velocity using consistent approximation spaces.
5 C
6 C-------------------------------------------------------------------------
7  include 'SIZE'
8  include 'INPUT'
9  include 'EIGEN'
10  include 'SOLN'
11  include 'TSTEP'
12 C
13  COMMON /scrhi/ h2inv(lx1,ly1,lz1,lelv)
14  COMMON /scrns/ resv1(lx1,ly1,lz1,lelv)
15  $ , resv2(lx1,ly1,lz1,lelv)
16  $ , resv3(lx1,ly1,lz1,lelv)
17  $ , dv1(lx1,ly1,lz1,lelv)
18  $ , dv2(lx1,ly1,lz1,lelv)
19  $ , dv3(lx1,ly1,lz1,lelv)
20  $ , wp(lx2,ly2,lz2,lelv)
21  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
22  $ , h2(lx1,ly1,lz1,lelv)
23  REAL G1 (LX1,LY1,LZ1,LELV)
24  REAL G2 (LX1,LY1,LZ1,LELV)
25  REAL G3 (LX1,LY1,LZ1,LELV)
26  equivalence(g1,resv1), (g2,resv2), (g3,resv3)
27  LOGICAL IFSTUZ
28 C
29  IF (igeom.EQ.1) THEN
30 C
31 C Old geometry
32 C
33  CALL makef
34  CALL lagvel
35 C
36  ELSE
37 C
38 C New geometry
39 C
40  CALL bcdirvc (vx,vy,vz,v1mask,v2mask,v3mask)
41  IF (ifstrs) CALL bcneutr
42 C
43 C Check if steady state
44 C
45  ifstuz = .false.
46  CALL convuz (ifstuz)
47 C... no steady state
48  ifstuz = .false.
49  IF (ifstuz) THEN
50  IF (nio.EQ.0) WRITE (6,*)
51  $ 'Steady state reached in the fluid solver'
52  return
53  ENDIF
54 C
55 C Uzawa decoupling: First, compute pressure.....
56 C
57  ntot1 = lx1*ly1*lz1*nelv
58  ntot2 = lx2*ly2*lz2*nelv
59  intype = 0
60  if (iftran) intype = -1
61  call sethlm (h1,h2,intype)
62  if (iftran) call invers2 (h2inv,h2,ntot1)
63  call makeg ( g1,g2,g3,h1,h2,intype)
64  call crespuz (wp,g1,g2,g3,h1,h2,h2inv,intype)
65  call uzawa (wp,h1,h2,h2inv,intype,icg)
66  if (icg.gt.0) call add2 (pr,wp,ntot2)
67 
68 C .... then, compute velocity:
69 
70  call cresvuz (resv1,resv2,resv3)
71  call ophinv (dv1,dv2,dv3,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
72  call opadd2 (vx,vy,vz,dv1,dv2,dv3)
73 
74  endif
75 
76  return
77  end
78 c-----------------------------------------------------------------------
79  subroutine crespuz (respr,g1,g2,g3,h1,h2,h2inv,intype)
80 
81 c Compute start-residual/right-hand-side in the pressure equation
82 
83  include 'SIZE'
84  include 'TOTAL'
85  REAL RESPR (LX2,LY2,LZ2,LELV)
86  REAL G1 (LX1,LY1,LZ1,LELV)
87  REAL G2 (LX1,LY1,LZ1,LELV)
88  REAL G3 (LX1,LY1,LZ1,LELV)
89  REAL H1 (LX1,LY1,LZ1,LELV)
90  REAL H2 (LX1,LY1,LZ1,LELV)
91  REAL H2INV (LX1,LY1,LZ1,LELV)
92  COMMON /scruz/ ta1(lx1,ly1,lz1,lelv)
93  $ , ta2(lx1,ly1,lz1,lelv)
94  $ , ta3(lx1,ly1,lz1,lelv)
95  COMMON /scrmg/ vbdry1(lx1,ly1,lz1,lelv)
96  $ , vbdry2(lx1,ly1,lz1,lelv)
97  $ , vbdry3(lx1,ly1,lz1,lelv)
98 
99  if ((intype.eq.0).or.(intype.eq.-1)) then
100  call ophinv (ta1,ta2,ta3,g1,g2,g3,h1,h2,tolhr,nmxp)
101  else
102  call opbinv (ta1,ta2,ta3,g1,g2,g3,h2inv)
103  endif
104  call opamask (vbdry1,vbdry2,vbdry3)
105  call opsub2 (ta1,ta2,ta3,vbdry1,vbdry2,vbdry3)
106  call opdiv (respr,ta1,ta2,ta3)
107  call ortho (respr)
108 
109  return
110  end
111 c-----------------------------------------------------------------------
112  subroutine cresvuz (resv1,resv2,resv3)
113 
114 c Compute the residual for the velocity - UZAWA SCHEME.
115 
116  include 'SIZE'
117  include 'GEOM'
118  include 'SOLN'
119  REAL RESV1 (LX1,LY1,LZ1,1)
120  REAL RESV2 (LX1,LY1,LZ1,1)
121  REAL RESV3 (LX1,LY1,LZ1,1)
122  COMMON /scrmg/ ta1(lx1,ly1,lz1,lelv)
123  $ , ta2(lx1,ly1,lz1,lelv)
124  $ , ta3(lx1,ly1,lz1,lelv)
125  COMMON /screv/ h1(lx1,ly1,lz1,lelv)
126  $ , h2(lx1,ly1,lz1,lelv)
127 C
128  inloc = -1
129  CALL sethlm (h1,h2,inloc)
130  CALL oprzero (resv1,resv2,resv3)
131  CALL opgradt (resv1,resv2,resv3,pr)
132  CALL opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
133  CALL ophx (ta1,ta2,ta3,vx,vy,vz,h1,h2)
134  CALL opsub2 (resv1,resv2,resv3,ta1,ta2,ta3)
135 C
136  return
137  END
138 C
139  subroutine makeg (out1,out2,out3,h1,h2,intype)
140 C----------------------------------------------------------------------
141 C
142 C Compute inhomogeneities for the elliptic solver in the pressure
143 C residual evaluation.
144 C INTYPE = 0 steady state
145 C INTYPE = 1 explicit, Euler forward
146 C INTYPE = -1 implicit, Euler backward
147 C
148 C-----------------------------------------------------------------------
149  include 'SIZE'
150  include 'TOTAL'
151  REAL OUT1 (LX1,LY1,LZ1,LELV)
152  REAL OUT2 (LX1,LY1,LZ1,LELV)
153  REAL OUT3 (LX1,LY1,LZ1,LELV)
154  REAL H1 (LX1,LY1,LZ1,LELV)
155  REAL H2 (LX1,LY1,LZ1,LELV)
156  COMMON /scrmg/ ta1(lx1,ly1,lz1,lelv)
157  $ ,ta2(lx1,ly1,lz1,lelv)
158  $ ,ta3(lx1,ly1,lz1,lelv)
159  COMMON /scruz/ tb1(lx1,ly1,lz1,lelv)
160  $ ,tb2(lx1,ly1,lz1,lelv)
161  $ ,tb3(lx1,ly1,lz1,lelv)
162  $ ,hzero(lx1,ly1,lz1,lelv)
163 C
164  ntot1 = lx1*ly1*lz1*nelv
165 C
166  CALL opgradt (out1,out2,out3,pr)
167  CALL opchsgn (out1,out2,out3)
168 C
169  IF ((intype.EQ.0.).OR.(intype.EQ.-1)) THEN
170 C
171 C Steady state or implicit scheme
172 C
173  CALL opamask (tb1,tb2,tb3)
174  CALL ophx (ta1,ta2,ta3,tb1,tb2,tb3,h1,h2)
175  CALL opadd2 (out1,out2,out3,ta1,ta2,ta3)
176  CALL opsub2 (out1,out2,out3,bfx,bfy,bfz)
177 C
178  ELSEIF (intype.EQ.1) THEN
179 C
180 C Explicit scheme
181 C
182  CALL rzero (hzero,ntot1)
183  CALL ophx (ta1,ta2,ta3,vx,vy,vz,h1,hzero)
184  CALL opadd2 (out1,out2,out3,ta1,ta2,ta3)
185  CALL opsub2 (out1,out2,out3,bfx,bfy,bfz)
186 C
187  ENDIF
188 C
189  return
190  END
191 
192 c-----------------------------------------------------------------------
193  subroutine ctolspl (tolspl,respr)
194 C
195 C Compute the pressure tolerance
196 C
197  include 'SIZE'
198  include 'MASS'
199  include 'TSTEP'
200  REAL RESPR (LX2,LY2,LZ2,LELV)
201  COMMON /scrmg/ work(lx1,ly1,lz1,lelv)
202 C
203  ntot1 = lx1*ly1*lz1*nelv
204  CALL invcol3 (work,respr,bm1,ntot1)
205  CALL col2 (work,respr,ntot1)
206  rinit = sqrt(glsum(work,ntot1)/volvm1)
207  IF (tolpdf.GT.0.) THEN
208  tolspl = tolpdf
209  tolmin = tolpdf
210  ELSE
211  tolspl = tolps/dt
212  tolmin = rinit*prelax
213  ENDIF
214  IF (tolspl.LT.tolmin) THEN
215  tolold = tolspl
216  tolspl = tolmin
217  IF (nio.EQ.0)
218  $ WRITE (6,*) 'Relax the pressure tolerance ',tolspl,tolold
219  ENDIF
220  return
221  end
222 c------------------------------------------------------------------------
223  subroutine ortho (respr)
224 
225 C Orthogonalize the residual in the pressure solver with respect
226 C to (1,1,...,1)T (only if all Dirichlet b.c.).
227 
228  include 'SIZE'
229  include 'GEOM'
230  include 'INPUT'
231  include 'PARALLEL'
232  include 'SOLN'
233  include 'TSTEP'
234  real respr (lx2,ly2,lz2,lelv)
235  integer*8 ntotg,nxyz2
236 
237  nxyz2 = lx2*ly2*lz2
238  ntot = nxyz2*nelv
239  ntotg = nxyz2*nelgv
240 
241  if (ifield.eq.1) then
242  if (ifvcor) then
243  rlam = glsum(respr,ntot)/ntotg
244  call cadd (respr,-rlam,ntot)
245  endif
246  elseif (ifield.eq.ifldmhd) then
247  if (ifbcor) then
248  rlam = glsum(respr,ntot)/ntotg
249  call cadd (respr,-rlam,ntot)
250  endif
251  else
252  call exitti('ortho: unaccounted ifield = $',ifield)
253  endif
254 
255  return
256  end
257 c------------------------------------------------------------------------
258  subroutine cdabdtp (ap,wp,h1,h2,h2inv,intype)
259 
260 C INTYPE= 0 Compute the matrix-vector product DA(-1)DT*p
261 C INTYPE= 1 Compute the matrix-vector product D(B/DT)(-1)DT*p
262 C INTYPE=-1 Compute the matrix-vector product D(A+B/DT)(-1)DT*p
263 
264  include 'SIZE'
265  include 'TOTAL'
266  REAL AP (LX2,LY2,LZ2,1)
267  REAL WP (LX2,LY2,LZ2,1)
268  REAL H1 (LX1,LY1,LZ1,1)
269  REAL H2 (LX1,LY1,LZ1,1)
270  REAL H2INV (LX1,LY1,LZ1,1)
271 C
272  COMMON /scrns/ ta1(lx1,ly1,lz1,lelv)
273  $ , ta2(lx1,ly1,lz1,lelv)
274  $ , ta3(lx1,ly1,lz1,lelv)
275  $ , tb1(lx1,ly1,lz1,lelv)
276  $ , tb2(lx1,ly1,lz1,lelv)
277  $ , tb3(lx1,ly1,lz1,lelv)
278 
279  call opgradt (ta1,ta2,ta3,wp)
280  if ((intype.eq.0).or.(intype.eq.-1)) then
281  tolhin=tolhs
282  call ophinv (tb1,tb2,tb3,ta1,ta2,ta3,h1,h2,tolhin,nmxv)
283  else
284  if (ifanls) then
285  dtbdi = dt/bd(1) ! scale by dt*backwd-diff coefficient
286  CALL opbinv1(tb1,tb2,tb3,ta1,ta2,ta3,dtbdi)
287  else
288  CALL opbinv (tb1,tb2,tb3,ta1,ta2,ta3,h2inv)
289  endif
290  endif
291  call opdiv (ap,tb1,tb2,tb3)
292 
293  return
294  end
295 C
296 C-----------------------------------------------------------------------
297  subroutine opgrad (out1,out2,out3,inp)
298 C---------------------------------------------------------------------
299 C
300 C Compute OUTi = Di*INP, i=1,2,3.
301 C the gradient of the scalar field INP.
302 C Note: OUTi is defined on the pressure mesh !!!
303 C
304 C---------------------------------------------------------------------
305  include 'SIZE'
306  include 'GEOM'
307  include 'INPUT'
308 
309  REAL OUT1 (LX2,LY2,LZ2,1)
310  REAL OUT2 (LX2,LY2,LZ2,1)
311  REAL OUT3 (LX2,LY2,LZ2,1)
312  REAL INP (LX1,LY1,LZ1,1)
313 
314  iflg = 0
315 
316  if (ifsplit .and. .not.ifaxis) then
317  call wgradm1(out1,out2,out3,inp,nelv) ! weak grad on FLUID mesh
318  return
319  endif
320 
321  ntot2 = lx2*ly2*lz2*nelv
322  CALL multd (out1,inp,rxm2,sxm2,txm2,1,iflg)
323  CALL multd (out2,inp,rym2,sym2,tym2,2,iflg)
324  IF (ldim.EQ.3)
325  $CALL multd (out3,inp,rzm2,szm2,tzm2,3,iflg)
326 C
327  return
328  END
329 c-----------------------------------------------------------------------
330  subroutine cdtp (dtx,x,rm2,sm2,tm2,isd)
331 C-------------------------------------------------------------
332 C
333 C Compute DT*X (entire field)
334 C
335 C-------------------------------------------------------------
336  include 'SIZE'
337  include 'WZ'
338  include 'DXYZ'
339  include 'IXYZ'
340  include 'GEOM'
341  include 'MASS'
342  include 'INPUT'
343  include 'ESOLV'
344 C
345  real dtx (lx1*ly1*lz1,lelv)
346  real x (lx2*ly2*lz2,lelv)
347  real rm2 (lx2*ly2*lz2,lelv)
348  real sm2 (lx2*ly2*lz2,lelv)
349  real tm2 (lx2*ly2*lz2,lelv)
350 C
351  common /ctmp1/ wx(lx1*ly1*lz1)
352  $ , ta1(lx1*ly1*lz1)
353  $ , ta2(lx1*ly1*lz1)
354  $ , ta3(lx1*ly1,lz1)
355 
356  REAL DUAX(LX1)
357 c
358  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
359  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
360  include 'CTIMER'
361 
362  integer e
363 C
364 #ifdef TIMER
365  if (icalld.eq.0) tcdtp=0.0
366  icalld=icalld+1
367  ncdtp=icalld
368  etime1=dnekclock()
369 #endif
370 
371  nxyz1 = lx1*ly1*lz1
372  nxyz2 = lx2*ly2*lz2
373  nyz2 = ly2*lz2
374  nxy1 = lx1*ly1
375 
376  n1 = lx1*ly1
377  n2 = lx1*ly2
378 
379  do e=1,nelv
380 
381 C Use the appropriate derivative- and interpolation operator in
382 C the y-direction (= radial direction if axisymmetric).
383  if (ifaxis) then
384  ly12 = ly1*ly2
385  if (ifrzer(e)) then
386  call copy (iym12,iam12,ly12)
387  call copy (dym12,dam12,ly12)
388  call copy (w3m2,w2am2,nxyz2)
389  else
390  call copy (iym12,icm12,ly12)
391  call copy (dym12,dcm12,ly12)
392  call copy (w3m2,w2cm2,nxyz2)
393  endif
394  endif
395 C
396 C Collocate with weights
397 C
398  if(ifsplit) then
399  call col3 (wx,bm1(1,1,1,e),x(1,e),nxyz1)
400  call invcol2(wx,jacm1(1,1,1,e),nxyz1)
401  else
402  if (.not.ifaxis) call col3 (wx,w3m2,x(1,e),nxyz2)
403 C
404  if (ifaxis) then
405  if (ifrzer(e)) then
406  call col3 (wx,x(1,e),bm2(1,1,1,e),nxyz2)
407  call invcol2 (wx,jacm2(1,1,1,e),nxyz2)
408  else
409  call col3 (wx,w3m2,x(1,e),nxyz2)
410  call col2 (wx,ym2(1,1,1,e),nxyz2)
411  endif
412  endif
413  endif
414 C
415  if (ldim.eq.2) then
416  if (.not.ifdfrm(e) .and. ifalgn(e)) then
417 C
418  if ( ifrsxy(e).and.isd.eq.1 .or.
419  $ .not.ifrsxy(e).and.isd.eq.2) then
420 C
421  call col3 (ta1,wx,rm2(1,e),nxyz2)
422  call mxm (dxtm12,lx1,ta1,lx2,ta2,nyz2)
423  call mxm (ta2,lx1,iym12,ly2,dtx(1,e),ly1)
424  else
425  call col3 (ta1,wx,sm2(1,e),nxyz2)
426  call mxm (ixtm12,lx1,ta1,lx2,ta2,nyz2)
427  call mxm (ta2,lx1,dym12,ly2,dtx(1,e),ly1)
428  endif
429  else
430  call col3 (ta1,wx,rm2(1,e),nxyz2)
431  call mxm (dxtm12,lx1,ta1,lx2,ta2,nyz2)
432  call mxm (ta2,lx1,iym12,ly2,dtx(1,e),ly1)
433 
434  call col3 (ta1,wx,sm2(1,e),nxyz2)
435  call mxm (ixtm12,lx1,ta1,lx2,ta2,nyz2)
436  call mxm (ta2,lx1,dym12,ly2,ta1,ly1)
437 
438  call add2 (dtx(1,e),ta1,nxyz1)
439  endif
440 
441  else
442  if (ifsplit) then
443 
444  call col3 (ta1,wx,rm2(1,e),nxyz2)
445  call mxm (dxtm12,lx1,ta1,lx2,dtx(1,e),nyz2)
446  call col3 (ta1,wx,sm2(1,e),nxyz2)
447  i1 = 1
448  i2 = 1
449  do iz=1,lz2
450  call mxm (ta1(i2),lx1,dym12,ly2,ta2(i1),ly1)
451  i1 = i1 + n1
452  i2 = i2 + n2
453  enddo
454  call add2 (dtx(1,e),ta2,nxyz1)
455  call col3 (ta1,wx,tm2(1,e),nxyz2)
456  call mxm (ta1,nxy1,dzm12,lz2,ta2,lz1)
457  call add2 (dtx(1,e),ta2,nxyz1)
458 
459  else
460 
461  call col3 (ta1,wx,rm2(1,e),nxyz2)
462  call mxm (dxtm12,lx1,ta1,lx2,ta2,nyz2)
463  i1 = 1
464  i2 = 1
465  do iz=1,lz2
466  call mxm (ta2(i2),lx1,iym12,ly2,ta1(i1),ly1)
467  i1 = i1 + n1
468  i2 = i2 + n2
469  enddo
470  call mxm (ta1,nxy1,izm12,lz2,dtx(1,e),lz1)
471 
472  call col3 (ta1,wx,sm2(1,e),nxyz2)
473  call mxm (ixtm12,lx1,ta1,lx2,ta2,nyz2)
474  i1 = 1
475  i2 = 1
476  do iz=1,lz2
477  call mxm (ta2(i2),lx1,dym12,ly2,ta1(i1),ly1)
478  i1 = i1 + n1
479  i2 = i2 + n2
480  enddo
481  call mxm (ta1,nxy1,izm12,lz2,ta2,lz1)
482  call add2 (dtx(1,e),ta2,nxyz1)
483 
484  call col3 (ta1,wx,tm2(1,e),nxyz2)
485  call mxm (ixtm12,lx1,ta1,lx2,ta2,nyz2)
486  i1 = 1
487  i2 = 1
488  do iz=1,lz2
489  call mxm (ta2(i2),lx1,iym12,ly2,ta1(i1),ly1)
490  i1 = i1 + n1
491  i2 = i2 + n2
492  enddo
493  call mxm (ta1,nxy1,dzm12,lz2,ta2,lz1)
494  call add2 (dtx(1,e),ta2,nxyz1)
495 
496  endif
497 
498  endif
499 C
500 C If axisymmetric, add an extra diagonal term in the radial
501 C direction (only if solving the momentum equations and ISD=2)
502 C NOTE: lz1=lz2=1
503 C
504 C
505  if(ifsplit) then
506 
507  if (ifaxis.and.(isd.eq.4)) then
508  call copy (ta1,x(1,e),nxyz1)
509  if (ifrzer(e)) THEN
510  call rzero(ta1, lx1)
511  call mxm (x(1,e),lx1,datm1,ly1,duax,1)
512  call copy (ta1,duax,lx1)
513  endif
514  call col2 (ta1,baxm1(1,1,1,e),nxyz1)
515  call add2 (dtx(1,e),ta1,nxyz1)
516  endif
517 
518  else
519 
520  if (ifaxis.and.(isd.eq.2)) then
521  call col3 (ta1,x(1,e),bm2(1,1,1,e),nxyz2)
522  call invcol2 (ta1,ym2(1,1,1,e),nxyz2)
523  call mxm (ixtm12,lx1,ta1,lx2,ta2,ly2)
524  call mxm (ta2,lx1,iym12,ly2,ta1,ly1)
525  call add2 (dtx(1,e),ta1,nxyz1)
526  endif
527 
528  endif
529 
530  enddo
531 C
532 #ifdef TIMER
533  tcdtp=tcdtp+(dnekclock()-etime1)
534 #endif
535  return
536  end
537 C
538  subroutine multd (dx,x,rm2,sm2,tm2,isd,iflg)
539 C---------------------------------------------------------------------
540 C
541 C Compute D*X
542 C X : input variable, defined on M1
543 C DX : output variable, defined on M2 (note: D is rectangular)
544 C RM2 : RXM2, RYM2 or RZM2
545 C SM2 : SXM2, SYM2 or SZM2
546 C TM2 : TXM2, TYM2 or TZM2
547 C ISD : spatial direction (x=1,y=2,z=3)
548 C IFLG: OPGRAD (iflg=0) or OPDIV (iflg=1)
549 C
550 C---------------------------------------------------------------------
551  include 'SIZE'
552  include 'WZ'
553  include 'DXYZ'
554  include 'IXYZ'
555  include 'GEOM'
556  include 'MASS'
557  include 'INPUT'
558  include 'ESOLV'
559 
560  real dx (lx2*ly2*lz2,lelv)
561  real x (lx1*ly1*lz1,lelv)
562  real rm2 (lx2*ly2*lz2,lelv)
563  real sm2 (lx2*ly2*lz2,lelv)
564  real tm2 (lx2*ly2*lz2,lelv)
565 
566  common /ctmp1/ ta1(lx1*ly1*lz1)
567  $ , ta2(lx1*ly1*lz1)
568  $ , ta3(lx1*ly1*lz1)
569 
570  real duax(lx1)
571 
572  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
573  logical ifdfrm, iffast, ifh2, ifsolv
574  include 'CTIMER'
575 
576  integer e
577 C
578 #ifdef TIMER
579  if (icalld.eq.0) tmltd=0.0
580  icalld=icalld+1
581  nmltd=icalld
582  etime1=dnekclock()
583 #endif
584 
585  nyz1 = ly1*lz1
586  nxy2 = lx2*ly2
587  nxyz1 = lx1*ly1*lz1
588  nxyz2 = lx2*ly2*lz2
589 
590  n1 = lx2*ly1
591  n2 = lx2*ly2
592 
593  do e=1,nelv
594 
595 c Use the appropriate derivative- and interpolation operator in
596 c the y-direction (= radial direction if axisymmetric).
597  if (ifaxis) then
598  ly12 = ly1*ly2
599  if (ifrzer(e)) then
600  call copy (iytm12,iatm12,ly12)
601  call copy (dytm12,datm12,ly12)
602  call copy (w3m2,w2am2,nxyz2)
603  else
604  call copy (iytm12,ictm12,ly12)
605  call copy (dytm12,dctm12,ly12)
606  call copy (w3m2,w2cm2,nxyz2)
607  endif
608  endif
609 
610  if (ldim.eq.2) then
611  if (.not.ifdfrm(e) .and. ifalgn(e)) then
612 c
613  if ( ifrsxy(e).and.isd.eq.1 .or.
614  $ .not.ifrsxy(e).and.isd.eq.2) then
615  call mxm (dxm12,lx2,x(1,e),lx1,ta1,nyz1)
616  call mxm (ta1,lx2,iytm12,ly1,dx(1,e),ly2)
617  call col2 (dx(1,e),rm2(1,e),nxyz2)
618  else
619  call mxm (ixm12,lx2,x(1,e),lx1,ta1,nyz1)
620  call mxm (ta1,lx2,dytm12,ly1,dx(1,e),ly2)
621  call col2 (dx(1,e),sm2(1,e),nxyz2)
622  endif
623  else
624  call mxm (dxm12,lx2,x(1,e),lx1,ta1,nyz1)
625  call mxm (ta1,lx2,iytm12,ly1,dx(1,e),ly2)
626  call col2 (dx(1,e),rm2(1,e),nxyz2)
627  call mxm (ixm12,lx2,x(1,e),lx1,ta1,nyz1)
628  call mxm (ta1,lx2,dytm12,ly1,ta3,ly2)
629  call addcol3 (dx(1,e),ta3,sm2(1,e),nxyz2)
630  endif
631 
632  else ! 3D
633 
634 c if (ifsplit) then
635 c
636 c call mxm (dxm12,lx2,x(1,e),lx1,dx(1,e),nyz1)
637 c call col2 (dx(1,e),rm2(1,e),nxyz2)
638 c i1=1
639 c i2=1
640 c do iz=1,lz1
641 c call mxm (x(1,e),lx2,dytm12,ly1,ta1(i2),ly2)
642 c i1=i1+n1
643 c i2=i2+n2
644 c enddo
645 c call addcol3 (dx(1,e),ta1,sm2(1,e),nxyz2)
646 c call mxm (x(1,e),nxy2,dztm12,lz1,ta1,lz2)
647 c call addcol3 (dx(1,e),ta1,tm2(1,e),nxyz2)
648 
649 c else ! PN - PN-2
650 
651  call mxm (dxm12,lx2,x(1,e),lx1,ta1,nyz1)
652  i1=1
653  i2=1
654  do iz=1,lz1
655  call mxm (ta1(i1),lx2,iytm12,ly1,ta2(i2),ly2)
656  i1=i1+n1
657  i2=i2+n2
658  enddo
659  call mxm (ta2,nxy2,iztm12,lz1,dx(1,e),lz2)
660  call col2 (dx(1,e),rm2(1,e),nxyz2)
661 
662  call mxm (ixm12,lx2,x(1,e),lx1,ta3,nyz1) ! reuse ta3 below
663  i1=1
664  i2=1
665  do iz=1,lz1
666  call mxm (ta3(i1),lx2,dytm12,ly1,ta2(i2),ly2)
667  i1=i1+n1
668  i2=i2+n2
669  enddo
670  call mxm (ta2,nxy2,iztm12,lz1,ta1,lz2)
671  call addcol3 (dx(1,e),ta1,sm2(1,e),nxyz2)
672 
673 c call mxm (ixm12,lx2,x(1,e),lx1,ta1,nyz1) ! reuse ta3 from above
674  i1=1
675  i2=1
676  do iz=1,lz1
677  call mxm (ta3(i1),lx2,iytm12,ly1,ta2(i2),ly2)
678  i1=i1+n1
679  i2=i2+n2
680  enddo
681  call mxm (ta2,nxy2,dztm12,lz1,ta3,lz2)
682  call addcol3 (dx(1,e),ta3,tm2(1,e),nxyz2)
683 c endif
684  endif
685 C
686 C Collocate with the weights on the pressure mesh
687 
688 
689  if(ifsplit) then
690  call col2 (dx(1,e),bm1(1,1,1,e),nxyz1)
691  call invcol2(dx(1,e),jacm1(1,1,1,e),nxyz1)
692  else
693  if (.not.ifaxis) call col2 (dx(1,e),w3m2,nxyz2)
694  if (ifaxis) then
695  if (ifrzer(e)) then
696  call col2 (dx(1,e),bm2(1,1,1,e),nxyz2)
697  call invcol2 (dx(1,e),jacm2(1,1,1,e),nxyz2)
698  else
699  call col2 (dx(1,e),w3m2,nxyz2)
700  call col2 (dx(1,e),ym2(1,1,1,e),nxyz2)
701  endif
702  endif
703  endif
704 
705 c If axisymmetric, add an extra diagonal term in the radial
706 c direction (ISD=2).
707 c NOTE: lz1=lz2=1
708 
709  if(ifsplit) then
710 
711  if (ifaxis.and.(isd.eq.2).and.iflg.eq.1) then
712  call copy (ta3,x(1,e),nxyz1)
713  if (ifrzer(e)) then
714  call rzero(ta3, lx1)
715  call mxm (x(1,e),lx1,datm1,ly1,duax,1)
716  call copy (ta3,duax,lx1)
717  endif
718  call col2 (ta3,baxm1(1,1,1,e),nxyz1)
719  call add2 (dx(1,e),ta3,nxyz2)
720  endif
721 
722  else
723 
724  if (ifaxis.and.(isd.eq.2)) then
725  call mxm (ixm12,lx2,x(1,e),lx1,ta1,ly1)
726  call mxm (ta1,lx2,iytm12,ly1,ta2,ly2)
727  call col3 (ta3,bm2(1,1,1,e),ta2,nxyz2)
728  call invcol2 (ta3,ym2(1,1,1,e),nxyz2)
729  call add2 (dx(1,e),ta3,nxyz2)
730  endif
731 
732  endif
733 C
734  enddo
735 C
736 #ifdef TIMER
737  tmltd=tmltd+(dnekclock()-etime1)
738 #endif
739  return
740  END
741 c-----------------------------------------------------------------------
742  subroutine ophx (out1,out2,out3,inp1,inp2,inp3,h1,h2)
743 C----------------------------------------------------------------------
744 C
745 C OUT = (H1*A+H2*B) * INP
746 C
747 C----------------------------------------------------------------------
748  include 'SIZE'
749  include 'INPUT'
750  include 'SOLN'
751  REAL OUT1 (LX1,LY1,LZ1,1)
752  REAL OUT2 (LX1,LY1,LZ1,1)
753  REAL OUT3 (LX1,LY1,LZ1,1)
754  REAL INP1 (LX1,LY1,LZ1,1)
755  REAL INP2 (LX1,LY1,LZ1,1)
756  REAL INP3 (LX1,LY1,LZ1,1)
757  REAL H1 (LX1,LY1,LZ1,1)
758  REAL H2 (LX1,LY1,LZ1,1)
759 C
760  imesh = 1
761 C
762  IF (ifstrs) THEN
763  matmod = 0
764  CALL axhmsf (out1,out2,out3,inp1,inp2,inp3,h1,h2,matmod)
765  ELSE
766  CALL axhelm (out1,inp1,h1,h2,imesh,1)
767  CALL axhelm (out2,inp2,h1,h2,imesh,2)
768  IF (ldim.EQ.3)
769  $ CALL axhelm (out3,inp3,h1,h2,imesh,3)
770  ENDIF
771 C
772  return
773  END
774 c-----------------------------------------------------------------------
775  subroutine opbinv (out1,out2,out3,inp1,inp2,inp3,h2inv)
776 C--------------------------------------------------------------------
777 C
778 C Compute OUT = (H2*B)-1 * INP (explicit)
779 C
780 C--------------------------------------------------------------------
781  include 'SIZE'
782  include 'INPUT'
783  include 'MASS'
784  include 'SOLN'
785 C
786  REAL OUT1 (1)
787  REAL OUT2 (1)
788  REAL OUT3 (1)
789  REAL INP1 (1)
790  REAL INP2 (1)
791  REAL INP3 (1)
792  REAL H2INV (1)
793 C
794 
795  include 'OPCTR'
796 C
797 #ifdef TIMER
798  if (isclld.eq.0) then
799  isclld=1
800  nrout=nrout+1
801  myrout=nrout
802  rname(myrout) = 'opbinv'
803  endif
804 #endif
805 C
806  call opmask (inp1,inp2,inp3)
807  call opdssum (inp1,inp2,inp3)
808 C
809  ntot=lx1*ly1*lz1*nelv
810 C
811 #ifdef TIMER
812  isbcnt = ntot*(1+ldim)
813  dct(myrout) = dct(myrout) + (isbcnt)
814  ncall(myrout) = ncall(myrout) + 1
815  dcount = dcount + (isbcnt)
816 #endif
817 
818  call invcol3 (out1,bm1,h2inv,ntot) ! this is expensive and should
819  call dssum (out1,lx1,ly1,lz1) ! be changed (pff, 3/18/09)
820  if (if3d) then
821  do i=1,ntot
822  tmp = 1./out1(i)
823  out1(i)=inp1(i)*tmp
824  out2(i)=inp2(i)*tmp
825  out3(i)=inp3(i)*tmp
826  enddo
827  else
828  do i=1,ntot
829  tmp = 1./out1(i)
830  out1(i)=inp1(i)*tmp
831  out2(i)=inp2(i)*tmp
832  enddo
833  endif
834 
835  return
836  end
837 c-----------------------------------------------------------------------
838  subroutine opbinv1(out1,out2,out3,inp1,inp2,inp3,SCALE)
839 C--------------------------------------------------------------------
840 C
841 C Compute OUT = (B)-1 * INP (explicit)
842 C
843 C--------------------------------------------------------------------
844  include 'SIZE'
845  include 'INPUT'
846  include 'MASS'
847  include 'SOLN'
848 C
849  REAL OUT1 (1)
850  REAL OUT2 (1)
851  REAL OUT3 (1)
852  REAL INP1 (1)
853  REAL INP2 (1)
854  REAL INP3 (1)
855 C
856 
857  include 'OPCTR'
858 C
859 #ifdef TIMER
860  if (isclld.eq.0) then
861  isclld=1
862  nrout=nrout+1
863  myrout=nrout
864  rname(myrout) = 'opbnv1'
865  endif
866 #endif
867 C
868  CALL opmask (inp1,inp2,inp3)
869  CALL opdssum (inp1,inp2,inp3)
870 C
871  ntot=lx1*ly1*lz1*nelv
872 C
873 #ifdef TIMER
874  isbcnt = ntot*(1+ldim)
875  dct(myrout) = dct(myrout) + (isbcnt)
876  ncall(myrout) = ncall(myrout) + 1
877  dcount = dcount + (isbcnt)
878 #endif
879 C
880  IF (if3d) THEN
881  DO 100 i=1,ntot
882  tmp =binvm1(i,1,1,1)*scale
883  out1(i)=inp1(i)*tmp
884  out2(i)=inp2(i)*tmp
885  out3(i)=inp3(i)*tmp
886  100 CONTINUE
887  ELSE
888  DO 200 i=1,ntot
889  tmp =binvm1(i,1,1,1)*scale
890  out1(i)=inp1(i)*tmp
891  out2(i)=inp2(i)*tmp
892  200 CONTINUE
893  ENDIF
894 C
895  return
896  END
897 C-----------------------------------------------------------------------
898  subroutine uzprec (rpcg,rcg,h1m1,h2m1,intype,wp)
899 C--------------------------------------------------------------------
900 C
901 C Uzawa preconditioner
902 C
903 C--------------------------------------------------------------------
904  include 'SIZE'
905  include 'GEOM'
906  include 'INPUT'
907  include 'MASS'
908  include 'TSTEP'
909  include 'PARALLEL'
910 c
911  REAL RPCG (LX2,LY2,LZ2,LELV)
912  REAL RCG (LX2,LY2,LZ2,LELV)
913  REAL WP (LX2,LY2,LZ2,LELV)
914  REAL H1M1 (LX1,LY1,LZ1,LELV)
915  REAL H2M1 (LX1,LY1,LZ1,LELV)
916  COMMON /scrch/ h1m2(lx2,ly2,lz2,lelv)
917  $ , h2m2(lx2,ly2,lz2,lelv)
918 C
919  integer kstep
920  save kstep
921  data kstep/-1/
922 c
923  integer*8 ntotg,nxyz2
924 c
925  ntot2 = lx2*ly2*lz2*nelv
926  if (istep.ne.kstep .and. .not.ifanls) then
927  kstep=istep
928  DO 100 ie=1,nelv
929  CALL map12 (h1m2(1,1,1,ie),h1m1(1,1,1,ie),ie)
930  CALL map12 (h2m2(1,1,1,ie),h2m1(1,1,1,ie),ie)
931  100 CONTINUE
932  endif
933 C
934  IF (intype.EQ.0) THEN
935  CALL col3 (wp,rcg,h1m2,ntot2)
936  CALL col3 (rpcg,wp,bm2inv,ntot2)
937  ELSEIF (intype.EQ.-1) THEN
938  CALL eprec (wp,rcg)
939  CALL col2 (wp,h2m2,ntot2)
940  CALL col3 (rpcg,rcg,bm2inv,ntot2)
941  CALL col2 (rpcg,h1m2,ntot2)
942  CALL add2 (rpcg,wp,ntot2)
943  ELSEIF (intype.EQ.1) THEN
944  if (ifanls) then
945  CALL eprec2 (rpcg,rcg)
946  dtbd = bd(1)/dt
947  CALL cmult (rpcg,dtbd,ntot2)
948  else
949  CALL eprec2 (rpcg,rcg)
950 c CALL COL2 (RPCG,H2M2,NTOT2)
951  endif
952  ELSE
953  CALL copy (rpcg,rcg,ntot2)
954  ENDIF
955 
956  call ortho (rpcg)
957 
958  return
959  end
960 C
961  subroutine eprec (z2,r2)
962 C----------------------------------------------------------------
963 C
964 C Precondition the explicit pressure operator (E) with
965 C a Neumann type (H1) Laplace operator: JT*A*J.
966 C Invert A by conjugate gradient iteration or multigrid.
967 C
968 C NOTE: SCRNS is used.
969 C
970 C----------------------------------------------------------------
971  include 'SIZE'
972  include 'INPUT'
973  include 'GEOM'
974  include 'SOLN'
975  include 'MASS'
976  include 'PARALLEL'
977  include 'TSTEP'
978  REAL Z2 (LX2,LY2,LZ2,LELV)
979  REAL R2 (LX2,LY2,LZ2,LELV)
980  COMMON /scrns/ mask(lx1,ly1,lz1,lelv)
981  $ ,r1(lx1,ly1,lz1,lelv)
982  $ ,x1(lx1,ly1,lz1,lelv)
983  $ ,w2(lx2,ly2,lz2,lelv)
984  $ ,h1(lx1,ly1,lz1,lelv)
985  $ ,h2(lx1,ly1,lz1,lelv)
986  REAL MASK
987  COMMON /cprint/ ifprint, ifhzpc
988  LOGICAL IFPRINT, IFHZPC
989  integer*8 ntotg,nxyz
990 
991  nxyz = lx1*ly1*lz1
992  ntotg = nxyz*nelgv
993  ntot1 = nxyz*nelv
994  ntot2 = lx2*ly2*lz2*nelv
995  nfaces = 2*ldim
996 C
997 C Set the tolerance for the preconditioner
998 C
999  CALL col3 (w2,r2,bm2inv,ntot2)
1000  rinit = sqrt(glsc2(w2,r2,ntot2)/volvm2)
1001  eps = 0.02
1002  tol = eps*rinit
1003 c if (param(91).gt.0) tol=param(91)*rinit
1004 c if (param(91).lt.0) tol=-param(91)
1005 C
1006  DO 100 iel=1,nelv
1007  CALL map21e (r1(1,1,1,iel),r2(1,1,1,iel),iel)
1008  100 CONTINUE
1009 C
1010  if (ifvcor) then
1011  otr1 = glsum(r1,ntot1)
1012  call rone (x1,ntot1)
1013  c2 = -otr1/ntotg
1014  call add2s2 (r1,x1,c2,ntot1)
1015 C
1016  otr1 = glsum(r1,ntot1)
1017  tolmin = 10.*abs(otr1)
1018  IF (tol .LT. tolmin) THEN
1019  if (nid.eq.0)
1020  $ write(6,*) 'Resetting tol in EPREC:(old,new)',tol,tolmin
1021  tol = tolmin
1022  ENDIF
1023  ENDIF
1024 C
1025  CALL rone (h1,ntot1)
1026  CALL rzero (h2,ntot1)
1027  ifhzpc = .true.
1028  CALL hmholtz ('PREC',x1,r1,h1,h2,pmask,vmult,imesh,tol,nmxp,1)
1029  ifhzpc = .false.
1030 C
1031  DO 200 iel=1,nelv
1032  CALL map12 (z2(1,1,1,iel),x1(1,1,1,iel),iel)
1033  200 CONTINUE
1034 C
1035  return
1036  END
1037 C
1038  subroutine convprn (iconv,rbnorm,rrpt,res,z,tol)
1039 C-----------------------------------------------------------------
1040 C T
1041 C Convergence test for the pressure step; r z
1042 C
1043 C-----------------------------------------------------------------
1044  include 'SIZE'
1045  include 'MASS'
1046  REAL RES (1)
1047  REAL Z (1)
1048  REAL wrk1(2),wrk2(2)
1049 c
1050  ntot2 = lx2*ly2*lz2*nelv
1051  wrk1(1) = vlsc21(res,bm2inv,ntot2) ! res*bm2inv*res
1052  wrk1(2) = vlsc2(res,z ,ntot2) ! res*z
1053  call gop(wrk1,wrk2,'+ ',2)
1054  rbnorm = sqrt(wrk1(1)/volvm2)
1055  rrpt = wrk1(2)
1056 c
1057 c CALL CONVPR (RCG,tolpss,ICONV,RNORM)
1058 c RRP1 = GLSC2 (RPCG,RCG,NTOT2)
1059 c RBNORM = SQRT (GLSC2 (BM2,TB,NTOT2)/VOLVM2)
1060 c
1061  iconv = 0
1062  IF (rbnorm.LT.tol) iconv=1
1063  return
1064  END
1065 C
1066  subroutine convpr (res,tol,iconv,rbnorm)
1067 C-----------------------------------------------------------------
1068 C
1069 C Convergence test for the pressure step
1070 C
1071 C-----------------------------------------------------------------
1072  include 'SIZE'
1073  include 'MASS'
1074  REAL RES (LX2,LY2,LZ2,LELV)
1075  COMMON /scrmg/ ta(lx2,ly2,lz2,lelv)
1076  $ , tb(lx2,ly2,lz2,lelv)
1077 C
1078  rbnorm = 0.
1079  ntot2 = lx2*ly2*lz2*nelv
1080  CALL col3 (ta,res,bm2inv,ntot2)
1081  CALL col3 (tb,ta,ta,ntot2)
1082  rbnorm = sqrt(glsc2(bm2,tb,ntot2)/volvm2)
1083 c
1084  iconv = 0
1085  IF (rbnorm.LT.tol) iconv=1
1086  return
1087  END
1088 C
1089  subroutine chktcg2 (tol,res,iconv)
1090 C-------------------------------------------------------------------
1091 C
1092 C Check that the tolerances are not too small for the CG-solver.
1093 C Important when calling the CG-solver (Gauss mesh) with
1094 C all Dirichlet velocity b.c. (zero Neumann for the pressure).
1095 C
1096 C-------------------------------------------------------------------
1097  include 'SIZE'
1098  include 'GEOM'
1099  include 'INPUT'
1100  include 'MASS'
1101  include 'TSTEP'
1102  REAL RES (LX2,LY2,LZ2,LELV)
1103  COMMON /ctmp0/ ta(lx2,ly2,lz2,lelv)
1104  $ , tb(lx2,ly2,lz2,lelv)
1105  COMMON /cprint/ ifprint
1106  LOGICAL IFPRINT
1107 C
1108  iconv = 0
1109 C
1110 C Single or double precision???
1111 C
1112  delta = 1.e-9
1113  x = 1.+delta
1114  y = 1.
1115  diff = abs(x-y)
1116  IF (diff.EQ.0.) eps = 1.e-5
1117  IF (diff.GT.0.) eps = 1.e-10
1118 C
1119 C Relaxed pressure iteration; maximum decrease in the residual (ER)
1120 C
1121  IF (prelax.NE.0.) eps = prelax
1122 C
1123  ntot2 = lx2*ly2*lz2*nelv
1124  CALL col3 (ta,res,bm2inv,ntot2)
1125  CALL col3 (tb,ta,ta,ntot2)
1126  CALL col2 (tb,bm2,ntot2)
1127  rinit = sqrt( glsum(tb,ntot2)/volvm2 )
1128  IF (rinit.LT.tol) THEN
1129  iconv = 1
1130  return
1131  ENDIF
1132  IF (tolpdf.GT.0.) THEN
1133  rmin = tolpdf
1134  ELSE
1135  rmin = eps*rinit
1136  ENDIF
1137  IF (tol.LT.rmin) THEN
1138  tolold = tol
1139  tol = rmin
1140  IF (nio.EQ.0 .AND. ifprint) WRITE (6,*)
1141  $ 'New CG2-tolerance (RINIT*10-5/10-10) = ',tol,tolold
1142  ENDIF
1143  IF (ifvcor) THEN
1144  otr = glsum(res,ntot2)
1145  tolmin = abs(otr)*100.
1146  IF (tol .LT. tolmin) THEN
1147  tolold = tol
1148  tol = tolmin
1149  IF (nio.EQ.0 .AND. ifprint)
1150  $ WRITE (6,*) 'New CG2-tolerance (OTR) = ',tolmin,tolold
1151  ENDIF
1152  ENDIF
1153  return
1154  END
1155 C
1156  subroutine dudxyz (du,u,rm1,sm1,tm1,jm1,imsh,isd)
1157 C--------------------------------------------------------------
1158 C
1159 C DU - dU/dx or dU/dy or dU/dz
1160 C U - a field variable defined on mesh 1
1161 C RM1 - dr/dx or dr/dy or dr/dz
1162 C SM1 - ds/dx or ds/dy or ds/dz
1163 C TM1 - dt/dx or dt/dy or dt/dz
1164 C JM1 - the Jacobian
1165 C IMESH - topology: velocity (1) or temperature (2) mesh
1166 C
1167 C--------------------------------------------------------------
1168  include 'SIZE'
1169  include 'DXYZ'
1170  include 'GEOM'
1171  include 'INPUT'
1172  include 'TSTEP'
1173 C
1174  REAL DU (LX1,LY1,LZ1,1)
1175  REAL U (LX1,LY1,LZ1,1)
1176  REAL RM1 (LX1,LY1,LZ1,1)
1177  REAL SM1 (LX1,LY1,LZ1,1)
1178  REAL TM1 (LX1,LY1,LZ1,1)
1179  REAL JM1 (LX1,LY1,LZ1,1)
1180 C
1181  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1182  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
1183 C
1184  REAL DRST(LX1,LY1,LZ1)
1185 C
1186  IF (imsh.EQ.1) nel = nelv
1187  IF (imsh.EQ.2) nel = nelt
1188  nxy1 = lx1*ly1
1189  nyz1 = ly1*lz1
1190  nxyz1 = lx1*ly1*lz1
1191  ntot = nxyz1*nel
1192 
1193  DO 1000 iel=1,nel
1194 C
1195  IF (ifaxis) CALL setaxdy (ifrzer(iel) )
1196 C
1197  IF (ldim.EQ.2) THEN
1198  CALL mxm (dxm1,lx1,u(1,1,1,iel),lx1,du(1,1,1,iel),nyz1)
1199  CALL col2 (du(1,1,1,iel),rm1(1,1,1,iel),nxyz1)
1200  CALL mxm (u(1,1,1,iel),lx1,dytm1,ly1,drst,ly1)
1201  CALL addcol3 (du(1,1,1,iel),drst,sm1(1,1,1,iel),nxyz1)
1202  ELSE
1203  CALL mxm (dxm1,lx1,u(1,1,1,iel),lx1,du(1,1,1,iel),nyz1)
1204  CALL col2 (du(1,1,1,iel),rm1(1,1,1,iel),nxyz1)
1205  DO 20 iz=1,lz1
1206  CALL mxm (u(1,1,iz,iel),lx1,dytm1,ly1,drst(1,1,iz),ly1)
1207  20 CONTINUE
1208  CALL addcol3 (du(1,1,1,iel),drst,sm1(1,1,1,iel),nxyz1)
1209  CALL mxm (u(1,1,1,iel),nxy1,dztm1,lz1,drst,lz1)
1210  CALL addcol3 (du(1,1,1,iel),drst,tm1(1,1,1,iel),nxyz1)
1211  ENDIF
1212 C
1213  1000 CONTINUE
1214 C
1215 c CALL INVCOL2 (DU,JM1,NTOT)
1216  CALL col2 (du,jacmi,ntot)
1217 C
1218  return
1219  END
1220 C
1221  subroutine convopo (conv,fi)
1222 C--------------------------------------------------------------------
1223 C
1224 C Compute the convective term CONV for a passive scalar field FI
1225 C using the skew-symmetric formulation.
1226 C The field variable FI is defined on mesh M1 (GLL) and
1227 C the velocity field is assumed given.
1228 C
1229 C IMPORTANT NOTE: Use the scratch-arrays carefully!!!!!
1230 C
1231 C The common-block SCRNS is used in CONV1 and CONV2.
1232 C The common-blocks CTMP0 and CTMP1 are also used as scratch-arrays
1233 C since there is no direct stiffness summation or Helmholtz-solves.
1234 C
1235 C--------------------------------------------------------------------
1236  include 'SIZE'
1237  include 'TOTAL'
1238 C
1239 C Use the common blocks CTMP0 and CTMP1 as work space.
1240 C
1241  COMMON /scrch/ cmask1(lx1,ly1,lz1,lelv)
1242  $ , cmask2(lx1,ly1,lz1,lelv)
1243  COMMON /ctmp1/ mfi(lx1,ly1,lz1,lelv)
1244  $ , dmfi(lx1,ly1,lz1,lelv)
1245  $ , mdmfi(lx1,ly1,lz1,lelv)
1246  REAL MFI,DMFI,MDMFI
1247 C
1248 C Arrays in parameter list
1249 C
1250  REAL CONV (LX1,LY1,LZ1,1)
1251  REAL FI (LX1,LY1,LZ1,1)
1252 C
1253 C
1254  nxyz1 = lx1*ly1*lz1
1255  ntot1 = lx1*ly1*lz1*nelv
1256  ntotz = lx1*ly1*lz1*nelt
1257 C
1258  CALL rzero (conv,ntotz)
1259 C
1260  IF (param(86).EQ.0.0) THEN
1261 C Always use the convective form !!! (ER)
1262  CALL conv1 (conv,fi)
1263  ELSE
1264 C Use skew-symmetric form
1265 C
1266 C Generate operator mask arrays CMASK1 and CMASK2
1267 C
1268  CALL cmask (cmask1,cmask2)
1269 C
1270  CALL col3 (mfi,fi,cmask1,ntot1)
1271  CALL conv1 (dmfi,mfi)
1272  CALL col3 (mdmfi,dmfi,cmask1,ntot1)
1273  CALL add2s2 (conv,mdmfi,0.5,ntot1)
1274 C
1275  CALL col3 (mdmfi,dmfi,cmask2,ntot1)
1276  CALL add2 (conv,mdmfi,ntot1)
1277 C
1278  CALL conv2 (dmfi,mfi)
1279  CALL col3 (mdmfi,dmfi,cmask1,ntot1)
1280  CALL add2s2 (conv,mdmfi,-0.5,ntot1)
1281 C
1282  CALL col3 (mfi,fi,cmask2,ntot1)
1283  CALL conv1 (dmfi,mfi)
1284  CALL col3 (mdmfi,dmfi,cmask2,ntot1)
1285  CALL add2s2 (conv,mdmfi,0.5,ntot1)
1286 C
1287  CALL conv2 (dmfi,mfi)
1288  CALL col3 (mdmfi,dmfi,cmask1,ntot1)
1289  CALL add2s2 (conv,mdmfi,-1.,ntot1)
1290 C
1291  CALL col3 (mdmfi,dmfi,cmask2,ntot1)
1292  CALL add2s2 (conv,mdmfi,0.5,ntot1)
1293  ENDIF
1294 C
1295  return
1296  END
1297 c-----------------------------------------------------------------------
1298  subroutine conv2 (dtfi,fi)
1299 C--------------------------------------------------------------------
1300 C
1301 C Compute DT*FI (part of the convection operator)
1302 C
1303 C--------------------------------------------------------------------
1304  include 'SIZE'
1305  include 'TOTAL'
1306  REAL DTFI (LX1,LY1,LZ1,1)
1307  REAL FI (LX1,LY1,LZ1,1)
1308  COMMON /scrns/ ta1(lx1,ly1,lz1,lelv)
1309  $ , ta2(lx1,ly1,lz1,lelv)
1310  $ , ta3(lx1,ly1,lz1,lelv)
1311  $ , tb1(lx1,ly1,lz1,lelv)
1312  $ , tb2(lx1,ly1,lz1,lelv)
1313  $ , tb3(lx1,ly1,lz1,lelv)
1314 C
1315  nxy1 = lx1*ly1
1316  nyz1 = ly1*lz1
1317  nxyz1 = lx1*ly1*lz1
1318  ntot1 = lx1*ly1*lz1*nelv
1319 C
1320  CALL rzero(dtfi,ntot1)
1321 C
1322  IF (ldim .EQ. 2) THEN
1323 C
1324 C 2-dimensional case
1325 C
1326  CALL col4 (ta1,rxm1,vx,fi,ntot1)
1327  CALL col4 (ta2,rym1,vy,fi,ntot1)
1328  CALL add2 (ta1,ta2,ntot1)
1329  DO 200 iel=1,nelv
1330  CALL col2 (ta1(1,1,1,iel),w3m1,nxyz1)
1331  CALL mxm (dxtm1,lx1,ta1(1,1,1,iel),lx1,tb1(1,1,1,iel),ly1)
1332  200 CONTINUE
1333  CALL copy(dtfi,tb1,ntot1)
1334 C
1335  CALL col4 (ta1,sxm1,vx,fi,ntot1)
1336  CALL col4 (ta2,sym1,vy,fi,ntot1)
1337  CALL add2 (ta1,ta2,ntot1)
1338  DO 300 iel=1,nelv
1339  CALL col2 (ta1(1,1,1,iel),w3m1,nxyz1)
1340  CALL mxm (ta1(1,1,1,iel),lx1,dym1,ly1,tb1(1,1,1,iel),ly1)
1341  300 CONTINUE
1342  CALL add2 (dtfi,tb1,ntot1)
1343  CALL invcol2 (dtfi,bm1,ntot1)
1344 C
1345  ELSE
1346 C
1347 C 3-dimensional case
1348 C
1349  CALL col4 (ta1,rxm1,vx,fi,ntot1)
1350  CALL col4 (ta2,rym1,vy,fi,ntot1)
1351  CALL col4 (ta3,rzm1,vz,fi,ntot1)
1352  CALL add2 (ta1,ta2,ntot1)
1353  CALL add2 (ta1,ta3,ntot1)
1354  DO 600 iel=1,nelv
1355  CALL col2 (ta1(1,1,1,iel),w3m1,nxyz1)
1356  CALL mxm (dxtm1,lx1,ta1(1,1,1,iel),lx1,tb1(1,1,1,iel),nyz1)
1357  600 CONTINUE
1358  CALL copy(dtfi,tb1,ntot1)
1359 C
1360  CALL col4 (ta1,sxm1,vx,fi,ntot1)
1361  CALL col4 (ta2,sym1,vy,fi,ntot1)
1362  CALL col4 (ta3,szm1,vz,fi,ntot1)
1363  CALL add2 (ta1,ta2,ntot1)
1364  CALL add2 (ta1,ta3,ntot1)
1365  DO 700 iel=1,nelv
1366  CALL col2 (ta1(1,1,1,iel),w3m1,nxyz1)
1367  DO 710 iz=1,lz1
1368  CALL mxm (ta1(1,1,iz,iel),lx1,dym1,ly1,tb1(1,1,iz,iel),ly1)
1369  710 CONTINUE
1370  700 CONTINUE
1371  CALL add2 (dtfi,tb1,ntot1)
1372 C
1373  CALL col4 (ta1,txm1,vx,fi,ntot1)
1374  CALL col4 (ta2,tym1,vy,fi,ntot1)
1375  CALL col4 (ta3,tzm1,vz,fi,ntot1)
1376  CALL add2 (ta1,ta2,ntot1)
1377  CALL add2 (ta1,ta3,ntot1)
1378  DO 800 iel=1,nelv
1379  CALL col2 (ta1(1,1,1,iel),w3m1,nxyz1)
1380  CALL mxm (ta1(1,1,1,iel),nxy1,dzm1,lz1,tb1(1,1,1,iel),lz1)
1381  800 CONTINUE
1382  CALL add2 (dtfi,tb1,ntot1)
1383  CALL invcol2 (dtfi,bm1,ntot1)
1384 C
1385  ENDIF
1386 C
1387  return
1388  END
1389 C
1390  subroutine cmask (cmask1,cmask2)
1391 C---------------------------------------------------------------
1392 C
1393 C Generate maskarrays for the convection operator
1394 C
1395 C CMASK1 = 1.0 in the interior
1396 C = 0.0 at outflow
1397 C CMASK2 = 0.0 in the interior
1398 C = 1.0 at outflow
1399 C
1400 C---------------------------------------------------------------
1401  include 'SIZE'
1402  include 'INPUT'
1403  include 'TSTEP'
1404 C
1405  REAL CMASK1 (LX1,LY1,LZ1,LELV)
1406  REAL CMASK2 (LX1,LY1,LZ1,LELV)
1407 C
1408  CHARACTER CB*1
1409 C
1410  ntot1 = lx1*ly1*lz1*nelv
1411  nfaces = 2*ldim
1412  CALL cfill (cmask1,1.,ntot1)
1413  CALL cfill (cmask2,0.,ntot1)
1414  DO 100 iel=1,nelv
1415  DO 100 iface=1,nfaces
1416  cb = cbc(iface,iel,ifield)
1417  IF (cb.EQ.'O' .OR. cb.EQ.'o') THEN
1418  CALL facev (cmask1,iel,iface,0.,lx1,ly1,lz1)
1419  CALL facev (cmask2,iel,iface,1.,lx1,ly1,lz1)
1420  ENDIF
1421  100 CONTINUE
1422  return
1423  END
1424 C
1425  subroutine makef
1426 C---------------------------------------------------------------------
1427 C
1428 C Compute and add: (1) user specified forcing function (FX,FY,FZ)
1429 C (2) driving force due to natural convection
1430 C (3) convection term
1431 C
1432 C !! NOTE: Do not change the arrays BFX, BFY, BFZ until the
1433 C current time step is completed.
1434 C
1435 C----------------------------------------------------------------------
1436  include 'SIZE'
1437  include 'SOLN'
1438  include 'MASS'
1439  include 'INPUT'
1440  include 'TSTEP'
1441  include 'CTIMER'
1442  include 'MVGEOM'
1443 
1444  etime1 = dnekclock()
1445  call makeuf
1446  if (filtertype.eq.2) call make_hpf
1447  if (ifexplvis.and.ifsplit) call makevis
1448  if (ifnav .and..not.ifchar) call advab
1449  if (ifmvbd.and..not.ifchar) call admeshv
1450  if (iftran) call makeabf
1451  if ((iftran.and..not.ifchar).or.
1452  $ (iftran.and..not.ifnav.and.ifchar)) call makebdf
1453  if (ifnav.and.ifchar) call advchar
1454 
1455 c Adding this call allows prescribed pressure bc for PnPn-2
1456 c if (.not.ifsplit.and..not.ifstrs) call bcneutr
1457 
1458  tmakf=tmakf+(dnekclock()-etime1)
1459 
1460  return
1461  end
1462 c
1463  subroutine makeuf
1464 C---------------------------------------------------------------------
1465 C
1466 C Compute and add: (1) user specified forcing function (FX,FY,FZ)
1467 C
1468 C----------------------------------------------------------------------
1469  include 'SIZE'
1470  include 'SOLN'
1471  include 'MASS'
1472  include 'TSTEP'
1473 C
1474  time = time-dt
1475  CALL nekuf (bfx,bfy,bfz)
1476  CALL opcolv (bfx,bfy,bfz,bm1)
1477  time = time+dt
1478 C
1479  return
1480  END
1481 C
1482  subroutine nekuf (f1,f2,f3)
1483  include 'SIZE'
1484  include 'PARALLEL'
1485  include 'NEKUSE'
1486  include 'CTIMER'
1487 
1488  REAL F1 (LX1,LY1,LZ1,LELV)
1489  REAL F2 (LX1,LY1,LZ1,LELV)
1490  REAL F3 (LX1,LY1,LZ1,LELV)
1491 
1492 #ifdef TIMER
1493  etime1=dnekclock_sync()
1494 #endif
1495 
1496  CALL oprzero (f1,f2,f3)
1497  DO 100 iel=1,nelv
1498  ielg = lglel(iel)
1499  DO 100 k=1,lz1
1500  DO 100 j=1,ly1
1501  DO 100 i=1,lx1
1502  if (optlevel.le.2) CALL nekasgn (i,j,k,iel)
1503  CALL userf (i,j,k,ielg)
1504  f1(i,j,k,iel) = ffx
1505  f2(i,j,k,iel) = ffy
1506  f3(i,j,k,iel) = ffz
1507  100 CONTINUE
1508 
1509 #ifdef TIMER
1510  tusfq=tusfq+(dnekclock()-etime1)
1511 #endif
1512 
1513  return
1514  END
1515 C
1516  subroutine advab
1517 C---------------------------------------------------------------
1518 C
1519 C Eulerian scheme, add convection term to forcing function
1520 C at current time step.
1521 C
1522 C---------------------------------------------------------------
1523  include 'SIZE'
1524  include 'SOLN'
1525  include 'MASS'
1526  include 'TSTEP'
1527 C
1528  COMMON /scruz/ ta1(lx1,ly1,lz1,lelv)
1529  $ , ta2(lx1,ly1,lz1,lelv)
1530  $ , ta3(lx1,ly1,lz1,lelv)
1531 C
1532  ntot1 = lx1*ly1*lz1*nelv
1533  CALL convop (ta1,vx)
1534  CALL convop (ta2,vy)
1535  CALL subcol3 (bfx,bm1,ta1,ntot1)
1536  CALL subcol3 (bfy,bm1,ta2,ntot1)
1537  IF (ldim.EQ.2) THEN
1538  CALL rzero (ta3,ntot1)
1539  ELSE
1540  CALL convop (ta3,vz)
1541  CALL subcol3 (bfz,bm1,ta3,ntot1)
1542  ENDIF
1543 C
1544 
1545  return
1546  END
1547 c-----------------------------------------------------------------------
1548  subroutine makebdf
1549 C
1550 C Add contributions to F from lagged BD terms.
1551 C
1552  include 'SIZE'
1553  include 'SOLN'
1554  include 'MASS'
1555  include 'GEOM'
1556  include 'INPUT'
1557  include 'TSTEP'
1558 C
1559  COMMON /scrns/ ta1(lx1,ly1,lz1,lelv)
1560  $ , ta2(lx1,ly1,lz1,lelv)
1561  $ , ta3(lx1,ly1,lz1,lelv)
1562  $ , tb1(lx1,ly1,lz1,lelv)
1563  $ , tb2(lx1,ly1,lz1,lelv)
1564  $ , tb3(lx1,ly1,lz1,lelv)
1565  $ , h2(lx1,ly1,lz1,lelv)
1566 C
1567  ntot1 = lx1*ly1*lz1*nelv
1568  const = 1./dt
1569 
1570  if(iflomach) then
1571  call cfill(h2,const,ntot1)
1572  else
1573  call cmult2(h2,vtrans(1,1,1,1,ifield),const,ntot1)
1574  endif
1575 
1576  CALL opcolv3c (tb1,tb2,tb3,vx,vy,vz,bm1,bd(2))
1577 C
1578  DO 100 ilag=2,nbd
1579  IF (ifgeom) THEN
1580  CALL opcolv3c(ta1,ta2,ta3,vxlag(1,1,1,1,ilag-1),
1581  $ vylag(1,1,1,1,ilag-1),
1582  $ vzlag(1,1,1,1,ilag-1),
1583  $ bm1lag(1,1,1,1,ilag-1),bd(ilag+1))
1584  ELSE
1585  CALL opcolv3c(ta1,ta2,ta3,vxlag(1,1,1,1,ilag-1),
1586  $ vylag(1,1,1,1,ilag-1),
1587  $ vzlag(1,1,1,1,ilag-1),
1588  $ bm1 ,bd(ilag+1))
1589  ENDIF
1590  CALL opadd2 (tb1,tb2,tb3,ta1,ta2,ta3)
1591  100 CONTINUE
1592  CALL opadd2col (bfx,bfy,bfz,tb1,tb2,tb3,h2)
1593 C
1594  return
1595  END
1596 c-----------------------------------------------------------------------
1597  subroutine makeabf
1598 C-----------------------------------------------------------------------
1599 C
1600 C Sum up contributions to kth order extrapolation scheme.
1601 C
1602 C-----------------------------------------------------------------------
1603  include 'SIZE'
1604  include 'INPUT'
1605  include 'SOLN'
1606  include 'TSTEP'
1607 C
1608  COMMON /scruz/ ta1(lx1,ly1,lz1,lelv)
1609  $ , ta2(lx1,ly1,lz1,lelv)
1610  $ , ta3(lx1,ly1,lz1,lelv)
1611 C
1612  ntot1 = lx1*ly1*lz1*nelv
1613 C
1614  ab0 = ab(1)
1615  ab1 = ab(2)
1616  ab2 = ab(3)
1617  CALL add3s2 (ta1,abx1,abx2,ab1,ab2,ntot1)
1618  CALL add3s2 (ta2,aby1,aby2,ab1,ab2,ntot1)
1619  CALL copy (abx2,abx1,ntot1)
1620  CALL copy (aby2,aby1,ntot1)
1621  CALL copy (abx1,bfx,ntot1)
1622  CALL copy (aby1,bfy,ntot1)
1623  CALL add2s1 (bfx,ta1,ab0,ntot1)
1624  CALL add2s1 (bfy,ta2,ab0,ntot1)
1625  if(.not.iflomach) CALL col2 (bfx,vtrans,ntot1) ! multiply by density
1626  if(.not.iflomach) CALL col2 (bfy,vtrans,ntot1)
1627  IF (ldim.EQ.3) THEN
1628  CALL add3s2 (ta3,abz1,abz2,ab1,ab2,ntot1)
1629  CALL copy (abz2,abz1,ntot1)
1630  CALL copy (abz1,bfz,ntot1)
1631  CALL add2s1 (bfz,ta3,ab0,ntot1)
1632  if(.not.iflomach) CALL col2 (bfz,vtrans,ntot1)
1633  ENDIF
1634 C
1635  return
1636  END
1637 C
1638 c-----------------------------------------------------------------------
1639  subroutine setab3 (ab0,ab1,ab2)
1640 C
1641 C Set coefficients for 3rd order Adams-Bashforth scheme
1642 C (variable time step).
1643 C
1644  include 'SIZE'
1645  include 'TSTEP'
1646 C
1647  IF (istep.LE.2) THEN
1648  ab0 = 1.
1649  ab1 = 0.
1650  ab2 = 0.
1651  ELSE
1652  dt0 = dtlag(1)
1653  dt1 = dtlag(2)
1654  dt2 = dtlag(3)
1655  ab2 = (dt0/dt2)*((dt0/3.+dt1/2.)/(dt1+dt2))
1656  ab1 = -(dt0/dt1)*(0.5+(dt0/3.+dt1/2.)/dt2)
1657  ab0 = 1.-ab1-ab2
1658  ENDIF
1659  return
1660  END
1661 C
1662  subroutine setabbd (ab,dtlag,nab,nbd)
1663 C-----------------------------------------------------------------------
1664 C
1665 C Compute Adams-Bashforth coefficients (order NAB, less or equal to 3)
1666 C
1667 C NBD .EQ. 1
1668 C Standard Adams-Bashforth coefficients
1669 C
1670 C NBD .GT. 1
1671 C Modified Adams-Bashforth coefficients to be used in con-
1672 C junction with Backward Differentiation schemes (order NBD)
1673 C
1674 C-----------------------------------------------------------------------
1675  REAL AB(NAB),DTLAG(NAB)
1676 C
1677  dt0 = dtlag(1)
1678  dt1 = dtlag(2)
1679  dt2 = dtlag(3)
1680 C
1681  IF ( nab.EQ.1 ) THEN
1682 C
1683  ab(1) = 1.0
1684 C
1685  ELSEIF ( nab.EQ.2 ) THEN
1686 C
1687  dta = dt0/dt1
1688 C
1689  IF ( nbd.EQ.1 ) THEN
1690 C
1691  ab(2) = -0.5*dta
1692  ab(1) = 1.0 - ab(2)
1693 C
1694  ELSEIF ( nbd.EQ.2 ) THEN
1695 C
1696  ab(2) = -dta
1697  ab(1) = 1.0 - ab(2)
1698 C
1699  ENDIF
1700 C
1701  ELSEIF ( nab.EQ.3 ) THEN
1702 C
1703  dts = dt1 + dt2
1704  dta = dt0 / dt1
1705  dtb = dt1 / dt2
1706  dtc = dt0 / dt2
1707  dtd = dts / dt1
1708  dte = dt0 / dts
1709 C
1710  IF ( nbd.EQ.1 ) THEN
1711 C
1712  ab(3) = dte*( 0.5*dtb + dtc/3. )
1713  ab(2) = -0.5*dta - ab(3)*dtd
1714  ab(1) = 1.0 - ab(2) - ab(3)
1715 C
1716  ELSEIF ( nbd.EQ.2 ) THEN
1717 C
1718  ab(3) = 2./3.*dtc*(1./dtd + dte)
1719  ab(2) = -dta - ab(3)*dtd
1720  ab(1) = 1.0 - ab(2) - ab(3)
1721 C
1722  ELSEIF ( nbd.EQ.3 ) THEN
1723 C
1724  ab(3) = dte*(dtb + dtc)
1725  ab(2) = -dta*(1.0 + dtb + dtc)
1726  ab(1) = 1.0 - ab(2) - ab(3)
1727 C
1728  ENDIF
1729 C
1730  ENDIF
1731 C
1732  return
1733  END
1734 C
1735  subroutine setbd (bd,dtbd,nbd)
1736 C-----------------------------------------------------------------------
1737 C
1738 C Compute bacward-differentiation coefficients of order NBD
1739 C
1740 C-----------------------------------------------------------------------
1741  parameter(ldim = 10)
1742  REAL BDMAT(ldim,ldim),BDRHS(ldim)
1743  INTEGER IR(ldim),IC(ldim)
1744  REAL BD(1),DTBD(1)
1745 C
1746  IF (nbd.EQ.1) THEN
1747  bd(1) = 1.
1748  bdf = 1.
1749  ELSEIF (nbd.GE.2) THEN
1750  nsys = nbd+1
1751  CALL bdsys (bdmat,bdrhs,dtbd,nbd,ldim)
1752  CALL lu (bdmat,nsys,ldim,ir,ic)
1753  CALL solve (bdrhs,bdmat,1,nsys,ldim,ir,ic)
1754  DO 30 i=1,nbd
1755  bd(i) = bdrhs(i)
1756  30 CONTINUE
1757  bdf = bdrhs(nbd+1)
1758  ENDIF
1759 C
1760 C Normalize
1761 C
1762  DO 100 ibd=nbd,1,-1
1763  bd(ibd+1) = bd(ibd)
1764  100 CONTINUE
1765  bd(1) = 1.
1766  DO 200 ibd=1,nbd+1
1767  bd(ibd) = bd(ibd)/bdf
1768  200 CONTINUE
1769 c write(6,1) (bd(k),k=1,nbd+1)
1770 c 1 format('bd:',1p8e13.5)
1771 C
1772  return
1773  END
1774 C
1775  subroutine bdsys (a,b,dt,nbd,ldim)
1776  REAL A(ldim,9),B(9),DT(9)
1777  CALL rzero (a,ldim**2)
1778  n = nbd+1
1779  DO 10 j=1,nbd
1780  a(1,j) = 1.
1781  10 CONTINUE
1782  a(1,nbd+1) = 0.
1783  b(1) = 1.
1784  DO 20 j=1,nbd
1785  sumdt = 0.
1786  DO 25 k=1,j
1787  sumdt = sumdt+dt(k)
1788  25 CONTINUE
1789  a(2,j) = sumdt
1790  20 CONTINUE
1791  a(2,nbd+1) = -dt(1)
1792  b(2) = 0.
1793  DO 40 i=3,nbd+1
1794  DO 30 j=1,nbd
1795  sumdt = 0.
1796  DO 35 k=1,j
1797  sumdt = sumdt+dt(k)
1798  35 CONTINUE
1799  a(i,j) = sumdt**(i-1)
1800  30 CONTINUE
1801  a(i,nbd+1) = 0.
1802  b(i) = 0.
1803  40 CONTINUE
1804  return
1805  END
1806 C
1807  subroutine tauinit (tau,ilag)
1808 C-------------------------------------------------------------------
1809 C
1810 C Set initial time for subintegration
1811 C
1812 C-------------------------------------------------------------------
1813  include 'SIZE'
1814  include 'TSTEP'
1815  tau = 0.
1816  DO 10 i=nbd,ilag+1,-1
1817  tau = tau+dtlag(i)
1818  10 CONTINUE
1819  return
1820  END
1821 C
1822  subroutine velinit (vel1,vel2,vel3,ilag)
1823 C-------------------------------------------------------------------
1824 C
1825 C Set initial conditions for subintegration
1826 C
1827 C-------------------------------------------------------------------
1828  include 'SIZE'
1829  include 'SOLN'
1830  REAL VEL1 (LX1,LY1,LZ1,LELV)
1831  REAL VEL2 (LX1,LY1,LZ1,LELV)
1832  REAL VEL3 (LX1,LY1,LZ1,LELV)
1833  IF (ilag.EQ.1) THEN
1834  CALL opcopy (vel1,vel2,vel3,vx,vy,vz)
1835  ELSE
1836  CALL opcopy (vel1,vel2,vel3,vxlag(1,1,1,1,ilag-1)
1837  $ ,vylag(1,1,1,1,ilag-1)
1838  $ ,vzlag(1,1,1,1,ilag-1) )
1839  ENDIF
1840  return
1841  END
1842 C
1843  subroutine velconv (vxn,vyn,vzn,tau)
1844 C--------------------------------------------------------------------
1845 C
1846 C Compute convecting velocity field (linearization)
1847 C
1848 C--------------------------------------------------------------------
1849  include 'SIZE'
1850  include 'SOLN'
1851  include 'TSTEP'
1852  REAL VXN (LX1,LY1,LZ1,LELV)
1853  REAL VYN (LX1,LY1,LZ1,LELV)
1854  REAL VZN (LX1,LY1,LZ1,LELV)
1855  CALL velchar (vx,vxn,vxlag,nbd,tau,dtlag)
1856  CALL velchar (vy,vyn,vylag,nbd,tau,dtlag)
1857  IF (ldim.EQ.3)
1858  $CALL velchar (vz,vzn,vzlag,nbd,tau,dtlag)
1859  return
1860  END
1861 C
1862  subroutine frkconv (y,x,mask)
1863 C--------------------------------------------------------------------
1864 C
1865 C Evaluate right-hand-side for Runge-Kutta scheme in the case of
1866 C pure convection.
1867 C
1868 C--------------------------------------------------------------------
1869  include 'SIZE'
1870  include 'MASS'
1871  include 'TSTEP'
1872  REAL Y (LX1,LY1,LZ1,1)
1873  REAL X (LX1,LY1,LZ1,1)
1874  REAL MASK (LX1,LY1,LZ1,1)
1875 C
1876  IF (imesh.EQ.1) nel=nelv
1877  IF (imesh.EQ.2) nel=nelt
1878  ntot1 = lx1*ly1*lz1*nel
1879  CALL convop (y,x)
1880  CALL col2 (y,bm1,ntot1)
1881  CALL dssum (y,lx1,ly1,lz1)
1882  IF (imesh.EQ.1) CALL col2 (y,binvm1,ntot1)
1883  IF (imesh.EQ.2) CALL col2 (y,bintm1,ntot1)
1884  CALL col2 (y,mask,ntot1)
1885 C
1886  return
1887  END
1888 C
1889  subroutine velchar (vel,vn,vlag,nbd,tau,dtbd)
1890 C-----------------------------------------------------------------------
1891 C
1892 C Compute linearized velocity field.
1893 C
1894 C-----------------------------------------------------------------------
1895  include 'SIZE'
1896  REAL VEL (LX1,LY1,LZ1,LELV)
1897  REAL VN (LX1,LY1,LZ1,LELV)
1898  REAL VLAG (LX1,LY1,LZ1,LELV,9)
1899  REAL DTBD (NBD)
1900 C
1901  ntot1 = lx1*ly1*lz1*nelv
1902  IF (nbd.EQ.1) THEN
1903  CALL copy (vel,vn,ntot1)
1904  return
1905  ELSEIF (nbd.EQ.2) THEN
1906  c1 = tau/dtbd(2)
1907  c2 = 1.-c1
1908  CALL add3s2 (vel,vn,vlag,c1,c2,ntot1)
1909  ELSEIF (nbd.EQ.3) THEN
1910  f1 = tau**2-dtbd(3)*tau
1911  f2 = tau**2-(dtbd(2)+dtbd(3))*tau
1912  f3 = dtbd(2)*dtbd(3)
1913  f4 = dtbd(2)*(dtbd(2)+dtbd(3))
1914  r1 = f2/f3
1915  r2 = f1/f4
1916  c1 = r2
1917  c2 = -r1
1918  c3 = 1+r1-r2
1919  CALL add3s2 (vel,vlag(1,1,1,1,1),vlag(1,1,1,1,2),c2,c3,ntot1)
1920  CALL add2s2 (vel,vn,c1,ntot1)
1921  ELSE
1922  WRITE (6,*) 'Need higher order expansion in VELCHAR'
1923  call exitt
1924  ENDIF
1925 C
1926  return
1927  END
1928 C
1929  subroutine lagvel
1930 C-----------------------------------------------------------------------
1931 C
1932 C Keep old velocity field(s)
1933 C
1934 C-----------------------------------------------------------------------
1935  include 'SIZE'
1936  include 'INPUT'
1937  include 'SOLN'
1938  include 'TSTEP'
1939 C
1940  ntot1 = lx1*ly1*lz1*nelv
1941 C
1942 c DO 100 ILAG=NBDINP-1,2,-1
1943  DO 100 ilag=3-1,2,-1
1944  CALL copy (vxlag(1,1,1,1,ilag),vxlag(1,1,1,1,ilag-1),ntot1)
1945  CALL copy (vylag(1,1,1,1,ilag),vylag(1,1,1,1,ilag-1),ntot1)
1946  IF (ldim.EQ.3)
1947  $ CALL copy (vzlag(1,1,1,1,ilag),vzlag(1,1,1,1,ilag-1),ntot1)
1948  100 CONTINUE
1949 C
1950  CALL opcopy (vxlag,vylag,vzlag,vx,vy,vz)
1951 C
1952  return
1953  END
1954 C
1955  subroutine hypmsk3 (hv1msk,hv2msk,hv3msk)
1956 C---------------------------------------------------------------------
1957 C
1958 C Generate mask-array for the hyperbolic system (velocity).
1959 C
1960 C---------------------------------------------------------------------
1961  include 'SIZE'
1962  include 'INPUT'
1963  include 'GEOM'
1964  include 'SOLN'
1965  include 'TSTEP'
1966  REAL HV1MSK (LX1,LY1,LZ1,1)
1967  REAL HV2MSK (LX1,LY1,LZ1,1)
1968  REAL HV3MSK (LX1,LY1,LZ1,1)
1969  CHARACTER CB*3
1970  parameter(lxyz1=lx1*ly1*lz1)
1971  COMMON /ctmp1/ work(lxyz1,lelt)
1972 C
1973  nfaces= 2*ldim
1974  ntot1 = lx1*ly1*lz1*nelv
1975  CALL rzero (work ,ntot1)
1976  CALL rone (hv1msk,ntot1)
1977 C
1978  IF (ifield.EQ.1) THEN
1979  DO 100 ie=1,nelv
1980  DO 100 iface=1,nfaces
1981  cb=cbc(iface,ie,ifield)
1982  IF (cb(1:1).EQ.'V' .OR. cb(1:1).EQ.'v') THEN
1983  CALL faccl3 (work(1,ie),vx(1,1,1,ie),unx(1,1,iface,ie),iface)
1984  CALL faddcl3(work(1,ie),vy(1,1,1,ie),uny(1,1,iface,ie),iface)
1985  IF (if3d)
1986  $ CALL faddcl3(work(1,ie),vz(1,1,1,ie),unz(1,1,iface,ie),iface)
1987  CALL fcaver (vaver,work,ie,iface)
1988 C
1989  IF (vaver.LT.0.) CALL facev (hv1msk,ie,iface,0.0,lx1,ly1,lz1)
1990  ENDIF
1991  IF (cb(1:2).EQ.'WS' .OR. cb(1:2).EQ.'ws')
1992  $ CALL facev (hv1msk,ie,iface,0.0,lx1,ly1,lz1)
1993  100 CONTINUE
1994  ENDIF
1995 C
1996  CALL copy (hv2msk,hv1msk,ntot1)
1997  CALL copy (hv3msk,hv1msk,ntot1)
1998 C
1999  return
2000  END
2001 C
2002  subroutine setordbd
2003 C----------------------------------------------------------------------
2004 C
2005 C Set up parameters for backward differentiation scheme.
2006 C
2007 C----------------------------------------------------------------------
2008  include 'SIZE'
2009  include 'INPUT'
2010  include 'TSTEP'
2011 C
2012 c IF (IFSPLIT .OR. NBDINP.EQ.0) THEN undid hardwire, 3/6/92 pff
2013  IF ( nbdinp.LT.1) THEN
2014  nbd = 1
2015  ELSE
2016  IF ((istep.EQ.0).OR.(istep.EQ.1)) nbd = 1
2017  IF ((istep.GT.1).AND.(istep.LE.nbdinp)) nbd = istep
2018  IF (istep.GT.nbdinp) nbd = nbdinp
2019  ENDIF
2020 C
2021  return
2022  END
2023 C
2024 c-----------------------------------------------------------------------
2025  subroutine normsc (h1,semi,l2,linf,x,imesh)
2026 C---------------------------------------------------------------
2027 C
2028 C Compute error norms of a (scalar) field variable X
2029 C defined on mesh 1 or mesh 2.
2030 C The error norms are normalized with respect to the volume
2031 C (except for Linf).
2032 C
2033 C---------------------------------------------------------------
2034  include 'SIZE'
2035  include 'MASS'
2036 C
2037  REAL X (LX1,LY1,LZ1,1)
2038  COMMON /scrnrm/ y(lx1,ly1,lz1,lelt)
2039  $ ,ta1(lx1,ly1,lz1,lelt)
2040  $ ,ta2(lx1,ly1,lz1,lelt)
2041  REAL H1,SEMI,L2,LINF
2042  REAL LENGTH
2043 C
2044  IF (imesh.EQ.1) THEN
2045  nel = nelv
2046  vol = volvm1
2047  ELSEIF (imesh.EQ.2) THEN
2048  nel = nelt
2049  vol = voltm1
2050  ENDIF
2051  length = vol**(1./(ldim))
2052  nxyz1 = lx1*ly1*lz1
2053  ntot1 = nxyz1*nel
2054 C
2055  h1 = 0.
2056  semi = 0.
2057  l2 = 0.
2058  linf = 0.
2059 C
2060  linf = glamax(x,ntot1)
2061 C
2062  CALL col3 (ta1,x,x,ntot1)
2063  CALL col2 (ta1,bm1,ntot1)
2064  l2 = glsum(ta1,ntot1)
2065  IF (l2.LT.0.0) l2 = 0.
2066 C
2067  CALL rone (ta1,ntot1)
2068  CALL rzero (ta2,ntot1)
2069  CALL axhelm (y,x,ta1,ta2,imesh,1)
2070  CALL col3 (ta1,y,x,ntot1)
2071  semi = glsum(ta1,ntot1)
2072  IF (semi.LT.0.0) semi = 0.
2073 C
2074  h1 = sqrt((semi*length**2+l2)/vol)
2075  semi = sqrt(semi/vol)
2076  l2 = sqrt(l2/vol)
2077  IF (h1.LT.0.) h1 = 0.
2078 C
2079  return
2080  END
2081 C
2082  subroutine normvc (h1,semi,l2,linf,x1,x2,x3)
2083 C---------------------------------------------------------------
2084 C
2085 C Compute error norms of a (vector) field variable (X1,X2,X3)
2086 C defined on mesh 1 (velocity mesh only !)
2087 C The error norms are normalized with respect to the volume
2088 C (except for Linf).
2089 C
2090 C---------------------------------------------------------------
2091  include 'SIZE'
2092  include 'MASS'
2093 C
2094  REAL X1 (LX1,LY1,LZ1,1)
2095  REAL X2 (LX1,LY1,LZ1,1)
2096  REAL X3 (LX1,LY1,LZ1,1)
2097  COMMON /scrmg/ y1(lx1,ly1,lz1,lelt)
2098  $ ,y2(lx1,ly1,lz1,lelt)
2099  $ ,y3(lx1,ly1,lz1,lelt)
2100  $ ,ta1(lx1,ly1,lz1,lelt)
2101  COMMON /scrch/ ta2(lx1,ly1,lz1,lelt)
2102  REAL H1,SEMI,L2,LINF
2103  REAL LENGTH
2104 C
2105  imesh = 1
2106  nel = nelv
2107  vol = volvm1
2108  length = vol**(1./(ldim))
2109  nxyz1 = lx1*ly1*lz1
2110  ntot1 = nxyz1*nel
2111 C
2112  h1 = 0.
2113  semi = 0.
2114  l2 = 0.
2115  linf = 0.
2116 C
2117  CALL col3 (ta1,x1,x1,ntot1)
2118  CALL col3 (ta2,x2,x2,ntot1)
2119  CALL add2 (ta1,ta2,ntot1)
2120  IF (ldim.EQ.3) THEN
2121  CALL col3 (ta2,x3,x3,ntot1)
2122  CALL add2 (ta1,ta2,ntot1)
2123  ENDIF
2124  linf = glamax(ta1,ntot1)
2125  linf = sqrt( linf )
2126 C
2127  CALL col3 (ta1,x1,x1,ntot1)
2128  CALL col3 (ta2,x2,x2,ntot1)
2129  CALL add2 (ta1,ta2,ntot1)
2130  IF (ldim.EQ.3) THEN
2131  CALL col3 (ta2,x3,x3,ntot1)
2132  CALL add2 (ta1,ta2,ntot1)
2133  ENDIF
2134  CALL col2 (ta1,bm1,ntot1)
2135  l2 = glsum(ta1,ntot1)
2136  IF (l2.LT.0.0) l2 = 0.
2137 C
2138  CALL rone (ta1,ntot1)
2139  CALL rzero (ta2,ntot1)
2140  CALL ophx (y1,y2,y3,x1,x2,x3,ta1,ta2)
2141  CALL col3 (ta1,y1,x1,ntot1)
2142  CALL col3 (ta2,y2,x2,ntot1)
2143  CALL add2 (ta1,ta2,ntot1)
2144  IF (ldim.EQ.3) THEN
2145  CALL col3 (ta2,y3,x3,ntot1)
2146  CALL add2 (ta1,ta2,ntot1)
2147  ENDIF
2148  semi = glsum(ta1,ntot1)
2149  IF (semi.LT.0.0) semi = 0.
2150 C
2151  h1 = sqrt((semi*length**2+l2)/vol)
2152  semi = sqrt(semi/vol)
2153  l2 = sqrt(l2 /vol)
2154  IF (h1.LT.0.) h1 = 0.
2155 C
2156  return
2157  END
2158 C
2159  subroutine genwp (wp,wm2,p)
2160 C------------------------------------------------------------------
2161 C
2162 C Collocate the weights on mesh M2 with the pressure or the
2163 C search direction in the cg-iteration.
2164 C
2165 C-----------------------------------------------------------------
2166  include 'SIZE'
2167 C
2168  REAL WP (LX2,LY2,LZ2,LELV)
2169  REAL P (LX2,LY2,LZ2,LELV)
2170  REAL WM2 (LX2,LY2,LZ2)
2171 C
2172  nxyz2 = lx2*ly2*lz2
2173  DO 100 iel=1,nelv
2174  CALL col3 (wp(1,1,1,iel),wm2(1,1,1),p(1,1,1,iel),nxyz2)
2175  100 CONTINUE
2176  return
2177  END
2178 C
2179  subroutine convuz (ifstuz)
2180 C--------------------------------------------------------------------
2181 C
2182 C Check convergence for the coupled form.
2183 C Consistent approximation spaces.
2184 C
2185 C--------------------------------------------------------------------
2186  include 'SIZE'
2187  include 'TOTAL'
2188  LOGICAL IFSTUZ
2189  COMMON /scrns/ resv1(lx1,ly1,lz1,lelv)
2190  $ , resv2(lx1,ly1,lz1,lelv)
2191  $ , resv3(lx1,ly1,lz1,lelv)
2192  $ , resp(lx2,ly2,lz2,lelv)
2193  $ , ta1(lx1,ly1,lz1,lelv)
2194  $ , ta2(lx1,ly1,lz1,lelv)
2195  $ , ta3(lx1,ly1,lz1,lelv)
2196  COMMON /scrmg/ tb1(lx1,ly1,lz1,lelv)
2197  $ , tb2(lx1,ly1,lz1,lelv)
2198  $ , tb3(lx1,ly1,lz1,lelv)
2199  $ , wp(lx2,ly2,lz2,lelv)
2200  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
2201  $ , h2(lx1,ly1,lz1,lelv)
2202 C
2203  ifstuz = .true.
2204  tcritv = tolhv*1.5
2205  tcritp = tolps*1.5
2206  ntot1 = lx1*ly1*lz1*nelv
2207  ntot2 = lx2*ly2*lz2*nelv
2208  intype = 0
2209  CALL sethlm (h1,h2,intype)
2210 C
2211 C Momentum
2212 C
2213  CALL opcopy (resv1,resv2,resv3,bfx,bfy,bfz)
2214  CALL opgradt (ta1,ta2,ta3,pr)
2215  CALL ophx (tb1,tb2,tb3,vx,vy,vz,h1,h2)
2216  CALL opadd2 (resv1,resv2,resv3,ta1,ta2,ta3)
2217  CALL opsub2 (resv1,resv2,resv3,tb1,tb2,tb3)
2218  CALL opmask (resv1,resv2,resv3)
2219  CALL opdssum (resv1,resv2,resv3)
2220  CALL opcolv3 (ta1,ta2,ta3,resv1 ,resv2 ,resv3,binvm1)
2221  rv1 = sqrt(glsc3(ta1,resv1,vmult,ntot1)/volvm1)
2222  rv2 = sqrt(glsc3(ta2,resv2,vmult,ntot1)/volvm1)
2223  IF (rv1 .GT. tcritv) ifstuz = .false.
2224  IF (rv2 .GT. tcritv) ifstuz = .false.
2225  IF (ldim.EQ.3) THEN
2226  rv3 = sqrt(glsc3(ta3,resv3,vmult,ntot1)/volvm1)
2227  IF (rv3 .GT. tcritv) ifstuz = .false.
2228  ENDIF
2229 C
2230 C Continuity
2231 C
2232  CALL opdiv (resp,vx,vy,vz)
2233  CALL col3 (wp,resp,bm2inv,ntot2)
2234  rp = sqrt(glsc2(wp,resp,ntot2)/volvm2)
2235  IF (rp .GT. tcritp) ifstuz = .false.
2236 C
2237  return
2238  END
2239 C
2240  subroutine convsp (ifstsp)
2241  LOGICAL IFSTSP
2242  ifstsp = .false.
2243  return
2244  END
2245 C
2246  subroutine antimsk (y,x,xmask,n)
2247 C------------------------------------------------------------------
2248 C
2249 C Return Dirichlet boundary values of X in the array Y
2250 C
2251 C-------------------------------------------------------------------
2252  REAL Y(1),X(1),XMASK(1)
2253  include 'OPCTR'
2254 C
2255 #ifdef TIMER
2256  if (isclld.eq.0) then
2257  isclld=1
2258  nrout=nrout+1
2259  myrout=nrout
2260  rname(myrout) = 'antmsk'
2261  endif
2262  isbcnt = 2*n
2263  dct(myrout) = dct(myrout) + (isbcnt)
2264  ncall(myrout) = ncall(myrout) + 1
2265  dcount = dcount + (isbcnt)
2266 #endif
2267 C
2268  DO 100 i=1,n
2269  y(i) = x(i)*(1.-xmask(i))
2270  100 CONTINUE
2271  return
2272  END
2273 C
2274  subroutine opamask (vbdry1,vbdry2,vbdry3)
2275 C----------------------------------------------------------------------
2276 C
2277 C Antimask the velocity arrays.
2278 C
2279 C----------------------------------------------------------------------
2280  include 'SIZE'
2281  include 'INPUT'
2282  include 'SOLN'
2283  include 'TSTEP'
2284  REAL VBDRY1 (LX1,LY1,LZ1,1)
2285  REAL VBDRY2 (LX1,LY1,LZ1,1)
2286  REAL VBDRY3 (LX1,LY1,LZ1,1)
2287 C
2288  ntot1 = lx1*ly1*lz1*nelv
2289 C
2290  IF (ifstrs) THEN
2291  if (ifield.eq.ifldmhd) then
2292  CALL amask (vbdry1,vbdry2,vbdry3,bx,by,bz,nelv)
2293  else
2294  CALL amask (vbdry1,vbdry2,vbdry3,vx,vy,vz,nelv)
2295  endif
2296  ELSE
2297  if (ifield.eq.ifldmhd) then
2298  CALL antimsk (vbdry1,bx,b1mask,ntot1)
2299  CALL antimsk (vbdry2,by,b2mask,ntot1)
2300  IF (ldim.EQ.3)
2301  $ CALL antimsk (vbdry3,bz,b3mask,ntot1)
2302  else
2303  CALL antimsk (vbdry1,vx,v1mask,ntot1)
2304  CALL antimsk (vbdry2,vy,v2mask,ntot1)
2305  IF (ldim.EQ.3)
2306  $ CALL antimsk (vbdry3,vz,v3mask,ntot1)
2307  endif
2308  ENDIF
2309 C
2310  return
2311  END
2312 C
2313  subroutine opmask (res1,res2,res3)
2314 C----------------------------------------------------------------------
2315 C
2316 C Mask the residual arrays.
2317 C
2318 C----------------------------------------------------------------------
2319  include 'SIZE'
2320  include 'INPUT'
2321  include 'SOLN'
2322  include 'TSTEP'
2323  REAL RES1(1),RES2(1),RES3(1)
2324 C
2325  ntot1 = lx1*ly1*lz1*nelv
2326 C
2327 c sv=glsum(v3mask,ntot1)
2328 c sb=glsum(b3mask,ntot1)
2329 c write(6,*) istep,' ifld:',ifield,intype,sv,sb
2330  IF (ifstrs) THEN
2331  CALL rmask (res1,res2,res3,nelv)
2332  ELSE
2333  if (ifield.eq.ifldmhd) then
2334  CALL col2 (res1,b1mask,ntot1)
2335  CALL col2 (res2,b2mask,ntot1)
2336  IF (ldim.EQ.3)
2337  $ CALL col2 (res3,b3mask,ntot1)
2338  else
2339  CALL col2 (res1,v1mask,ntot1)
2340  CALL col2 (res2,v2mask,ntot1)
2341  IF (ldim.EQ.3)
2342  $ CALL col2 (res3,v3mask,ntot1)
2343  endif
2344  ENDIF
2345 C
2346  return
2347  END
2348 C
2349  subroutine opadd2 (a1,a2,a3,b1,b2,b3)
2350  include 'SIZE'
2351  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1)
2352  ntot1=lx1*ly1*lz1*nelv
2353  CALL add2(a1,b1,ntot1)
2354  CALL add2(a2,b2,ntot1)
2355  IF(ldim.EQ.3)CALL add2(a3,b3,ntot1)
2356  return
2357  END
2358 C
2359  subroutine opsub2 (a1,a2,a3,b1,b2,b3)
2360  include 'SIZE'
2361  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1)
2362  ntot1=lx1*ly1*lz1*nelv
2363  CALL sub2(a1,b1,ntot1)
2364  CALL sub2(a2,b2,ntot1)
2365  IF(ldim.EQ.3)CALL sub2(a3,b3,ntot1)
2366  return
2367  END
2368 C
2369  subroutine opsub3 (a1,a2,a3,b1,b2,b3,c1,c2,c3)
2370  include 'SIZE'
2371  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1),C1(1),C2(1),C3(1)
2372  ntot1=lx1*ly1*lz1*nelv
2373  CALL sub3(a1,b1,c1,ntot1)
2374  CALL sub3(a2,b2,c2,ntot1)
2375  IF(ldim.EQ.3)CALL sub3(a3,b3,c3,ntot1)
2376  return
2377  END
2378 C
2379  subroutine opcolv3(a1,a2,a3,b1,b2,b3,c)
2380  include 'SIZE'
2381  REAL A1(1),A2(1),A3(1)
2382  REAL B1(1),B2(1),B3(1)
2383  REAL C (1)
2384  include 'OPCTR'
2385 C
2386  ntot1=lx1*ly1*lz1*nelv
2387 
2388 #ifdef TIMER
2389  if (isclld.eq.0) then
2390  isclld=1
2391  nrout=nrout+1
2392  myrout=nrout
2393  rname(myrout) = 'opcolv'
2394  endif
2395 C
2396  isbcnt = ntot1*ldim
2397  dct(myrout) = dct(myrout) + (isbcnt)
2398  ncall(myrout) = ncall(myrout) + 1
2399  dcount = dcount + (isbcnt)
2400 #endif
2401 C
2402  IF (ldim.EQ.3) THEN
2403  DO 100 i=1,ntot1
2404  a1(i)=b1(i)*c(i)
2405  a2(i)=b2(i)*c(i)
2406  a3(i)=b3(i)*c(i)
2407  100 CONTINUE
2408  ELSE
2409  DO 200 i=1,ntot1
2410  a1(i)=b1(i)*c(i)
2411  a2(i)=b2(i)*c(i)
2412  200 CONTINUE
2413  ENDIF
2414  return
2415  END
2416 C
2417  subroutine opcolv (a1,a2,a3,c)
2418  include 'SIZE'
2419  REAL A1(1),A2(1),A3(1),C(1)
2420  include 'OPCTR'
2421 C
2422  ntot1=lx1*ly1*lz1*nelv
2423 
2424 #ifdef TIMER
2425  if (isclld.eq.0) then
2426  isclld=1
2427  nrout=nrout+1
2428  myrout=nrout
2429  rname(myrout) = 'opcolv'
2430  endif
2431 C
2432  isbcnt = ntot1*ldim
2433  dct(myrout) = dct(myrout) + (isbcnt)
2434  ncall(myrout) = ncall(myrout) + 1
2435  dcount = dcount + (isbcnt)
2436 #endif
2437 C
2438  IF (ldim.EQ.3) THEN
2439  DO 100 i=1,ntot1
2440  a1(i)=a1(i)*c(i)
2441  a2(i)=a2(i)*c(i)
2442  a3(i)=a3(i)*c(i)
2443  100 CONTINUE
2444  ELSE
2445  DO 200 i=1,ntot1
2446  a1(i)=a1(i)*c(i)
2447  a2(i)=a2(i)*c(i)
2448  200 CONTINUE
2449  ENDIF
2450  return
2451  END
2452 C
2453  subroutine opcol2 (a1,a2,a3,b1,b2,b3)
2454  include 'SIZE'
2455  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1)
2456  ntot1=lx1*ly1*lz1*nelv
2457  CALL col2(a1,b1,ntot1)
2458  CALL col2(a2,b2,ntot1)
2459  IF(ldim.EQ.3)CALL col2(a3,b3,ntot1)
2460  return
2461  END
2462 C
2463  subroutine opchsgn (a,b,c)
2464  include 'SIZE'
2465  REAL A(1),B(1),C(1)
2466  ntot1=lx1*ly1*lz1*nelv
2467  CALL chsign(a,ntot1)
2468  CALL chsign(b,ntot1)
2469  IF(ldim.EQ.3)CALL chsign(c,ntot1)
2470  return
2471  END
2472 c
2473  subroutine opcopy (a1,a2,a3,b1,b2,b3)
2474  include 'SIZE'
2475  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1)
2476  ntot1=lx1*ly1*lz1*nelv
2477  CALL copy(a1,b1,ntot1)
2478  CALL copy(a2,b2,ntot1)
2479  IF(ldim.EQ.3)CALL copy(a3,b3,ntot1)
2480  return
2481  END
2482 
2483 c-----------------------------------------------------------------------
2484  subroutine rotate_cyc(r1,r2,r3,idir)
2485 
2486  include 'SIZE'
2487  include 'GEOM'
2488  include 'INPUT'
2489  include 'PARALLEL'
2490  include 'TSTEP'
2491 
2492  real r1(lx1,ly1,lz1,1)
2493  $ , r2(lx1,ly1,lz1,1)
2494  $ , r3(lx1,ly1,lz1,1)
2495 
2496  integer e,f
2497  logical ifxy
2498  real length
2499 
2500 c (1) Face n-t transformation
2501 
2502 
2503  nface = 2*ldim
2504  do e=1,nelfld(ifield)
2505  do f=1,nface
2506 
2507  if(cbc(f,e,ifield) .eq. 'P '.or.cbc(f,e,ifield).eq.'p ')then
2508 
2509  call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f)
2510  if (idir.eq.1) then
2511  k=0
2512  do j2=js2,jf2,jskip2
2513  do j1=js1,jf1,jskip1
2514  k=k+1
2515 
2516  dotprod = unx(k,1,f,e)*ym1(j1,j2,1,e)
2517  $ -uny(k,1,f,e)*xm1(j1,j2,1,e)
2518 
2519  ifxy = .true.
2520 
2521  length = unx(k,1,f,e)*unx(k,1,f,e)
2522  $ + uny(k,1,f,e)*uny(k,1,f,e)
2523  length = sqrt(length)
2524 
2525  cost = unx(k,1,f,e)/length
2526  sint = uny(k,1,f,e)/length
2527 
2528  rnor = ( r1(j1,j2,1,e)*cost + r2(j1,j2,1,e)*sint )
2529  rtn1 = (-r1(j1,j2,1,e)*sint + r2(j1,j2,1,e)*cost )
2530 
2531  if (ifxy.and.dotprod .ge. 0.0) then
2532  r1(j1,j2,1,e) = rnor
2533  r2(j1,j2,1,e) = rtn1
2534  elseif (ifxy) then
2535  r1(j1,j2,1,e) =-rnor
2536  r2(j1,j2,1,e) =-rtn1
2537  endif
2538  enddo
2539  enddo
2540 
2541  else ! reverse rotate
2542 
2543  k=0
2544  do j2=js2,jf2,jskip2
2545  do j1=js1,jf1,jskip1
2546  k=k+1
2547 
2548  dotprod = unx(k,1,f,e)*ym1(j1,j2,1,e)
2549  $ -uny(k,1,f,e)*xm1(j1,j2,1,e)
2550  ifxy = .true.
2551 
2552  length = unx(k,1,f,e)*unx(k,1,f,e)
2553  $ + uny(k,1,f,e)*uny(k,1,f,e)
2554  length = sqrt(length)
2555 
2556  cost = unx(k,1,f,e)/length
2557  sint = uny(k,1,f,e)/length
2558 
2559  rnor = ( r1(j1,j2,1,e)*cost - r2(j1,j2,1,e)*sint )
2560  rtn1 = ( r1(j1,j2,1,e)*sint + r2(j1,j2,1,e)*cost )
2561 
2562  if(ifxy.and.dotprod .ge. 0.0) then
2563  r1(j1,j2,1,e) = rnor
2564  r2(j1,j2,1,e) = rtn1
2565  elseif (ifxy) then
2566  r1(j1,j2,1,e) =-rnor
2567  r2(j1,j2,1,e) =-rtn1
2568  endif
2569  enddo
2570  enddo
2571  endif
2572 
2573  endif
2574 
2575  enddo
2576  enddo
2577 
2578  return
2579  end
2580 c-----------------------------------------------------------------------
2581  subroutine opdssum (a,b,c)! NOTE: opdssum works on FLUID/MHD arrays only!
2582 
2583  include 'SIZE'
2584  include 'INPUT'
2585  include 'PARALLEL'
2586  include 'TSTEP'
2587  include 'GEOM'
2588 
2589  real a(1),b(1),c(1)
2590 
2591  if (ifcyclic) then
2592  call rotate_cyc (a,b,c,1)
2593  call vec_dssum (a,b,c,lx1,ly1,lz1)
2594  call rotate_cyc (a,b,c,0)
2595  else
2596  call vec_dssum (a,b,c,lx1,ly1,lz1)
2597  endif
2598 
2599  return
2600  end
2601 c-----------------------------------------------------------------------
2602  subroutine opdsop (a,b,c,op)! opdsop works on FLUID/MHD arrays only!
2603 
2604  include 'SIZE'
2605  include 'INPUT'
2606  include 'PARALLEL'
2607  include 'TSTEP'
2608  include 'GEOM'
2609 
2610  real a(1),b(1),c(1)
2611  character*3 op
2612 
2613  if (ifcyclic) then
2614 
2615  if (op.eq.'* ' .or. op.eq.'mul' .or. op.eq.'MUL') then
2616  call vec_dsop (a,b,c,lx1,ly1,lz1,op)
2617  else
2618  call rotate_cyc (a,b,c,1)
2619  call vec_dsop (a,b,c,lx1,ly1,lz1,op)
2620  call rotate_cyc (a,b,c,0)
2621  endif
2622 
2623  else
2624 
2625  call vec_dsop (a,b,c,lx1,ly1,lz1,op)
2626 
2627  endif
2628 
2629  return
2630  end
2631 c-----------------------------------------------------------------------
2632  subroutine opicol2 (a1,a2,a3,b1,b2,b3)
2633  include 'SIZE'
2634  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1)
2635  ntot1=lx1*ly1*lz1*nelv
2636  CALL invcol2(a1,b1,ntot1)
2637  CALL invcol2(a2,b2,ntot1)
2638  IF(ldim.EQ.3)CALL invcol2(a3,b3,ntot1)
2639  return
2640  END
2641 C
2642  subroutine oprzero (a,b,c)
2643  include 'SIZE'
2644  REAL A(1),B(1),C(1)
2645  ntot1=lx1*ly1*lz1*nelv
2646  CALL rzero(a,ntot1)
2647  CALL rzero(b,ntot1)
2648  IF(ldim.EQ.3) CALL rzero(c,ntot1)
2649  return
2650  END
2651 C
2652  subroutine oprone (a,b,c)
2653  include 'SIZE'
2654  REAL A(1),B(1),C(1)
2655  ntot1=lx1*ly1*lz1*nelv
2656  CALL rone(a,ntot1)
2657  CALL rone(b,ntot1)
2658  IF(ldim.EQ.3) CALL rone(c,ntot1)
2659  return
2660  END
2661 C
2662  subroutine opcmult (a,b,c,const)
2663  include 'SIZE'
2664  REAL A(1),B(1),C(1)
2665  ntot1=lx1*ly1*lz1*nelv
2666  CALL cmult(a,const,ntot1)
2667  CALL cmult(b,const,ntot1)
2668  IF(ldim.EQ.3) CALL cmult(c,const,ntot1)
2669  return
2670  END
2671 c-----------------------------------------------------------------------
2672  subroutine opcolv2c(a1,a2,a3,b1,b2,c)
2673  include 'SIZE'
2674  REAL A1(1),A2(1),A3(1)
2675  REAL B1(1),B2(1)
2676  include 'OPCTR'
2677 C
2678  ntot1=lx1*ly1*lz1*nelv
2679 
2680 #ifdef TIMER
2681  if (isclld.eq.0) then
2682  isclld=1
2683  nrout=nrout+1
2684  myrout=nrout
2685  rname(myrout) = 'opcv2c'
2686  endif
2687 C
2688  isbcnt = ntot1*(ldim+2)
2689  dct(myrout) = dct(myrout) + (isbcnt)
2690  ncall(myrout) = ncall(myrout) + 1
2691  dcount = dcount + (isbcnt)
2692 #endif
2693 C
2694  IF (ldim.EQ.3) THEN
2695  DO 100 i=1,ntot1
2696  tmp = c*b1(i)*b2(i)
2697  a1(i)=a1(i)*tmp
2698  a2(i)=a2(i)*tmp
2699  a3(i)=a3(i)*tmp
2700  100 CONTINUE
2701  ELSE
2702  DO 200 i=1,ntot1
2703  tmp = c*b1(i)*b2(i)
2704  a1(i)=a1(i)*tmp
2705  a2(i)=a2(i)*tmp
2706  200 CONTINUE
2707  ENDIF
2708  return
2709  END
2710 c-----------------------------------------------------------------------
2711  subroutine opcolv2(a1,a2,a3,b1,b2)
2712  include 'SIZE'
2713  REAL A1(1),A2(1),A3(1)
2714  REAL B1(1),B2(1)
2715  include 'OPCTR'
2716 C
2717  ntot1=lx1*ly1*lz1*nelv
2718 
2719 #ifdef TIMER
2720  if (isclld.eq.0) then
2721  isclld=1
2722  nrout=nrout+1
2723  myrout=nrout
2724  rname(myrout) = 'opclv2'
2725  endif
2726 C
2727  isbcnt = ntot1*(ldim+1)
2728  dct(myrout) = dct(myrout) + (isbcnt)
2729  ncall(myrout) = ncall(myrout) + 1
2730  dcount = dcount + (isbcnt)
2731 #endif
2732 C
2733  IF (ldim.EQ.3) THEN
2734  DO 100 i=1,ntot1
2735  tmp = b1(i)*b2(i)
2736  a1(i)=a1(i)*tmp
2737  a2(i)=a2(i)*tmp
2738  a3(i)=a3(i)*tmp
2739  100 CONTINUE
2740  ELSE
2741  DO 200 i=1,ntot1
2742  tmp = b1(i)*b2(i)
2743  a1(i)=a1(i)*tmp
2744  a2(i)=a2(i)*tmp
2745  200 CONTINUE
2746  ENDIF
2747  return
2748  END
2749 c-----------------------------------------------------------------------
2750  subroutine opadd2col(a1,a2,a3,b1,b2,b3,c)
2751  include 'SIZE'
2752  REAL A1(1),A2(1),A3(1)
2753  REAL B1(1),B2(1),B3(1),C(1)
2754  include 'OPCTR'
2755 C
2756  ntot1=lx1*ly1*lz1*nelv
2757 
2758 #ifdef TIMER
2759  if (isclld.eq.0) then
2760  isclld=1
2761  nrout=nrout+1
2762  myrout=nrout
2763  rname(myrout) = 'opa2cl'
2764  endif
2765 C
2766  isbcnt = ntot1*(ldim*2)
2767  dct(myrout) = dct(myrout) + (isbcnt)
2768  ncall(myrout) = ncall(myrout) + 1
2769  dcount = dcount + (isbcnt)
2770 #endif
2771 C
2772  IF (ldim.EQ.3) THEN
2773  DO 100 i=1,ntot1
2774  a1(i)=a1(i)+b1(i)*c(i)
2775  a2(i)=a2(i)+b2(i)*c(i)
2776  a3(i)=a3(i)+b3(i)*c(i)
2777  100 CONTINUE
2778  ELSE
2779  DO 200 i=1,ntot1
2780  a1(i)=a1(i)+b1(i)*c(i)
2781  a2(i)=a2(i)+b2(i)*c(i)
2782  200 CONTINUE
2783  ENDIF
2784  return
2785  END
2786 c-----------------------------------------------------------------------
2787  subroutine opcolv3c(a1,a2,a3,b1,b2,b3,c,d)
2788  include 'SIZE'
2789  REAL A1(1),A2(1),A3(1)
2790  REAL B1(1),B2(1),B3(1)
2791  REAL C (1)
2792  include 'OPCTR'
2793 C
2794  ntot1=lx1*ly1*lz1*nelv
2795 
2796 #ifdef TIMER
2797  if (isclld.eq.0) then
2798  isclld=1
2799  nrout=nrout+1
2800  myrout=nrout
2801  rname(myrout) = 'opcv3c'
2802  endif
2803 C
2804  isbcnt = ntot1*ldim*2
2805  dct(myrout) = dct(myrout) + (isbcnt)
2806  ncall(myrout) = ncall(myrout) + 1
2807  dcount = dcount + (isbcnt)
2808 #endif
2809 C
2810  IF (ldim.EQ.3) THEN
2811  DO 100 i=1,ntot1
2812  a1(i)=b1(i)*c(i)*d
2813  a2(i)=b2(i)*c(i)*d
2814  a3(i)=b3(i)*c(i)*d
2815  100 CONTINUE
2816  ELSE
2817  DO 200 i=1,ntot1
2818  a1(i)=b1(i)*c(i)*d
2819  a2(i)=b2(i)*c(i)*d
2820  200 CONTINUE
2821  ENDIF
2822  return
2823  END
2824 c-----------------------------------------------------------------------
2825 C
2826  subroutine uzawa (rcg,h1,h2,h2inv,intype,iter)
2827 C-----------------------------------------------------------------------
2828 C
2829 C Solve the pressure equation by (nested) preconditioned
2830 C conjugate gradient iteration.
2831 C INTYPE = 0 (steady)
2832 C INTYPE = 1 (explicit)
2833 C INTYPE = -1 (implicit)
2834 C
2835 C-----------------------------------------------------------------------
2836  include 'SIZE'
2837  include 'TOTAL'
2838  COMMON /ctolpr/ divex
2839  COMMON /cprint/ ifprint
2840  LOGICAL IFPRINT
2841  REAL RCG (LX2,LY2,LZ2,LELV)
2842  REAL H1 (LX1,LY1,LZ1,LELV)
2843  REAL H2 (LX1,LY1,LZ1,LELV)
2844  REAL H2INV(LX1,LY1,LZ1,LELV)
2845  COMMON /scruz/ wp(lx2,ly2,lz2,lelv)
2846  $ , xcg(lx2,ly2,lz2,lelv)
2847  $ , pcg(lx2,ly2,lz2,lelv)
2848  $ , rpcg(lx2,ly2,lz2,lelv)
2849 
2850  real*8 etime1,dnekclock
2851  integer*8 ntotg,nxyz2
2852 
2853 
2854  etime1 = dnekclock()
2855  divex = 0.
2856  iter = 0
2857 c
2858  CALL chktcg2 (tolps,rcg,iconv)
2859  if (param(21).gt.0.and.tolps.gt.abs(param(21)))
2860  $ tolps = abs(param(21))
2861 C
2862 c IF (ICONV.EQ.1) THEN
2863 c IF (NID.EQ.0) WRITE(6,9999) ITER,DIVEX,TOLPS
2864 c return
2865 c ENDIF
2866 
2867  nxyz2 = lx2*ly2*lz2
2868  ntot2 = nxyz2*nelv
2869  ntotg = nxyz2*nelgv
2870 
2871  CALL uzprec (rpcg,rcg,h1,h2,intype,wp)
2872  rrp1 = glsc2(rpcg,rcg,ntot2)
2873  CALL copy (pcg,rpcg,ntot2)
2874  CALL rzero (xcg,ntot2)
2875  if (rrp1.eq.0) return
2876  beta = 0.
2877  div0=0.
2878 C
2879  tolpss = tolps
2880  DO 1000 iter=1,nmxp
2881 C
2882 C CALL CONVPR (RCG,tolpss,ICONV,RNORM)
2883  call convprn (iconv,rnorm,rrp1,rcg,rpcg,tolpss)
2884 
2885  if (iter.eq.1) div0 = rnorm
2886  if (param(21).lt.0) tolpss = abs(param(21))*div0
2887 
2888  ratio = rnorm/div0
2889  IF (ifprint.AND.nio.EQ.0)
2890  $ WRITE (6,66) iter,tolpss,rnorm,div0,ratio,istep
2891  66 format(i5,1p4e12.5,i8,' Divergence')
2892 c
2893  IF (iconv.EQ.1.and.iter.gt.1) GOTO 9000
2894 c IF (ICONV.EQ.1.and.(iter.gt.1.or.istep.le.2)) GOTO 9000
2895 c IF (ICONV.EQ.1) GOTO 9000
2896 c if (ratio.le.1.e-5) goto 9000
2897 
2898 
2899  IF (iter .NE. 1) THEN
2900  beta = rrp1/rrp2
2901  CALL add2s1 (pcg,rpcg,beta,ntot2)
2902  ENDIF
2903 
2904  CALL cdabdtp (wp,pcg,h1,h2,h2inv,intype)
2905  pap = glsc2(pcg,wp,ntot2)
2906 
2907  IF (pap.NE.0.) THEN
2908  alpha = rrp1/pap
2909  ELSE
2910  pcgmx = glamax(pcg,ntot2)
2911  wp_mx = glamax(wp ,ntot2)
2912  ntot1 = lx1*ly1*lz1*nelv
2913  h1_mx = glamax(h1 ,ntot1)
2914  h2_mx = glamax(h2 ,ntot1)
2915  if (nid.eq.0) write(6,*) 'ERROR: pap=0 in uzawa.'
2916  $ ,iter,pcgmx,wp_mx,h1_mx,h2_mx
2917  call exitt
2918  ENDIF
2919  CALL add2s2 (xcg,pcg,alpha,ntot2)
2920  CALL add2s2 (rcg,wp,-alpha,ntot2)
2921 
2922  if (iter.eq.-1) then
2923  call convprn (iconv,rnrm1,rrpx,rcg,rpcg,tolpss)
2924  if (iconv.eq.1) then
2925  rnorm = rnrm1
2926  ratio = rnrm1/div0
2927  if (nio.eq.0)
2928  $ write (6,66) iter,tolpss,rnrm1,div0,ratio,istep
2929  goto 9000
2930  endif
2931  endif
2932 
2933  call ortho(rcg)
2934 
2935  rrp2 = rrp1
2936  CALL uzprec (rpcg,rcg,h1,h2,intype,wp)
2937 c RRP1 = GLSC2 (RPCG,RCG,NTOT2)
2938 
2939  1000 CONTINUE
2940  if (nid.eq.0) WRITE (6,3001) iter,rnorm,tolpss
2941 c if (istep.gt.20) CALL EMERXIT
2942  3001 FORMAT(i6,' **ERROR**: Failed to converge in UZAWA:',6e13.4)
2943  9000 CONTINUE
2944 
2945  divex = rnorm
2946  iter = iter-1
2947 
2948  if (iter.gt.0) call copy (rcg,xcg,ntot2)
2949  call ortho(rcg)
2950 
2951  etime1 = dnekclock()-etime1
2952  IF (nio.EQ.0) WRITE(6,9999) istep, ' U-Press std. ',
2953  & iter,divex,div0,tolpss,etime1
2954  9999 FORMAT(i11,a,i7,1p4e13.4)
2955 19999 FORMAT(i11,' U-Press 1.e-5: ',i7,1p4e13.4)
2956 C
2957 C
2958  return
2959  END
2960 c-----------------------------------------------------------------------
2961  subroutine mapw(md,nd,m1,n1,mflg)
2962 c
2963 c Interpolate from mesh "1" to "d" if mflg = 1
2964 c
2965  include 'SIZE'
2966  include 'DEALIAS'
2967 c
2968  real w(lxd*lxd*lx1)
2969  real md(lxd,lyd,lzd,lelv),m1(lx1,ly1,lz1,lelv)
2970 c
2971  integer icalld
2972  save icalld
2973  data icalld /0/
2974 c
2975  if (icalld .eq. 0) then
2976  call setmap(n1,nd)
2977  icalld = icalld + 1
2978  endif
2979 c
2980  if (mflg .eq.1) then
2981  do ie = 1,nelv
2982  call specmp(md(1,1,1,ie),nd,m1(1,1,1,ie),n1,im1d,im1dt,w)
2983  enddo
2984  else
2985  do ie = 1,nelv
2986  call specmp(m1(1,1,1,ie),n1,md(1,1,1,ie),nd,imd1,imd1t,w)
2987  enddo
2988  endif
2989 c
2990  return
2991  end
2992 c-----------------------------------------------------------------------
2993  subroutine mapwp(md,nd,m1,n1,mflg)
2994 c
2995 c Project from "d" to 1
2996 c
2997  include 'SIZE'
2998  include 'DEALIAS'
2999 c
3000  real w(lxd*lxd*lx1)
3001  real md(lxd,lyd,lzd,lelv),m1(lx1,ly1,lz1,lelv)
3002 c
3003  integer icalld
3004  save icalld
3005  data icalld /0/
3006 c
3007  if (icalld .eq. 0) then
3008  call setproj(n1,nd)
3009  icalld = icalld + 1
3010  endif
3011 c
3012  do ie = 1,nelv
3013  call specmp(m1(1,1,1,ie),n1,md(1,1,1,ie),nd,pmd1,pmd1t,w)
3014  enddo
3015 c
3016  return
3017  end
3018 c-----------------------------------------------------------------------
3019  subroutine specmp(b,nb,a,na,ba,ab,w)
3020 C
3021 C - Spectral interpolation from A to B via tensor products
3022 C - scratch arrays: w(na*na*nb)
3023 C
3024 C
3025  include 'SIZE'
3026  include 'INPUT'
3027  real b(nb,nb,nb),a(na,na,na)
3028  real w(1)
3029 C
3030 C
3031  if (if3d) then
3032  nab = na*nb
3033  nbb = nb*nb
3034  call mxm(ba,nb,a,na,b,na*na)
3035  k=1
3036  l=1
3037  do iz=1,na
3038  call mxm(b(k,1,1),nb,ab,na,w(l),nb)
3039  k=k+nab
3040  l=l+nbb
3041  enddo
3042  call mxm(w,nbb,ab,na,b,nb)
3043  else
3044  call mxm(ba,nb,a,na,w,na)
3045  call mxm(w,nb,ab,na,b,nb)
3046  endif
3047  return
3048  end
3049 c-----------------------------------------------------------------------
3050  subroutine setmap(n1,nd)
3051 c
3052  include 'SIZE'
3053  include 'DEALIAS'
3054 c
3055  parameter(lx=80)
3056  real z1(lx),zd(lx),w(lx)
3057 c
3058  if (n1.gt.lx.or.nd.gt.lx) then
3059  write(6,*)'ERROR: increase lx in setmap to max:',n1,nd
3060  call exitt
3061  endif
3062 c
3063  call zwgll(z1,w,n1)
3064  call zwgll(zd,w,nd)
3065  call igllm(im1d,im1dt,z1,zd,n1,nd,n1,nd)
3066  call igllm(imd1,imd1t,zd,z1,nd,n1,nd,n1)
3067 c
3068  return
3069  end
3070 c-----------------------------------------------------------------------
3071  subroutine set_pnd(P,LkD,LkNt,N,D)
3072 c
3073  integer N,D
3074  real P(N,D),LkD(0:N-1,D),LkNt(N,0:N-1)
3075 c
3076  parameter(lx=80)
3077  real zN(lx),zD(lx),w(lx)
3078 c
3079 c Compute Lagrangian interpolant points
3080 c
3081  call zwgll(zn,w,n)
3082  call zwgll(zd,w,d)
3083 c
3084 c
3085 c D D
3086 c Compute L (x ) and L (x )
3087 c k i k j
3088 c
3089  do k=0,n-1
3090 c
3091  do j=1,d
3092  lkd(k,j) = pnleg(zd(j),k)
3093  enddo
3094 c
3095  do j=1,n
3096  lknt(j,k) = pnleg(zn(j),k)
3097  enddo
3098 c
3099  enddo
3100 c
3101 c Find scale factors to normalize the first N-1 Legendre polynomials
3102 c such that (L_i,L_j) = delta_ij
3103 c
3104  do k=0,n-1
3105  s = 0
3106  do j=1,d
3107  s = s + lkd(k,j)*lkd(k,j)*w(j)
3108  enddo
3109  s = 1./sqrt(s)
3110 c
3111 c Normalize polynomials
3112 c
3113  do j=1,d
3114  lkd(k,j) = s * lkd(k,j)
3115  enddo
3116 c
3117  do j=1,n
3118  lknt(j,k) = s * lknt(j,k)
3119  enddo
3120 c
3121  enddo
3122 c
3123 c Scale columns of LkD by w_j
3124 c
3125  do j=1,d
3126  do k=0,n-1
3127  lkd(k,j) = lkd(k,j)*w(j)
3128  enddo
3129  enddo
3130 c
3131 c Compute P = LkNt * LkD
3132 c
3133  call mxm(lknt,n,lkd,n,p,d)
3134 c
3135  return
3136  end
3137 c-----------------------------------------------------------------------
3138  subroutine convop(conv,fi)
3139 C
3140 C Compute the convective term CONV for a passive scalar field FI
3141 C using the skew-symmetric formulation.
3142 C The field variable FI is defined on mesh M1 (GLL) and
3143 C the velocity field is assumed given.
3144 C
3145 C IMPORTANT NOTE: Use the scratch-arrays carefully!!!!!
3146 C
3147 C The common-block SCRNS is used in CONV1 and CONV2.
3148 C The common-blocks CTMP0 and CTMP1 are also used as scratch-arrays
3149 C since there is no direct stiffness summation or Helmholtz-solves.
3150 C
3151  include 'SIZE'
3152  include 'TOTAL'
3153  include 'CTIMER'
3154 C
3155 C Use the common blocks CTMP0 and CTMP1 as work space.
3156 C
3157  COMMON /scrch/ cmask1(lx1,ly1,lz1,lelv)
3158  $ , cmask2(lx1,ly1,lz1,lelv)
3159  COMMON /ctmp1/ mfi(lx1,ly1,lz1,lelv)
3160  $ , dmfi(lx1,ly1,lz1,lelv)
3161  $ , mdmfi(lx1,ly1,lz1,lelv)
3162  REAL MFI,DMFI,MDMFI
3163 C
3164 C Arrays in parameter list
3165 C
3166  REAL CONV (LX1,LY1,LZ1,1)
3167  REAL FI (LX1,LY1,LZ1,1)
3168 
3169  if (nio.eq.0.and.loglevel.gt.2)
3170  $ write(6,*) 'convop', ifield, ifdeal(ifield)
3171 
3172 #ifdef TIMER
3173  if (icalld.eq.0) tadvc=0.0
3174  icalld=icalld+1
3175  nadvc=icalld
3176  etime1=dnekclock()
3177 #endif
3178 
3179  nxyz1 = lx1*ly1*lz1
3180  ntot1 = lx1*ly1*lz1*nelv
3181  ntotz = lx1*ly1*lz1*nelfld(ifield)
3182  ntott = lx1*ly1*lz1*nelt
3183 
3184  call rzero (conv,ntott)
3185 
3186  if (ifdgfld(ifield)) then
3187  call convect_dg (conv,fi,.false.,vxd,vyd,vzd,.true.)
3188  goto 100
3189  elseif (param(86).ne.0.0) then ! skew-symmetric form
3190  call convopo(conv,fi)
3191  goto 100
3192  endif
3193 
3194  if (.not. ifdeal(ifield)) then
3195  call conv1 (conv,fi)
3196  elseif (param(99).eq.2.or.param(99).eq.3) then
3197  call conv1d(conv,fi)
3198  elseif (param(99).eq.4) then
3199  if (ifpert) then
3200  call convect_new (conv,fi,.false.,vx,vy,vz,.false.)
3201  else
3202  call convect_new (conv,fi,.false.,vxd,vyd,vzd,.true.)
3203  endif
3204  call invcol2 (conv,bm1,ntot1) ! local mass inverse
3205  elseif (param(99).eq.5) then
3206  call convect_cons(conv,fi,.false.,vx,vy,vz,.false.)
3207  call invcol2 (conv,bm1,ntot1) ! local mass inverse
3208  else
3209  call conv1 (conv,fi)
3210  endif
3211 
3212  100 continue
3213 
3214 #ifdef TIMER
3215  tadvc=tadvc+(dnekclock()-etime1)
3216 #endif
3217 
3218  return
3219  END
3220 c-----------------------------------------------------------------------
3221  subroutine conv1d (dfi,fi)
3222 C--------------------------------------------------------------------
3223 C
3224 C Compute D*FI (part of the convection operator)
3225 C De-aliased version 3/11/97
3226 C
3227 C--------------------------------------------------------------------
3228  include 'SIZE'
3229  include 'TOTAL'
3230  REAL DFI (LX1,LY1,LZ1,1)
3231  REAL FI (LX1,LY1,LZ1,1)
3232 c
3233  COMMON /ctmp0/ ta1(lx1,ly1,lz1,lelv)
3234  $ , dfid(lxd,lyd,lzd,lelv)
3235  $ , ta1d(lxd,lyd,lzd,lelv)
3236 C
3237  integer icalld
3238  save icalld
3239  data icalld /0/
3240 c
3241  ntotd = lxd*lyd*lzd*nelv
3242 c
3243 c
3244 c interpolate ta1 and vx onto larger mesh
3245 c
3246  CALL dudxyz (ta1,fi,rxm1,sxm1,txm1,jacm1,imesh,1)
3247  call mapw (ta1d,lxd,ta1,lx1,1)
3248  call mapw (vxd ,lxd,vx ,lx1,1)
3249  CALL col3 (dfid,ta1d,vxd,ntotd)
3250 c
3251 c
3252 c interpolate ta1 and vy onto larger mesh
3253 c
3254  CALL dudxyz (ta1,fi,rym1,sym1,tym1,jacm1,imesh,2)
3255  call mapw (ta1d,lxd,ta1,lx1,1)
3256  call mapw (vyd ,lxd,vy ,lx1,1)
3257  CALL addcol3 (dfid,ta1d,vyd,ntotd)
3258 c
3259  IF (if3d) THEN
3260 c
3261 c interpolate ta1 and vy onto larger mesh
3262 c
3263  CALL dudxyz (ta1,fi,rzm1,szm1,tzm1,jacm1,imesh,3)
3264  call mapw (ta1d,lxd,ta1,lx1,1)
3265  call mapw (vzd ,lxd,vz ,lx1,1)
3266  CALL addcol3 (dfid,ta1d,vzd,ntotd)
3267 c
3268  ENDIF
3269 c
3270 c Now, *project* DFID onto mesh 1 using L2 projection
3271 c
3272  call mapwp(dfid,lxd,dfi,lx1,-1)
3273  return
3274  END
3275 C------------------------------------------------------------------------
3276  subroutine conv1(du,u) ! used to be conv1n
3277 c
3278  include 'SIZE'
3279  include 'DXYZ'
3280  include 'INPUT'
3281  include 'GEOM'
3282  include 'SOLN'
3283  include 'TSTEP'
3284 c
3285  real du (lx1*ly1*lz1,1)
3286  real u (lx1,ly1,lz1,1)
3287 c
3288  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
3289  logical ifdfrm, iffast, ifh2, ifsolv
3290 C
3291 C Store the inverse jacobian to speed this operation up
3292 C
3293  common /ctmp0/ dudr(lx1,ly1,lz1)
3294  $ , duds(lx1,ly1,lz1)
3295  $ , dudt(lx1,ly1,lz1)
3296 
3297  nel = nelv
3298  if (imesh.eq.2) nel = nelt
3299  nxy1 = lx1*ly1
3300  nyz1 = ly1*lz1
3301  nxyz1 = lx1*ly1*lz1
3302  ntot = nxyz1*nel
3303 C
3304 C Compute vel.grad(u)
3305 C
3306  do ie=1,nel
3307 C
3308  if (if3d) then
3309 c
3310  call mxm (dxm1,lx1,u(1,1,1,ie),lx1,dudr,nyz1)
3311  do iz=1,lz1
3312  call mxm (u(1,1,iz,ie),lx1,dytm1,ly1,duds(1,1,iz),ly1)
3313  enddo
3314  call mxm (u(1,1,1,ie),nxy1,dztm1,lz1,dudt,lz1)
3315 c
3316  do i=1,nxyz1
3317  du(i,ie) = jacmi(i,ie)*(
3318  $ vx(i,1,1,ie)*(
3319  $ rxm1(i,1,1,ie)*dudr(i,1,1)
3320  $ + sxm1(i,1,1,ie)*duds(i,1,1)
3321  $ + txm1(i,1,1,ie)*dudt(i,1,1) )
3322  $ + vy(i,1,1,ie)*(
3323  $ rym1(i,1,1,ie)*dudr(i,1,1)
3324  $ + sym1(i,1,1,ie)*duds(i,1,1)
3325  $ + tym1(i,1,1,ie)*dudt(i,1,1) )
3326  $ + vz(i,1,1,ie)*(
3327  $ rzm1(i,1,1,ie)*dudr(i,1,1)
3328  $ + szm1(i,1,1,ie)*duds(i,1,1)
3329  $ + tzm1(i,1,1,ie)*dudt(i,1,1) ) )
3330  enddo
3331 c
3332  else
3333 c
3334 c 2D
3335  call mxm (dxm1,lx1,u(1,1,1,ie),lx1,dudr,nyz1)
3336  call mxm (u(1,1,1,ie),lx1,dytm1,ly1,duds,ly1)
3337  do i=1,nxyz1
3338  du(i,ie) = jacmi(i,ie)*(
3339  $ vx(i,1,1,ie)*(
3340  $ rxm1(i,1,1,ie)*dudr(i,1,1)
3341  $ + sxm1(i,1,1,ie)*duds(i,1,1) )
3342  $ + vy(i,1,1,ie)*(
3343  $ rym1(i,1,1,ie)*dudr(i,1,1)
3344  $ + sym1(i,1,1,ie)*duds(i,1,1) ) )
3345  enddo
3346  endif
3347 
3348  enddo
3349 c
3350  return
3351  end
3352 c-----------------------------------------------------------------------
3353  subroutine conv1no(du,u)
3354 c
3355  include 'SIZE'
3356  include 'DXYZ'
3357  include 'INPUT'
3358  include 'GEOM'
3359  include 'SOLN'
3360  include 'TSTEP'
3361 c
3362  real du (lx1*ly1*lz1,1)
3363  real u (lx1,ly1,lz1,1)
3364 c
3365  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
3366  logical ifdfrm, iffast, ifh2, ifsolv
3367 C
3368 C Store the inverse jacobian to speed this operation up
3369 C
3370 C
3371  common /ctmp0/ dudr(lx1,ly1,lz1)
3372  $ , duds(lx1,ly1,lz1)
3373  $ , dudt(lx1,ly1,lz1)
3374 C
3375  nel = nelv
3376  if (imesh.eq.2) nel = nelt
3377  nxy1 = lx1*ly1
3378  nyz1 = ly1*lz1
3379  nxyz1 = lx1*ly1*lz1
3380  ntot = nxyz1*nel
3381 C
3382 C Compute vel.grad(u)
3383 C
3384  do ie=1,nel
3385 C
3386  if (if3d) then
3387 c
3388  call mxm (dxm1,lx1,u(1,1,1,ie),lx1,dudr,nyz1)
3389  do iz=1,lz1
3390  call mxm (u(1,1,iz,ie),lx1,dytm1,ly1,duds(1,1,iz),ly1)
3391  enddo
3392  call mxm (u(1,1,1,ie),nxy1,dztm1,lz1,dudt,lz1)
3393 c
3394  do i=1,nxyz1
3395  du(i,ie) = jacmi(i,ie)*(
3396  $ vx(i,1,1,ie)*(
3397  $ rxm1(i,1,1,ie)*dudr(i,1,1)
3398  $ + sxm1(i,1,1,ie)*duds(i,1,1)
3399  $ + txm1(i,1,1,ie)*dudt(i,1,1) )
3400  $ + vy(i,1,1,ie)*(
3401  $ rym1(i,1,1,ie)*dudr(i,1,1)
3402  $ + sym1(i,1,1,ie)*duds(i,1,1)
3403  $ + tym1(i,1,1,ie)*dudt(i,1,1) )
3404  $ + vz(i,1,1,ie)*(
3405  $ rzm1(i,1,1,ie)*dudr(i,1,1)
3406  $ + szm1(i,1,1,ie)*duds(i,1,1)
3407  $ + tzm1(i,1,1,ie)*dudt(i,1,1) ) )
3408  enddo
3409 c
3410  else
3411 c
3412 c 2D
3413  call mxm (dxm1,lx1,u(1,1,1,ie),lx1,dudr,nyz1)
3414  call mxm (u(1,1,1,ie),lx1,dytm1,ly1,duds,ly1)
3415  do i=1,nxyz1
3416  du(i,ie) = jacmi(i,ie)*(
3417  $ vx(i,1,1,ie)*(
3418  $ rxm1(i,1,1,ie)*dudr(i,1,1)
3419  $ + sxm1(i,1,1,ie)*duds(i,1,1) )
3420  $ + vy(i,1,1,ie)*(
3421  $ rym1(i,1,1,ie)*dudr(i,1,1)
3422  $ + sym1(i,1,1,ie)*duds(i,1,1) ) )
3423  enddo
3424  endif
3425 
3426  enddo
3427 c
3428  return
3429  end
3430 c-----------------------------------------------------------------------
3431  subroutine conv1rk(du,dv,dw,u,v,w)
3432 c
3433  include 'SIZE'
3434  include 'DXYZ'
3435  include 'INPUT'
3436  include 'GEOM'
3437  include 'SOLN'
3438  include 'TSTEP'
3439 c
3440  real du(lx1*ly1*lz1,1),dv(lx1*ly1*lz1,1),dw(lx1*ly1*lz1,1)
3441  real u (lx1,ly1,lz1,1),v (lx1,ly1,lz1,1),w (lx1,ly1,lz1,1)
3442 c
3443  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
3444  logical ifdfrm, iffast, ifh2, ifsolv
3445 C
3446  common /ctmp0/ duds(lx1,ly1,lz1)
3447  $ , dvds(lx1,ly1,lz1)
3448  $ , dwds(lx1,ly1,lz1)
3449 C
3450  nel = nelv
3451  if (imesh.eq.2) nel = nelt
3452  nxy1 = lx1*ly1
3453  nyz1 = ly1*lz1
3454  nxyz1 = lx1*ly1*lz1
3455 C
3456 C Compute vel.grad(u)
3457 C
3458  do ie=1,nel
3459 C
3460  if (if3d) then
3461  call mxm (dxm1,lx1,u(1,1,1,ie),lx1,du(1,ie),nyz1)
3462  call mxm (dxm1,lx1,v(1,1,1,ie),lx1,dv(1,ie),nyz1)
3463  call mxm (dxm1,lx1,w(1,1,1,ie),lx1,dw(1,ie),nyz1)
3464  do i=1,nxyz1
3465  du(i,ie) = du(i,ie)*vx(i,1,1,ie)
3466  dv(i,ie) = dv(i,ie)*vx(i,1,1,ie)
3467  dw(i,ie) = dw(i,ie)*vx(i,1,1,ie)
3468  enddo
3469 c
3470  do iz=1,lz1
3471  call mxm (u(1,1,iz,ie),lx1,dytm1,ly1,duds(1,1,iz),ly1)
3472  enddo
3473  do iz=1,lz1
3474  call mxm (v(1,1,iz,ie),lx1,dytm1,ly1,dvds(1,1,iz),ly1)
3475  enddo
3476  do iz=1,lz1
3477  call mxm (w(1,1,iz,ie),lx1,dytm1,ly1,dwds(1,1,iz),ly1)
3478  enddo
3479  do i=1,nxyz1
3480  du(i,ie) = du(i,ie) + duds(i,1,1)*vy(i,1,1,ie)
3481  dv(i,ie) = dv(i,ie) + dvds(i,1,1)*vy(i,1,1,ie)
3482  dw(i,ie) = dw(i,ie) + dwds(i,1,1)*vy(i,1,1,ie)
3483  enddo
3484 c
3485  call mxm (u(1,1,1,ie),nxy1,dztm1,lz1,duds,lz1)
3486  call mxm (v(1,1,1,ie),nxy1,dztm1,lz1,dvds,lz1)
3487  call mxm (w(1,1,1,ie),nxy1,dztm1,lz1,dwds,lz1)
3488  do i=1,nxyz1
3489  du(i,ie) = du(i,ie) + duds(i,1,1)*vz(i,1,1,ie)
3490  dv(i,ie) = dv(i,ie) + dvds(i,1,1)*vz(i,1,1,ie)
3491  dw(i,ie) = dw(i,ie) + dwds(i,1,1)*vz(i,1,1,ie)
3492  enddo
3493  else
3494 c 2D
3495  call mxm (dxm1,lx1,u(1,1,1,ie),lx1,du(1,ie),nyz1)
3496  call mxm (dxm1,lx1,v(1,1,1,ie),lx1,dv(1,ie),nyz1)
3497  do i=1,nxyz1
3498  du(i,ie) = du(i,ie)*vx(i,1,1,ie)
3499  dv(i,ie) = dv(i,ie)*vx(i,1,1,ie)
3500  enddo
3501 c
3502  call mxm (u(1,1,1,ie),lx1,dytm1,ly1,duds,ly1)
3503  call mxm (v(1,1,1,ie),lx1,dytm1,ly1,dvds,ly1)
3504  do i=1,nxyz1
3505  du(i,ie) = du(i,ie) + duds(i,1,1)*vy(i,1,1,ie)
3506  dv(i,ie) = dv(i,ie) + dvds(i,1,1)*vy(i,1,1,ie)
3507  enddo
3508 c
3509  endif
3510 c
3511  enddo
3512 c
3513  return
3514  end
3515 c-----------------------------------------------------------------------
3516  subroutine velconvv(vxn,vyn,vzn,tau)
3517 c
3518 c Compute convecting velocity field (linearization)
3519 c
3520  include 'SIZE'
3521  include 'GEOM'
3522  include 'MASS'
3523  include 'SOLN'
3524  include 'TSTEP'
3525  real vxn(1),vyn(1),vzn(1)
3526 c
3527  include 'OPCTR'
3528 
3529 c Operation count
3530 #ifdef TIMER
3531  integer opct
3532  if (isclld.eq.0) then
3533  isclld=1
3534  nrout=nrout+1
3535  myrout=nrout
3536  rname(myrout) = 'velcvv'
3537  endif
3538  ncall(myrout) = ncall(myrout) + 1
3539  opct = 0
3540 #endif
3541 
3542  call velchar (vx,vxn,vxlag,nbd,tau,dtlag)
3543  call velchar (vy,vyn,vylag,nbd,tau,dtlag)
3544  if (ldim.eq.3)
3545  $call velchar (vz,vzn,vzlag,nbd,tau,dtlag)
3546 c
3547  ntot = lx1*ly1*lz1*nelv
3548  if (ldim.eq.3) then
3549  do i=1,ntot
3550 c
3551  tx = vx(i,1,1,1)*bm1(i,1,1,1)*jacmi(i,1)
3552  ty = vy(i,1,1,1)*bm1(i,1,1,1)*jacmi(i,1)
3553  tz = vz(i,1,1,1)*bm1(i,1,1,1)*jacmi(i,1)
3554 c
3555  vx(i,1,1,1) = tx*rxm1(i,1,1,1)
3556  $ + ty*rym1(i,1,1,1)
3557  $ + tz*rzm1(i,1,1,1)
3558  vy(i,1,1,1) = tx*sxm1(i,1,1,1)
3559  $ + ty*sym1(i,1,1,1)
3560  $ + tz*szm1(i,1,1,1)
3561  vz(i,1,1,1) = tx*txm1(i,1,1,1)
3562  $ + ty*tym1(i,1,1,1)
3563  $ + tz*tzm1(i,1,1,1)
3564  enddo
3565  else
3566  do i=1,ntot
3567 c
3568  tx = vx(i,1,1,1)*bm1(i,1,1,1)*jacmi(i,1)
3569  ty = vy(i,1,1,1)*bm1(i,1,1,1)*jacmi(i,1)
3570 c
3571  vx(i,1,1,1) = tx*rxm1(i,1,1,1)
3572  $ + ty*rym1(i,1,1,1)
3573  vy(i,1,1,1) = tx*sxm1(i,1,1,1)
3574  $ + ty*sym1(i,1,1,1)
3575  enddo
3576  endif
3577 c
3578 #ifdef TIMER
3579  if (ldim.eq.3) then
3580  opct = ntot*21
3581  else
3582  opct = ntot*10
3583  endif
3584  dct(myrout) = dct(myrout) + opct
3585  dcount = dcount + opct
3586 #endif
3587 C
3588 c
3589  return
3590  end
3591 c-----------------------------------------------------------------------
3592  subroutine frkconvv (du,dv,dw,u,v,w,mu)
3593 c
3594  include 'SIZE'
3595  include 'DXYZ'
3596  include 'INPUT'
3597  include 'GEOM'
3598  include 'SOLN'
3599  include 'MASS'
3600  include 'TSTEP'
3601 c
3602  real du(1),dv(1),dw(1)
3603  real u (1),v (1),w (1)
3604  integer mu(0:1)
3605 c
3606  include 'OPCTR'
3607 c
3608 c Operation count
3609 c
3610 #ifdef TIMER
3611  integer opct
3612  if (isclld.eq.0) then
3613  isclld=1
3614  nrout=nrout+1
3615  myrout=nrout
3616  rname(myrout) = 'frkcvv'
3617  endif
3618  ncall(myrout) = ncall(myrout) + 1
3619 #endif
3620 c
3621 c
3622 c Evaluate right-hand-side for Runge-Kutta scheme in the case of
3623 c pure convection.
3624 c
3625  ntot = lx1*ly1*lz1*nelv
3626  call conv1rk (du,dv,dw,u,v,w)
3627  CALL opdssum (du,dv,dw)
3628 c
3629  if (ldim.eq.3) then
3630  do i=1,ntot
3631  du(i) = du(i)*binvm1(i,1,1,1)
3632  dv(i) = dv(i)*binvm1(i,1,1,1)
3633  dw(i) = dw(i)*binvm1(i,1,1,1)
3634  enddo
3635  else
3636  do i=1,ntot
3637  du(i) = du(i)*binvm1(i,1,1,1)
3638  dv(i) = dv(i)*binvm1(i,1,1,1)
3639  enddo
3640  endif
3641 c
3642 c Mask
3643 c
3644  nu = mu(0)
3645  if (ldim.eq.3) then
3646  do i=1,nu
3647  du(mu(i)) = 0.
3648  dv(mu(i)) = 0.
3649  dw(mu(i)) = 0.
3650  enddo
3651  else
3652  do i=1,nu
3653  du(mu(i)) = 0.
3654  dv(mu(i)) = 0.
3655  enddo
3656  endif
3657 c
3658 #ifdef TIMER
3659  opct = ldim*ntot
3660  dct(myrout) = dct(myrout) + opct
3661  dcount = dcount + opct
3662 #endif
3663 c
3664  return
3665  end
3666 c-----------------------------------------------------------------------
3667  subroutine conv1rk2(du,dv,dw,u,v,w,cu,cv,cw,beta,wk)
3668 c
3669  include 'SIZE'
3670  include 'DXYZ'
3671  include 'INPUT'
3672  include 'GEOM'
3673  include 'SOLN'
3674  include 'TSTEP'
3675 c
3676  real du(lx1*ly1*lz1,1),dv(lx1*ly1*lz1,1),dw(lx1*ly1*lz1,1)
3677  real u (lx1,ly1,lz1,1),v (lx1,ly1,lz1,1),w (lx1,ly1,lz1,1)
3678  real cu(lx1,ly1,lz1,1),cv(lx1,ly1,lz1,1),cw(lx1,ly1,lz1,1)
3679  real wk(lx1,ly1,lz1,3)
3680 c
3681  common /ctmp0/ duds(lx1,ly1,lz1)
3682  $ , dvds(lx1,ly1,lz1)
3683  $ , dwds(lx1,ly1,lz1)
3684 C
3685  include 'OPCTR'
3686 c
3687 c Operation count
3688 c
3689 #ifdef TIMER
3690  integer opct
3691  if (isclld.eq.0) then
3692  isclld=1
3693  nrout=nrout+1
3694  myrout=nrout
3695  rname(myrout) = 'cv1rk2'
3696  endif
3697  ncall(myrout) = ncall(myrout) + 1
3698  opct = 0
3699 #endif
3700 c
3701  nel = nelv
3702  if (imesh.eq.2) nel = nelt
3703  nxy1 = lx1*ly1
3704  nyz1 = ly1*lz1
3705  nxyz1 = lx1*ly1*lz1
3706  ntot = nxyz1*nel
3707 C
3708 C Compute vel.grad(u)
3709 C
3710  do ie=1,nel
3711 C
3712  if (if3d) then
3713  do i=1,nxyz1
3714  wk(i,1,1,1)=u(i,1,1,ie)+beta*cu(i,1,1,ie)
3715  wk(i,1,1,2)=v(i,1,1,ie)+beta*cv(i,1,1,ie)
3716  wk(i,1,1,3)=w(i,1,1,ie)+beta*cw(i,1,1,ie)
3717  enddo
3718 c
3719  call mxm (dxm1,lx1,wk(1,1,1,1),lx1,du(1,ie),nyz1)
3720  call mxm (dxm1,lx1,wk(1,1,1,2),lx1,dv(1,ie),nyz1)
3721  call mxm (dxm1,lx1,wk(1,1,1,3),lx1,dw(1,ie),nyz1)
3722  do i=1,nxyz1
3723  du(i,ie) = du(i,ie)*vx(i,1,1,ie)
3724  dv(i,ie) = dv(i,ie)*vx(i,1,1,ie)
3725  dw(i,ie) = dw(i,ie)*vx(i,1,1,ie)
3726  enddo
3727 c
3728  do iz=1,lz1
3729  call mxm (wk(1,1,iz,1),lx1,dytm1,ly1,duds(1,1,iz),ly1)
3730  enddo
3731  do iz=1,lz1
3732  call mxm (wk(1,1,iz,2),lx1,dytm1,ly1,dvds(1,1,iz),ly1)
3733  enddo
3734  do iz=1,lz1
3735  call mxm (wk(1,1,iz,3),lx1,dytm1,ly1,dwds(1,1,iz),ly1)
3736  enddo
3737  do i=1,nxyz1
3738  du(i,ie) = du(i,ie) + duds(i,1,1)*vy(i,1,1,ie)
3739  dv(i,ie) = dv(i,ie) + dvds(i,1,1)*vy(i,1,1,ie)
3740  dw(i,ie) = dw(i,ie) + dwds(i,1,1)*vy(i,1,1,ie)
3741  enddo
3742 c
3743  call mxm (wk(1,1,1,1),nxy1,dztm1,lz1,duds,lz1)
3744  call mxm (wk(1,1,1,2),nxy1,dztm1,lz1,dvds,lz1)
3745  call mxm (wk(1,1,1,3),nxy1,dztm1,lz1,dwds,lz1)
3746  do i=1,nxyz1
3747  du(i,ie) = du(i,ie) + duds(i,1,1)*vz(i,1,1,ie)
3748  dv(i,ie) = dv(i,ie) + dvds(i,1,1)*vz(i,1,1,ie)
3749  dw(i,ie) = dw(i,ie) + dwds(i,1,1)*vz(i,1,1,ie)
3750  enddo
3751  else
3752 c 2D
3753  do i=1,nxyz1
3754  wk(i,1,1,1)=u(i,1,1,ie)+beta*cu(i,1,1,ie)
3755  wk(i,1,1,2)=v(i,1,1,ie)+beta*cv(i,1,1,ie)
3756  enddo
3757 c
3758  call mxm (dxm1,lx1,wk(1,1,1,1),lx1,du(1,ie),nyz1)
3759  call mxm (dxm1,lx1,wk(1,1,1,2),lx1,dv(1,ie),nyz1)
3760  do i=1,nxyz1
3761  du(i,ie) = du(i,ie)*vx(i,1,1,ie)
3762  dv(i,ie) = dv(i,ie)*vx(i,1,1,ie)
3763  enddo
3764 c
3765  call mxm (wk(1,1,1,1),lx1,dytm1,ly1,duds,ly1)
3766  call mxm (wk(1,1,1,2),lx1,dytm1,ly1,dvds,ly1)
3767  do i=1,nxyz1
3768  du(i,ie) = du(i,ie) + duds(i,1,1)*vy(i,1,1,ie)
3769  dv(i,ie) = dv(i,ie) + dvds(i,1,1)*vy(i,1,1,ie)
3770  enddo
3771  endif
3772 c
3773  enddo
3774 c
3775 #ifdef TIMER
3776  if (if3d) then
3777  opct = 21*ntot
3778  else
3779  opct = 10*ntot
3780  endif
3781 c
3782  dct(myrout) = dct(myrout) + opct
3783  dcount = dcount + opct
3784 #endif
3785 c
3786  return
3787  end
3788 c-----------------------------------------------------------------------
3789  subroutine frkconvv2(du,dv,dw,u,v,w,cu,cv,cw,beta,mu,wk)
3790 c
3791  include 'SIZE'
3792  include 'DXYZ'
3793  include 'INPUT'
3794  include 'GEOM'
3795  include 'SOLN'
3796  include 'MASS'
3797  include 'TSTEP'
3798 c
3799  real du(1),dv(1),dw(1)
3800  real u (1),v (1),w (1)
3801  real cu(1),cv(1),cw(1)
3802  real wk(lx1*ly1*lz1,3)
3803  integer mu(0:1)
3804 c
3805  include 'OPCTR'
3806 c
3807 c Operation count
3808 c
3809 #ifdef TIMER
3810  integer opct
3811  if (isclld.eq.0) then
3812  isclld=1
3813  nrout=nrout+1
3814  myrout=nrout
3815  rname(myrout) = 'frkcv2'
3816  endif
3817  ncall(myrout) = ncall(myrout) + 1
3818 #endif
3819 c
3820 c
3821 c Evaluate right-hand-side for Runge-Kutta scheme in the case of
3822 c pure convection.
3823 c
3824 C
3825  ntot = lx1*ly1*lz1*nelv
3826  call conv1rk2 (du,dv,dw,u,v,w,cu,cv,cw,beta,wk)
3827  CALL opdssum (du,dv,dw)
3828 c
3829  if (ldim.eq.3) then
3830  do i=1,ntot
3831  du(i) = du(i)*binvm1(i,1,1,1)
3832  dv(i) = dv(i)*binvm1(i,1,1,1)
3833  dw(i) = dw(i)*binvm1(i,1,1,1)
3834  enddo
3835  else
3836  do i=1,ntot
3837  du(i) = du(i)*binvm1(i,1,1,1)
3838  dv(i) = dv(i)*binvm1(i,1,1,1)
3839  enddo
3840  endif
3841 c
3842 c Mask
3843 c
3844  nu = mu(0)
3845  if (ldim.eq.3) then
3846  do i=1,nu
3847  du(mu(i)) = 0.
3848  dv(mu(i)) = 0.
3849  dw(mu(i)) = 0.
3850  enddo
3851  else
3852  do i=1,nu
3853  du(mu(i)) = 0.
3854  dv(mu(i)) = 0.
3855  enddo
3856  endif
3857 c
3858 #ifdef TIMER
3859  opct = ldim*ntot
3860  dct(myrout) = dct(myrout) + opct
3861  dcount = dcount + opct
3862 #endif
3863 C
3864  return
3865  end
3866 c-----------------------------------------------------------------------
3867  subroutine hypmsk3v(msk,mask)
3868 C---------------------------------------------------------------------
3869 C
3870 C Generate mask-array for the hyperbolic system (velocity).
3871 C
3872 C---------------------------------------------------------------------
3873  include 'SIZE'
3874  include 'INPUT'
3875  include 'GEOM'
3876  include 'SOLN'
3877  include 'TSTEP'
3878  integer msk(0:1)
3879  CHARACTER CB*3
3880  parameter(lxyz1=lx1*ly1*lz1)
3881  COMMON /ctmp1/ work(lxyz1,lelt)
3882  real mask(lxyz1,lelt)
3883 C
3884  nfaces= 2*ldim
3885  ntot1 = lx1*ly1*lz1*nelv
3886  CALL rzero (work ,ntot1)
3887  CALL rone (mask,ntot1)
3888 C
3889  IF (ifield.EQ.1) THEN
3890  DO 100 ie=1,nelv
3891  DO 100 iface=1,nfaces
3892  cb=cbc(iface,ie,ifield)
3893  IF (cb(1:1).EQ.'V' .OR. cb(1:1).EQ.'v') THEN
3894  CALL faccl3 (work(1,ie),vx(1,1,1,ie),unx(1,1,iface,ie),iface)
3895  CALL faddcl3(work(1,ie),vy(1,1,1,ie),uny(1,1,iface,ie),iface)
3896  IF (if3d)
3897  $ CALL faddcl3(work(1,ie),vz(1,1,1,ie),unz(1,1,iface,ie),iface)
3898  CALL fcaver (vaver,work,ie,iface)
3899 C
3900  IF (vaver.LT.0.) CALL facev (mask,ie,iface,0.0,lx1,ly1,lz1)
3901  ENDIF
3902  IF (cb(1:2).EQ.'WS' .OR. cb(1:2).EQ.'ws')
3903  $ CALL facev (mask,ie,iface,0.0,lx1,ly1,lz1)
3904  100 CONTINUE
3905  ENDIF
3906 C
3907  nm = 0
3908  ntot = lx1*ly1*lz1*nelv
3909  do i=1,ntot
3910  if (mask(i,1).eq.0) then
3911  nm = nm+1
3912  msk(nm) = i
3913  endif
3914  enddo
3915  msk(0) = nm
3916 C
3917  return
3918  END
3919 c-----------------------------------------------------------------------
3920  subroutine ophyprk (vel1,vel2,vel3,ilag)
3921 C---------------------------------------------------------------------------
3922 C
3923 C Convection of all velocity components.
3924 C Runge-Kutta scheme.
3925 C
3926 C--------------------------------------------------------------------------
3927  include 'SIZE'
3928  include 'MASS'
3929  include 'SOLN'
3930  include 'TSTEP'
3931 C
3932  REAL VEL1 (LX1,LY1,LZ1,1)
3933  REAL VEL2 (LX1,LY1,LZ1,1)
3934  REAL VEL3 (LX1,LY1,LZ1,1)
3935  COMMON /scrns/ vxn(lx1,ly1,lz1,lelv)
3936  $ , vyn(lx1,ly1,lz1,lelv)
3937  $ , vzn(lx1,ly1,lz1,lelv)
3938  $ , hv1msk(lx1,ly1,lz1,lelv)
3939  $ , hv2msk(lx1,ly1,lz1,lelv)
3940  $ , hv3msk(lx1,ly1,lz1,lelv)
3941  $ , work(lx1,ly1,lz1,lelv)
3942  COMMON /ctmp1/ rkx1(lx1,ly1,lz1,lelv)
3943  $ , rkx2(lx1,ly1,lz1,lelv)
3944  $ , rkx3(lx1,ly1,lz1,lelv)
3945  $ , rkx4(lx1,ly1,lz1,lelv)
3946  COMMON /scrmg/ rky1(lx1,ly1,lz1,lelv)
3947  $ , rky2(lx1,ly1,lz1,lelv)
3948  $ , rky3(lx1,ly1,lz1,lelv)
3949  $ , rky4(lx1,ly1,lz1,lelv)
3950  COMMON /screv/ rkz1(lx1,ly1,lz1,lelv)
3951  $ , rkz2(lx1,ly1,lz1,lelv)
3952  COMMON /scrch/ rkz3(lx1,ly1,lz1,lelv)
3953  $ , rkz4(lx1,ly1,lz1,lelv)
3954 C
3955  ntot1 = lx1*ly1*lz1*nelv
3956 C
3957  CALL opcopy (vxn,vyn,vzn,vx,vy,vz)
3958  CALL hypmsk3 (hv1msk,hv2msk,hv3msk)
3959  CALL tauinit (tau,ilag)
3960  CALL velinit (vel1,vel2,vel3,ilag)
3961  CALL velconv (vxn,vyn,vzn,tau)
3962 C
3963  DO 1000 jlag=ilag,1,-1
3964 C
3965  dtau = dtlag(jlag)/(ntaubd)
3966  dthalf = dtau/2.
3967  crk1 = dtau/6.
3968  crk2 = dtau/3.
3969 C
3970  DO 900 itau=1,ntaubd
3971 C
3972 C Stage 1.
3973 C
3974  CALL frkconv (rkx1,vel1,hv1msk)
3975  CALL frkconv (rky1,vel2,hv2msk)
3976  IF (ldim.EQ.3)
3977  $ CALL frkconv (rkz1,vel3,hv3msk)
3978 C
3979 C
3980 C Stage 2.
3981 C
3982  tau = tau + dthalf
3983  CALL velconv (vxn,vyn,vzn,tau)
3984 C
3985  CALL copy (work,vel1,ntot1)
3986  CALL add2s2 (work,rkx1,-dthalf,ntot1)
3987  CALL frkconv (rkx2,work,hv1msk)
3988 C
3989  CALL copy (work,vel2,ntot1)
3990  CALL add2s2 (work,rky1,-dthalf,ntot1)
3991  CALL frkconv (rky2,work,hv2msk)
3992 C
3993  IF (ldim.EQ.3) THEN
3994  CALL copy (work,vel3,ntot1)
3995  CALL add2s2 (work,rkz1,-dthalf,ntot1)
3996  CALL frkconv (rkz2,work,hv3msk)
3997  ENDIF
3998 C
3999 C
4000 C Stage 3.
4001 C
4002  CALL copy (work,vel1,ntot1)
4003  CALL add2s2 (work,rkx2,-dthalf,ntot1)
4004  CALL frkconv (rkx3,work,hv1msk)
4005 
4006  CALL copy (work,vel2,ntot1)
4007  CALL add2s2 (work,rky2,-dthalf,ntot1)
4008  CALL frkconv (rky3,work,hv2msk)
4009 C
4010  IF (ldim.EQ.3) THEN
4011  CALL copy (work,vel3,ntot1)
4012  CALL add2s2 (work,rkz2,-dthalf,ntot1)
4013  CALL frkconv (rkz3,work,hv3msk)
4014  ENDIF
4015 C
4016 C
4017 C Stage 4.
4018 C
4019  tau = tau + dthalf
4020  CALL velconv (vxn,vyn,vzn,tau)
4021 C
4022  CALL copy (work,vel1,ntot1)
4023  CALL add2s2 (work,rkx3,-dtau,ntot1)
4024  CALL frkconv (rkx4,work,hv1msk)
4025 
4026  CALL copy (work,vel2,ntot1)
4027  CALL add2s2 (work,rky3,-dtau,ntot1)
4028  CALL frkconv (rky4,work,hv2msk)
4029 C
4030  IF (ldim.EQ.3) THEN
4031  CALL copy (work,vel3,ntot1)
4032  CALL add2s2 (work,rkz3,-dtau,ntot1)
4033  CALL frkconv (rkz4,work,hv3msk)
4034  ENDIF
4035 C
4036 C
4037 C Sum up contributions from 4 stages.
4038 C
4039  CALL add2s2 (vel1,rkx1,-crk1,ntot1)
4040  CALL add2s2 (vel1,rkx2,-crk2,ntot1)
4041  CALL add2s2 (vel1,rkx3,-crk2,ntot1)
4042  CALL add2s2 (vel1,rkx4,-crk1,ntot1)
4043 C
4044  CALL add2s2 (vel2,rky1,-crk1,ntot1)
4045  CALL add2s2 (vel2,rky2,-crk2,ntot1)
4046  CALL add2s2 (vel2,rky3,-crk2,ntot1)
4047  CALL add2s2 (vel2,rky4,-crk1,ntot1)
4048 C
4049  IF (ldim.EQ.3) THEN
4050  CALL add2s2 (vel3,rkz1,-crk1,ntot1)
4051  CALL add2s2 (vel3,rkz2,-crk2,ntot1)
4052  CALL add2s2 (vel3,rkz3,-crk2,ntot1)
4053  CALL add2s2 (vel3,rkz4,-crk1,ntot1)
4054  ENDIF
4055 C
4056  900 CONTINUE
4057  1000 CONTINUE
4058 C
4059  CALL opcopy (vx,vy,vz,vxn,vyn,vzn)
4060 C
4061  return
4062  END
4063 c-----------------------------------------------------------------------
4064  subroutine opdiv(outfld,inpx,inpy,inpz)
4065 C---------------------------------------------------------------------
4066 C
4067 C Compute OUTFLD = SUMi Di*INPi,
4068 C the divergence of the vector field (INPX,INPY,INPZ)
4069 C
4070 C---------------------------------------------------------------------
4071  include 'SIZE'
4072  include 'GEOM'
4073  real outfld (lx2,ly2,lz2,1)
4074  real inpx (lx1,ly1,lz1,1)
4075  real inpy (lx1,ly1,lz1,1)
4076  real inpz (lx1,ly1,lz1,1)
4077  common /ctmp0/ work(lx2,ly2,lz2,lelv)
4078 C
4079  iflg = 1
4080 
4081  ntot2 = lx2*ly2*lz2*nelv
4082  call multd (work,inpx,rxm2,sxm2,txm2,1,iflg)
4083  call copy (outfld,work,ntot2)
4084  call multd (work,inpy,rym2,sym2,tym2,2,iflg)
4085  call add2 (outfld,work,ntot2)
4086  if (ldim.eq.3) then
4087  call multd (work,inpz,rzm2,szm2,tzm2,3,iflg)
4088  call add2 (outfld,work,ntot2)
4089  endif
4090 C
4091  return
4092  end
4093 C
4094 c-----------------------------------------------------------------------
4095  subroutine opgradt(outx,outy,outz,inpfld)
4096 C------------------------------------------------------------------------
4097 C
4098 C Compute DTx, DTy, DTz of an input field INPFLD
4099 C
4100 C-----------------------------------------------------------------------
4101  include 'SIZE'
4102  include 'TOTAL'
4103  real outx (lx1,ly1,lz1,1)
4104  real outy (lx1,ly1,lz1,1)
4105  real outz (lx1,ly1,lz1,1)
4106  real inpfld (lx2,ly2,lz2,1)
4107 C
4108  call cdtp (outx,inpfld,rxm2,sxm2,txm2,1)
4109  call cdtp (outy,inpfld,rym2,sym2,tym2,2)
4110  if (ldim.eq.3)
4111  $ call cdtp (outz,inpfld,rzm2,szm2,tzm2,3)
4112 C
4113  return
4114  end
4115 c-----------------------------------------------------------------------
4116  subroutine setproj(n1,nd)
4117 c
4118  include 'SIZE'
4119  include 'DEALIAS'
4120  include 'INPUT'
4121 c
4122  parameter(lx=80)
4123  real LkN(lx,lx),LkD(lx,lx),LkNt(lx,lx)
4124 c
4125  if (n1.gt.lx.or.nd.gt.lx) then
4126  write(6,*)'ERROR: increase lx in setmap to max:',n1,nd
4127  call exitt
4128  endif
4129 c
4130 c
4131  if (param(99).eq.2) then
4132  call set_pnd (pmd1 ,lkd,lknt,n1,nd)
4133  elseif (param(99).eq.3) then
4134  call set_pndoi(pmd1 ,lkd,lknt,n1,nd)
4135  endif
4136  call transpose(pmd1t,nd,pmd1,n1)
4137 c
4138  return
4139  end
4140 c-----------------------------------------------------------------------
4141  subroutine set_pndoi(Pt,P,LkNt,N,D)
4142 
4143  include 'SIZE' ! for write stmt
4144 
4145 c
4146 c Set up operators for overintegration and interpolation
4147 c
4148  integer N,D
4149  real Pt(N,D),P(D,N),LkNt(N,0:N-1)
4150 c
4151  parameter(lx=80)
4152  real zN(lx),zD(lx),wN(lx),wD(lx)
4153 c
4154 c Compute Lagrangian interpolant points
4155 c
4156  call zwgll(zn,wn,n)
4157  call zwgll(zd,wd,d)
4158 c
4159  if (nio.eq.0) write(6,*) 'dealias, pndoi:',n,d
4160  call igllm (p,pt,zn,zd,n,d,n,d)
4161 c
4162  do j=1,d
4163  do i=1,n
4164  pt(i,j) = wd(j)*pt(i,j)/wn(i)
4165  enddo
4166  enddo
4167  return
4168  end
4169 c-----------------------------------------------------------------------
4170  subroutine wgradm1(ux,uy,uz,u,nel) ! weak form of grad
4171 c
4172 c Compute gradient of T -- mesh 1 to mesh 1 (vel. to vel.)
4173 c
4174  include 'SIZE'
4175  include 'DXYZ'
4176  include 'GEOM'
4177  include 'INPUT'
4178  include 'TSTEP'
4179  include 'WZ'
4180 c
4181  parameter(lxyz=lx1*ly1*lz1)
4182  real ux(lxyz,1),uy(lxyz,1),uz(lxyz,1),u(lxyz,1)
4183 c
4184  common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
4185 
4186  integer e
4187 
4188  n = lx1-1
4189  do e=1,nel
4190  if (if3d) then
4191  call local_grad3(ur,us,ut,u,n,e,dxm1,dxtm1)
4192  do i=1,lxyz
4193  ux(i,e) = w3m1(i,1,1)*(ur(i)*rxm1(i,1,1,e)
4194  $ + us(i)*sxm1(i,1,1,e)
4195  $ + ut(i)*txm1(i,1,1,e) )
4196  uy(i,e) = w3m1(i,1,1)*(ur(i)*rym1(i,1,1,e)
4197  $ + us(i)*sym1(i,1,1,e)
4198  $ + ut(i)*tym1(i,1,1,e) )
4199  uz(i,e) = w3m1(i,1,1)*(ur(i)*rzm1(i,1,1,e)
4200  $ + us(i)*szm1(i,1,1,e)
4201  $ + ut(i)*tzm1(i,1,1,e) )
4202  enddo
4203  else
4204 
4205  if (ifaxis) then
4206  call setaxdy (ifrzer(e)) ! reset dytm1
4207  call setaxw1 (ifrzer(e)) ! reset w3m1
4208  endif
4209 
4210  call local_grad2(ur,us,u,n,e,dxm1,dytm1)
4211 
4212  do i=1,lxyz
4213  ux(i,e) =w3m1(i,1,1)*(ur(i)*rxm1(i,1,1,e)
4214  $ + us(i)*sxm1(i,1,1,e) )
4215  uy(i,e) =w3m1(i,1,1)*(ur(i)*rym1(i,1,1,e)
4216  $ + us(i)*sym1(i,1,1,e) )
4217  enddo
4218  endif
4219 
4220  enddo
4221 c
4222  return
4223  end
4224 c-----------------------------------------------------------------------
4225  SUBROUTINE makevis
4226 C----------------------------------------------------------------------
4227 C
4228 C construct viscous term:
4229 c
4230 c DEL*[mue*(DEL V + (DEL V)^T - 2/3*DEL*V *I)]
4231 c =
4232 c 2*DEL*[mue*(S - 1/3*QTL*I)]
4233 c
4234 c where mue is the viscosity, S the strain rate tensor,
4235 c tr(S) the trace of S and I the indentitiy matrix.
4236 c
4237 c NOTE: for now only for incompressible flows
4238 c In the compressible case we need to compute S using
4239 c a kth-order extrapolation scheme because we cannot
4240 c use the existing MAKEABF extrapolater and we need
4241 c to use mue/QTL from the thermo-chemical subsystem.
4242 c CAUTION: we cannot use BFX,BFY,BFZ anymore. Some
4243 c extra handling in plan4 is needed to add the
4244 c viscous term to the RHS!
4245 
4246 C----------------------------------------------------------------------
4247  include 'SIZE'
4248  include 'SOLN'
4249  include 'MASS'
4250  include 'TSTEP'
4251 C
4252  COMMON /scruz/ w1(lx1,ly1,lz1,lelv)
4253  $ , w2(lx1,ly1,lz1,lelv)
4254  $ , w3(lx1,ly1,lz1,lelv)
4255 C
4256  COMMON /scrns/ sxz(lx1,ly1,lz1,lelt)
4257  $ , syz(lx1,ly1,lz1,lelt)
4258  $ , sxx(lx1,ly1,lz1,lelt)
4259  $ , sxy(lx1,ly1,lz1,lelt)
4260  $ , syy(lx1,ly1,lz1,lelt)
4261  $ , szz(lx1,ly1,lz1,lelt)
4262 
4263  REAL fac
4264 C----------------------------------------------------------------------
4265 
4266  ntot = lx1*ly1*lz1*nelv
4267 
4268  ! CONSTRUCT strain rate tensor S (SXX, ..., SZZ)
4269  ! CALL MAKEABS
4270  CALL comp_siej(vx,vy,vz)
4271 
4272  ! substract trace of S
4273 c CALL ADD4 (W1,SXX,SYY,SZZ,NTOT)
4274  CALL copy (w1,qtl,ntot)
4275  fac = -1./3.
4276  CALL cmult (w1,fac,ntot)
4277  CALL add2 (sxx,w1,ntot)
4278  CALL add2 (syy,w1,ntot)
4279  CALL add2 (szz,w1,ntot)
4280 
4281  CALL opcolv(sxx,syy,szz,vdiff_e)
4282  CALL opcolv(sxy,sxz,syz,vdiff_e)
4283 
4284  ! not sure if that is really needed
4285  CALL opcolv (sxx,syy,szz,bm1)
4286  CALL opcolv (sxy,sxz,syz,bm1)
4287  CALL opdssum(sxx,syy,szz)
4288  CALL opdssum(sxy,sxz,syz)
4289  CALL opcolv (sxx,syy,szz,binvm1)
4290  CALL opcolv (sxy,sxz,syz,binvm1)
4291 
4292  CALL rone(w2,ntot)
4293  fac = 2.0
4294  CALL cmult (w2,fac,ntot)
4295 
4296 c add to RHS (BFX,BFY,BFZ)
4297  CALL opdiv (w1,sxx,sxy,sxz)
4298  CALL col2 (w1,w2,ntot)
4299  CALL add2 (bfx,w1,ntot)
4300 
4301  CALL opdiv (w1,sxy,syy,syz)
4302  CALL col2 (w1,w2,ntot)
4303  CALL add2 (bfy,w1,ntot)
4304 
4305  IF (ldim.EQ.3) THEN
4306  CALL opdiv (w1,sxz,syz,szz)
4307  CALL col2 (w1,w2,ntot)
4308  CALL add2 (bfz,w1,ntot)
4309  ENDIF
4310 
4311 
4312  RETURN
4313  END
4314 c-----------------------------------------------------------------------
4315 
4316  SUBROUTINE comp_siej (U1,U2,U3)
4317 C
4318 C Compute strainrates
4319 C
4320 C CAUTION : CB SCRNS is used for data change
4321 C
4322  include 'SIZE'
4323  include 'INPUT'
4324  include 'GEOM'
4325 
4326  COMMON /scrns/ exz(lx1,ly1,lz1,lelt)
4327  $ , eyz(lx1,ly1,lz1,lelt)
4328  $ , exx(lx1,ly1,lz1,lelt)
4329  $ , exy(lx1,ly1,lz1,lelt)
4330  $ , eyy(lx1,ly1,lz1,lelt)
4331  $ , ezz(lx1,ly1,lz1,lelt)
4332 
4333 C
4334  dimension u1(lx1,ly1,lz1,1)
4335  $ , u2(lx1,ly1,lz1,1)
4336  $ , u3(lx1,ly1,lz1,1)
4337 
4338 
4339  nel = nelv
4340  ntot1 = lx1*ly1*lz1*nel
4341 
4342  CALL rzero3 (exx,eyy,ezz,ntot1)
4343  CALL rzero3 (exy,exz,eyz,ntot1)
4344 C
4345  CALL uxyz (u1,exx,exy,exz,nel)
4346  CALL uxyz (u2,exy,eyy,eyz,nel)
4347  IF (ldim.EQ.3) CALL uxyz (u3,exz,eyz,ezz,nel)
4348 C
4349  CALL col2 (exx,jacmi,ntot1)
4350  CALL col2 (exy,jacmi,ntot1)
4351  CALL col2 (eyy,jacmi,ntot1)
4352 C
4353  IF (ifaxis) CALL axiezz (u2,eyy,ezz,nel)
4354 C
4355  IF (ldim.EQ.3) THEN
4356  CALL col2 (exz,jacmi,ntot1)
4357  CALL col2 (eyz,jacmi,ntot1)
4358  CALL col2 (ezz,jacmi,ntot1)
4359  ENDIF
4360 C
4361  fac = 0.5
4362  CALL cmult (exy,fac,ntot1)
4363  CALL cmult (exz,fac,ntot1)
4364  CALL cmult (eyz,fac,ntot1)
4365 
4366  RETURN
4367  END
4368 c-----------------------------------------------------------------------
4369  subroutine wlaplacian(out,a,diff,ifld)
4370 c
4371 c compute weak form of the laplacian operator including the boundary
4372 c contribution
4373 c
4374  include 'SIZE'
4375  include 'TOTAL'
4376 
4377  real out(1),a(1),diff(1)
4378  real wrk(lx1,ly1,lz1,lelt)
4379  real h2(lx1,ly1,lz1,lelt)
4380 
4381  ntot = lx1*ly1*lz1*nelfld(ifld)
4382  if (.not.iftmsh(ifld)) imesh = 1
4383  if ( iftmsh(ifld)) imesh = 2
4384 
4385  call rzero(h2,ntot)
4386 
4387  ifield_ = ifield
4388  ifield = ifld
4389 
4390  call bcneusc(out,1)
4391  call axhelm(wrk,a,diff,h2,imesh,1)
4392  call sub2 (out,wrk,ntot)
4393 
4394  ifield = ifield_
4395 
4396  return
4397  end
4398 c-----------------------------------------------------------------------
4399  subroutine explstrs ! Explicit stress tensor w/ variable viscosity
4400 
4401  include 'SIZE'
4402  include 'TOTAL'
4403 
4404  common /scruz/ u(lx1*ly1*lz1),v(lx1*ly1*lz1),w(lx1*ly1*lz1)
4405 
4406  integer e
4407 
4408  nxyz = lx1*ly1*lz1
4409 
4410 
4411  do e=1,nelv
4412 
4413  call expl_strs_e(u,v,w,vx(1,1,1,e),vy(1,1,1,e),vz(1,1,1,e),e)
4414 
4415  do i=1,nxyz
4416  bfx(i,1,1,e) = bfx(i,1,1,e) - u(i)
4417  bfy(i,1,1,e) = bfy(i,1,1,e) - v(i)
4418  bfz(i,1,1,e) = bfz(i,1,1,e) - w(i)
4419  enddo
4420 
4421  enddo
4422 
4423  return
4424  end
4425 c-----------------------------------------------------------------------
4426  subroutine expl_strs(w1,w2,w3,u1,u2,u3)
4427  include 'SIZE'
4428  include 'TOTAL'
4429 
4430  real w1(1),w2(1),w3(1),u1(1),u2(1),u3(1)
4431 
4432  integer e
4433 
4434  nxyz = lx1*ly1*lz1
4435  k = 1
4436  do e=1,nelv
4437 
4438  call expl_strs_e(w1(k),w2(k),w3(k),u1(k),u2(k),u3(k),e)
4439  k = k+nxyz
4440 
4441  enddo
4442 
4443  return
4444  end
4445 c-----------------------------------------------------------------------
4446  subroutine expl_strs_e(w1,w2,w3,u1,u2,u3,e)
4447  include 'SIZE'
4448  include 'INPUT' ! if3d
4449  include 'SOLN' ! nu_star
4450 
4451  real w1(1),w2(1),w3(1),u1(1),u2(1),u3(1)
4452  integer e
4453 
4454  integer icalld
4455  save icalld
4456  data icalld /0/
4457 
4458  if (nio.eq.0.and.icalld.eq.0) write(6,*) 'nu_star:',nu_star
4459  icalld=1
4460 
4461  if (if3d) then
4462  call expl_strs_e_3d (w1,w2,w3,u1,u2,u3,e)
4463  else
4464  call expl_strs_e_2d (w1,w2,u1,u2,e)
4465  endif
4466 
4467  return
4468  end
4469 c-----------------------------------------------------------------------
4470  subroutine expl_strs_e_3d(w1,w2,w3,u1,u2,u3,e)
4471 
4472 c Evaluate, for element e,
4473 c
4474 c /dvi\T / dui duj \
4475 c |---| | dnu --- + nu --- | (no boundary terms at present)
4476 c \dxj/ \ dxj dxi /
4477 
4478 
4479  real w1(1),w2(1),w3(1),u1(1),u2(1),u3(1)
4480  integer e
4481 
4482  include 'SIZE'
4483  include 'GEOM' ! jacmi,rxm1, etc.
4484  include 'INPUT' ! if3d
4485  include 'MASS' ! bm1
4486  include 'SOLN' ! vtrans,vdiff,nu_star
4487  include 'TSTEP' ! dt
4488  include 'WZ' ! w3m1
4489 
4490  real nu
4491 
4492  parameter(lxyz=lx1*ly1*lz1)
4493  common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
4494  $ , vr(lxyz),vs(lxyz),vt(lxyz)
4495  $ , wr(lxyz),ws(lxyz),wt(lxyz)
4496 
4497  call gradl_rst(ur,us,ut,u1,lx1,if3d) ! Grad on GLL
4498  call gradl_rst(vr,vs,vt,u2,lx1,if3d)
4499  call gradl_rst(wr,ws,wt,u3,lx1,if3d)
4500 
4501  do i=1,lxyz
4502 
4503  nu = vdiff(i,1,1,e,1)
4504  dnu = nu - nu_star ! nu_star is the constant implicit part
4505 
4506 c uij := jac*( du_i / dx_j )
4507 
4508  u11=ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e)
4509  u21=vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e)
4510  u31=wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e)
4511  u12=ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e)
4512  u22=vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e)+vt(i)*tym1(i,1,1,e)
4513  u32=wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e)
4514  u13=ur(i)*rzm1(i,1,1,e)+us(i)*szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e)
4515  u23=vr(i)*rzm1(i,1,1,e)+vs(i)*szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e)
4516  u33=wr(i)*rzm1(i,1,1,e)+ws(i)*szm1(i,1,1,e)+wt(i)*tzm1(i,1,1,e)
4517 
4518  w11 = dnu*u11 + nu*u11
4519  w12 = dnu*u12 + nu*u21
4520  w13 = dnu*u13 + nu*u31
4521  w21 = dnu*u21 + nu*u12
4522  w22 = dnu*u22 + nu*u22
4523  w23 = dnu*u23 + nu*u32
4524  w31 = dnu*u31 + nu*u13
4525  w32 = dnu*u32 + nu*u23
4526  w33 = dnu*u33 + nu*u33
4527 
4528  w = w3m1(i,1,1)*jacmi(i,e) ! note, ry has jac in it.
4529 
4530  ur(i)=(w11*rxm1(i,1,1,e)+w12*rym1(i,1,1,e)+w13*rzm1(i,1,1,e))*w
4531  us(i)=(w11*sxm1(i,1,1,e)+w12*sym1(i,1,1,e)+w13*szm1(i,1,1,e))*w
4532  ut(i)=(w11*txm1(i,1,1,e)+w12*tym1(i,1,1,e)+w13*tzm1(i,1,1,e))*w
4533  vr(i)=(w21*rxm1(i,1,1,e)+w22*rym1(i,1,1,e)+w23*rzm1(i,1,1,e))*w
4534  vs(i)=(w21*sxm1(i,1,1,e)+w22*sym1(i,1,1,e)+w23*szm1(i,1,1,e))*w
4535  vt(i)=(w21*txm1(i,1,1,e)+w22*tym1(i,1,1,e)+w23*tzm1(i,1,1,e))*w
4536  wr(i)=(w31*rxm1(i,1,1,e)+w32*rym1(i,1,1,e)+w33*rzm1(i,1,1,e))*w
4537  ws(i)=(w31*sxm1(i,1,1,e)+w32*sym1(i,1,1,e)+w33*szm1(i,1,1,e))*w
4538  wt(i)=(w31*txm1(i,1,1,e)+w32*tym1(i,1,1,e)+w33*tzm1(i,1,1,e))*w
4539 
4540  enddo
4541 
4542  call gradl_rst_t(w1,ur,us,ut,lx1,if3d)
4543  call gradl_rst_t(w2,vr,vs,vt,lx1,if3d)
4544  call gradl_rst_t(w3,wr,ws,wt,lx1,if3d)
4545 
4546  return
4547  end
4548 c-----------------------------------------------------------------------
4549  subroutine expl_strs_e_2d(w1,w2,u1,u2,e)
4550 
4551 c
4552 c Evaluate, for element e,
4553 c
4554 c /dvi\T / dui duj \
4555 c |---| | dnu --- + nu --- | (no boundary terms at present)
4556 c \dxj/ \ dxj dxi /
4557 c
4558 
4559 
4560  real w1(1),w2(1),u1(1),u2(1)
4561  integer e
4562 
4563  include 'SIZE'
4564  include 'GEOM' ! jacmi,rxm1, etc.
4565  include 'INPUT' ! if3d
4566  include 'MASS' ! bm1
4567  include 'SOLN' ! vtrans,vdiff,nu_star
4568  include 'TSTEP' ! dt
4569  include 'WZ' ! w3m1
4570 
4571  real nu
4572 
4573  parameter(lxyz=lx1*ly1*lz1)
4574  common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
4575  $ , vr(lxyz),vs(lxyz),vt(lxyz)
4576  $ , wr(lxyz),ws(lxyz),wt(lxyz)
4577 
4578  call gradl_rst(ur,us,ut,u1,lx1,if3d) ! Grad on GLL
4579  call gradl_rst(vr,vs,vt,u2,lx1,if3d)
4580 
4581  do i=1,lxyz
4582 
4583  nu = vdiff(i,1,1,e,1)
4584  dnu = nu - nu_star ! nu_star is the constant implicit part
4585 
4586 c uij := jac*( du_i / dx_j )
4587 
4588  u11=ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e)
4589  u21=vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e)
4590  u12=ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)
4591  u22=vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e)
4592 
4593  w11 = dnu*u11 + nu*u11
4594  w12 = dnu*u12 + nu*u21
4595  w21 = dnu*u21 + nu*u12
4596  w22 = dnu*u22 + nu*u22
4597 
4598  w = w3m1(i,1,1)*jacmi(i,e) ! note, ry has jac in it.
4599 
4600  ur(i)=(w11*rxm1(i,1,1,e)+w12*rym1(i,1,1,e))*w
4601  us(i)=(w11*sxm1(i,1,1,e)+w12*sym1(i,1,1,e))*w
4602  vr(i)=(w21*rxm1(i,1,1,e)+w22*rym1(i,1,1,e))*w
4603  vs(i)=(w21*sxm1(i,1,1,e)+w22*sym1(i,1,1,e))*w
4604 
4605  enddo
4606 
4607  call gradl_rst_t(w1,ur,us,ut,lx1,if3d)
4608  call gradl_rst_t(w2,vr,vs,vt,lx1,if3d)
4609 
4610  return
4611  end
4612 c-----------------------------------------------------------------------
4613  subroutine gradl_rst_t(u,ur,us,ut,md,if3d) ! GLL grad-transpose
4614 c
4615 c Thus routine originally from fsi file: u5.usr (May 2010)
4616 c
4617 c
4618  include 'SIZE'
4619  include 'DXYZ'
4620 
4621  real u(1),ur(1),us(1),ut(1)
4622  logical if3d
4623 
4624 c dgradl holds GLL-based derivative / interpolation operators
4625 
4626  parameter(ldg=lxd**3,lwkd=2*ldg)
4627  common /dgradl/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
4628  $ , wkd(lwkd)
4629  real jgl,jgt
4630 
4631  m0 = md-1
4632  call get_dgll_ptr (ip,md,md)
4633  if (if3d) then
4634  call local_grad3_t(u,ur,us,ut,m0,1,d(ip),dt(ip),wkd)
4635  else
4636  call local_grad2_t(u,ur,us ,m0,1,d(ip),dt(ip),wkd)
4637  endif
4638 
4639  return
4640  end
4641 c-----------------------------------------------------------------------
4642  subroutine gradl_rst(ur,us,ut,u,md,if3d) ! GLL-based gradient
4643 c
4644  include 'SIZE'
4645  include 'DXYZ'
4646 
4647  real ur(1),us(1),ut(1),u(1)
4648  logical if3d
4649 
4650 c dgradl holds GLL-based derivative / interpolation operators
4651 
4652  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
4653  common /dgradl/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
4654  $ , wkd(lwkd)
4655  real jgl,jgt
4656 
4657  m0 = md-1
4658  call get_dgll_ptr (ip,md,md)
4659  if (if3d) then
4660  call local_grad3(ur,us,ut,u,m0,1,d(ip),dt(ip))
4661  else
4662  call local_grad2(ur,us ,u,m0,1,d(ip),dt(ip))
4663  endif
4664 
4665  return
4666  end
4667 c-----------------------------------------------------------------------
4668  subroutine local_grad3_t(u,ur,us,ut,N,e,D,Dt,w)
4669 c Output: ur,us,ut Input:u,N,e,D,Dt
4670  real u (0:N,0:N,0:N,1)
4671  real ur(0:N,0:N,0:N),us(0:N,0:N,0:N),ut(0:N,0:N,0:N)
4672  real D (0:N,0:N),Dt(0:N,0:N)
4673  real w (0:N,0:N,0:N)
4674  integer e
4675 
4676  m1 = n+1
4677  m2 = m1*m1
4678  m3 = m1*m1*m1
4679 
4680  call mxm(dt,m1,ur,m1,u(0,0,0,e),m2)
4681 
4682  do k=0,n
4683  call mxm(us(0,0,k),m1,d ,m1,w(0,0,k),m1)
4684  enddo
4685  call add2(u(0,0,0,e),w,m3)
4686 
4687  call mxm(ut,m2,d ,m1,w,m1)
4688  call add2(u(0,0,0,e),w,m3)
4689 
4690  return
4691  end
4692 c-----------------------------------------------------------------------
4693  subroutine local_grad2_t(u,ur,us,N,e,D,Dt,w)
4694 c Output: ur,us Input:u,N,e,D,Dt
4695  real u (0:N,0:N,1)
4696  real ur(0:N,0:N),us(0:N,0:N)
4697  real D (0:N,0:N),Dt(0:N,0:N)
4698  real w (0:N,0:N)
4699  integer e
4700 
4701  m1 = n+1
4702  m2 = m1*m1
4703 
4704  call mxm(dt,m1,ur,m1,u(0,0,e),m1)
4705  call mxm(us,m1,d ,m1,w ,m1)
4706  call add2(u(0,0,e),w,m2)
4707 
4708  return
4709  end
4710 c-----------------------------------------------------------------------
4711  subroutine get_dgll_ptr (ip,mx,md)
4712 c
4713 c Get pointer to GLL-GLL interpolation dgl() for pair (mx,md)
4714 c
4715  include 'SIZE'
4716 
4717 c dgradl holds GLL-based derivative / interpolation operators
4718 
4719  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
4720  common /dgradl/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
4721  $ , wkd(lwkd)
4722  real jgl,jgt
4723 
4724 c Pointers into GLL-based derivative / interpolation operators
4725 
4726  parameter(ld=2*lxd)
4727  common /jgradl/ pd(0:ld*ld)
4728  $ , pdg(0:ld*ld)
4729  $ , pjgl(0:ld*ld)
4730  integer pd , pdg , pjgl
4731 
4732  ij = md + ld*(mx-1)
4733  ip = pdg(ij)
4734 
4735  if (ip.eq.0) then
4736 
4737  nstore = pdg(0)
4738  pdg(ij) = nstore+1
4739  nstore = nstore + md*mx
4740  pdg(0) = nstore
4741  ip = pdg(ij)
4742 
4743  if (nid.eq.985) write(6,*) nstore,ldg,ij,md,mx,ip,' NSTOR'
4744 c
4745  nwrkd = mx + md
4746  call lim_chk(nstore,ldg ,'dg ','ldg ','get_dgl_pt')
4747  call lim_chk(nwrkd ,lwkd,'wkd ','lwkd ','get_dgl_pt')
4748 c
4749  call gen_dgll(d(ip),dt(ip),md,mx,wkd)
4750  endif
4751 c
4752  return
4753  end
4754 c-----------------------------------------------------------------------
4755  subroutine gen_dgll(dgl,dgt,mp,np,w)
4756 c
4757 c Generate derivative from np GL points onto mp GL points
4758 c
4759 c dgl = derivative matrix, mapping from velocity nodes to pressure
4760 c dgt = transpose of derivative matrix
4761 c w = work array of size (3*np+mp)
4762 c
4763 c np = number of points on GLL grid
4764 c mp = number of points on GL grid
4765 c
4766 c
4767 c
4768  real dgl(mp,np),dgt(np*mp),w(1)
4769 c
4770 c
4771  iz = 1
4772  id = iz + np
4773 c
4774  call zwgll (w(iz),dgt,np) ! GL points
4775  call zwgll (w(id),dgt,mp) ! GL points
4776 c
4777  ndgt = 2*np
4778  ldgt = mp*np
4779  call lim_chk(ndgt,ldgt,'ldgt ','dgt ','gen_dgl ')
4780 c
4781  n = np-1
4782  do i=1,mp
4783  call fd_weights_full(w(id+i-1),w(iz),n,1,dgt) ! 1=1st deriv.
4784  do j=1,np
4785  dgl(i,j) = dgt(np+j) ! Derivative matrix
4786  enddo
4787  enddo
4788 c
4789  call transpose(dgt,np,dgl,mp)
4790 c
4791  return
4792  end
4793 c-----------------------------------------------------------------------
bcneusc
subroutine bcneusc(S, ITYPE)
Definition: bdry.f:786
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Definition: bdry.f:574
nekasgn
subroutine nekasgn(ix, iy, iz, e)
Definition: bdry.f:1121
rzero3
subroutine rzero3(A, B, C, N)
Definition: bdry.f:1821
facind2
subroutine facind2(JS1, JF1, JSKIP1, JS2, JF2, JSKIP2, IFC)
Definition: bdry.f:2098
bcneutr
subroutine bcneutr
Definition: bdry.f:1200
exitt
void exitt()
Definition: comm_mpi.f:604
map21e
subroutine map21e(y, x, iel)
Definition: coef.f:1609
map12
subroutine map12(y, x, iel)
Definition: coef.f:1530
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dnekclock_sync
real *8 function dnekclock_sync()
Definition: comm_mpi.f:401
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
convect_dg
subroutine convect_dg(du, u, ifuf, cr, cs, ct, ifcf)
Definition: convect.f:1667
convect_cons
subroutine convect_cons(bdu, u, ifuf, cx, cy, cz, ifcf)
Definition: convect.f:732
lim_chk
subroutine lim_chk(n, m, avar5, lvar5, sub_name10)
Definition: convect.f:454
convect_new
subroutine convect_new(bdu, u, ifuf, cx, cy, cz, ifcf)
Definition: convect.f:639
advchar
subroutine advchar
Definition: convect.f:904
vec_dssum
subroutine vec_dssum(u, v, w, nx, ny, nz)
Definition: dssum.f:164
vec_dsop
subroutine vec_dsop(u, v, w, nx, ny, nz, op)
Definition: dssum.f:212
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
fd_weights_full
subroutine fd_weights_full(xx, x, n, m, c)
Definition: fast3d.f:1293
solve
subroutine solve(F, A, K, N, ldim, IR, IC)
Definition: gauss.f:63
lu
subroutine lu(A, N, ldim, IR, IC)
Definition: gauss.f:2
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
make_hpf
subroutine make_hpf
Definition: hpf.f:2
mfi
subroutine mfi(fname_in, ifile)
Definition: ic.f:2355
ophinv
subroutine ophinv(o1, o2, o3, i1, i2, i3, h1, h2, tolh, nmxhi)
Definition: induct.f:1023
cadd
subroutine cadd(a, const, n)
Definition: math.f:326
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
glsc2
function glsc2(x, y, n)
Definition: math.f:794
vlsc21
real function vlsc21(x, y, n)
Definition: math.f:753
glsc3
function glsc3(a, b, mult, n)
Definition: math.f:776
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
copy
subroutine copy(a, b, n)
Definition: math.f:260
col4
subroutine col4(a, b, c, d, n)
Definition: math.f:120
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
subcol3
subroutine subcol3(a, b, c, n)
Definition: math.f:186
glsum
function glsum(x, n)
Definition: math.f:861
add3s2
subroutine add3s2(a, b, c, c1, c2, n)
Definition: math.f:716
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
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rone
subroutine rone(a, n)
Definition: math.f:230
invcol3
subroutine invcol3(a, b, c, n)
Definition: math.f:98
add2s1
subroutine add2s1(a, b, c1, n)
Definition: math.f:678
rzero
subroutine rzero(a, n)
Definition: math.f:208
chsign
subroutine chsign(a, n)
Definition: math.f:305
admeshv
subroutine admeshv
Definition: mvmesh.f:82
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
frkconv
subroutine frkconv(y, x, mask)
Definition: navier1.f:1863
opgradt
subroutine opgradt(outx, outy, outz, inpfld)
Definition: navier1.f:4096
cdtp
subroutine cdtp(dtx, x, rm2, sm2, tm2, isd)
Definition: navier1.f:331
makeuf
subroutine makeuf
Definition: navier1.f:1464
makebdf
subroutine makebdf
Definition: navier1.f:1549
wlaplacian
subroutine wlaplacian(out, a, diff, ifld)
Definition: navier1.f:4370
velinit
subroutine velinit(vel1, vel2, vel3, ilag)
Definition: navier1.f:1823
opbinv1
subroutine opbinv1(out1, out2, out3, inp1, inp2, inp3, SCALE)
Definition: navier1.f:839
gradl_rst
subroutine gradl_rst(ur, us, ut, u, md, if3d)
Definition: navier1.f:4643
uzawa
subroutine uzawa(rcg, h1, h2, h2inv, intype, iter)
Definition: navier1.f:2827
convprn
subroutine convprn(iconv, rbnorm, rrpt, res, z, tol)
Definition: navier1.f:1039
opcolv2c
subroutine opcolv2c(a1, a2, a3, b1, b2, c)
Definition: navier1.f:2673
conv1d
subroutine conv1d(dfi, fi)
Definition: navier1.f:3222
tauinit
subroutine tauinit(tau, ilag)
Definition: navier1.f:1808
specmp
subroutine specmp(b, nb, a, na, ba, ab, w)
Definition: navier1.f:3020
get_dgll_ptr
subroutine get_dgll_ptr(ip, mx, md)
Definition: navier1.f:4712
plan1
subroutine plan1(igeom)
Definition: navier1.f:2
frkconvv
subroutine frkconvv(du, dv, dw, u, v, w, mu)
Definition: navier1.f:3593
opamask
subroutine opamask(vbdry1, vbdry2, vbdry3)
Definition: navier1.f:2275
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
opcolv3c
subroutine opcolv3c(a1, a2, a3, b1, b2, b3, c, d)
Definition: navier1.f:2788
opcol2
subroutine opcol2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2454
advab
subroutine advab
Definition: navier1.f:1517
convuz
subroutine convuz(ifstuz)
Definition: navier1.f:2180
convop
subroutine convop(conv, fi)
Definition: navier1.f:3139
opcmult
subroutine opcmult(a, b, c, const)
Definition: navier1.f:2663
uzprec
subroutine uzprec(rpcg, rcg, h1m1, h2m1, intype, wp)
Definition: navier1.f:899
ortho
subroutine ortho(respr)
Definition: navier1.f:224
convpr
subroutine convpr(res, tol, iconv, rbnorm)
Definition: navier1.f:1067
expl_strs_e
subroutine expl_strs_e(w1, w2, w3, u1, u2, u3, e)
Definition: navier1.f:4447
ophyprk
subroutine ophyprk(vel1, vel2, vel3, ilag)
Definition: navier1.f:3921
opbinv
subroutine opbinv(out1, out2, out3, inp1, inp2, inp3, h2inv)
Definition: navier1.f:776
local_grad3_t
subroutine local_grad3_t(u, ur, us, ut, N, e, D, Dt, w)
Definition: navier1.f:4669
opadd2col
subroutine opadd2col(a1, a2, a3, b1, b2, b3, c)
Definition: navier1.f:2751
opcopy
subroutine opcopy(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2474
gradl_rst_t
subroutine gradl_rst_t(u, ur, us, ut, md, if3d)
Definition: navier1.f:4614
velchar
subroutine velchar(vel, vn, vlag, nbd, tau, dtbd)
Definition: navier1.f:1890
opdiv
subroutine opdiv(outfld, inpx, inpy, inpz)
Definition: navier1.f:4065
genwp
subroutine genwp(wp, wm2, p)
Definition: navier1.f:2160
conv1rk
subroutine conv1rk(du, dv, dw, u, v, w)
Definition: navier1.f:3432
cdabdtp
subroutine cdabdtp(ap, wp, h1, h2, h2inv, intype)
Definition: navier1.f:259
wgradm1
subroutine wgradm1(ux, uy, uz, u, nel)
Definition: navier1.f:4171
opcolv3
subroutine opcolv3(a1, a2, a3, b1, b2, b3, c)
Definition: navier1.f:2380
opcolv2
subroutine opcolv2(a1, a2, a3, b1, b2)
Definition: navier1.f:2712
expl_strs
subroutine expl_strs(w1, w2, w3, u1, u2, u3)
Definition: navier1.f:4427
opmask
subroutine opmask(res1, res2, res3)
Definition: navier1.f:2314
rotate_cyc
subroutine rotate_cyc(r1, r2, r3, idir)
Definition: navier1.f:2485
makef
subroutine makef
Definition: navier1.f:1426
bdsys
subroutine bdsys(a, b, dt, nbd, ldim)
Definition: navier1.f:1776
lagvel
subroutine lagvel
Definition: navier1.f:1930
setab3
subroutine setab3(ab0, ab1, ab2)
Definition: navier1.f:1640
conv1
subroutine conv1(du, u)
Definition: navier1.f:3277
cmask
subroutine cmask(cmask1, cmask2)
Definition: navier1.f:1391
convopo
subroutine convopo(conv, fi)
Definition: navier1.f:1222
setordbd
subroutine setordbd
Definition: navier1.f:2003
frkconvv2
subroutine frkconvv2(du, dv, dw, u, v, w, cu, cv, cw, beta, mu, wk)
Definition: navier1.f:3790
set_pndoi
subroutine set_pndoi(Pt, P, LkNt, N, D)
Definition: navier1.f:4142
opcolv
subroutine opcolv(a1, a2, a3, c)
Definition: navier1.f:2418
set_pnd
subroutine set_pnd(P, LkD, LkNt, N, D)
Definition: navier1.f:3072
normvc
subroutine normvc(h1, semi, l2, linf, x1, x2, x3)
Definition: navier1.f:2083
crespuz
subroutine crespuz(respr, g1, g2, g3, h1, h2, h2inv, intype)
Definition: navier1.f:80
setmap
subroutine setmap(n1, nd)
Definition: navier1.f:3051
setbd
subroutine setbd(bd, dtbd, nbd)
Definition: navier1.f:1736
mapw
subroutine mapw(md, nd, m1, n1, mflg)
Definition: navier1.f:2962
oprone
subroutine oprone(a, b, c)
Definition: navier1.f:2653
setabbd
subroutine setabbd(ab, dtlag, nab, nbd)
Definition: navier1.f:1663
antimsk
subroutine antimsk(y, x, xmask, n)
Definition: navier1.f:2247
expl_strs_e_3d
subroutine expl_strs_e_3d(w1, w2, w3, u1, u2, u3, e)
Definition: navier1.f:4471
expl_strs_e_2d
subroutine expl_strs_e_2d(w1, w2, u1, u2, e)
Definition: navier1.f:4550
makeg
subroutine makeg(out1, out2, out3, h1, h2, intype)
Definition: navier1.f:140
ophx
subroutine ophx(out1, out2, out3, inp1, inp2, inp3, h1, h2)
Definition: navier1.f:743
conv1no
subroutine conv1no(du, u)
Definition: navier1.f:3354
comp_siej
subroutine comp_siej(U1, U2, U3)
Definition: navier1.f:4317
setproj
subroutine setproj(n1, nd)
Definition: navier1.f:4117
velconvv
subroutine velconvv(vxn, vyn, vzn, tau)
Definition: navier1.f:3517
makeabf
subroutine makeabf
Definition: navier1.f:1598
velconv
subroutine velconv(vxn, vyn, vzn, tau)
Definition: navier1.f:1844
cresvuz
subroutine cresvuz(resv1, resv2, resv3)
Definition: navier1.f:113
chktcg2
subroutine chktcg2(tol, res, iconv)
Definition: navier1.f:1090
opdsop
subroutine opdsop(a, b, c, op)
Definition: navier1.f:2603
conv1rk2
subroutine conv1rk2(du, dv, dw, u, v, w, cu, cv, cw, beta, wk)
Definition: navier1.f:3668
nekuf
subroutine nekuf(f1, f2, f3)
Definition: navier1.f:1483
ctolspl
subroutine ctolspl(tolspl, respr)
Definition: navier1.f:194
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
hypmsk3
subroutine hypmsk3(hv1msk, hv2msk, hv3msk)
Definition: navier1.f:1956
convsp
subroutine convsp(ifstsp)
Definition: navier1.f:2241
hypmsk3v
subroutine hypmsk3v(msk, mask)
Definition: navier1.f:3868
opsub3
subroutine opsub3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
Definition: navier1.f:2370
explstrs
subroutine explstrs
Definition: navier1.f:4400
opicol2
subroutine opicol2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2633
multd
subroutine multd(dx, x, rm2, sm2, tm2, isd, iflg)
Definition: navier1.f:539
makevis
subroutine makevis
Definition: navier1.f:4226
local_grad2_t
subroutine local_grad2_t(u, ur, us, N, e, D, Dt, w)
Definition: navier1.f:4694
eprec
subroutine eprec(z2, r2)
Definition: navier1.f:962
dudxyz
subroutine dudxyz(du, u, rm1, sm1, tm1, jm1, imsh, isd)
Definition: navier1.f:1157
opchsgn
subroutine opchsgn(a, b, c)
Definition: navier1.f:2464
opsub2
subroutine opsub2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2360
gen_dgll
subroutine gen_dgll(dgl, dgt, mp, np, w)
Definition: navier1.f:4756
mapwp
subroutine mapwp(md, nd, m1, n1, mflg)
Definition: navier1.f:2994
oprzero
subroutine oprzero(a, b, c)
Definition: navier1.f:2643
normsc
subroutine normsc(h1, semi, l2, linf, x, imesh)
Definition: navier1.f:2026
opgrad
subroutine opgrad(out1, out2, out3, inp)
Definition: navier1.f:298
conv2
subroutine conv2(dtfi, fi)
Definition: navier1.f:1299
eprec2
subroutine eprec2(Z2, R2)
Definition: navier3.f:2
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
zwgll
subroutine zwgll(Z, W, NP)
Definition: speclib.f:108
pnleg
real function pnleg(Z, N)
Definition: speclib.f:876
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
axhmsf
subroutine axhmsf(au1, au2, au3, u1, u2, u3, h1, h2, matmod)
Definition: subs1.f:1241
fcaver
subroutine fcaver(xaver, a, iel, iface1)
Definition: subs1.f:870
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342
faccl3
subroutine faccl3(a, b, c, iface1)
Definition: subs1.f:944
uxyz
subroutine uxyz(u, ex, ey, ez, nel)
Definition: subs1.f:1458
axiezz
subroutine axiezz(u2, eyy, ezz, nel)
Definition: subs1.f:1534
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
amask
subroutine amask(VB1, VB2, VB3, V1, V2, V3, NEL)
Definition: subs2.f:1757
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341
setaxw1
subroutine setaxw1(IFAXWG)
Definition: subs2.f:2
rmask
subroutine rmask(R1, R2, R3, NEL)
Definition: subs2.f:1801