Context Navigation

navier1.f

main

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

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

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

Property mode set to 100755

File size: 131.7 KB

Line
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 = lx1ly1lz1*nelv
58	ntot2 = lx2ly2lz2*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 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx2ly2lz2
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 = lx2ly2lz2*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 (lx1ly1lz1,lelv)
346	real x (lx2ly2lz2,lelv)
347	real rm2 (lx2ly2lz2,lelv)
348	real sm2 (lx2ly2lz2,lelv)
349	real tm2 (lx2ly2lz2,lelv)
350	C
351	common /ctmp1/ wx (lx1ly1lz1)
352	$ , ta1 (lx1ly1lz1)
353	$ , ta2 (lx1ly1lz1)
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 = lx1ly1lz1
372	nxyz2 = lx2ly2lz2
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 (lx2ly2lz2,lelv)
561	real x (lx1ly1lz1,lelv)
562	real rm2 (lx2ly2lz2,lelv)
563	real sm2 (lx2ly2lz2,lelv)
564	real tm2 (lx2ly2lz2,lelv)
565
566	common /ctmp1/ ta1 (lx1ly1lz1)
567	$ , ta2 (lx1ly1lz1)
568	$ , ta3 (lx1ly1lz1)
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 = lx1ly1lz1
588	nxyz2 = lx2ly2lz2
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 = (H1A+H2B) * 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 = (H2B)-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=lx1ly1lz1*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=lx1ly1lz1*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 = lx2ly2lz2*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: JTAJ.
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 = lx1ly1lz1
992	ntotg = nxyz*nelgv
993	ntot1 = nxyz*nelv
994	ntot2 = lx2ly2lz2*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 = lx2ly2lz2*nelv
1051	wrk1(1) = vlsc21 (res,bm2inv,ntot2) ! resbm2invres
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 = lx2ly2lz2*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 = lx2ly2lz2*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 = lx1ly1lz1
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 = lx1ly1lz1
1255	NTOT1 = lx1ly1lz1*NELV
1256	NTOTZ = lx1ly1lz1*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 = lx1ly1lz1
1318	NTOT1 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx1ly1lz1*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.5DTB + DTC/3. )
1713	AB(2) = -0.5DTA - 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 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx1ly1lz1*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=LX1LY1LZ1)
1971	COMMON /CTMP1/ WORK (LXYZ1,LELT)
1972	C
1973	NFACES= 2*ldim
1974	NTOT1 = lx1ly1lz1*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 = lx1ly1lz1
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((SEMILENGTH*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 = lx1ly1lz1
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((SEMILENGTH*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 = lx2ly2lz2
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 = lx1ly1lz1*NELV
2207	NTOT2 = lx2ly2lz2*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 = lx1ly1lz1*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 = lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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=lx1ly1lz1*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 = cb1(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 = cb1(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=lx1ly1lz1*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=lx1ly1lz1*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(ldim2)
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=lx1ly1lz1*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 = ntot1ldim2
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 = lx2ly2lz2
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 = lx1ly1lz1*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(lxdlxdlx1)
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(lxdlxdlx1)
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(nananb)
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 = lx1ly1lz1
3180	ntot1 = lx1ly1lz1*nelv
3181	ntotz = lx1ly1lz1*nelfld(ifield)
3182	ntott = lx1ly1lz1*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 = lxdlydlzd*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 (lx1ly1lz1,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 = lx1ly1lz1
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 (lx1ly1lz1,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 = lx1ly1lz1
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(lx1ly1lz1,1),dv(lx1ly1lz1,1),dw(lx1ly1lz1,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 = lx1ly1lz1
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 = lx1ly1lz1*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 = lx1ly1lz1*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(lx1ly1lz1,1),dv(lx1ly1lz1,1),dw(lx1ly1lz1,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 = lx1ly1lz1
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(lx1ly1lz1,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 = lx1ly1lz1*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=LX1LY1LZ1)
3881	COMMON /CTMP1/ WORK(LXYZ1,LELT)
3882	real mask(lxyz1,lelt)
3883	C
3884	NFACES= 2*ldim
3885	NTOT1 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx2ly2lz2*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=lx1ly1lz1)
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/3DELV *I)]
4231	c =
4232	c 2DEL[mue(S - 1/3QTL*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 = lx1ly1lz1*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 = lx1ly1lz1*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 = lx1ly1lz1*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(lx1ly1lz1),v(lx1ly1lz1),w(lx1ly1lz1)
4405
4406	integer e
4407
4408	nxyz = lx1ly1lz1
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 = lx1ly1lz1
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=lx1ly1lz1)
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 = dnuu11 + nuu11
4519	w12 = dnuu12 + nuu21
4520	w13 = dnuu13 + nuu31
4521	w21 = dnuu21 + nuu12
4522	w22 = dnuu22 + nuu22
4523	w23 = dnuu23 + nuu32
4524	w31 = dnuu31 + nuu13
4525	w32 = dnuu32 + nuu23
4526	w33 = dnuu33 + nuu33
4527
4528	w = w3m1(i,1,1)*jacmi(i,e) ! note, ry has jac in it.
4529
4530	ur(i)=(w11rxm1(i,1,1,e)+w12rym1(i,1,1,e)+w13rzm1(i,1,1,e))w
4531	us(i)=(w11sxm1(i,1,1,e)+w12sym1(i,1,1,e)+w13szm1(i,1,1,e))w
4532	ut(i)=(w11txm1(i,1,1,e)+w12tym1(i,1,1,e)+w13tzm1(i,1,1,e))w
4533	vr(i)=(w21rxm1(i,1,1,e)+w22rym1(i,1,1,e)+w23rzm1(i,1,1,e))w
4534	vs(i)=(w21sxm1(i,1,1,e)+w22sym1(i,1,1,e)+w23szm1(i,1,1,e))w
4535	vt(i)=(w21txm1(i,1,1,e)+w22tym1(i,1,1,e)+w23tzm1(i,1,1,e))w
4536	wr(i)=(w31rxm1(i,1,1,e)+w32rym1(i,1,1,e)+w33rzm1(i,1,1,e))w
4537	ws(i)=(w31sxm1(i,1,1,e)+w32sym1(i,1,1,e)+w33szm1(i,1,1,e))w
4538	wt(i)=(w31txm1(i,1,1,e)+w32tym1(i,1,1,e)+w33tzm1(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=lx1ly1lz1)
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 = dnuu11 + nuu11
4594	w12 = dnuu12 + nuu21
4595	w21 = dnuu21 + nuu12
4596	w22 = dnuu22 + nuu22
4597
4598	w = w3m1(i,1,1)*jacmi(i,e) ! note, ry has jac in it.
4599
4600	ur(i)=(w11rxm1(i,1,1,e)+w12rym1(i,1,1,e))*w
4601	us(i)=(w11sxm1(i,1,1,e)+w12sym1(i,1,1,e))*w
4602	vr(i)=(w21rxm1(i,1,1,e)+w22rym1(i,1,1,e))*w
4603	vs(i)=(w21sxm1(i,1,1,e)+w22sym1(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=2ldg)
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=4lxd*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 = m1m1m1
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=4lxd*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-----------------------------------------------------------------------

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

Context Navigation

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

Download in other formats: