KTH framework for Nek5000 toolboxes; testing version  0.0.1
bdry.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  SUBROUTINE setlog(ifecho)
3 C
4 C Subroutine to initialize logical flags
5 C
6  include 'SIZE'
7  include 'GEOM'
8  include 'INPUT'
9  include 'TSTEP'
10  include 'CTIMER'
11  include 'ADJOINT'
12 
13  logical ifecho
14 
15  COMMON /cprint/ ifprint
16 C
17  common /nekcb/ cb
18  CHARACTER CB*3
19  LOGICAL IFALGN,IFNORX,IFNORY,IFNORZ,IFPRINT
20 C
21  nface = 2*ldim
22 C
23  ifprint = .true.
24  ifvcor = .true.
25  ifgeom = .false.
26  ifintq = .false.
27  ifsurt = .false.
28  ifwcno = .false.
29  DO 10 ifield=1,nfield
30  ifnonl(ifield) = .false.
31  10 CONTINUE
32 C
33  CALL lfalse (ifeppm,nface*nelv)
34  CALL lfalse (ifqinp,nface*nelv)
35 C
36  IF (ifmvbd) THEN
37  ifgeom = .true.
38  IF ( ifflow .AND. .NOT.ifnav ) ifwcno = .true.
39  IF ( ifmelt .AND. .NOT.ifflow ) ifwcno = .true.
40  ENDIF
41 C
42  IF (ifflow) THEN
43  ierr = 0
44  ifield = 1
45  DO 100 iel=1,nelv
46  DO 100 ifc=1,nface
47  cb = cbc(ifc,iel,ifield)
48  CALL chknord (ifalgn,ifnorx,ifnory,ifnorz,ifc,iel)
49  CALL chkcbc (cb,iel,ifc,ifalgn,ierr)
50  IF (cb.EQ.'O ' .OR. cb.EQ.'o ' .OR.
51  $ cb.EQ.'ON ' .OR. cb.EQ.'on ' .OR.
52  $ cb.EQ.'S ' .OR. cb.EQ.'s ' .OR.
53  $ cb.EQ.'SL ' .OR. cb.EQ.'sl ' .OR.
54  $ cb.EQ.'MM ' .OR. cb.EQ.'mm ' .OR.
55  $ cb.EQ.'MS ' .OR. cb.EQ.'ms ') THEN
56  ifvcor = .false.
57  ifeppm(ifc,iel) = .true.
58  ENDIF
59  IF (cb.EQ.'VL ' .OR. cb.EQ.'vl ' .OR.
60  $ cb.EQ.'WSL' .OR. cb.EQ.'wsl' .OR.
61  $ cb.EQ.'SL ' .OR. cb.EQ.'sl ' .OR.
62  $ cb.EQ.'SHL' .OR. cb.EQ.'shl' .OR.
63  $ cb.EQ.'MM ' .OR. cb.EQ.'mm ' .OR.
64  $ cb.EQ.'MS ' .OR. cb.EQ.'ms ' .OR.
65  $ cb.EQ.'O ' .OR. cb.EQ.'o ' .OR.
66  $ cb.EQ.'ON ' .OR. cb.EQ.'on ') THEN
67  ifqinp(ifc,iel) = .true.
68  ENDIF
69  IF (cb.EQ.'MS ' .OR. cb.EQ.'ms ' .OR.
70  $ cb.EQ.'MM ' .OR. cb.EQ.'mm ' .OR.
71  $ cb.EQ.'MSI' .OR. cb.EQ.'msi' ) THEN
72  ifsurt = .true.
73  ENDIF
74  100 CONTINUE
75 
76  ierr = iglsum(ierr,1)
77  if (ierr.gt.0) call exitt
78  ENDIF
79 C
80  IF (ifheat) THEN
81 C
82  DO 250 ifield=2,nfield
83  DO 250 iel=1,nelfld(ifield)
84  DO 250 ifc=1,nface
85  cb=cbc(ifc,iel,ifield)
86  IF (cb.EQ.'r ' .OR. cb.EQ.'R ') THEN
87  ifnonl(ifield) = .true.
88  ENDIF
89  250 CONTINUE
90 C
91  ENDIF
92 
93  if (ifmhd) call set_ifbcor
94 C
95 C Establish global consistency of LOGICALS amongst all processors.
96 C
97  CALL gllog(ifvcor , .false.)
98  CALL gllog(ifsurt , .true. )
99  CALL gllog(ifwcno , .true. )
100  DO 400 ifield=2,nfield
101  CALL gllog(ifnonl(ifield),.true.)
102  400 CONTINUE
103 C
104  IF (nio.EQ.0 .AND. ifecho) THEN
105  WRITE (6,*) 'IFTRAN =',iftran
106  WRITE (6,*) 'IFFLOW =',ifflow
107  WRITE (6,*) 'IFHEAT =',ifheat
108  WRITE (6,*) 'IFSPLIT =',ifsplit
109  WRITE (6,*) 'IFLOMACH =',iflomach
110  WRITE (6,*) 'IFUSERVP =',ifuservp
111  WRITE (6,*) 'IFUSERMV =',ifusermv
112  WRITE (6,*) 'IFPERT =',ifpert
113  WRITE (6,*) 'IFADJ =',ifadj
114  WRITE (6,*) 'IFSTRS =',ifstrs
115  WRITE (6,*) 'IFCHAR =',ifchar
116  WRITE (6,*) 'IFCYCLIC =',ifcyclic
117  WRITE (6,*) 'IFAXIS =',ifaxis
118  WRITE (6,*) 'IFMVBD =',ifmvbd
119  WRITE (6,*) 'IFMELT =',ifmelt
120  WRITE (6,*) 'IFNEKNEK =',ifneknek
121  WRITE (6,*) 'IFNEKNEKC =',ifneknekc
122  WRITE (6,*) 'IFSYNC =',ifsync
123  WRITE (6,*) ' '
124  WRITE (6,*) 'IFVCOR =',ifvcor
125  WRITE (6,*) 'IFINTQ =',ifintq
126  WRITE (6,*) 'IFGEOM =',ifgeom
127  WRITE (6,*) 'IFSURT =',ifsurt
128  WRITE (6,*) 'IFWCNO =',ifwcno
129 
130  DO 500 ifield=1,nfield
131  WRITE (6,*) ' '
132  WRITE (6,*) 'IFTMSH for field',ifield,' = ',iftmsh(ifield)
133  WRITE (6,*) 'IFADVC for field',ifield,' = ',ifadvc(ifield)
134  WRITE (6,*) 'IFNONL for field',ifield,' = ',ifnonl(ifield)
135  500 CONTINUE
136  WRITE (6,*) ' '
137  if (param(99).gt.0) write(6,*) 'Dealiasing enabled, nxd=', nxd
138  ENDIF
139 C
140  RETURN
141  END
142 C
143 c-----------------------------------------------------------------------
144  SUBROUTINE setrzer
145 C-------------------------------------------------------------------
146 C
147 C Check for axisymmetric case.
148 C Are some of the elements close to the axis?
149 C
150 C-------------------------------------------------------------------
151  include 'SIZE'
152  include 'GEOM'
153  include 'INPUT'
154 C
155 C Single or double precision???
156 C
157  delta = 1.e-9
158  x = 1.+delta
159  y = 1.
160  diff = abs(x-y)
161  IF (diff.EQ.0.) eps = 1.e-7
162  IF (diff.GT.0.) eps = 1.e-14
163  eps1 = 1.e-6 ! for prenek mesh in real*4
164 C
165  DO 100 iel=1,nelt
166  ifrzer(iel) = .false.
167  IF (ifaxis) THEN
168  nvert = 0
169  DO 10 ic=1,4
170  IF(abs(yc(ic,iel)).LT.eps1) THEN
171  nvert = nvert+1
172  yc(ic,iel) = 0.0 ! exactly on the axis
173  ENDIF
174  10 CONTINUE
175  ENDIF
176  iedge = 1
177  IF ((nvert.EQ.2).AND.(ccurve(iedge,iel).EQ.' '))
178  $ ifrzer(iel) = .true.
179  100 CONTINUE
180  RETURN
181  END
182 C
183 c-----------------------------------------------------------------------
184  SUBROUTINE chknord (IFALGN,IFNORX,IFNORY,IFNORZ,IFC,IEL)
185 C
186 C Check direction of normal of an element face for
187 C alignment with the X, Y, or Z axis.
188 C
189  include 'SIZE'
190  include 'GEOM'
191 C
192  LOGICAL IFALGN,IFNORX,IFNORY,IFNORZ
193 C
194  sumx = 0.0
195  sumy = 0.0
196  sumz = 0.0
197  tolnor = 1.0e-3
198  ifalgn = .false.
199  ifnorx = .false.
200  ifnory = .false.
201  ifnorz = .false.
202 C
203  IF (ldim.EQ.2) THEN
204 C
205  ncpf = lx1
206  DO 100 ix=1,lx1
207  sumx = sumx + abs( abs(unx(ix,1,ifc,iel)) - 1.0 )
208  sumy = sumy + abs( abs(uny(ix,1,ifc,iel)) - 1.0 )
209  100 CONTINUE
210  sumx = sumx / ncpf
211  sumy = sumy / ncpf
212  IF ( sumx.LT.tolnor ) THEN
213  ifnorx = .true.
214  ifalgn = .true.
215  ENDIF
216  IF ( sumy.LT.tolnor ) THEN
217  ifnory = .true.
218  ifalgn = .true.
219  ENDIF
220 C
221  ELSE
222 C
223  ncpf = lx1*lx1
224  DO 200 ix=1,lx1
225  DO 200 iy=1,ly1
226  sumx = sumx + abs( abs(unx(ix,iy,ifc,iel)) - 1.0 )
227  sumy = sumy + abs( abs(uny(ix,iy,ifc,iel)) - 1.0 )
228  sumz = sumz + abs( abs(unz(ix,iy,ifc,iel)) - 1.0 )
229  200 CONTINUE
230  sumx = sumx / ncpf
231  sumy = sumy / ncpf
232  sumz = sumz / ncpf
233  IF ( sumx.LT.tolnor ) THEN
234  ifnorx = .true.
235  ifalgn = .true.
236  ENDIF
237  IF ( sumy.LT.tolnor ) THEN
238  ifnory = .true.
239  ifalgn = .true.
240  ENDIF
241  IF ( sumz.LT.tolnor ) THEN
242  ifnorz = .true.
243  ifalgn = .true.
244  ENDIF
245 C
246  ENDIF
247 C
248  RETURN
249  END
250 c-----------------------------------------------------------------------
251  SUBROUTINE chkaxcb
252 C
253  include 'SIZE'
254  include 'INPUT'
255  CHARACTER CB*3
256 C
257  ifld = 1
258  nface = 2*ldim
259 C
260  DO 100 iel=1,nelv
261  DO 100 ifc=1,nface
262  cb = cbc(ifc,iel,ifld)
263  IF (cb.EQ.'A ' .AND. ifc.NE.1) GOTO 9000
264  100 CONTINUE
265 C
266  RETURN
267 C
268  9000 WRITE (6,*) ' Element face on the axis of symmetry must be FACE 1'
269  WRITE (6,*) ' Element',iel,' face',ifc,' is on the axis.'
270  call exitt
271 C
272  END
273 c-----------------------------------------------------------------------
274  SUBROUTINE chkcbc (CB,IEL,IFC,IFALGN,IERR)
275  include 'SIZE'
276  include 'PARALLEL'
277  include 'INPUT'
278 C
279 C Check for illegal boundary conditions
280 C
281  CHARACTER CB*3
282  LOGICAL IFALGN
283 
284  if (ifstrs) return
285 
286  ieg = lglel(iel)
287 
288 C Laplacian formulation only
289  IF (cb.EQ.'SH ' .OR. cb.EQ.'sh ' .OR.
290  $ cb.EQ.'SHL' .OR. cb.EQ.'shl' .OR.
291  $ cb.EQ.'S ' .OR. cb.EQ.'s ' .OR.
292  $ cb.EQ.'SL ' .OR. cb.EQ.'sl ' .OR.
293  $ cb.EQ.'MM ' .OR. cb.EQ.'mm ' .OR.
294  $ cb.EQ.'MS ' .OR. cb.EQ.'ms ' .OR.
295  $ cb.EQ.'MSI' .OR. cb.EQ.'msi' ) GOTO 9001
296 
297  IF ( .NOT.ifalgn .AND.
298  $ (cb.EQ.'ON ' .OR. cb.EQ.'on ' .OR. cb.EQ.'SYM') ) GOTO 9010
299 
300  RETURN
301 
302  9001 WRITE (6,*) ' Illegal traction boundary conditions detected for'
303  GOTO 9999
304 
305  9010 WRITE (6,*) ' Mixed B.C. on a side nonaligned with either the X,Y,
306  $ or Z axis detected for'
307 
308  9999 WRITE (6,*) ' Element',ieg,' side',ifc
309  WRITE (6,*) ' Requires PN/PN-2 STRESS FORMULATION'
310 
311  ierr = err + 1
312 
313  END
314 c-----------------------------------------------------------------------
315  SUBROUTINE bcmask
316 C
317 C Zero out masks corresponding to Dirichlet boundary points.
318 C
319  include 'SIZE'
320  include 'TSTEP'
321  include 'INPUT'
322  include 'MVGEOM'
323  include 'SOLN'
324  include 'TOPOL'
325 
326  common /nekcb/ cb
327  character*3 cb
328  character*1 cb1(3)
329  equivalence(cb1,cb)
330 
331  logical ifalgn,ifnorx,ifnory,ifnorz
332  integer e,f
333 
334  nfaces=2*ldim
335  nxyz =lx1*ly1*lz1
336 
337 C
338 C Masks for moving mesh
339 C
340  IF (ifmvbd) THEN
341  ifield = 0
342  CALL stsmask (w1mask,w2mask,w3mask)
343  do e=1,nelv
344  do f=1,nfaces
345  if (cbc(f,e,1).eq.'msi'.or.cbc(f,e,1).eq.'msi') then
346  call facev(w1mask,e,f,0.0,lx1,ly1,lz1)
347  call facev(w2mask,e,f,0.0,lx1,ly1,lz1)
348  call facev(w3mask,e,f,0.0,lx1,ly1,lz1)
349  endif
350  enddo
351  enddo
352  ENDIF
353 C
354 C Masks for flow variables
355 C
356  IF (ifflow) THEN
357  ifield = 1
358  nel = nelfld(ifield)
359  ntot = nxyz*nel
360 C
361 C Pressure mask
362 C
363  call rone(pmask,ntot)
364  do 50 iel=1,nelt
365  do 50 iface=1,nfaces
366  cb=cbc(iface,iel,ifield)
367  if (cb.eq.'O ' .or. cb.eq.'ON ' .or.
368  $ cb.eq.'o ' .or. cb.eq.'on ')
369  $ call facev(pmask,iel,iface,0.0,lx1,ly1,lz1)
370  50 continue
371  if (nelt.gt.nelv) then
372  nn=lx1*ly1*lz1*(nelt-nelv)
373  call rzero(pmask(1,1,1,nelv+1),nn)
374  endif
375 C
376 C Zero out mask at Neumann-Dirichlet interfaces
377 C
378  CALL dsop(pmask,'MUL',lx1,ly1,lz1)
379 C
380 C Velocity masks
381 C
382  IF (ifstrs) THEN
383  CALL stsmask (v1mask,v2mask,v3mask)
384  ELSE
385 C
386  CALL rone(v1mask,ntot)
387  CALL rone(v2mask,ntot)
388  CALL rone(v3mask,ntot)
389 C
390  DO 100 iel=1,nelv
391  DO 100 iface=1,nfaces
392  cb =cbc(iface,iel,ifield)
393  CALL chknord (ifalgn,ifnorx,ifnory,ifnorz,iface,iel)
394 C
395 C All-Dirichlet boundary conditions
396 C
397  IF (cb.EQ.'v ' .OR. cb.EQ.'V ' .OR. cb.EQ.'vl ' .OR.
398  $ cb.eq.'MV ' .or. cb.eq.'mv ' .or.
399  $ cb.EQ.'VL ' .OR. cb.EQ.'W ') THEN
400  CALL facev (v1mask,iel,iface,0.0,lx1,ly1,lz1)
401  CALL facev (v2mask,iel,iface,0.0,lx1,ly1,lz1)
402  CALL facev (v3mask,iel,iface,0.0,lx1,ly1,lz1)
403  GOTO 100
404  ENDIF
405 C
406 C Mixed-Dirichlet-Neumann boundary conditions
407 C
408  IF (cb.EQ.'SYM') THEN
409  IF ( .NOT.ifalgn .OR. ifnorx )
410  $ CALL facev (v1mask,iel,iface,0.0,lx1,ly1,lz1)
411  IF ( ifnory )
412  $ CALL facev (v2mask,iel,iface,0.0,lx1,ly1,lz1)
413  IF ( ifnorz )
414  $ CALL facev (v3mask,iel,iface,0.0,lx1,ly1,lz1)
415  GOTO 100
416  ENDIF
417 
418  IF (cb.EQ.'ON ' .OR. cb.EQ.'on ') THEN
419  IF ( ifnory .OR. ifnorz )
420  $ CALL facev (v1mask,iel,iface,0.0,lx1,ly1,lz1)
421  IF ( .NOT.ifalgn .OR. ifnorx .OR. ifnorz )
422  $ CALL facev (v2mask,iel,iface,0.0,lx1,ly1,lz1)
423  IF ( .NOT.ifalgn .OR. ifnorx .OR. ifnory )
424  $ CALL facev (v3mask,iel,iface,0.0,lx1,ly1,lz1)
425  GOTO 100
426  ENDIF
427  IF (cb.EQ.'A ') THEN
428  CALL facev (v2mask,iel,iface,0.0,lx1,ly1,lz1)
429  ENDIF
430  100 CONTINUE
431 
432  call opdsop(v1mask,v2mask,v3mask,'MUL') ! no rotation for mul
433 
434  ENDIF
435 
436  CALL rone(omask,ntot)
437  DO 200 iel=1,nelv
438  DO 200 iface=1,nfaces
439  cb =cbc(iface,iel,ifield)
440  IF (cb.EQ.'A ') THEN
441  CALL facev (omask,iel,iface,0.0,lx1,ly1,lz1)
442  ENDIF
443  200 CONTINUE
444  CALL dsop(omask,'MUL',lx1,ly1,lz1)
445 C
446  ENDIF
447 C
448 C Masks for passive scalars
449 C
450  IF (ifheat) THEN
451 C
452  DO 1200 ifield=2,nfield
453  ipscal = ifield-1
454  nel = nelfld(ifield)
455  ntot = nxyz*nel
456  CALL rone (tmask(1,1,1,1,ipscal),ntot)
457  DO 1100 iel=1,nel
458  DO 1100 iface=1,nfaces
459  cb =cbc(iface,iel,ifield)
460 C
461 C Assign mask values.
462 C
463  IF (cb.EQ.'T ' .OR. cb.EQ.'t ' .OR.
464  $ (cb.EQ.'A ' .AND. ifaziv) .OR.
465  $ cb.EQ.'MCI' .OR. cb.EQ.'MLI' .OR.
466  $ cb.EQ.'KD ' .OR. cb.EQ.'kd ' .OR.
467  $ cb.EQ.'ED ' .OR. cb.EQ.'ed ' .OR.
468  $ cb.EQ.'KW ' .OR. cb.EQ.'KWS' .OR. cb.EQ.'EWS')
469  $ CALL facev (tmask(1,1,1,1,ipscal),
470  $ iel,iface,0.0,lx1,ly1,lz1)
471  1100 CONTINUE
472  CALL dsop (tmask(1,1,1,1,ipscal),'MUL',lx1,ly1,lz1)
473  1200 CONTINUE
474 C
475  ENDIF
476 C
477 C Masks for B-field
478 C
479  if (ifmhd) then
480  ifield = ifldmhd
481  nel = nelfld(ifield)
482  ntot = nxyz*nel
483 C
484 C B-field pressure mask
485 C
486  call rone(bpmask,ntot)
487  do iel=1,nelv
488  do iface=1,nfaces
489  cb=cbc(iface,iel,ifield)
490  if (cb.eq.'O ' .or. cb.eq.'ON ')
491  $ call facev(bpmask,iel,iface,0.0,lx1,ly1,lz1)
492  enddo
493  enddo
494 C
495 C Zero out mask at Neumann-Dirichlet interfaces
496 C
497  call dsop(bpmask,'MUL',lx1,ly1,lz1)
498 C
499 C B-field masks
500 C
501  if (ifstrs) then
502  call stsmask (b1mask,b2mask,b3mask)
503  else
504 C
505  call rone(b1mask,ntot)
506  call rone(b2mask,ntot)
507  call rone(b3mask,ntot)
508 C
509  do iel=1,nelv
510  do iface=1,nfaces
511  cb =cbc(iface,iel,ifield)
512  call chknord (ifalgn,ifnorx,ifnory,ifnorz,iface,iel)
513 c
514  if (cb.eq.'v ' .or. cb.eq.'V ' .or. cb.eq.'vl ' .or.
515  $ cb.eq.'VL ' .or. cb.eq.'W ') then
516 c
517 c All-Dirichlet boundary conditions
518 c
519  call facev (b1mask,iel,iface,0.0,lx1,ly1,lz1)
520  call facev (b2mask,iel,iface,0.0,lx1,ly1,lz1)
521  call facev (b3mask,iel,iface,0.0,lx1,ly1,lz1)
522 c
523  elseif (cb.eq.'SYM') then
524 c
525 c Mixed-Dirichlet-Neumann boundary conditions
526 c
527  if ( .not.ifalgn .or. ifnorx )
528  $ call facev (b1mask,iel,iface,0.0,lx1,ly1,lz1)
529  if ( ifnory )
530  $ call facev (b2mask,iel,iface,0.0,lx1,ly1,lz1)
531  if ( ifnorz )
532  $ call facev (b3mask,iel,iface,0.0,lx1,ly1,lz1)
533 c
534  elseif (cb.eq.'ON ') then
535 c
536 c Mixed-Dirichlet-Neumann boundary conditions
537 c
538  if ( ifnory .or. ifnorz )
539  $ call facev (b1mask,iel,iface,0.0,lx1,ly1,lz1)
540  if ( .not.ifalgn .or. ifnorx .or. ifnorz )
541  $ call facev (b2mask,iel,iface,0.0,lx1,ly1,lz1)
542  if ( .not.ifalgn .or. ifnorx .or. ifnory )
543  $ call facev (b3mask,iel,iface,0.0,lx1,ly1,lz1)
544 c
545  elseif (cb.eq.'A ') then
546 c
547 c axisymmetric centerline
548 c
549  call facev (b2mask,iel,iface,0.0,lx1,ly1,lz1)
550 c
551  else
552 c
553  if ( cb1(1).eq.'d' )
554  $ call facev (b1mask,iel,iface,0.0,lx1,ly1,lz1)
555  if ( cb1(2).eq.'d' )
556  $ call facev (b2mask,iel,iface,0.0,lx1,ly1,lz1)
557  if ( cb1(3).eq.'d' .and. if3d )
558  $ call facev (b3mask,iel,iface,0.0,lx1,ly1,lz1)
559 c
560  endif
561  enddo
562  enddo
563 c
564  call dsop(b1mask,'MUL',lx1,ly1,lz1)
565  call dsop(b2mask,'MUL',lx1,ly1,lz1)
566  if (ldim.eq.3) call dsop(b3mask,'MUL',lx1,ly1,lz1)
567  endif
568  endif
569 C
570  RETURN
571  END
572 c-----------------------------------------------------------------------
573  SUBROUTINE bcdirvc(V1,V2,V3,mask1,mask2,mask3)
574 C
575 C Apply Dirichlet boundary conditions to surface of vector (V1,V2,V3).
576 C Use IFIELD as a guide to which boundary conditions are to be applied.
577 C
578  include 'SIZE'
579  include 'TSTEP'
580  include 'INPUT'
581  include 'GEOM'
582  include 'SOLN'
583  include 'TOPOL'
584  include 'CTIMER'
585  COMMON /scruz/ tmp1(lx1,ly1,lz1,lelv)
586  $ , tmp2(lx1,ly1,lz1,lelv)
587  $ , tmp3(lx1,ly1,lz1,lelv)
588  COMMON /scrmg/ tmq1(lx1,ly1,lz1,lelv)
589  $ , tmq2(lx1,ly1,lz1,lelv)
590  $ , tmq3(lx1,ly1,lz1,lelv)
591 C
592  REAL V1(lx1,ly1,lz1,LELV),V2(lx1,ly1,lz1,LELV)
593  $ ,V3(lx1,ly1,lz1,LELV)
594  real mask1(lx1,ly1,lz1,lelv),mask2(lx1,ly1,lz1,lelv)
595  $ ,mask3(lx1,ly1,lz1,lelv)
596 c
597  common /nekcb/ cb
598  character cb*3
599  character*1 cb1(3)
600  equivalence(cb1,cb)
601 c
602  logical ifonbc
603 c
604  ifonbc = .false.
605 c
606  if (icalld.eq.0) then
607  tusbc=0.0
608  nusbc=0
609  icalld=icalld+1
610  endif
611  nusbc=nusbc+1
612  etime1=dnekclock()
613 C
614 C
615  nfaces=2*ldim
616  nxyz =lx1*ly1*lz1
617  nel =nelfld(ifield)
618  ntot =nxyz*nel
619 C
620  CALL rzero(tmp1,ntot)
621  CALL rzero(tmp2,ntot)
622  IF (if3d) CALL rzero(tmp3,ntot)
623 C
624 C Velocity boundary conditions
625 C
626 c write(6,*) 'BCDIRV: ifield',ifield
627  DO 2100 isweep=1,2
628  DO 2000 ie=1,nel
629  DO 2000 iface=1,nfaces
630  cb = cbc(iface,ie,ifield)
631  bc1 = bc(1,iface,ie,ifield)
632  bc2 = bc(2,iface,ie,ifield)
633  bc3 = bc(3,iface,ie,ifield)
634 
635  IF (cb.EQ.'V ' .OR. cb.EQ.'VL ' .OR.
636  $ cb.EQ.'WS ' .OR. cb.EQ.'WSL') THEN
637  CALL facev (tmp1,ie,iface,bc1,lx1,ly1,lz1)
638  CALL facev (tmp2,ie,iface,bc2,lx1,ly1,lz1)
639  IF (if3d) CALL facev (tmp3,ie,iface,bc3,lx1,ly1,lz1)
640  IF ( ifqinp(iface,ie) )
641  $ CALL globrot (tmp1(1,1,1,ie),tmp2(1,1,1,ie),
642  $ tmp3(1,1,1,ie),ie,iface)
643  ENDIF
644 
645  IF (cb.EQ.'v ' .OR. cb.EQ.'vl ' .OR.
646  $ cb.EQ.'ws ' .OR. cb.EQ.'wsl' .OR.
647  $ cb.EQ.'mv ' .OR. cb.EQ.'mvn' .OR.
648  $ cb1(1).eq.'d'.or.cb1(2).eq.'d'.or.cb1(3).eq.'d') then
649 
650  call faceiv (cb,tmp1(1,1,1,ie),tmp2(1,1,1,ie),
651  $ tmp3(1,1,1,ie),ie,iface,lx1,ly1,lz1)
652 
653  IF ( ifqinp(iface,ie) )
654  $ CALL globrot (tmp1(1,1,1,ie),tmp2(1,1,1,ie),
655  $ tmp3(1,1,1,ie),ie,iface)
656  ENDIF
657 
658  IF (cb.EQ.'ON ' .OR. cb.EQ.'on ') then ! 5/21/01 pff
659  ifonbc =.true.
660  CALL faceiv ('v ',tmp1(1,1,1,ie),tmp2(1,1,1,ie),
661  $ tmp3(1,1,1,ie),ie,iface,lx1,ly1,lz1)
662  ENDIF
663 
664  2000 CONTINUE
665  DO 2010 ie=1,nel
666  DO 2010 iface=1,nfaces
667  IF (cbc(iface,ie,ifield).EQ.'W ') THEN
668  CALL facev (tmp1,ie,iface,0.0,lx1,ly1,lz1)
669  CALL facev (tmp2,ie,iface,0.0,lx1,ly1,lz1)
670  IF (if3d) CALL facev (tmp3,ie,iface,0.0,lx1,ly1,lz1)
671  ENDIF
672  2010 CONTINUE
673 C
674 C Take care of Neumann-Dirichlet shared edges...
675 C
676  if (isweep.eq.1) then
677  call opdsop(tmp1,tmp2,tmp3,'MXA')
678  else
679  call opdsop(tmp1,tmp2,tmp3,'MNA')
680  endif
681  2100 CONTINUE
682 C
683 C Copy temporary array to velocity arrays.
684 C
685  IF ( .NOT.ifstrs ) THEN
686  CALL col2(v1,mask1,ntot)
687  CALL col2(v2,mask2,ntot)
688  IF (if3d) CALL col2(v3,mask3,ntot)
689  if (ifonbc) then
690  call antimsk1(tmp1,mask1,ntot)
691  call antimsk1(tmp2,mask2,ntot)
692  if (if3d) call antimsk1(tmp3,mask3,ntot)
693  endif
694  ELSE
695  CALL rmask (v1,v2,v3,nelv)
696  ENDIF
697 
698  CALL add2(v1,tmp1,ntot)
699  CALL add2(v2,tmp2,ntot)
700  IF (if3d) CALL add2(v3,tmp3,ntot)
701 
702  if (ifneknekc) call fix_surface_flux
703 
704  tusbc=tusbc+(dnekclock()-etime1)
705 
706  RETURN
707  END
708 c-----------------------------------------------------------------------
709  SUBROUTINE bcdirsc(S)
710 C
711 C Apply Dirichlet boundary conditions to surface of scalar, S.
712 C Use IFIELD as a guide to which boundary conditions are to be applied.
713 C
714  include 'SIZE'
715  include 'TSTEP'
716  include 'INPUT'
717  include 'SOLN'
718  include 'TOPOL'
719  include 'CTIMER'
720 C
721  dimension s(lx1,ly1,lz1,lelt)
722  COMMON /scrsf/ tmp(lx1,ly1,lz1,lelt)
723  $ , tma(lx1,ly1,lz1,lelt)
724  $ , smu(lx1,ly1,lz1,lelt)
725  common /nekcb/ cb
726  CHARACTER CB*3
727 
728  if (icalld.eq.0) then
729  tusbc=0.0
730  nusbc=0
731  icalld=icalld+1
732  endif
733  nusbc=nusbc+1
734  etime1=dnekclock()
735 C
736  ifld = 1
737  nfaces = 2*ldim
738  nxyz = lx1*ly1*lz1
739  nel = nelfld(ifield)
740  ntot = nxyz*nel
741  nfldt = nfield - 1
742 C
743  CALL rzero(tmp,ntot)
744 C
745 C Temperature boundary condition
746 C
747  DO 2100 isweep=1,2
748 C
749  DO 2010 ie=1,nel
750  DO 2010 iface=1,nfaces
751  cb=cbc(iface,ie,ifield)
752  bc1=bc(1,iface,ie,ifield)
753  bc2=bc(2,iface,ie,ifield)
754  bc3=bc(3,iface,ie,ifield)
755  bc4=bc(4,iface,ie,ifield)
756  bck=bc(4,iface,ie,ifld)
757  bce=bc(5,iface,ie,ifld)
758  IF (cb.EQ.'T ') CALL facev (tmp,ie,iface,bc1,lx1,ly1,lz1)
759  IF (cb.EQ.'MCI') CALL facev (tmp,ie,iface,bc4,lx1,ly1,lz1)
760  IF (cb.EQ.'MLI') CALL facev (tmp,ie,iface,bc4,lx1,ly1,lz1)
761  IF (cb.EQ.'KD ') CALL facev (tmp,ie,iface,bck,lx1,ly1,lz1)
762  IF (cb.EQ.'ED ') CALL facev (tmp,ie,iface,bce,lx1,ly1,lz1)
763  IF (cb.EQ.'t ' .OR. cb.EQ.'kd ' .or.
764  $ cb.EQ.'ed ' .or. cb.eq.'o ' .or. cb.eq.'on ')
765  $ CALL faceis (cb,tmp(1,1,1,ie),ie,iface,lx1,ly1,lz1)
766  2010 CONTINUE
767 C
768 C Take care of Neumann-Dirichlet shared edges...
769 C
770  IF (isweep.EQ.1) CALL dsop(tmp,'MXA',lx1,ly1,lz1)
771  IF (isweep.EQ.2) CALL dsop(tmp,'MNA',lx1,ly1,lz1)
772  2100 CONTINUE
773 C
774 C Copy temporary array to temperature array.
775 C
776  CALL col2(s,tmask(1,1,1,1,ifield-1),ntot)
777  CALL add2(s,tmp,ntot)
778 
779  tusbc=tusbc+(dnekclock()-etime1)
780 
781  RETURN
782  END
783 C
784 c-----------------------------------------------------------------------
785  SUBROUTINE bcneusc(S,ITYPE)
786 C
787 C Apply Neumann boundary conditions to surface of scalar, S.
788 C Use IFIELD as a guide to which boundary conditions are to be applied.
789 C
790 C If ITYPE = 1, then S is returned as the rhs contribution to the
791 C volumetric flux.
792 C
793 C If ITYPE =-1, then S is returned as the lhs contribution to the
794 C diagonal of A.
795 C
796 C
797  include 'SIZE'
798  include 'TOTAL'
799  include 'CTIMER'
800  include 'NEKUSE'
801 C
802  dimension s(lx1,ly1,lz1,lelt)
803  common /nekcb/ cb
804  CHARACTER CB*3
805 C
806  if (icalld.eq.0) then
807  tusbc=0.0
808  nusbc=0
809  icalld=icalld+1
810  endif
811  nusbc=nusbc+1
812  etime1=dnekclock()
813 C
814  nfaces=2*ldim
815  nxyz =lx1*ly1*lz1
816  nel =nelfld(ifield)
817  ntot =nxyz*nel
818  CALL rzero(s,ntot)
819 C
820  IF (itype.EQ.-1) THEN
821 C
822 C Compute diagonal contributions to accomodate Robin boundary conditions
823 C
824  DO 1000 ie=1,nel
825  DO 1000 iface=1,nfaces
826  ieg=lglel(ie)
827  cb =cbc(iface,ie,ifield)
828  IF (cb.EQ.'C ' .OR. cb.EQ.'c ' .OR.
829  $ cb.EQ.'R ' .OR. cb.EQ.'r ') THEN
830 C
831  IF (cb.EQ.'C ') hc = bc(2,iface,ie,ifield)
832  IF (cb.EQ.'R ') THEN
833  tinf = bc(1,iface,ie,ifield)
834  hrad = bc(2,iface,ie,ifield)
835  ENDIF
836  ia=0
837 C
838 C IA is areal counter, assumes advancing fastest index first. (IX...IY...IZ)
839 C
840  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,iface)
841  DO 100 iz=kz1,kz2
842  DO 100 iy=ky1,ky2
843  DO 100 ix=kx1,kx2
844  ia = ia + 1
845  ts = t(ix,iy,iz,ie,ifield-1)
846  IF (cb.EQ.'c ' .OR. cb.EQ.'r ') THEN
847  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,ie)
848  CALL userbc (ix,iy,iz,iface,ieg)
849  ENDIF
850  IF (cb.EQ.'r ' .OR. cb.EQ.'R ')
851  $ hc = hrad * (tinf**2 + ts**2) * (tinf + ts)
852  s(ix,iy,iz,ie) = s(ix,iy,iz,ie) +
853  $ hc*area(ia,1,iface,ie)/bm1(ix,iy,iz,ie)
854  100 CONTINUE
855  ENDIF
856  1000 CONTINUE
857  ENDIF
858  IF (itype.EQ.1) THEN
859 C
860 C Add passive scalar fluxes to rhs
861 C
862  DO 2000 ie=1,nel
863  DO 2000 iface=1,nfaces
864  ieg=lglel(ie)
865  cb =cbc(iface,ie,ifield)
866  IF (cb.EQ.'F ' .OR. cb.EQ.'f ' .OR.
867  $ cb.EQ.'C ' .OR. cb.EQ.'c ' .OR.
868  $ cb.EQ.'R ' .OR. cb.EQ.'r ' ) THEN
869 C
870  IF (cb.EQ.'F ') flux=bc(1,iface,ie,ifield)
871  IF (cb.EQ.'C ') flux=bc(1,iface,ie,ifield)
872  $ *bc(2,iface,ie,ifield)
873  IF (cb.EQ.'R ') THEN
874  tinf=bc(1,iface,ie,ifield)
875  hrad=bc(2,iface,ie,ifield)
876  ENDIF
877 C
878 C Add local weighted flux values to rhs, S.
879 C
880 C IA is areal counter, assumes advancing fastest index first. (IX...IY...IZ)
881  ia=0
882  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,iface)
883  DO 200 iz=kz1,kz2
884  DO 200 iy=ky1,ky2
885  DO 200 ix=kx1,kx2
886  ia = ia + 1
887  ts = t(ix,iy,iz,ie,ifield-1)
888  IF (cb.EQ.'f ') THEN
889  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,ie)
890  CALL userbc (ix,iy,iz,iface,ieg)
891  ENDIF
892  IF (cb.EQ.'c ') THEN
893  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,ie)
894  CALL userbc (ix,iy,iz,iface,ieg)
895  flux = tinf*hc
896  ENDIF
897  IF (cb.EQ.'r ') THEN
898  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,ie)
899  CALL userbc (ix,iy,iz,iface,ieg)
900  ENDIF
901  IF (cb.EQ.'R ' .OR. cb.EQ.'r ')
902  $ flux = hrad*(tinf**2 + ts**2)*(tinf + ts) * tinf
903 C
904 C Add computed fluxes to boundary surfaces:
905 C
906  s(ix,iy,iz,ie) = s(ix,iy,iz,ie)
907  $ + flux*area(ia,1,iface,ie)
908  200 CONTINUE
909  ENDIF
910  2000 CONTINUE
911  ENDIF
912 C
913  tusbc=tusbc+(dnekclock()-etime1)
914 C
915  RETURN
916  END
917 c-----------------------------------------------------------------------
918  SUBROUTINE faceis (CB,S,IEL,IFACE,NX,NY,NZ)
919 C
920 C Assign inflow boundary conditions to face(IE,IFACE)
921 C for scalar S.
922 C
923  include 'SIZE'
924  include 'PARALLEL'
925  include 'NEKUSE'
926  include 'TSTEP' ! ifield 11/19/2010
927  include 'SOLN' ! tmask() 11/19/2010
928 C
929  dimension s(lx1,ly1,lz1)
930  CHARACTER CB*3
931 c
932  common /nekcb/ cb3
933  character*3 cb3
934  cb3 = cb
935 
936  ifld1 = ifield-1
937 
938 
939 C Passive scalar term
940 
941  ieg = lglel(iel)
942  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
943 
944  if (cb.eq.'t ') then
945  DO 100 iz=kz1,kz2 ! 11/19/2010: The tmask() screen
946  DO 100 iy=ky1,ky2 ! added here so users can leave
947  DO 100 ix=kx1,kx2 ! certain points to be Neumann,
948  if (tmask(ix,iy,iz,iel,ifld1).eq.0) then ! if desired.
949  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
950  CALL userbc (ix,iy,iz,iface,ieg)
951  s(ix,iy,iz) = temp
952  endif
953  100 CONTINUE
954  RETURN
955 
956  elseif (cb.eq.'o ' .or. cb.eq.'on ') then
957  DO 101 iz=kz1,kz2 ! 11/19/2010: The tmask() screen
958  DO 101 iy=ky1,ky2 ! added here so users can leave
959  DO 101 ix=kx1,kx2 ! certain points to be Neumann,
960  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
961  CALL userbc (ix,iy,iz,iface,ieg)
962  s(ix,iy,iz) = pa
963  101 CONTINUE
964  RETURN
965 
966  ELSEIF (cb.EQ.'ms ' .OR. cb.EQ.'msi') THEN
967 
968  DO 200 iz=kz1,kz2
969  DO 200 iy=ky1,ky2
970  DO 200 ix=kx1,kx2
971  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
972  CALL userbc (ix,iy,iz,iface,ieg)
973  s(ix,iy,iz) = sigma
974  200 CONTINUE
975 C
976  ELSEIF (cb.EQ.'kd ') THEN
977 C
978  DO 300 iz=kz1,kz2
979  DO 300 iy=ky1,ky2
980  DO 300 ix=kx1,kx2
981  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
982  CALL userbc (ix,iy,iz,iface,ieg)
983  s(ix,iy,iz) = turbk
984  300 CONTINUE
985 C
986  ELSEIF (cb.EQ.'ed ') THEN
987 C
988  DO 400 iz=kz1,kz2
989  DO 400 iy=ky1,ky2
990  DO 400 ix=kx1,kx2
991  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
992  CALL userbc (ix,iy,iz,iface,ieg)
993  s(ix,iy,iz) = turbe
994  400 CONTINUE
995 C
996  ENDIF
997 C
998  RETURN
999  END
1000 c-----------------------------------------------------------------------
1001  SUBROUTINE faceiv (CB,V1,V2,V3,IEL,IFACE,NX,NY,NZ)
1002 C
1003 C Assign fortran function boundary conditions to
1004 C face IFACE of element IEL for vector (V1,V2,V3).
1005 C
1006  include 'SIZE'
1007  include 'NEKUSE'
1008  include 'PARALLEL'
1009 C
1010  dimension v1(nx,ny,nz),v2(nx,ny,nz),v3(nx,ny,nz)
1011  character cb*3
1012 c
1013  character*1 cb1(3)
1014 c
1015  common /nekcb/ cb3
1016  character*3 cb3
1017  cb3 = cb
1018 c
1019  call chcopy(cb1,cb,3)
1020 c
1021  ieg = lglel(iel)
1022  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
1023 C
1024  IF (cb.EQ.'v ' .OR. cb.EQ.'ws ' .OR. cb.EQ.'mv '.OR.
1025  $ cb.EQ.'mvn') THEN
1026 C
1027  DO 100 iz=kz1,kz2
1028  DO 100 iy=ky1,ky2
1029  DO 100 ix=kx1,kx2
1030  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
1031  CALL userbc (ix,iy,iz,iface,ieg)
1032  v1(ix,iy,iz) = ux
1033  v2(ix,iy,iz) = uy
1034  v3(ix,iy,iz) = uz
1035  100 CONTINUE
1036  RETURN
1037 C
1038  elseif (cb1(1).eq.'d'.or.cb1(2).eq.'d'.or.cb1(3).eq.'d') then
1039 C
1040  do iz=kz1,kz2
1041  do iy=ky1,ky2
1042  do ix=kx1,kx2
1043  if (optlevel.le.2) call nekasgn (ix,iy,iz,iel)
1044  call userbc (ix,iy,iz,iface,ieg)
1045  if (cb1(1).eq.'d') v1(ix,iy,iz) = ux
1046  if (cb1(2).eq.'d') v2(ix,iy,iz) = uy
1047  if (cb1(3).eq.'d') v3(ix,iy,iz) = uz
1048  enddo
1049  enddo
1050  enddo
1051  return
1052 C
1053  ELSEIF (cb.EQ.'vl ' .OR. cb.EQ.'wsl') THEN
1054 C
1055  DO 120 iz=kz1,kz2
1056  DO 120 iy=ky1,ky2
1057  DO 120 ix=kx1,kx2
1058  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
1059  CALL userbc (ix,iy,iz,iface,ieg)
1060  v1(ix,iy,iz) = un
1061  v2(ix,iy,iz) = u1
1062  v3(ix,iy,iz) = u2
1063  120 CONTINUE
1064  RETURN
1065 C
1066  ELSEIF (cb.EQ.'s ' .OR. cb.EQ.'sh ') THEN
1067 C
1068  DO 200 iz=kz1,kz2
1069  DO 200 iy=ky1,ky2
1070  DO 200 ix=kx1,kx2
1071  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
1072  CALL userbc (ix,iy,iz,iface,ieg)
1073  v1(ix,iy,iz) = trx
1074  v2(ix,iy,iz) = try
1075  v3(ix,iy,iz) = trz
1076  200 CONTINUE
1077  RETURN
1078 C
1079  ELSEIF (cb.EQ.'sl ' .OR. cb.EQ.'shl') THEN
1080 C
1081  DO 220 iz=kz1,kz2
1082  DO 220 iy=ky1,ky2
1083  DO 220 ix=kx1,kx2
1084  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
1085  CALL userbc (ix,iy,iz,iface,ieg)
1086  v1(ix,iy,iz) = trn
1087  v2(ix,iy,iz) = tr1
1088  v3(ix,iy,iz) = tr2
1089  220 CONTINUE
1090 C
1091  ELSEIF (cb.EQ.'ms ') THEN
1092 C
1093  DO 240 iz=kz1,kz2
1094  DO 240 iy=ky1,ky2
1095  DO 240 ix=kx1,kx2
1096  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
1097  CALL userbc (ix,iy,iz,iface,ieg)
1098  v1(ix,iy,iz) = -pa
1099  v2(ix,iy,iz) = tr1
1100  v3(ix,iy,iz) = tr2
1101  240 CONTINUE
1102 C
1103  ELSEIF (cb.EQ.'on ' .OR. cb.EQ.'o ') THEN
1104 C
1105  DO 270 iz=kz1,kz2
1106  DO 270 iy=ky1,ky2
1107  DO 270 ix=kx1,kx2
1108  if (optlevel.le.2) CALL nekasgn (ix,iy,iz,iel)
1109  CALL userbc (ix,iy,iz,iface,ieg)
1110  v1(ix,iy,iz) = -pa
1111  v2(ix,iy,iz) = 0.0
1112  v3(ix,iy,iz) = 0.0
1113  270 CONTINUE
1114 C
1115  ENDIF
1116 C
1117  RETURN
1118  END
1119 c-----------------------------------------------------------------------
1120  subroutine nekasgn (ix,iy,iz,e)
1121 C
1122 C Assign NEKTON variables for definition (by user) of
1123 C boundary conditions at collocation point (IX,IY,IZ)
1124 C of element IEL.
1125 C
1126 C X X-coordinate
1127 C Y Y-coordinate
1128 C Z Z-coordinate
1129 C UX X-velocity
1130 C UY Y-velocity
1131 C UZ Z-velocity
1132 C TEMP Temperature
1133 C PS1 Passive scalar No. 1
1134 C PS2 Passive scalar No. 2
1135 C . .
1136 C . .
1137 C PS9 Passive scalar No. 9
1138 C SI2 Strainrate invariant II
1139 C SI3 Strainrate invariant III
1140 C
1141 C Variables to be defined by user for imposition of
1142 C boundary conditions :
1143 C
1144 C SH1 Shear component No. 1
1145 C SH2 Shear component No. 2
1146 C TRX X-traction
1147 C TRY Y-traction
1148 C TRZ Z-traction
1149 C SIGMA Surface-tension coefficient
1150 C FLUX Flux
1151 C HC Convection heat transfer coefficient
1152 C HRAD Radiation heat transfer coefficient
1153 C TINF Temperature at infinity
1154 C
1155  include 'SIZE'
1156  include 'GEOM'
1157  include 'SOLN'
1158  include 'INPUT'
1159  include 'TSTEP'
1160  include 'NEKUSE'
1161 
1162  integer e
1163 
1164  common /nekcb/ cb
1165  character cb*3
1166 
1167  COMMON /screv / sii(lx1,ly1,lz1,lelt)
1168  $ , siii(lx1,ly1,lz1,lelt)
1169 
1170  x = xm1(ix,iy,iz,e)
1171  y = ym1(ix,iy,iz,e)
1172  z = zm1(ix,iy,iz,e)
1173 
1174  r = x**2+y**2
1175  if (r.gt.0.0) r=sqrt(r)
1176  if (x.ne.0.0 .or. y.ne.0.0) theta = atan2(y,x)
1177 
1178  ux = vx(ix,iy,iz,e)
1179  uy = vy(ix,iy,iz,e)
1180  uz = vz(ix,iy,iz,e)
1181  temp = t(ix,iy,iz,e,1)
1182  do ips=1,npscal
1183  ps(ips) = t(ix,iy,iz,e,ips+1)
1184  enddo
1185 
1186  pa = pr(ix,iy,iz,e)
1187  p0 = p0th
1188 
1189  si2 = sii(ix,iy,iz,e)
1190  si3 = siii(ix,iy,iz,e)
1191  udiff = vdiff(ix,iy,iz,e,ifield)
1192  utrans= vtrans(ix,iy,iz,e,ifield)
1193 
1194  cbu = cb
1195 
1196  return
1197  end
1198 c-----------------------------------------------------------------------
1199  SUBROUTINE bcneutr
1200 C
1201  include 'SIZE'
1202  include 'SOLN'
1203  include 'GEOM'
1204  include 'INPUT'
1205  COMMON /scrsf/ trx(lx1,ly1,lz1)
1206  $ , try(lx1,ly1,lz1)
1207  $ , trz(lx1,ly1,lz1)
1208  COMMON /ctmp0/ stc(lx1,ly1,lz1)
1209  REAL SIGST(LX1,LY1)
1210 C
1211  LOGICAL IFALGN,IFNORX,IFNORY,IFNORZ
1212  common /nekcb/ cb
1213  CHARACTER CB*3
1214 C
1215  ifld = 1
1216  nface = 2*ldim
1217  nxy1 = lx1*ly1
1218  nxyz1 = lx1*ly1*lz1
1219 C
1220  DO 100 iel=1,nelv
1221  DO 100 ifc=1,nface
1222 C
1223  cb = cbc(ifc,iel,ifld)
1224  bc1 = bc(1,ifc,iel,ifld)
1225  bc2 = bc(2,ifc,iel,ifld)
1226  bc3 = bc(3,ifc,iel,ifld)
1227  bc4 = bc(4,ifc,iel,ifld)
1228  CALL rzero3 (trx,try,trz,nxyz1)
1229 C
1230 C Prescribed tractions and shear tractions
1231 C
1232  IF (cb.EQ.'S ' .OR. cb.EQ.'SL ' .OR.
1233  $ cb.EQ.'SH ' .OR. cb.EQ.'SHL' ) THEN
1234  CALL trcon (trx,try,trz,bc1,bc2,bc3,iel,ifc)
1235  IF (ifqinp(ifc,iel)) CALL globrot (trx,try,trz,iel,ifc)
1236  GOTO 120
1237  ENDIF
1238  IF (cb.EQ.'s ' .OR. cb.EQ.'sl ' .OR.
1239  $ cb.EQ.'sh ' .OR. cb.EQ.'shl' ) THEN
1240  CALL faceiv (cb,trx,try,trz,iel,ifc,lx1,ly1,lz1)
1241  CALL faccvs (trx,try,trz,area(1,1,ifc,iel),ifc)
1242  IF (ifqinp(ifc,iel)) CALL globrot (trx,try,trz,iel,ifc)
1243  GOTO 120
1244  ENDIF
1245 C
1246 C Prescribed outflow ambient pressure
1247 C
1248  IF (cb.EQ.'ON ' .OR. cb.EQ.'O ') THEN
1249  bcn = -bc1
1250  bc2 = 0.0
1251  bc3 = 0.0
1252  CALL trcon (trx,try,trz,bcn,bc2,bc3,iel,ifc)
1253  CALL globrot (trx,try,trz,iel,ifc)
1254  GOTO 120
1255  ENDIF
1256  IF (cb.EQ.'on ' .OR. cb.EQ.'o ') THEN
1257  CALL faceiv (cb,trx,try,trz,iel,ifc,lx1,ly1,lz1)
1258  CALL faccvs (trx,try,trz,area(1,1,ifc,iel),ifc)
1259  CALL globrot (trx,try,trz,iel,ifc)
1260  GOTO 120
1261  ENDIF
1262 C
1263 C Surface-tension
1264 C
1265  IF (cb.EQ.'MS ' .OR. cb.EQ.'MSI' .OR.
1266  $ cb.EQ.'MM ' .OR. cb.EQ.'mm ' .OR.
1267  $ cb.EQ.'ms ' .OR. cb.EQ.'msi') THEN
1268  IF (cb.EQ.'MS '.or.cb.eq.'MM ') THEN
1269  bcn = -bc1
1270  CALL trcon (trx,try,trz,bcn,bc2,bc3,iel,ifc)
1271  CALL globrot (trx,try,trz,iel,ifc)
1272  ENDIF
1273 c IF (CB.EQ.'ms '.or.cb.eq.'mm ') THEN
1274  IF (cb.EQ.'ms '.or.cb.eq.'msi') THEN
1275  CALL faceiv (cb,trx,try,trz,iel,ifc,lx1,ly1,lz1)
1276  CALL faccvs (trx,try,trz,area(1,1,ifc,iel),ifc)
1277  CALL globrot (trx,try,trz,iel,ifc)
1278  ENDIF
1279  IF (cb(1:1).EQ.'M') THEN
1280  CALL cfill (sigst,bc4,nxy1)
1281  ELSE
1282  CALL faceis (cb,stc,iel,ifc,lx1,ly1,lz1)
1283  CALL facexs (sigst,stc,ifc,0)
1284  ENDIF
1285  IF (ifaxis) THEN
1286  CALL trstax (trx,try,sigst,iel,ifc)
1287  ELSEIF (ldim.EQ.2) THEN
1288  CALL trst2d (trx,try,sigst,iel,ifc)
1289  ELSE
1290  CALL trst3d (trx,try,trz,sigst,iel,ifc)
1291  ENDIF
1292  ENDIF
1293 C
1294  120 CALL add2 (bfx(1,1,1,iel),trx,nxyz1)
1295  CALL add2 (bfy(1,1,1,iel),try,nxyz1)
1296  IF (ldim.EQ.3) CALL add2 (bfz(1,1,1,iel),trz,nxyz1)
1297 C
1298  100 CONTINUE
1299 C
1300  RETURN
1301  END
1302 c-----------------------------------------------------------------------
1303  SUBROUTINE trcon (TRX,TRY,TRZ,TR1,TR2,TR3,IEL,IFC)
1304 C
1305  include 'SIZE'
1306  include 'GEOM'
1307  include 'TOPOL'
1308 C
1309  dimension trx(lx1,ly1,lz1)
1310  $ , try(lx1,ly1,lz1)
1311  $ , trz(lx1,ly1,lz1)
1312 C
1313  CALL dsset(lx1,ly1,lz1)
1314  iface = eface1(ifc)
1315  js1 = skpdat(1,iface)
1316  jf1 = skpdat(2,iface)
1317  jskip1 = skpdat(3,iface)
1318  js2 = skpdat(4,iface)
1319  jf2 = skpdat(5,iface)
1320  jskip2 = skpdat(6,iface)
1321  i = 0
1322 C
1323  IF (ldim.EQ.2) THEN
1324  DO 100 j2=js2,jf2,jskip2
1325  DO 100 j1=js1,jf1,jskip1
1326  i = i + 1
1327  trx(j1,j2,1) = tr1*area(i,1,ifc,iel)
1328  try(j1,j2,1) = tr2*area(i,1,ifc,iel)
1329  100 CONTINUE
1330  ELSE
1331  DO 200 j2=js2,jf2,jskip2
1332  DO 200 j1=js1,jf1,jskip1
1333  i = i + 1
1334  trx(j1,j2,1) = tr1*area(i,1,ifc,iel)
1335  try(j1,j2,1) = tr2*area(i,1,ifc,iel)
1336  trz(j1,j2,1) = tr3*area(i,1,ifc,iel)
1337  200 CONTINUE
1338  ENDIF
1339 C
1340  RETURN
1341  END
1342 c-----------------------------------------------------------------------
1343  SUBROUTINE trst2d (TRX,TRY,SIGST,IEL,IFC)
1344 C
1345 C Compute taction due to surface tension (2D)
1346 C
1347  include 'SIZE'
1348  include 'GEOM'
1349  include 'DXYZ'
1350  include 'TOPOL'
1351  include 'WZ'
1352  COMMON /ctmp1/ a1x(lx1),a1y(lx1),stx(lx1),sty(lx1)
1353 C
1354  dimension trx(lx1,ly1,lz1),try(lx1,ly1,lz1),sigst(lx1,1)
1355  dimension cang(2),sang(2)
1356  dimension ixn(2),iyn(2),ian(2)
1357 C
1358  DO 100 ix=1,lx1
1359  aa = sigst(ix,1) * wxm1(ix)
1360  stx(ix) = t1x(ix,1,ifc,iel) * aa
1361  sty(ix) = t1y(ix,1,ifc,iel) * aa
1362  100 CONTINUE
1363 C
1364  IF (ifc.EQ.3 .OR. ifc.EQ.4) THEN
1365  CALL chsign (stx,lx1)
1366  CALL chsign (sty,lx1)
1367  ENDIF
1368 C
1369  IF (ifc.EQ.1 .OR. ifc.EQ.3) THEN
1370  CALL mxm (dxtm1,lx1,stx,lx1,a1x,1)
1371  CALL mxm (dxtm1,lx1,sty,lx1,a1y,1)
1372  ELSE
1373  CALL mxm (dytm1,ly1,stx,ly1,a1x,1)
1374  CALL mxm (dytm1,ly1,sty,ly1,a1y,1)
1375  ENDIF
1376 C
1377  CALL dsset (lx1,ly1,lz1)
1378  iface = eface1(ifc)
1379  js1 = skpdat(1,iface)
1380  jf1 = skpdat(2,iface)
1381  jskip1 = skpdat(3,iface)
1382  js2 = skpdat(4,iface)
1383  jf2 = skpdat(5,iface)
1384  jskip2 = skpdat(6,iface)
1385  i = 0
1386 C
1387  DO 200 j2=js2,jf2,jskip2
1388  DO 200 j1=js1,jf1,jskip1
1389  i = i + 1
1390  trx(j1,j2,1) = trx(j1,j2,1) - a1x(i)
1391  try(j1,j2,1) = try(j1,j2,1) - a1y(i)
1392  200 CONTINUE
1393 C
1394 C Contact angle corrections
1395 C
1396  CALL ctang2d (cang,sang,ixn,iyn,ian,ifc,iel)
1397  DO 500 i=1,2
1398  ix = ixn(i)
1399  iy = iyn(i)
1400  ia = ian(i)
1401  trx(ix,iy,1)=trx(ix,iy,1) + sigst(ia,1)*cang(i)
1402  try(ix,iy,1)=try(ix,iy,1) + sigst(ia,1)*sang(i)
1403  500 CONTINUE
1404 C
1405  RETURN
1406  END
1407 c-----------------------------------------------------------------------
1408  SUBROUTINE trstax (TRX,TRY,SIGST,IEL,IFC)
1409 C
1410 C Compute taction due to surface tension (axisymmetric)
1411 C
1412  include 'SIZE'
1413  include 'GEOM'
1414  include 'DXYZ'
1415  include 'TOPOL'
1416  include 'WZ'
1417  COMMON /ctmp1/ a1x(lx1),a1y(lx1),a2x(lx1),a2y(lx1)
1418  $ , stx(lx1),sty(lx1),xjm1(lx1)
1419  COMMON /ctmp0/ xfm1(lx1),yfm1(lx1),t1xf(lx1),t1yf(lx1)
1420 C
1421  dimension trx(lx1,ly1,lz1),try(lx1,ly1,lz1),sigst(lx1,ly1)
1422  dimension cang(2),sang(2)
1423  dimension ixn(2),iyn(2),ian(2)
1424  LOGICAL IFGLJ
1425 C
1426  ifglj = .false.
1427  IF ( ifrzer(iel) .AND. (ifc.EQ.2 .OR. ifc.EQ.4) ) ifglj = .true.
1428  CALL facec2 (xfm1,yfm1,xm1(1,1,1,iel),ym1(1,1,1,iel),ifc)
1429 C
1430  IF (ifglj) THEN
1431  CALL mxm (dam1,ly1,xfm1,ly1,t1xf,1)
1432  CALL mxm (dam1,ly1,yfm1,ly1,t1yf,1)
1433  ys0 = t1yf(1)
1434  ELSE
1435  CALL mxm (dxm1,lx1,xfm1,lx1,t1xf,1)
1436  CALL mxm (dxm1,lx1,yfm1,lx1,t1yf,1)
1437  ENDIF
1438 C
1439  DO 10 ix=1,lx1
1440  xjm1(ix)=sqrt( t1xf(ix)**2 + t1yf(ix)**2 )
1441  t1xf(ix)=t1xf(ix) / xjm1(ix)
1442  t1yf(ix)=t1yf(ix) / xjm1(ix)
1443  10 CONTINUE
1444 C
1445  IF ( ifglj ) THEN
1446  CALL mxm (dam1,1,t1xf,ly1,t1xs0,1)
1447  CALL mxm (dam1,1,uny(1,1,ifc,iel),ly1,unys0,1)
1448  ddx = wam1(1)*sigst(1,1)*t1xs0*ys0
1449  ddy = wam1(1)*sigst(1,1)*t1yf(1)*ys0*2.0
1450  a2x(1) = wam1(1)*sigst(1,1)*xjm1(1)*unx(1,1,ifc,iel)*unys0
1451  a2y(1) = 0.0
1452  stx(1) = 0.0
1453  sty(1) = 0.0
1454  DO 100 iy=2,ly1
1455  aa = wam1(iy) * sigst(iy,1) / (1.0 + zam1(iy))
1456  stx(iy) = t1xf(iy) * aa
1457  sty(iy) = t1yf(iy) * aa
1458  aa = aa * xjm1(iy) * uny(iy,1,ifc,iel)
1459  a2x(iy) = unx(iy,1,ifc,iel) * aa
1460  a2y(iy) = uny(iy,1,ifc,iel) * aa
1461  100 CONTINUE
1462  ELSE
1463  DO 200 ix=1,lx1
1464  aa = sigst(ix,1) * wxm1(ix)
1465  stx(ix) = t1xf(ix) * aa
1466  sty(ix) = t1yf(ix) * aa
1467  aa = aa * xjm1(ix) * uny(ix,1,ifc,iel)
1468  a2x(ix) = unx(ix,1,ifc,iel) * aa
1469  a2y(ix) = uny(ix,1,ifc,iel) * aa
1470  200 CONTINUE
1471  ENDIF
1472 C
1473  IF (ifglj) THEN
1474  DO 220 iy=1,ly1
1475  ysiy = t1yf(iy)*xjm1(iy)
1476  dtx1 = 0.0
1477  dty1 = datm1(iy,1)*ddy
1478  dtx2 = ysiy*stx(iy)
1479  dty2 = ysiy*sty(iy)
1480  dty3 = 0.0
1481  DO 240 j=2,ly1
1482  dtys = datm1(iy,j)*yfm1(j)
1483  dtx1 = dtx1 + dtys*stx(j)
1484  dty3 = dty3 + dtys*sty(j)
1485  240 CONTINUE
1486  a1x(iy) = dtx1 + dtx2
1487  a1y(iy) = dty1 + dty2 + dty3
1488  220 CONTINUE
1489  a1x(1) = a1x(1) + ddx
1490  ELSE
1491  CALL mxm (dxtm1,lx1,stx,lx1,a1x,1)
1492  CALL mxm (dxtm1,lx1,sty,lx1,a1y,1)
1493  CALL col2 (a1x,yfm1,lx1)
1494  CALL col2 (a1y,yfm1,lx1)
1495  ENDIF
1496 C
1497  CALL dsset (lx1,ly1,lz1)
1498  iface = eface1(ifc)
1499  js1 = skpdat(1,iface)
1500  jf1 = skpdat(2,iface)
1501  jskip1 = skpdat(3,iface)
1502  js2 = skpdat(4,iface)
1503  jf2 = skpdat(5,iface)
1504  jskip2 = skpdat(6,iface)
1505  i = 0
1506 C
1507  DO 300 j2=js2,jf2,jskip2
1508  DO 300 j1=js1,jf1,jskip1
1509  i = i + 1
1510  trx(j1,j2,1) = trx(j1,j2,1) - a2x(i) - a1x(i)
1511  try(j1,j2,1) = try(j1,j2,1) - a2y(i) - a1y(i)
1512  300 CONTINUE
1513 C
1514 C Contact angle corrections
1515 C
1516  CALL ctang2d (cang,sang,ixn,iyn,ian,ifc,iel)
1517  DO 500 i=1,2
1518  ix = ixn(i)
1519  iy = iyn(i)
1520  ia = ian(i)
1521  aa = sigst(ia,1)*ym1(ix,iy,1,iel)
1522  trx(ix,iy,1)=trx(ix,iy,1) + aa*cang(i)
1523  try(ix,iy,1)=try(ix,iy,1) + aa*sang(i)
1524  500 CONTINUE
1525 C
1526  RETURN
1527  END
1528 c-----------------------------------------------------------------------
1529  SUBROUTINE ctang2d (CANG,SANG,IXN,IYN,IAN,IFC,IEL)
1530 C
1531  include 'SIZE'
1532  include 'GEOM'
1533  include 'SOLN'
1534  include 'INPUT'
1535 C
1536  dimension cang(2),sang(2)
1537  dimension ixn(2),iyn(2),ian(2),isn(2),nebpt(4,2)
1538  CHARACTER CBN*3
1539 C
1540  DATA nebpt /4,1,2,3, 2,3,4,1/
1541  ifld = 1
1542  eps = 1.e-6
1543 C
1544  DO 100 i=1,2
1545  ifcn = nebpt(ifc,i)
1546  cbn = cbc(ifcn,iel,ifld)
1547  ixn(i) = 1
1548  iyn(i) = 1
1549  ian(i) = 1
1550  isn(i) = 1
1551  cang(i) = 0.0
1552  sang(i) = 0.0
1553  IF (cbn.EQ.'E '.OR.cbn.EQ.'P '.OR.cbn.eq.'p '.or.
1554  $ cbn(1:1).EQ.'M' .OR. cbn(1:1).EQ.'m') GOTO 100
1555  nc = ifc
1556  IF (i.EQ.2) nc=ifcn
1557  IF (nc .EQ.2 .OR. nc .EQ.3) ixn(i) = lx1
1558  IF (nc .EQ.3 .OR. nc .EQ.4) iyn(i) = ly1
1559  IF (ifc .EQ.2 .OR. ifc .EQ.3) isn(i) = lx1
1560  IF (ifcn.EQ.2 .OR. ifcn.EQ.3) ian(i) = lx1
1561  ix = ixn(i)
1562  iy = iyn(i)
1563  ia = ian(i)
1564  is = isn(i)
1565  IF (cbn(1:1).EQ.'V' .OR. cbn(1:1).EQ.'v' .OR.
1566  $ cbn .EQ.'S ' .OR. cbn .EQ.'s ' .OR.
1567  $ cbn .EQ.'SL ' .OR. cbn .EQ.'sl ' .OR.
1568  $ cbn(1:1).EQ.'O' .OR. cbn(1:1).EQ.'o' ) THEN
1569  ux=vx(ix,iy,1,iel)
1570  uy=vy(ix,iy,1,iel)
1571  um=ux**2 + uy**2
1572  IF (um.GT.eps) THEN
1573  unlx=unx(is,1,ifcn,iel)
1574  unly=uny(is,1,ifcn,iel)
1575  um=sqrt(um)
1576  dot =ux*unlx + uy*unly
1577  IF (dot.LT.0.0) um=-um
1578  cang(i)=ux/um
1579  sang(i)=uy/um
1580  GOTO 100
1581  ENDIF
1582  ENDIF
1583  cang(i)=unx(is,1,ifcn,iel)
1584  sang(i)=uny(is,1,ifcn,iel)
1585  100 CONTINUE
1586 C
1587  RETURN
1588  END
1589 c-----------------------------------------------------------------------
1590  SUBROUTINE trst3d (TRX,TRY,TRZ,SIGST,IEL,IFC)
1591 C
1592 C Compute taction due to surface tension (3D)
1593 C
1594  include 'SIZE'
1595  include 'GEOM'
1596  include 'WZ'
1597  COMMON /ctmp0/ xfm1(lx1,ly1),yfm1(lx1,ly1),zfm1(lx1,ly1)
1598  COMMON /ctmp1/ drm1(lx1,lx1),drtm1(lx1,ly1)
1599  $ , dsm1(lx1,lx1),dstm1(lx1,ly1)
1600  $ , wgs(lx1,ly1)
1601  COMMON /scrmg/ xrm1(lx1,ly1),yrm1(lx1,ly1),zrm1(lx1,ly1)
1602  $ , xsm1(lx1,ly1),ysm1(lx1,ly1),zsm1(lx1,ly1)
1603  COMMON /scruz/ s1x(lx1,ly1),s1y(lx1,ly1),s1z(lx1,ly1)
1604  $ , s2x(lx1,ly1),s2y(lx1,ly1),s2z(lx1,ly1)
1605  COMMON /scrns/ g1x(lx1,ly1),g1y(lx1,ly1),g1z(lx1,ly1)
1606  $ , g2x(lx1,ly1),g2y(lx1,ly1),g2z(lx1,ly1)
1607  $ , gbs(lx1,ly1),gb1l(lx1,ly1),gb2l(lx1,ly1)
1608 C
1609  dimension trx(lx1,ly1,lz1),try(lx1,ly1,lz1),trz(lx1,ly1,lz1)
1610  dimension sigst(lx1,ly1)
1611 C
1612  nxy1 = lx1*ly1
1613 C
1614  CALL rzero3 (s1x,s1y,s1z,nxy1)
1615  CALL rzero3 (s2x,s2y,s2z,nxy1)
1616  CALL facexv (xfm1,yfm1,zfm1,xm1(1,1,1,iel),ym1(1,1,1,iel),
1617  $ zm1(1,1,1,iel),ifc,0)
1618  CALL setdrs (drm1,drtm1,dsm1,dstm1,ifc)
1619 C
1620  CALL mxm (drm1,lx1, xfm1,lx1,xrm1,ly1)
1621  CALL mxm (drm1,lx1, yfm1,lx1,yrm1,ly1)
1622  CALL mxm (drm1,lx1, zfm1,lx1,zrm1,ly1)
1623  CALL mxm (xfm1,lx1,dstm1,ly1,xsm1,ly1)
1624  CALL mxm (yfm1,lx1,dstm1,ly1,ysm1,ly1)
1625  CALL mxm (zfm1,lx1,dstm1,ly1,zsm1,ly1)
1626 C
1627  DO 100 ix=1,lx1
1628  DO 100 iy=1,ly1
1629  gb1x=xrm1(ix,iy)
1630  gb1y=yrm1(ix,iy)
1631  gb1z=zrm1(ix,iy)
1632  gb2x=xsm1(ix,iy)
1633  gb2y=ysm1(ix,iy)
1634  gb2z=zsm1(ix,iy)
1635  gb11=gb1x*gb1x + gb1y*gb1y + gb1z*gb1z
1636  gb12=gb1x*gb2x + gb1y*gb2y + gb1z*gb2z
1637  gb22=gb2x*gb2x + gb2y*gb2y + gb2z*gb2z
1638  gdet=gb11*gb22 - gb12*gb12
1639  IF (gdet .LT. 1.e-20) GO TO 9001
1640  gt11= gb22/gdet
1641  gt12=-gb12/gdet
1642  gt22= gb11/gdet
1643  gb1l(ix,iy)=sqrt(gb11)
1644  gb2l(ix,iy)=sqrt(gb22)
1645  gbs(ix,iy)=sqrt(gdet)
1646  wgs(ix,iy)=wxm1(ix)*wym1(iy)*sigst(ix,iy)
1647  bb = gbs(ix,iy) * wgs(ix,iy)
1648  g1x(ix,iy) = bb * ( gt11*gb1x + gt12*gb2x )
1649  g1y(ix,iy) = bb * ( gt11*gb1y + gt12*gb2y )
1650  g1z(ix,iy) = bb * ( gt11*gb1z + gt12*gb2z )
1651  g2x(ix,iy) = bb * ( gt12*gb1x + gt22*gb2x )
1652  g2y(ix,iy) = bb * ( gt12*gb1y + gt22*gb2y )
1653  g2z(ix,iy) = bb * ( gt12*gb1z + gt22*gb2z )
1654  100 CONTINUE
1655 C
1656  CALL mxm (drtm1,lx1,g1x,lx1,s1x,ly1)
1657  CALL mxm (drtm1,lx1,g1y,lx1,s1y,ly1)
1658  CALL mxm (drtm1,lx1,g1z,lx1,s1z,ly1)
1659 C
1660  CALL mxm (g2x,lx1,dsm1,ly1,s2x,ly1)
1661  CALL mxm (g2y,lx1,dsm1,ly1,s2y,ly1)
1662  CALL mxm (g2z,lx1,dsm1,ly1,s2z,ly1)
1663 C
1664  CALL add2 (s1x,s2x,nxy1)
1665  CALL add2 (s1y,s2y,nxy1)
1666  CALL add2 (s1z,s2z,nxy1)
1667 C
1668 C Contact angle option on hold
1669 C
1670 C ICONTAC=INT(BC2)
1671 C IF (ICONTAC.NE.0) THEN
1672 C IX=1
1673 C IY=1
1674 C IF (ICONTAC.GE.3) IY=ly1
1675 C IF (ICONTAC.EQ.2 .OR. ICONTAC.EQ.3) IX=lx1
1676 C ANG = BC3 * PI / 180.00
1677 C RR = YM1(IX,IY,IZ,IEL)
1678 C TRX(IX,IY,IZ)=TRX(IX,IY,IZ) + RR*SIGST*COS( ANG )
1679 C TRY(IX,IY,IZ)=TRY(IX,IY,IZ) + RR*SIGST*SIN( ANG )
1680 C ENDIF
1681 C
1682  CALL facsub2 (trx,try,trz,s1x,s1y,s1z,ifc)
1683 C
1684  RETURN
1685 C
1686  9001 WRITE ( 6,*) 'Zero area for Element=',iel,' Face=',ifc
1687  call exitt
1688 C
1689  END
1690 c-----------------------------------------------------------------------
1691  SUBROUTINE setdrs (DRM1,DRTM1,DSM1,DSTM1,IFC)
1692 C
1693  include 'SIZE'
1694  include 'DXYZ'
1695 C
1696  dimension drm1(lx1,lx1),drtm1(lx1,lx1)
1697  $ , dsm1(ly1,ly1),dstm1(ly1,ly1)
1698 C
1699  nxy1=lx1*ly1
1700 C
1701  IF (ifc.EQ.5 .OR. ifc.EQ.6) THEN
1702  CALL copy (drm1 ,dxm1 ,nxy1)
1703  CALL copy (dsm1 ,dym1 ,nxy1)
1704  CALL copy (drtm1,dxtm1,nxy1)
1705  CALL copy (dstm1,dytm1,nxy1)
1706  ELSEIF (ifc.EQ.2 .OR. ifc.EQ.4) THEN
1707  CALL copy (drm1 ,dym1 ,nxy1)
1708  CALL copy (dsm1 ,dzm1 ,nxy1)
1709  CALL copy (drtm1,dytm1,nxy1)
1710  CALL copy (dstm1,dztm1 ,nxy1)
1711  ELSE
1712  CALL copy (drm1 ,dzm1 ,nxy1)
1713  CALL copy (dsm1 ,dxm1 ,nxy1)
1714  CALL copy (drtm1,dztm1,nxy1)
1715  CALL copy (dstm1,dxtm1,nxy1)
1716  ENDIF
1717 C
1718  RETURN
1719  END
1720 c-----------------------------------------------------------------------
1721  SUBROUTINE globrot (R1,R2,R3,IEL,IFC)
1722 C
1723 C Rotate vector components R1,R2,R3 at face IFC
1724 C of element IEL from local to global system.
1725 C
1726 C R1, R2, R3 have the (NX,NY,NZ) data structure
1727 C IFACE1 is in the preprocessor notation
1728 C IFACE is the dssum notation.
1729 C
1730  include 'SIZE'
1731  include 'GEOM'
1732  include 'TOPOL'
1733 C
1734  dimension r1(lx1,ly1,lz1)
1735  $ , r2(lx1,ly1,lz1)
1736  $ , r3(lx1,ly1,lz1)
1737 C
1738  CALL dsset (lx1,ly1,lz1)
1739  iface = eface1(ifc)
1740  js1 = skpdat(1,iface)
1741  jf1 = skpdat(2,iface)
1742  jskip1 = skpdat(3,iface)
1743  js2 = skpdat(4,iface)
1744  jf2 = skpdat(5,iface)
1745  jskip2 = skpdat(6,iface)
1746  i = 0
1747 C
1748  IF (ldim.EQ.2) THEN
1749  DO 200 j2=js2,jf2,jskip2
1750  DO 200 j1=js1,jf1,jskip1
1751  i = i+1
1752  rnorl = r1(j1,j2,1)
1753  rtan1 = r2(j1,j2,1)
1754  r1(j1,j2,1) = rnorl*unx(i,1,ifc,iel) +
1755  $ rtan1*t1x(i,1,ifc,iel)
1756  r2(j1,j2,1) = rnorl*uny(i,1,ifc,iel) +
1757  $ rtan1*t1y(i,1,ifc,iel)
1758  200 CONTINUE
1759  ELSE
1760  DO 300 j2=js2,jf2,jskip2
1761  DO 300 j1=js1,jf1,jskip1
1762  i = i+1
1763  rnorl = r1(j1,j2,1)
1764  rtan1 = r2(j1,j2,1)
1765  rtan2 = r3(j1,j2,1)
1766  r1(j1,j2,1) = rnorl*unx(i,1,ifc,iel) +
1767  $ rtan1*t1x(i,1,ifc,iel) +
1768  $ rtan2*t2x(i,1,ifc,iel)
1769  r2(j1,j2,1) = rnorl*uny(i,1,ifc,iel) +
1770  $ rtan1*t1y(i,1,ifc,iel) +
1771  $ rtan2*t2y(i,1,ifc,iel)
1772  r3(j1,j2,1) = rnorl*unz(i,1,ifc,iel) +
1773  $ rtan1*t1z(i,1,ifc,iel) +
1774  $ rtan2*t2z(i,1,ifc,iel)
1775  300 CONTINUE
1776  ENDIF
1777 C
1778  RETURN
1779  END
1780 c-----------------------------------------------------------------------
1781  SUBROUTINE facec2 (A1,A2,B1,B2,IFC)
1782 C
1783 C 2-D Geometry only
1784 C Extract A1,A2 from B1,B2 on surface IFC.
1785 C
1786 C A1, A2 have the (lx1, 1,NFACE) data structure
1787 C B1, B2 have the (lx1,ly1, 1) data structure
1788 C
1789  include 'SIZE'
1790 C
1791  dimension a1(lx1),a2(lx1),b1(lx1,ly1),b2(lx1,ly1)
1792 C
1793  ix=1
1794  iy=1
1795  IF (ifc.EQ.1 .OR. ifc.EQ.3) THEN
1796  IF (ifc.EQ.3) iy = ly1
1797  DO 10 ix=1,lx1
1798  a1(ix)=b1(ix,iy)
1799  a2(ix)=b2(ix,iy)
1800  10 CONTINUE
1801  ELSE
1802  IF (ifc.EQ.2) ix = lx1
1803  DO 20 iy=1,ly1
1804  a1(iy)=b1(ix,iy)
1805  a2(iy)=b2(ix,iy)
1806  20 CONTINUE
1807  ENDIF
1808 C
1809  RETURN
1810  END
1811 c-----------------------------------------------------------------------
1812  SUBROUTINE lfalse (IFA,N)
1813  LOGICAL IFA(1)
1814  DO 100 i=1,n
1815  ifa(i)=.false.
1816  100 CONTINUE
1817  RETURN
1818  END
1819 c-----------------------------------------------------------------------
1820  SUBROUTINE rzero3 (A,B,C,N)
1821  dimension a(1),b(1),c(1)
1822  DO 100 i=1,n
1823  a(i)=0.0
1824  b(i)=0.0
1825  c(i)=0.0
1826  100 CONTINUE
1827  RETURN
1828  END
1829 c-----------------------------------------------------------------------
1830  SUBROUTINE unitvec (X,Y,Z,N)
1831  dimension x(1),y(1),z(1)
1832  DO 100 i=1,n
1833  xlngth = sqrt( x(i)**2 + y(i)**2 + z(i)**2 )
1834  IF (xlngth.NE.0.0) THEN
1835  x(i) = x(i)/xlngth
1836  y(i) = y(i)/xlngth
1837  z(i) = z(i)/xlngth
1838  ENDIF
1839  100 CONTINUE
1840  RETURN
1841  END
1842 c-----------------------------------------------------------------------
1843  SUBROUTINE chkzvn (VMAX,IEL,IFC,IVNORL)
1844 C
1845  include 'SIZE'
1846  include 'GEOM'
1847  include 'SOLN'
1848  COMMON /scrmg/ v1(lx1,ly1,lz1,lelv)
1849  $ , v2(lx1,ly1,lz1,lelv)
1850  $ , v3(lx1,ly1,lz1,lelv)
1851  $ , vv(lx1,ly1,lz1,lelv)
1852 C
1853  nxz1 = lx1*lz1
1854  tolv = 0.01*vmax
1855 C
1856  vnor1 = facdot(v1(1,1,1,iel),unx(1,1,ifc,iel),ifc)
1857  vnor2 = facdot(v2(1,1,1,iel),uny(1,1,ifc,iel),ifc)
1858  vnor = vnor1 + vnor2
1859  IF (ldim.EQ.3) THEN
1860  vnor3 = facdot(v3(1,1,1,iel),unz(1,1,ifc,iel),ifc)
1861  vnor = vnor + vnor3
1862  ENDIF
1863  vnor = abs(vnor) / nxz1
1864 C
1865  ivnorl = 1
1866  IF (vnor .LT. tolv) ivnorl = 0
1867 C
1868  RETURN
1869  END
1870 c-----------------------------------------------------------------------
1871  SUBROUTINE bctwall (TMP1,TMP2,TMP3)
1872 C
1873 C Apply Dirichlet boundary conditions to surface of vector (V1,V2,V3)
1874 C (No antimask operation is applied).
1875 C
1876  include 'SIZE'
1877  include 'GEOM'
1878  include 'INPUT'
1879  include 'TSTEP'
1880 C
1881  dimension tmp1(lx1,ly1,lz1,1)
1882  $ , tmp2(lx1,ly1,lz1,1)
1883  $ , tmp3(lx1,ly1,lz1,1)
1884  common /nekcb/ cb
1885  CHARACTER CB*3
1886 C
1887  nface = 2*ldim
1888  ntot1 = lx1*ly1*lz1*nelv
1889 C
1890  CALL rzero (tmp1,ntot1)
1891  CALL rzero (tmp2,ntot1)
1892  IF (if3d) CALL rzero (tmp3,ntot1)
1893 C
1894  DO 2000 iel=1,nelv
1895  DO 2000 ifc=1,nface
1896  cb = cbc(ifc,iel,ifield)
1897  bc1 = bc(1,ifc,iel,ifield)
1898  bc2 = bc(2,ifc,iel,ifield)
1899  bc3 = bc(3,ifc,iel,ifield)
1900  IF (cb.EQ.'V ' .OR. cb.EQ.'VL ' .OR.
1901  $ cb.EQ.'WS ' .OR. cb.EQ.'WSL') THEN
1902  CALL facev (tmp1,iel,ifc,bc1,lx1,ly1,lz1)
1903  CALL facev (tmp2,iel,ifc,bc2,lx1,ly1,lz1)
1904  IF (ldim.EQ.3) CALL facev (tmp3,iel,ifc,bc3,lx1,ly1,lz1)
1905  IF (cb.EQ.'VL ' .OR. cb.EQ.'WSL')
1906  $ CALL globrot (tmp1(1,1,1,iel),tmp2(1,1,1,iel),
1907  $ tmp3(1,1,1,iel),iel,ifc)
1908  ENDIF
1909  IF (cb.EQ.'v ' .OR. cb.EQ.'vl ' .OR.
1910  $ cb.EQ.'ws ' .OR. cb.EQ.'wsl' .OR.
1911  $ cb.EQ.'mv ' .OR. cb.EQ.'mvn') THEN
1912  CALL faceiv (cb,tmp1(1,1,1,iel),tmp2(1,1,1,iel),
1913  $ tmp3(1,1,1,iel),iel,ifc,lx1,ly1,lz1)
1914  IF (cb.EQ.'vl ' .OR. cb.EQ.'wsl')
1915  $ CALL globrot (tmp1(1,1,1,iel),tmp2(1,1,1,iel),
1916  $ tmp3(1,1,1,iel),iel,ifc)
1917  ENDIF
1918  2000 CONTINUE
1919 C
1920  RETURN
1921  END
1922 c-----------------------------------------------------------------------
1923  SUBROUTINE antimsk1(X,XMASK,N)
1924 C------------------------------------------------------------------
1925 C
1926 C Return only Dirichlet boundary values of X
1927 C
1928 C-------------------------------------------------------------------
1929  REAL X(1),XMASK(1)
1930  include 'OPCTR'
1931 C
1932  DO 100 i=1,n
1933  x(i) = x(i)*(1.-xmask(i))
1934  100 CONTINUE
1935  RETURN
1936  END
1937 c-----------------------------------------------------------------------
1938  subroutine check_cyclic ! check for cyclic bcs
1939  include 'SIZE'
1940  include 'TOTAL'
1941 
1942  common /scrmg/ v1(lx1,ly1,lz1,lelt)
1943  $ , v2(lx1,ly1,lz1,lelt)
1944  $ , v3(lx1,ly1,lz1,lelt)
1945 
1946  integer e,f
1947 
1948  nface = 2*ldim
1949 
1950  n = lx1*ly1*lz1*nelt
1951  call rzero(v1,n)
1952  call rzero(v2,n)
1953  call rzero(v3,n)
1954 
1955  ifield = 1
1956  do e=1,nelt ! possibly U or B field
1957  do f=1,nface
1958 
1959  if (cbc(f,e,ifield).eq.'P '.or.cbc(f,e,ifield).eq.'p ') then
1960  call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f)
1961  k = 0
1962  do j2=js2,jf2,jskip2
1963  do j1=js1,jf1,jskip1
1964  k = k+1
1965  v1(j1,j2,1,e) = unx(j1,j2,1,e)
1966  v2(j1,j2,1,e) = uny(j1,j2,1,e)
1967  v3(j1,j2,1,e) = unz(j1,j2,1,e)
1968  enddo
1969  enddo
1970  endif
1971 
1972  enddo
1973  enddo
1974 
1975  ifcyclic = .false.
1976  call opdssum(v1,v2,v3)
1977 
1978  eps = 1.e-4
1979  if (ldim.eq.2) call rzero(v3,n)
1980 
1981  do e=1,nelt ! Check for turning angle
1982  do f=1,nface
1983 
1984  if (cbc(f,e,ifield).eq.'P '.or.cbc(f,e,ifield).eq.'p ') then
1985 
1986  call facindr(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,f) ! restricted indx
1987  snorm = 0.
1988  dnorm = 0.
1989  do k=k0,k1
1990  do j=j0,j1
1991  do i=i0,i1
1992  snorm = abs(v1(i,j,k,e))
1993  $ + abs(v2(i,j,k,e))
1994  $ + abs(v3(i,j,k,e))
1995  enddo
1996  enddo
1997  enddo
1998  if (snorm.gt.eps) ifcyclic = .true.
1999 
2000  endif
2001 
2002  enddo
2003  enddo
2004 
2005  itest = 0
2006  if (ifcyclic) itest = 1
2007  itest = iglmax(itest,1)
2008 
2009  if (itest.gt.0) ifcyclic = .true.
2010 
2011  return
2012  end
2013 c-----------------------------------------------------------------------
2014  real function glcflux(tx,ty,tz)
2015 c
2016  include 'SIZE'
2017  include 'TOTAL'
2018 
2019  real tx(lx1,ly1,lz1,lelv)
2020  real ty(lx1,ly1,lz1,lelv)
2021  real tz(lx1,ly1,lz1,lelv)
2022 
2023  character cb*3
2024 
2025  nxyz1= lx1*ly1*lz1
2026  ntot1= nxyz1*nelv
2027  nfaces = 2*ldim
2028 
2029  terma = 0.0
2030  termvl= 0.0
2031 
2032  do 100 iel=1,nelv
2033  do 100 iface=1,nfaces
2034  cb = cbc(iface,iel,1)
2035  if (cb.eq.'v ' .or. cb.eq.'V ' .or. cb.eq.'mv ') then
2036  call facind(kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,iface)
2037  ia = 0
2038  do 10 iz=kz1,kz2
2039  do 10 iy=ky1,ky2
2040  do 10 ix=kx1,kx2
2041  ia =ia + 1
2042  termxyz = tx(ix,iy,iz,iel)*unx(ia,1,iface,iel)
2043  $ + ty(ix,iy,iz,iel)*uny(ia,1,iface,iel)
2044  $ + tz(ix,iy,iz,iel)*unz(ia,1,iface,iel)
2045  terma = terma + area(ia,1,iface,iel)
2046  termvl = termvl+ termxyz * area(ia,1,iface,iel)
2047  10 continue
2048  endif
2049  100 continue
2050 
2051  glcflux = glsum(termvl,1) ! sum across processors
2052 
2053  return
2054  end
2055 c-----------------------------------------------------------------------
2056  subroutine local_bflux(flux,tx,ty,tz,ifld)
2057 c
2058  include 'SIZE'
2059  include 'TOTAL'
2060 
2061  real tx(lx1,ly1,lz1,1),
2062  $ ty(lx1,ly1,lz1,1),
2063  $ tz(lx1,ly1,lz1,1),
2064  $ flux(lx1,ly1,lz1,1)
2065 
2066  character cb*3
2067 
2068  nel = nelfld(ifld)
2069  nxyz = lx1*ly1*lz1
2070  ntot = nxyz*nel
2071  nfaces = 2*ldim
2072 
2073  call rzero(flux,ntot)
2074 
2075  do 100 iel=1,nel
2076  do 100 iface=1,nfaces
2077  cb = cbc(iface,iel,ifld)
2078  if (cb.ne.'E ') then
2079  call facind(kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,iface)
2080  ia = 0
2081  do 10 iz=kz1,kz2
2082  do 10 iy=ky1,ky2
2083  do 10 ix=kx1,kx2
2084  ia =ia + 1
2085  dtmp = tx(ix,iy,iz,iel)*unx(ia,1,iface,iel)
2086  $ + ty(ix,iy,iz,iel)*uny(ia,1,iface,iel)
2087  $ + tz(ix,iy,iz,iel)*unz(ia,1,iface,iel)
2088  flux(ix,iy,iz,iel) = flux(ix,iy,iz,iel)
2089  $ + dtmp*area(ia,1,iface,iel)
2090  10 continue
2091  endif
2092  100 continue
2093 
2094  return
2095  end
2096 c-----------------------------------------------------------------------
2097  SUBROUTINE facind2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
2098 C
2099  include 'SIZE'
2100  include 'TOPOL'
2101 C
2102  CALL dsset (lx1,ly1,lz1)
2103  iface = eface1(ifc)
2104  js1 = skpdat(1,iface)
2105  jf1 = skpdat(2,iface)
2106  jskip1 = skpdat(3,iface)
2107  js2 = skpdat(4,iface)
2108  jf2 = skpdat(5,iface)
2109  jskip2 = skpdat(6,iface)
2110 C
2111  RETURN
2112  END
2113 c-----------------------------------------------------------------------
2114  subroutine create_obj(iobjo,sid_list,n)
2115 c
2116 c defines an object for a given list of surface ids
2117 c
2118  include 'SIZE'
2119  include 'TOTAL'
2120 
2121  integer sid_list(n)
2122 
2123  integer e,f
2124 
2125  nobj = nobj + 1
2126  iobj = nobj
2127 
2128  if (maxobj.lt.nobj)
2129  $ call exitti('maxobj too small, increate in SIZE.$',ierr)
2130 
2131  do e=1,nelv
2132  do f=1,2*ndim
2133  do i=1,n
2134  if (boundaryid(f,e) .eq. sid_list(i)) then
2135  nmember(iobj) = nmember(iobj) + 1
2136  mem = nmember(iobj)
2137  ieg = lglel(e)
2138  object(iobj,mem,1) = ieg
2139  object(iobj,mem,2) = f
2140 c write(6,1) iobj,mem,f,ieg,e,nid,' OBJ'
2141  1 format(6i9,a4)
2142  endif
2143  enddo
2144  enddo
2145  enddo
2146 
2147  iobjo = iobj
2148 
2149  return
2150  end
2151 c-----------------------------------------------------------------------
2152  subroutine setbc(bid,ifld,cbci)
2153 c
2154 c sets boundary condition for a given surface id and field
2155 c
2156  include 'SIZE'
2157  include 'INPUT'
2158  include 'GEOM'
2159 
2160  character*3 cbci
2161  integer bid
2162 
2163  if (bid.lt.1 .or. bid.gt.lbid)
2164  $ call exitti('invalid boundary id!$',bid)
2165 
2166  cbc_bmap(bid,ifld) = cbci
2167 
2168  if (iftmsh(ifld)) then
2169  do iel = 1,nelt
2170  do ifc = 1,2*ndim
2171  if (boundaryidt(ifc,iel).eq.bid)
2172  $ cbc(ifc,iel,ifld) = cbc_bmap(bid,ifld)
2173  enddo
2174  enddo
2175  else
2176  do iel = 1,nelv
2177  do ifc = 1,2*ndim
2178  if (boundaryid(ifc,iel).eq.bid)
2179  $ cbc(ifc,iel,ifld) = cbc_bmap(bid,ifld)
2180  enddo
2181  enddo
2182  endif
2183 
2184  return
2185  end
unitvec
subroutine unitvec(X, Y, Z, N)
Definition: bdry.f:1831
create_obj
subroutine create_obj(iobjo, sid_list, n)
Definition: bdry.f:2115
bcdirsc
subroutine bcdirsc(S)
Definition: bdry.f:710
chknord
subroutine chknord(IFALGN, IFNORX, IFNORY, IFNORZ, IFC, IEL)
Definition: bdry.f:185
setbc
subroutine setbc(bid, ifld, cbci)
Definition: bdry.f:2153
check_cyclic
subroutine check_cyclic
Definition: bdry.f:1939
chkzvn
subroutine chkzvn(VMAX, IEL, IFC, IVNORL)
Definition: bdry.f:1844
globrot
subroutine globrot(R1, R2, R3, IEL, IFC)
Definition: bdry.f:1722
setrzer
subroutine setrzer
Definition: bdry.f:145
setdrs
subroutine setdrs(DRM1, DRTM1, DSM1, DSTM1, IFC)
Definition: bdry.f:1692
local_bflux
subroutine local_bflux(flux, tx, ty, tz, ifld)
Definition: bdry.f:2057
ctang2d
subroutine ctang2d(CANG, SANG, IXN, IYN, IAN, IFC, IEL)
Definition: bdry.f:1530
bcneusc
subroutine bcneusc(S, ITYPE)
Definition: bdry.f:786
trstax
subroutine trstax(TRX, TRY, SIGST, IEL, IFC)
Definition: bdry.f:1409
faceis
subroutine faceis(CB, S, IEL, IFACE, NX, NY, NZ)
Definition: bdry.f:919
setlog
subroutine setlog(ifecho)
Definition: bdry.f:3
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Definition: bdry.f:574
glcflux
real function glcflux(tx, ty, tz)
Definition: bdry.f:2015
trst2d
subroutine trst2d(TRX, TRY, SIGST, IEL, IFC)
Definition: bdry.f:1344
trcon
subroutine trcon(TRX, TRY, TRZ, TR1, TR2, TR3, IEL, IFC)
Definition: bdry.f:1304
nekasgn
subroutine nekasgn(ix, iy, iz, e)
Definition: bdry.f:1121
lfalse
subroutine lfalse(IFA, N)
Definition: bdry.f:1813
rzero3
subroutine rzero3(A, B, C, N)
Definition: bdry.f:1821
antimsk1
subroutine antimsk1(X, XMASK, N)
Definition: bdry.f:1924
facind2
subroutine facind2(JS1, JF1, JSKIP1, JS2, JF2, JSKIP2, IFC)
Definition: bdry.f:2098
chkcbc
subroutine chkcbc(CB, IEL, IFC, IFALGN, IERR)
Definition: bdry.f:275
bcneutr
subroutine bcneutr
Definition: bdry.f:1200
trst3d
subroutine trst3d(TRX, TRY, TRZ, SIGST, IEL, IFC)
Definition: bdry.f:1591
faceiv
subroutine faceiv(CB, V1, V2, V3, IEL, IFACE, NX, NY, NZ)
Definition: bdry.f:1002
facec2
subroutine facec2(A1, A2, B1, B2, IFC)
Definition: bdry.f:1782
chkaxcb
subroutine chkaxcb
Definition: bdry.f:252
bctwall
subroutine bctwall(TMP1, TMP2, TMP3)
Definition: bdry.f:1872
bcmask
subroutine bcmask
Definition: bdry.f:316
exitt
void exitt()
Definition: comm_mpi.f:604
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
facind
subroutine facind(kx1, kx2, ky1, ky2, kz1, kz2, nx, ny, nz, iface)
Definition: connect1.f:1028
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
facindr
subroutine facindr(kx1, kx2, ky1, ky2, kz1, kz2, nx, ny, nz, iface)
Definition: connect1.f:1045
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
dot
real function dot(V1, V2, N)
Definition: genxyz.f:885
set_ifbcor
subroutine set_ifbcor
Definition: induct.f:1092
iglsum
function iglsum(a, n)
Definition: math.f:926
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
copy
subroutine copy(a, b, n)
Definition: math.f:260
iglmax
function iglmax(a, n)
Definition: math.f:913
gllog
subroutine gllog(la, lb)
Definition: math.f:986
glsum
function glsum(x, n)
Definition: math.f:861
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rone
subroutine rone(a, n)
Definition: math.f:230
rzero
subroutine rzero(a, n)
Definition: math.f:208
chsign
subroutine chsign(a, n)
Definition: math.f:305
fix_surface_flux
subroutine fix_surface_flux
Definition: multimesh.f:530
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
opdsop
subroutine opdsop(a, b, c, op)
Definition: navier1.f:2603
facdot
function facdot(A, B, IFACE1)
Definition: subs1.f:834
facsub2
subroutine facsub2(A1, A2, A3, B1, B2, B3, IFACE1)
Definition: subs2.f:220
facexv
subroutine facexv(A1, A2, A3, B1, B2, B3, IFACE1, IOP)
Definition: subs2.f:170
stsmask
subroutine stsmask(C1MASK, C2MASK, C3MASK)
Definition: subs2.f:676
faccvs
subroutine faccvs(A1, A2, A3, B, IFACE1)
Definition: subs2.f:392
facexs
subroutine facexs(A, B, IFACE1, IOP)
Definition: subs2.f:125
rmask
subroutine rmask(R1, R2, R3, NEL)
Definition: subs2.f:1801