source: CIVL/examples/fortran/nek5000/core/mvmesh.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: 34.1 KB
Line 
1c-----------------------------------------------------------------------
2 subroutine cbcmesh
3C
4C Generate boundary conditions (CBC arrays) for mesh solver
5C
6 include 'SIZE'
7 include 'GEOM'
8 include 'INPUT'
9 include 'TSTEP'
10C
11 CHARACTER CBM*1,CBF*3,CBT*3,CB*3
12C
13 IFLD = 0
14 NFACE = 2*ldim
15 IFMELT = .FALSE.
16 IF (IFTMSH(IFLD)) IFMELT=.TRUE.
17C
18 DO 100 IEL=1,NELT
19 DO 100 IFC=1,NFACE
20 CBM = CBC(IFC,IEL,0)
21 CBF = CBC(IFC,IEL,1)
22 CBT = CBC(IFC,IEL,2)
23C
24 IF (CBT(1:1).EQ.'M') THEN
25 IFLD = 2
26 CB = CBT
27 GOTO 200
28 ENDIF
29 IF (CBF(1:1).EQ.'M' .OR. CBF(1:1).EQ.'m') THEN
30 IFLD = 1
31 CB = CBF
32C IF (CBF.EQ.'mv ' .AND. CBT.EQ.'E ') THEN
33C IFTMSH(0)=.TRUE.
34C ENDIF
35 IF (CBF.EQ.'mv ' .AND. CBM.EQ.'+' ) THEN
36 CB = 'mvn'
37 CBC(IFC,IEL,1) = CB
38 ENDIF
39 GOTO 200
40 ENDIF
41 IF (CBF(1:1).EQ.'V' .OR. CBF(1:1).EQ.'v' .OR.
42 $ CBF(1:1).EQ.'W' ) THEN
43 IFLD = 1
44 CB = 'FIX'
45 IF (IFMELT .OR. CBM.EQ.'+') CB='SYM'
46 GOTO 200
47 ENDIF
48 IF (CBT.EQ.'T ' .OR. CBT.EQ.'t ') THEN
49 IFLD = 2
50 CB = 'FIX'
51 IF (CBM.EQ.'+') CB='SYM'
52 GOTO 200
53 ENDIF
54 IF (CBF.EQ.'P ' .OR. CBF.EQ.'E ') THEN
55 IFLD = 1
56 CB = CBF
57 IF (CBM.EQ.'-') CB='FIX'
58 IF (CBM.EQ.'+') CB='SYM'
59 GOTO 200
60 ENDIF
61 IF (CBT.EQ.'P ' .OR. CBT.EQ.'E ') THEN
62 IFLD = 2
63 CB = CBT
64 IF (CBM.EQ.'-') CB='FIX'
65 IF (CBM.EQ.'+') CB='SYM'
66 GOTO 200
67 ENDIF
68 IFLD = 1
69 IF (CBF.EQ.' ') IFLD = 2
70 CB = 'SYM'
71 IF (CBM.EQ.'-') CB = 'FIX'
72C
73 200 CBC(IFC,IEL,0) = CB
74 DO 250 I=1,5
75 250 BC(I,IFC,IEL,0)=BC(I,IFC,IEL,IFLD)
76 100 CONTINUE
77C
78 return
79 end
80c-----------------------------------------------------------------------
81 subroutine admeshv
82C
83 include 'SIZE'
84 include 'SOLN'
85 include 'TSTEP'
86C
87 COMMON /SCRUZ/ FM1(LX1,LY1,LZ1,LELT)
88 $ , FM2(LX1,LY1,LZ1,LELT)
89 $ , FM3(LX1,LY1,LZ1,LELT)
90 $ , PHI(LX1,LY1,LZ1,LELT)
91C
92 NTOT1=lx1*ly1*lz1*NELV
93C
94 CALL RZERO (FM1,NTOT1)
95 CALL RZERO (FM2,NTOT1)
96 CALL RZERO (FM3,NTOT1)
97C
98 CALL DIVWS (FM1,VX,PHI,NELV,1)
99 CALL DIVWS (FM2,VY,PHI,NELV,2)
100 CALL ADD2 (BFX,FM1,NTOT1)
101 CALL ADD2 (BFY,FM2,NTOT1)
102 IF (ldim.EQ.3) THEN
103 CALL DIVWS (FM3,VZ,PHI,NELV,3)
104 CALL ADD2 (BFZ,FM3,NTOT1)
105 ENDIF
106C
107 return
108 end
109c-----------------------------------------------------------------------
110 subroutine admesht
111C
112 include 'SIZE'
113 include 'SOLN'
114 include 'TSTEP'
115C
116 COMMON /SCRUZ/ FMT(LX1,LY1,LZ1,LELT)
117 $ , PHI(LX1,LY1,LZ1,LELT)
118C
119 IFLD = 0
120 NEL = NELFLD(IFLD)
121 NTOT1= lx1*ly1*lz1*NEL
122C
123 CALL RZERO (FMT,NTOT1)
124 CALL DIVWS (FMT,T(1,1,1,1,IFIELD-1),PHI,NEL,1)
125 CALL ADDCOL3 (BQ(1,1,1,1,IFIELD-1),FMT,VTRANS(1,1,1,1,IFIELD),
126 & NTOT1)
127C
128 return
129 end
130c-----------------------------------------------------------------------
131 subroutine divws (fms,sfv,phi,nel,idir)
132C
133 include 'SIZE'
134 include 'GEOM'
135 include 'MASS'
136 include 'MVGEOM'
137 include 'WZ'
138 include 'INPUT'
139C
140 COMMON /SCRSF/ PHR(LX1,LY1,LZ1,LELT)
141 $ , PHS(LX1,LY1,LZ1,LELT)
142 $ , PHT(LX1,LY1,LZ1,LELT)
143
144 DIMENSION FMS(LX1,LY1,LZ1,1)
145 $ , SFV(LX1,LY1,LZ1,1)
146 $ , PHI(LX1,LY1,LZ1,1)
147C
148 NXYZ1 = lx1*ly1*lz1
149 NTOT1 = NXYZ1*NEL
150C
151 CALL COL3 (PHI,SFV,WX,NTOT1)
152 CALL URST (PHI,PHR,PHS,PHT,NEL)
153 CALL ADDCOL3 (FMS,RXM1,PHR,NTOT1)
154 CALL ADDCOL3 (FMS,SXM1,PHS,NTOT1)
155 IF (ldim.EQ.3) CALL ADDCOL3 (FMS,TXM1,PHT,NTOT1)
156C
157 CALL COL3 (PHI,SFV,WY,NTOT1)
158 CALL URST (PHI,PHR,PHS,PHT,NEL)
159 CALL ADDCOL3 (FMS,RYM1,PHR,NTOT1)
160 CALL ADDCOL3 (FMS,SYM1,PHS,NTOT1)
161 IF (ldim.EQ.3) CALL ADDCOL3 (FMS,TYM1,PHT,NTOT1)
162C
163 IF (ldim.EQ.3) THEN
164 CALL COL3 (PHI,SFV,WZ,NTOT1)
165 CALL URST (PHI,PHR,PHS,PHT,NEL)
166 CALL ADDCOL3 (FMS,RZM1,PHR,NTOT1)
167 CALL ADDCOL3 (FMS,SZM1,PHS,NTOT1)
168 CALL ADDCOL3 (FMS,TZM1,PHT,NTOT1)
169 ENDIF
170C
171 CALL COL2 (FMS,BM1,NTOT1)
172 CALL INVCOL2 (FMS,JACM1,NTOT1)
173C
174 IF (IFAXIS) CALL AXIFMS (FMS,SFV,PHI,NEL,IDIR)
175C
176 return
177 end
178c-----------------------------------------------------------------------
179 subroutine axifms (fms,sfv,phi,nel,idir)
180C
181 include 'SIZE'
182 include 'DXYZ'
183 include 'GEOM'
184 include 'INPUT'
185 include 'MASS'
186 include 'MVGEOM'
187 include 'WZ'
188 COMMON /SCRSF/ PHR(LX1,LY1,LZ1,LELT)
189 $ , PHS(LX1,LY1,LZ1,LELT)
190 $ , PHT(LX1,LY1,LZ1,LELT)
191C
192 DIMENSION FMS(LX1,LY1,LZ1,1)
193 $ , PHI(LX1,LY1,LZ1,1)
194 $ , SFV(LX1,LY1,LZ1,1)
195 $ , WYS(LX1)
196 EQUIVALENCE (WYS(1),PHT(1,1,1,1))
197C
198 NXYZ1 = lx1*ly1*lz1
199 NTOT1 = NXYZ1*NEL
200 CALL COL3 (PHI,SFV,WY,NTOT1)
201C
202 DO 100 IEL=1,NEL
203 IF ( IFRZER(IEL) ) THEN
204 IF (IDIR.EQ.1) THEN
205 CALL MXM (WY(1,1,1,IEL),lx1,DATM1,ly1,WYS,1)
206 DO 220 IX=1,lx1
207 FMS(IX,1,1,IEL)= FMS(IX,1,1,IEL) + WXM1(IX)*WAM1(1)*
208 $ WYS(IX)*SFV(IX,1,1,IEL)*JACM1(IX,1,1,IEL)
209 220 CONTINUE
210 ENDIF
211 DO 320 IX=1,lx1
212 DO 320 IY=2,ly1
213 FMS(IX,IY,1,IEL)=FMS(IX,IY,1,IEL) + PHI(IX,IY,1,IEL) *
214 $ BM1(IX,IY,1,IEL) / YM1(IX,IY,1,IEL)
215 320 CONTINUE
216 ELSE
217 CALL ADDCOL4 (FMS(1,1,1,IEL),PHI(1,1,1,IEL),JACM1(1,1,1,IEL),
218 $ W2CM1,NXYZ1)
219 ENDIF
220 100 CONTINUE
221C
222 return
223 end
224c-----------------------------------------------------------------------
225 subroutine updcoor
226C-----------------------------------------------------------------------
227C
228C Subroutine to update geometry for moving boundary problems
229C
230C-----------------------------------------------------------------------
231 include 'SIZE'
232 include 'TSTEP'
233 include 'INPUT'
234C
235 IFIELD = 0
236 NEL = NELFLD(IFIELD)
237C
238C Update collocation points coordinates
239C
240 CALL UPDXYZ (NEL)
241C
242C Shift lagged mesh velocity
243C
244 if (.not.ifrich) call lagmshv (nel)
245C
246 return
247 end
248c-----------------------------------------------------------------------
249 subroutine mvbdry (nel)
250
251c Evaluate mesh velocities at all moving boundaries
252
253 include 'SIZE'
254 include 'GEOM'
255 include 'INPUT'
256 include 'MVGEOM'
257 include 'SOLN'
258 include 'TSTEP'
259
260 common /scrsf/ wvx(lx1*ly1*lz1,lelt)
261 $ , wvy(lx1*ly1*lz1,lelt)
262 $ , wvz(lx1*ly1*lz1,lelt)
263 common /scrch/ wtx(lx1*ly1*lz1,lelt)
264 $ , wty(lx1*ly1*lz1,lelt)
265 common /scrmg/ wtz(lx1*ly1*lz1,lelt)
266 $ , rnx(lx1*ly1*lz1,lelt)
267 $ , rny(lx1*ly1*lz1,lelt)
268 $ , rnz(lx1*ly1*lz1,lelt)
269 common /scruz/ dsa(lx1*ly1*lz1,lelt)
270 $ , qni(lx1*ly1*lz1,lelt)
271 $ , smt(lx1*ly1*lz1,lelt)
272 $ , ta (lx1*ly1*lz1,lelt)
273
274 logical ifalgn,ifnorx,ifnory,ifnorz,ifdsmv,ifregw
275 character cb*3
276 integer e,f
277
278 ifield = 0
279 nxyz1 = lx1*ly1*lz1
280 n = lx1*ly1*lz1*nel
281 nface = 2*ldim
282 call rzero3 (rnx,rny,rnz,n)
283
284 do 100 e=1,nel
285 do 100 f=1,nface
286 cb = cbc(f,e,ifield)
287 if (cb.eq.'ms ' .or. cb.eq.'MS ' .or.
288 $ cb.eq.'msi' .or. cb.eq.'MSI' .or.
289 $ cb.eq.'mm ' .or. cb.eq.'MM ' .or.
290 $ cb.eq.'mv ' .OR. cb.eq.'mvn' .or.
291 $ cb.eq.'MLI') then
292 call facexv (unx(1,1,f,e),uny(1,1,f,e),
293 $ unz(1,1,f,e),rnx(1,e),
294 $ rny(1,e),rnz(1,e),f,1)
295 endif
296 100 continue
297
298 call opdssum (rnx,rny,rnz)
299 call unitvec (rnx,rny,rnz,n)
300
301 call rzero3 (wvx,wvy,wvz,n)
302 call rzero3 (wtx,wty,wtz,n)
303 do 1000 isweep=1,2
304
305 ifregw = .false.
306 ifdsmv = .false.
307 call rzero (dsa,n)
308 call rzero (ta,n)
309
310 if (ifflow) then
311 ifield = 1
312 call rzero (smt,n)
313 do 210 e=1,nelv
314 do 210 f=1,nface
315 cb = cbc(f,e,ifield)
316 if (cb.eq.'mv ' .or. cb.eq.'mvn' .or.
317 $ cb.eq.'mm ' .or. cb.eq.'MM ' .or.
318 $ cb.eq.'msi' .or. cb.eq.'MSI' .or.
319 $ cb.eq.'ms ' .or. cb.eq.'MS ') then
320 ifregw = .true.
321 call facec3 (wvx(1,e),wvy(1,e),wvz(1,e),
322 $ vx (1,1,1,e),vy (1,1,1,e),vz (1,1,1,e),f)
323 if (cb.ne.'mv ')
324 $ call norcmp2(wvx(1,e),wvy(1,e),wvz(1,e),e,f)
325 endif
326c if (cb.eq.'msi' .or. cb.eq.'MSI') then
327c ifdsmv = .true.
328c call facsmt (smt(1,e),f)
329c endif
330 210 continue
331
332 call dsavg(wvx)
333 call dsavg(wvy)
334 if (if3d) call dsavg(wvz)
335
336 if (istep.eq.0) call opcopy(wx,wy,wz,wvx,wvy,wvz)
337c if (istep.eq.4) then
338c call opcopy(wx,wy,wz,wvx,wvy,wvz)
339c call outpost(vx,vy,vz,wtx,wtx,' ')
340c call outpost(wx,wy,wz,wtx,wtx,' ')
341c write(6,*) 'quit1',istep
342c stop
343c endif
344
345 iregw = 0 ! Global handshake on logicals ifregw, ifdsmv
346 if (ifregw) iregw = 1 ! pff 9/11/07
347 iregw = iglmax(iregw,1)
348 if (iregw.eq.1) ifregw = .true.
349
350 idsmv = 0
351 if (ifdsmv) idsmv = 1
352 idsmv = iglmax(idsmv,1)
353 if (idsmv.eq.1) ifdsmv = .true.
354
355 ifdsmv = .false.
356 ifregw = .true.
357
358c if (ifdsmv) then
359c call dssum (smt,lx1,ly1,lz1)
360c do 215 e=1,nelv
361c do 215 f=1,nface
362c cb = cbc(f,e,ifield)
363c if (cb.eq.'msi' .or. cb.eq.'MSI') then
364c call facec3 (wtx(1,e),wty(1,e),wtz(1,e),
365c $ vx(1,1,1,e),vy(1,1,1,e),vz(1,1,1,e),f)
366c call facemv (wtx(1,e),wty(1,e),wtz(1,e),
367c $ rnx(1,e),rny(1,e),rnz(1,e),
368c $ smt(1,e),f)
369c endif
370c 215 continue
371c call opdssum(wtx,wty,wtz)
372c endif
373
374 endif
375C
376 if (ifmelt .and. istep.gt.0) then
377 ifield = 2
378 call rzero (smt,n)
379 call cqnet (qni,ta,nel)
380 do 220 e=1,nelt
381 do 220 f=1,nface
382 cb = cbc(f,e,ifield)
383 if (cb.eq.'MLI') call facsmt (smt(1,e),f)
384 if (cb.eq.'MLI' .or. cb.eq.'MCI') then
385 call facexs (area(1,1,f,e),ta,f,1)
386 call add2 (dsa(1,e),ta,nxyz1)
387 endif
388 220 continue
389 call dssum (smt,lx1,ly1,lz1)
390 call dssum (dsa,lx1,ly1,lz1)
391 do 280 e=1,nelt
392 do 280 f=1,nface
393 cb = cbc(f,e,ifield)
394 if (cb.eq.'MLI') then
395 rhola = -0.5 * bc(5,f,e,ifield)
396 call facemt (wtx(1,e),wty(1,e),wtz(1,e),
397 $ rnx(1,e),rny(1,e),rnz(1,e),
398 $ qni(1,e),dsa(1,e),smt(1,e),
399 $ rhola,f)
400 endif
401 280 continue
402 call opdssum (wtx,wty,wtz)
403 endif
404
405 ifield = 0
406 do 330 e=1,nel
407 do 330 f=1,nface
408 cb = cbc(f,e,ifield)
409 if (cb.eq.'SYM') then
410 call chknord (ifalgn,ifnorx,ifnory,ifnorz,f,e)
411 if (ifregw) then
412 if (ifnorx) call facev (wvx,e,f,0.0,lx1,ly1,lz1)
413 if (ifnory) call facev (wvy,e,f,0.0,lx1,ly1,lz1)
414 if (ifnorz) call facev (wvz,e,f,0.0,lx1,ly1,lz1)
415 if (.not.ifalgn) call faczqn (wvx(1,e),wvy(1,e),
416 $ wvz(1,e),f,e)
417 endif
418 if (ifdsmv .or. ifmelt) then
419 if (ifnorx) call facev (wtx,e,f,0.0,lx1,ly1,lz1)
420 if (ifnory) call facev (wty,e,f,0.0,lx1,ly1,lz1)
421 if (ifnorz) call facev (wtz,e,f,0.0,lx1,ly1,lz1)
422 if (.not.ifalgn) call faczqn (wtx(1,e),wty(1,e),
423 $ wtz(1,e),f,e)
424 endif
425 endif
426 330 continue
427
428 do 350 e=1,nel
429 do 350 f=1,nface
430 cb = cbc(f,e,ifield)
431 if (cb.eq.'FIX') then
432 if (ifregw) then
433 call facev (wvx,e,f,0.0,lx1,ly1,lz1)
434 call facev (wvy,e,f,0.0,lx1,ly1,lz1)
435 if (ldim.eq.3) call facev (wvz,e,f,0.0,lx1,ly1,lz1)
436 endif
437 if (ifdsmv .or. ifmelt) then
438 call facev (wtx,e,f,0.0,lx1,ly1,lz1)
439 call facev (wty,e,f,0.0,lx1,ly1,lz1)
440 if (ldim.eq.3) call facev (wtz,e,f,0.0,lx1,ly1,lz1)
441 endif
442 endif
443 350 continue
444
445 if (isweep.eq.1) then
446 if (ifregw) then
447 call dsop (wvx,'MXA',lx1,ly1,lz1)
448 call dsop (wvy,'MXA',lx1,ly1,lz1)
449 if (ldim.eq.3) call dsop (wvz,'MXA',lx1,ly1,lz1)
450 endif
451 if (ifdsmv .or. ifmelt) then
452 call dsop (wtx,'MXA',lx1,ly1,lz1)
453 call dsop (wty,'MXA',lx1,ly1,lz1)
454 if (ldim.eq.3) call dsop (wtz,'MXA',lx1,ly1,lz1)
455 endif
456 else
457 if (ifregw) then
458 call dsop (wvx,'MNA',lx1,ly1,lz1)
459 call dsop (wvy,'MNA',lx1,ly1,lz1)
460 if (ldim.eq.3) call dsop (wvz,'MNA',lx1,ly1,lz1)
461 endif
462 if (ifdsmv .or. ifmelt) then
463 call dsop (wtx,'MNA',lx1,ly1,lz1)
464 call dsop (wty,'MNA',lx1,ly1,lz1)
465 if (ldim.eq.3) call dsop (wtz,'MNA',lx1,ly1,lz1)
466 endif
467 endif
468
469 1000 continue
470
471 call rmask (wx,wy,wz,nel)
472c if (istep.eq.2) then
473c call outpost(wx,wy,wz,wtx,wtx,' ')
474c write(6,*) 'quit2'
475c stop
476c endif
477
478 if (ifregw) then
479 call add2 (wx,wvx,n)
480 call add2 (wy,wvy,n)
481 if (ldim.eq.3) call add2 (wz,wvz,n)
482 endif
483 if (ifdsmv .or. ifmelt) then
484 call add2 (wx,wtx,n)
485 call add2 (wy,wty,n)
486 if (ldim.eq.3) call add2 (wz,wtz,n)
487 endif
488
489c if (istep.gt.4) then
490c call outpost(wx,wy,wz,wtx,wtx,' ')
491c write(6,*) 'quit3'
492c stop
493c endif
494
495 return
496 end
497c-----------------------------------------------------------------------
498 subroutine norcmp2(wvx,wvy,wvz,e,f)
499 include 'SIZE'
500 include 'GEOM'
501 include 'INPUT'
502
503
504 real wvx(lx1,ly1,lz1),wvy(lx1,ly1,lz1),wvz(lx1,ly1,lz1)
505
506 integer e,f
507
508 common /scruz/ r1(lx1,ly1,lz1),r2(lx1,ly1,lz1),r3(lx1,ly1,lz1)
509
510 call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,f)
511
512 l=0
513 do k=k0,k1
514 do j=j0,j1
515 do i=i0,i1
516 l=l+1
517 scale=wvx(i,j,k)*unx(l,1,f,e)
518 $ +wvy(i,j,k)*uny(l,1,f,e)
519 $ +wvz(i,j,k)*unz(l,1,f,e)
520 wvx(i,j,k) = scale*unx(l,1,f,e)
521 wvy(i,j,k) = scale*uny(l,1,f,e)
522 wvz(i,j,k) = scale*unz(l,1,f,e)
523 enddo
524 enddo
525 enddo
526
527c wvxm = vlamax(wvx,lx1*ly1)
528c write(6,*) f,e,wvxm,' w-max'
529
530 return
531 end
532c-----------------------------------------------------------------------
533 subroutine norcmp (wt1,wt2,wt3,rnx,rny,rnz,ifc)
534C
535 include 'SIZE'
536 COMMON /SCRUZ/ R1(LX1,LY1,LZ1),R2(LX1,LY1,LZ1),R3(LX1,LY1,LZ1)
537C
538 DIMENSION WT1(LX1,LY1,LZ1),WT2(LX1,LY1,LZ1),WT3(LX1,LY1,LZ1)
539 $ , RNX(LX1,LY1,LZ1),RNY(LX1,LY1,LZ1),RNZ(LX1,LY1,LZ1)
540C
541 NXYZ1 = lx1*ly1*lz1
542C
543 CALL COPY (R1,WT1,NXYZ1)
544 CALL COPY (R2,WT2,NXYZ1)
545 IF (ldim.EQ.3) CALL COPY (R3,WT3,NXYZ1)
546 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
547C
548 IF (ldim.EQ.2) THEN
549 DO 200 J2=JS2,JF2,JSKIP2
550 DO 200 J1=JS1,JF1,JSKIP1
551 WN = R1(J1,J2,1)*RNX(J1,J2,1) +
552 $ R2(J1,J2,1)*RNY(J1,J2,1)
553 WT1(J1,J2,1) = WN *RNX(J1,J2,1)
554 WT2(J1,J2,1) = WN *RNY(J1,J2,1)
555 200 CONTINUE
556 ELSE
557 DO 300 J2=JS2,JF2,JSKIP2
558 DO 300 J1=JS1,JF1,JSKIP1
559 WN = R1(J1,J2,1)*RNX(J1,J2,1) +
560 $ R2(J1,J2,1)*RNY(J1,J2,1) +
561 $ R3(J1,J2,1)*RNZ(J1,J2,1)
562 WT1(J1,J2,1) = WN *RNX(J1,J2,1)
563 WT2(J1,J2,1) = WN *RNY(J1,J2,1)
564 WT3(J1,J2,1) = WN *RNZ(J1,J2,1)
565 300 CONTINUE
566 ENDIF
567C
568 return
569 end
570c-----------------------------------------------------------------------
571 subroutine facemv (wt1,wt2,wt3,rnx,rny,rnz,smt,ifc)
572C
573 include 'SIZE'
574 COMMON /SCRUZ/ R1(LX1,LY1,LZ1),R2(LX1,LY1,LZ1),R3(LX1,LY1,LZ1)
575C
576 DIMENSION WT1(LX1,LY1,LZ1),WT2(LX1,LY1,LZ1),WT3(LX1,LY1,LZ1)
577 $ , RNX(LX1,LY1,LZ1),RNY(LX1,LY1,LZ1),RNZ(LX1,LY1,LZ1)
578 $ , SMT(LX1,LY1,LZ1)
579C
580 NXYZ1 = lx1*ly1*lz1
581C
582 CALL COPY (R1,WT1,NXYZ1)
583 CALL COPY (R2,WT2,NXYZ1)
584 IF (ldim.EQ.3) CALL COPY (R3,WT3,NXYZ1)
585 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
586C
587 IF (ldim.EQ.2) THEN
588 DO 200 J2=JS2,JF2,JSKIP2
589 DO 200 J1=JS1,JF1,JSKIP1
590 WN = ( R1(J1,J2,1)*RNX(J1,J2,1) +
591 $ R2(J1,J2,1)*RNY(J1,J2,1) ) / SMT(J1,J2,1)
592 WT1(J1,J2,1) = WN *RNX(J1,J2,1)
593 WT2(J1,J2,1) = WN *RNY(J1,J2,1)
594 200 CONTINUE
595 ELSE
596 DO 300 J2=JS2,JF2,JSKIP2
597 DO 300 J1=JS1,JF1,JSKIP1
598 WN = ( R1(J1,J2,1)*RNX(J1,J2,1) +
599 $ R2(J1,J2,1)*RNY(J1,J2,1) +
600 $ R3(J1,J2,1)*RNZ(J1,J2,1) ) / SMT(J1,J2,1)
601 WT1(J1,J2,1) = WN *RNX(J1,J2,1)
602 WT2(J1,J2,1) = WN *RNY(J1,J2,1)
603 WT3(J1,J2,1) = WN *RNZ(J1,J2,1)
604 300 CONTINUE
605 ENDIF
606C
607 return
608 end
609c-----------------------------------------------------------------------
610 subroutine faczqn (wt1,wt2,wt3,ifc,iel)
611C
612 include 'SIZE'
613 include 'GEOM'
614 include 'TOPOL'
615 COMMON /SCRUZ/ R1(LX1,LY1,LZ1),R2(LX1,LY1,LZ1),R3(LX1,LY1,LZ1)
616C
617 DIMENSION WT1(LX1,LY1,LZ1),WT2(LX1,LY1,LZ1),WT3(LX1,LY1,LZ1)
618C
619 NXYZ1 = lx1*ly1*lz1
620 CALL COPY (R1,WT1,NXYZ1)
621 CALL COPY (R2,WT2,NXYZ1)
622 IF (ldim.EQ.3) CALL COPY (R3,WT3,NXYZ1)
623C
624 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
625 I = 0
626C
627 IF (ldim.EQ.2) THEN
628 DO 200 J2=JS2,JF2,JSKIP2
629 DO 200 J1=JS1,JF1,JSKIP1
630 I = I+1
631 W1 = R1(J1,J2,1)*T1X(I,1,IFC,IEL) +
632 $ R2(J1,J2,1)*T1Y(I,1,IFC,IEL)
633 WT1(J1,J2,1) = W1 *T1X(I,1,IFC,IEL)
634 WT2(J1,J2,1) = W1 *T1Y(I,1,IFC,IEL)
635 200 CONTINUE
636 ELSE
637 DO 300 J2=JS2,JF2,JSKIP2
638 DO 300 J1=JS1,JF1,JSKIP1
639 I = I+1
640 W1 = R1(J1,J2,1)*T1X(I,1,IFC,IEL) +
641 $ R2(J1,J2,1)*T1Y(I,1,IFC,IEL) +
642 $ R3(J1,J2,1)*T1Z(I,1,IFC,IEL)
643 WT1(J1,J2,1) = W1 *T1X(I,1,IFC,IEL)
644 WT2(J1,J2,1) = W1 *T1Y(I,1,IFC,IEL)
645 WT3(J1,J2,1) = W1 *T1Z(I,1,IFC,IEL)
646 300 CONTINUE
647 ENDIF
648C
649 return
650 end
651c-----------------------------------------------------------------------
652 subroutine facsmt (smt,ifc)
653C
654 include 'SIZE'
655 DIMENSION SMT(LX1,LY1,LZ1)
656C
657 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
658C
659 DO 100 J2=JS2,JF2,JSKIP2
660 DO 100 J1=JS1,JF1,JSKIP1
661 SMT(J1,J2,1)=SMT(J1,J2,1) + 1.0
662 100 CONTINUE
663C
664 return
665 end
666c-----------------------------------------------------------------------
667 subroutine cqnet (qni,ta,nel)
668C
669 include 'SIZE'
670 include 'INPUT'
671 include 'SOLN'
672 include 'TSTEP'
673 COMMON /SCRVH/ H1(LX1,LY1,LZ1,LELT)
674 $ , H2(LX1,LY1,LZ1,LELT)
675C
676 DIMENSION QNI(LX1,LY1,LZ1,1)
677 $ , TA (LX1,LY1,LZ1,1)
678C
679 INTLOC = -1
680 IMSHL = 2
681 NTOT1 = lx1*ly1*lz1*NEL
682C
683 CALL SETHLM (H1,H2,INTLOC)
684 CALL AXHELM (TA,T(1,1,1,1,IFIELD-1),H1,H2,IMSHL,1)
685 CALL SUB3 (QNI,TA,BQ(1,1,1,1,IFIELD-1),NTOT1)
686 CALL DSSUM (QNI,lx1,ly1,lz1)
687C
688 return
689 end
690c-----------------------------------------------------------------------
691 subroutine facemt (w1,w2,w3,rnx,rny,rnz,qni,dsa,smt,rhola,ifc)
692C
693 include 'SIZE'
694 include 'GEOM'
695C
696 DIMENSION W1 (LX1,LY1,LZ1)
697 $ , W2 (LX1,LY1,LZ1)
698 $ , W3 (LX1,LY1,LZ1)
699 $ , RNX(LX1,LY1,LZ1)
700 $ , RNY(LX1,LY1,LZ1)
701 $ , RNZ(LX1,LY1,LZ1)
702 $ , QNI(LX1,LY1,LZ1)
703 $ , DSA(LX1,LY1,LZ1)
704 $ , SMT(LX1,LY1,LZ1)
705C
706 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
707C
708 IF (ldim.EQ.2) THEN
709 DO 200 J2=JS2,JF2,JSKIP2
710 DO 200 J1=JS1,JF1,JSKIP1
711 AA = QNI(J1,J2,1) / ( DSA(J1,J2,1)*SMT(J1,J2,1)*RHOLA )
712 W1(J1,J2,1) = RNX(J1,J2,1) * AA
713 W2(J1,J2,1) = RNY(J1,J2,1) * AA
714 200 CONTINUE
715 ELSE
716 DO 300 J2=JS2,JF2,JSKIP2
717 DO 300 J1=JS1,JF1,JSKIP1
718 AA = QNI(J1,J2,1) / ( DSA(J1,J2,1)*SMT(J1,J2,1)*RHOLA )
719 W1(J1,J2,1) = RNX(J1,J2,1) * AA
720 W2(J1,J2,1) = RNY(J1,J2,1) * AA
721 W3(J1,J2,1) = RNZ(J1,J2,1) * AA
722 300 CONTINUE
723 ENDIF
724C
725 return
726 end
727c-----------------------------------------------------------------------
728 subroutine elasolv (nel)
729C
730C Elastostatic solver for mesh deformation
731C
732 include 'SIZE'
733 include 'GEOM'
734 include 'INPUT'
735 include 'MVGEOM'
736 include 'SOLN'
737 include 'TSTEP'
738C
739 COMMON /SCRNS/ DW1 (LX1,LY1,LZ1,LELT)
740 $ , DW2 (LX1,LY1,LZ1,LELT)
741 $ , DW3 (LX1,LY1,LZ1,LELT)
742 $ , AW1 (LX1,LY1,LZ1,LELT)
743 $ , AW2 (LX1,LY1,LZ1,LELT)
744 $ , AW3 (LX1,LY1,LZ1,LELT)
745 COMMON /SCRVH/ H1 (LX1,LY1,LZ1,LELT)
746 $ , H2 (LX1,LY1,LZ1,LELT)
747 common /scruz/ prt (lx1,ly1,lz1,lelt)
748 COMMON /FASTMD/ IFDFRM(LELT), IFFAST(LELT), IFH2, IFSOLV
749 LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
750
751 if (ifusermv) return ! Compute wx,wy,wz in userchk.
752c
753 NTOT1 = lx1*ly1*lz1*NEL
754 MAXIT = 1000
755 MATMOD = -1
756 IFH2 = .FALSE.
757 IFSOLV = .TRUE.
758 IMSOLV = 0
759 VNU = 0.0
760 VNU = param(47)
761 if (vnu.eq.0) VNU = 0.4
762 vnu = max(0.00,vnu)
763 vnu = min(0.499,vnu)
764
765C Set up elastic material constants
766
767 CE = 1./(1. + VNU)
768 C2 = VNU * CE / (1. - 2.*VNU)
769 C3 = 0.5 * CE
770 CALL CFILL (H1,C2,NTOT1)
771 CALL CFILL (H2,C3,NTOT1)
772
773C Solve for interior mesh velocities
774
775 CALL MESHTOL (AW1,TOLMSH,NEL,IMSOLV)
776 IF (IMSOLV.EQ.1) return
777
778 CALL AXHMSF (AW1,AW2,AW3,WX,WY,WZ,H1,H2,MATMOD)
779
780c if (istep.eq.2) then
781c call outpost(wx,wy,wz,h1,h2,' ')
782c call outpost(aw1,aw2,aw3,h1,h2,' ')
783c write(6,*) 'quit elas 1'
784c stop
785c endif
786
787 CALL CHSIGN (AW1,NTOT1)
788 CALL CHSIGN (AW2,NTOT1)
789 IF (ldim.EQ.3) CALL CHSIGN (AW3,NTOT1)
790 CALL HMHZSF ('NOMG',DW1,DW2,DW3,AW1,AW2,AW3,H1,H2,
791 $ W1MASK,W2MASK,W3MASK,WMULT,TOLMSH,
792 $ MAXIT,MATMOD)
793C
794C Update mesh velocities
795C
796 CALL ADD2 (WX,DW1,NTOT1)
797 CALL ADD2 (WY,DW2,NTOT1)
798 IF (ldim.EQ.3) CALL ADD2 (WZ,DW3,NTOT1)
799
800 ifldt = ifield
801 ifield=1
802 if (ifheat) ifield=2
803 call dsavg(wx)
804 call dsavg(wy)
805 call dsavg(wz)
806 ifield = ifldt
807
808c if (istep.gt.1) then
809c ifldx = ifield
810c ifield = 1
811c call incomprn (wx,wy,wz,prt) ! project U onto div-free space
812c ifield = ifldx
813c return
814c endif
815
816 return
817 end
818c-----------------------------------------------------------------------
819 subroutine meshtol (ta,tolmsh,nel,imsolv)
820C
821 include 'SIZE'
822 include 'EIGEN'
823 include 'MVGEOM'
824 include 'TSTEP'
825 DIMENSION TA(LX1,LY1,LZ1,1)
826C
827 NTOT1 = lx1*ly1*lz1*NEL
828 TOLAB = TOLREL
829C
830 DELTA = 1.0E-9
831 X = 1.0 + DELTA
832 Y = 1.0
833 DIFF = ABS(X - Y)
834 IF (DIFF .EQ. 0.0) EPS = 1.0E-05
835 IF (DIFF .GT. 0.0) EPS = 1.0E-12
836C
837 CALL OPDOT (TA,WX,WY,WZ,WX,WY,WZ,NTOT1)
838C
839 WDOT = GLMAX(TA,NTOT1)
840 WMAX = SQRT(WDOT)
841 IF (WMAX .LT. EPS) THEN
842 IMSOLV = 1
843 return
844 ENDIF
845C
846 TOLMSH = TOLAB * WMAX * SQRT(EIGAA)
847C
848 return
849 end
850c-----------------------------------------------------------------------
851 subroutine updxyz (nel)
852C
853 include 'SIZE'
854 include 'TSTEP'
855 include 'MVGEOM'
856 include 'GEOM'
857 include 'INPUT'
858 COMMON /SCRSF/ UX(LX1,LY1,LZ1,LELT)
859 $ , UY(LX1,LY1,LZ1,LELT)
860 $ , UZ(LX1,LY1,LZ1,LELT)
861 DIMENSION ABM(3)
862C
863 NTOT1 = lx1*ly1*lz1*NEL
864C
865 DO 10 I=1,NBD
866 10 ABM(I) = DT*ABMSH(I)
867C
868 IF (ISTEP.EQ.0) THEN
869 CALL COPY (UX,WX,NTOT1)
870 CALL COPY (UY,WY,NTOT1)
871 IF (ldim.EQ.3) CALL COPY (UZ,WZ,NTOT1)
872 else
873 if (ifrich) then
874 call cmult2(ux,wx,dt,ntot1)
875 call cmult2(uy,wy,dt,ntot1)
876 if (ldim.eq.3) call cmult2(uz,wz,dt,ntot1)
877 else
878 CALL CMULT2 (UX,WX,ABM(1),NTOT1)
879 CALL CMULT2 (UY,WY,ABM(1),NTOT1)
880 IF (ldim.EQ.3) CALL CMULT2 (UZ,WZ,ABM(1),NTOT1)
881 DO 100 ILAG=2,NBD
882 CALL ADD2S2 (UX,WXLAG(1,1,1,1,ILAG-1),ABM(ILAG),NTOT1)
883 CALL ADD2S2 (UY,WYLAG(1,1,1,1,ILAG-1),ABM(ILAG),NTOT1)
884 IF (ldim.EQ.3)
885 $ CALL ADD2S2 (UZ,WZLAG(1,1,1,1,ILAG-1),ABM(ILAG),NTOT1)
886 100 CONTINUE
887 endif
888 ENDIF
889C
890 CALL ADD2 (XM1,UX,NTOT1)
891 CALL ADD2 (YM1,UY,NTOT1)
892 IF (ldim.EQ.3) CALL ADD2 (ZM1,UZ,NTOT1)
893C
894 return
895 end
896c-----------------------------------------------------------------------
897 subroutine lagmshv (nel)
898C-----------------------------------------------------------------------
899C
900C Keep old mesh velocity
901C
902C-----------------------------------------------------------------------
903 include 'SIZE'
904 include 'INPUT'
905 include 'MVGEOM'
906 include 'TSTEP'
907C
908 NTOT1 = lx1*ly1*lz1*NEL
909C
910 DO 100 ILAG=NBDINP-1,2,-1
911 CALL COPY (WXLAG(1,1,1,1,ILAG),WXLAG(1,1,1,1,ILAG-1),NTOT1)
912 CALL COPY (WYLAG(1,1,1,1,ILAG),WYLAG(1,1,1,1,ILAG-1),NTOT1)
913 IF (ldim.EQ.3)
914 $ CALL COPY (WZLAG(1,1,1,1,ILAG),WZLAG(1,1,1,1,ILAG-1),NTOT1)
915 100 CONTINUE
916C
917 CALL COPY (WXLAG(1,1,1,1,1),WX,NTOT1)
918 CALL COPY (WYLAG(1,1,1,1,1),WY,NTOT1)
919 IF (ldim.EQ.3)
920 $ CALL COPY (WZLAG(1,1,1,1,1),WZ,NTOT1)
921C
922 return
923 end
924c-----------------------------------------------------------------------
925 subroutine facec3 (a1,a2,a3,b1,b2,b3,ifc)
926C
927C Copy the face (IFC) of B1,B2,B3 to A1,A2,A3.
928C IFACE is the input in the pre-processor ordering scheme.
929C
930 include 'SIZE'
931 DIMENSION A1(LX1,LY1,LZ1)
932 $ , A2(LX1,LY1,LZ1)
933 $ , A3(LX1,LY1,LZ1)
934 $ , B1(LX1,LY1,LZ1)
935 $ , B2(LX1,LY1,LZ1)
936 $ , B3(LX1,LY1,LZ1)
937C
938 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
939C
940 DO 100 J2=JS2,JF2,JSKIP2
941 DO 100 J1=JS1,JF1,JSKIP1
942 A1(J1,J2,1)=B1(J1,J2,1)
943 A2(J1,J2,1)=B2(J1,J2,1)
944 A3(J1,J2,1)=B3(J1,J2,1)
945 100 CONTINUE
946 return
947 end
948c-----------------------------------------------------------------------
949 subroutine ptbgeom
950C-----------------------------------------------------------------------
951C
952C Subroutine to impose perturbation to geometry before solution
953C for moving boundary problems
954C
955C-----------------------------------------------------------------------
956 include 'SIZE'
957 include 'GEOM'
958 include 'MVGEOM'
959 include 'SOLN'
960 include 'TSTEP'
961 include 'INPUT'
962 COMMON /SCRUZ/ XM3 (LX3,LY3,LZ3,LELT)
963 $ , YM3 (LX3,LY3,LZ3,LELT)
964 $ , ZM3 (LX3,LY3,LZ3,LELT)
965C
966 IF (ISTEP .EQ. 0) return
967 IFIELD = 0
968 NEL = NELFLD(IFIELD)
969 NTOT1 = lx1*ly1*lz1*NEL
970 IMESH = 1
971 IF ( IFTMSH(IFIELD) ) IMESH = 2
972C
973 CALL IBDGEOM (NEL)
974 CALL ELASOLV (NEL)
975 CALL UPDXYZ (NEL)
976 CALL GEOM1 (XM3,YM3,ZM3)
977 CALL GEOM2
978 CALL UPDMSYS (0)
979 CALL VOLUME
980 CALL SETINVM
981 CALL LAGMASS
982C
983 CALL RZERO (WX,NTOT1)
984 CALL RZERO (WY,NTOT1)
985 IF (ldim.EQ.3) CALL RZERO (WZ,NTOT1)
986C
987 return
988 end
989c-----------------------------------------------------------------------
990 subroutine ibdgeom (nel)
991C
992C Routine to evaluate mesh velocities at all moving boundaries
993C
994 include 'SIZE'
995 include 'GEOM'
996 include 'INPUT'
997 include 'PARALLEL'
998 include 'MVGEOM'
999 include 'TSTEP'
1000C
1001 CHARACTER CB*1
1002C
1003 NFACE = 2*ldim
1004 NTOT1 = lx1*ly1*lz1*NEL
1005C
1006 CALL RZERO (WX,NTOT1)
1007 CALL RZERO (WY,NTOT1)
1008 CALL RZERO (WZ,NTOT1)
1009C
1010 DO 1000 ISWEEP=1,2
1011C
1012 IFLD = 0
1013 DO 110 IEL=1,NEL
1014 DO 110 IFC=1,NFACE
1015 ieg = lglel(iel)
1016 CB = CBC(IFC,IEL,IFLD)
1017 IF (CB.EQ.'M' .OR. CB.EQ.'m') THEN
1018 CALL FACIND (KX1,KX2,KY1,KY2,KZ1,KZ2,lx1,ly1,lz1,IFC)
1019 DO 140 IZ=KZ1,KZ2
1020 DO 140 IY=KY1,KY2
1021 DO 140 IX=KX1,KX2
1022 CALL INIGEOM (WX(IX,IY,IZ,IEL),WY(IX,IY,IZ,IEL),
1023 $ WZ(IX,IY,IZ,IEL),XM1(IX,IY,IZ,IEL),
1024 $ YM1(IX,IY,IZ,IEL),ZM1(IX,IY,IZ,IEL),
1025 $ IFC,IEG)
1026 140 CONTINUE
1027 ENDIF
1028 110 CONTINUE
1029C
1030 IF (ISWEEP.EQ.1) THEN
1031 CALL DSOP (WX,'MXA',lx1,ly1,lz1)
1032 CALL DSOP (WY,'MXA',lx1,ly1,lz1)
1033 IF (ldim.EQ.3) CALL DSOP (WZ,'MXA',lx1,ly1,lz1)
1034 ELSE
1035 CALL DSOP (WX,'MNA',lx1,ly1,lz1)
1036 CALL DSOP (WY,'MNA',lx1,ly1,lz1)
1037 IF (ldim.EQ.3) CALL DSOP (WZ,'MNA',lx1,ly1,lz1)
1038 ENDIF
1039C
1040 1000 CONTINUE
1041C
1042 return
1043 end
1044c-----------------------------------------------------------------------
1045 subroutine inigeom (ux,uy,uz,x,y,z,iside,iel)
1046C
1047 include 'SIZE'
1048 include 'TSTEP'
1049C
1050 UX = 0.0
1051 UY = 0.0
1052 UZ = 0.0
1053C
1054 return
1055 end
1056c-----------------------------------------------------------------------
1057 subroutine quickmv
1058 include 'SIZE'
1059 include 'TOTAL'
1060 include 'ZPER'
1061c
1062 if (if3d) then
1063 call quickmv3d
1064 else
1065 call quickmv2d
1066 endif
1067 return
1068 end
1069c-----------------------------------------------------------------------
1070 subroutine quickmv2d
1071 include 'SIZE'
1072 include 'TOTAL'
1073 include 'ZPER'
1074c
1075 integer e,ex,ey,ez,eg
1076 common /surfa/ zsurf(lx1,lz1,lelx,lely)
1077 $ , wsurf(lx1,lz1,lelx,lely)
1078 real nxs,nys,nzs
1079c
1080 icount = 0
1081 do ex=1,nelx
1082 do ix=1,lx1
1083 zsurf(ix,1,ex,1) = -1.e20
1084 wsurf(ix,1,ex,1) = -1.e20
1085 ey=nely
1086 eg = ex + nelx*(ey-1)
1087 mid = gllnid(eg)
1088 e = gllel (eg)
1089 if (mid.eq.nid) then
1090 zsurf(ix,1,ex,1) = ym1(ix,ly1,1,e)
1091 vxs = vx (ix,ly1,1,e)
1092 vys = vy (ix,ly1,1,e)
1093 nxs = unx(ix,1,3,e) ! Face 3 is on top in 2D
1094 nys = uny(ix,1,3,e)
1095 gamma_s = (vxs*nxs + vys*nys)/(nys)
1096 wsurf(ix,1,ex,1) = gamma_s ! vertical component of wsurf
1097 endif
1098 zsurf(ix,1,ex,1) = glmax(zsurf(ix,1,ex,1),1)
1099 wsurf(ix,1,ex,1) = glmax(wsurf(ix,1,ex,1),1)
1100 icount = icount+1
1101
1102c write(6,6) ex,e,ix,xm1(ix,ly1,1,e),ym1(ix,ly1,1,e)
1103c $ ,vxs,vys,nxs,nys,gamma_s,wsurf(ix,1,ex,1),zsurf(ix,1,ex,1)
1104c 6 format(3i3,1p9e12.4,' srf')
1105
1106 enddo
1107 enddo
1108 zmin = glmin(ym1,lx1*ly1*lz1*nelv)
1109c
1110 do ex=1,nelx
1111 do ix=1,lx1
1112 do ey=1,nely
1113 eg = ex + nelx*(ey-1)
1114 mid = gllnid(eg)
1115 e = gllel (eg)
1116 if (mid.eq.nid) then
1117 do iy=1,ly1
1118 wy (ix,iy,1,e) = wsurf(ix,1,ex,1)
1119 $ * (ym1(ix,iy,1,e)-zmin)/(zsurf(ix,1,ex,1)-zmin)
1120 enddo
1121 endif
1122 enddo
1123 enddo
1124 enddo
1125c
1126 n = lx1*ly1*lz1*nelv
1127 call rzero(wx,n)
1128 call dsavg(wy)
1129
1130c call opcopy(vx,vy,vz,wx,wy,wz)
1131c call outpost(vx,vy,vz,pr,t,' ')
1132c call exitt
1133
1134 return
1135 end
1136c-----------------------------------------------------------------------
1137 subroutine quickmv3d
1138 include 'SIZE'
1139 include 'TOTAL'
1140 include 'ZPER'
1141c
1142 integer e,ex,ey,ez,eg
1143 common /surfa/ zsurf(lx1,lz1,lelx,lely)
1144 $ , wsurf(lx1,lz1,lelx,lely)
1145 real nxs,nys,nzs
1146c
1147 icount = 0
1148 do ey=1,nely
1149 do ex=1,nelx
1150 do iy=1,ly1
1151 do ix=1,lx1
1152 zsurf(ix,iy,ex,ey) = -1.e20
1153 wsurf(ix,iy,ex,ey) = -1.e20
1154 ez = nelz
1155 eg = ex + nelx*(ey-1) + nelx*nely*(ez-1)
1156 mid = gllnid(eg)
1157 e = gllel (eg)
1158 if (mid.eq.nid) then
1159 zsurf(ix,iy,ex,ey) = zm1(ix,iy,lz1,e)
1160 vxs = vx (ix,iy,lz1,e)
1161 vys = vy (ix,iy,lz1,e)
1162 vzs = vz (ix,iy,lz1,e)
1163 nxs = unx(ix,iy,6,e) ! Face 6 is on top in 3D
1164 nys = uny(ix,iy,6,e)
1165 nzs = unz(ix,iy,6,e)
1166 gamma_s = (vxs*nxs+vys*nys+vzs*nzs)/(nzs)
1167 wsurf(ix,iy,ex,ey) = gamma_s ! vertical component of wsurf
1168 endif
1169 zsurf(ix,iy,ex,ey) = glmax(zsurf(ix,iy,ex,ey),1)
1170 wsurf(ix,iy,ex,ey) = glmax(wsurf(ix,iy,ex,ey),1)
1171 icount = icount+1
1172 enddo
1173 enddo
1174 enddo
1175 enddo
1176
1177 n = lx1*ly1*lz1*nelv
1178 zmin = glmin(zm1,n)
1179c
1180 do ey=1,nely
1181 do ex=1,nelx
1182 do iy=1,ly1
1183 do ix=1,lx1
1184 do ez=1,nelz
1185 eg = ex + nelx*(ey-1) + nelx*nely*(ez-1)
1186 mid = gllnid(eg)
1187 e = gllel (eg)
1188 if (mid.eq.nid) then
1189 do iz=1,lz1
1190 wz (ix,iy,iz,e) = wsurf(ix,iy,ex,ey)
1191 $ * (zm1(ix,iy,iz,e)-zmin)/(zsurf(ix,iy,ex,ey)-zmin)
1192 enddo
1193 endif
1194 enddo
1195 enddo
1196 enddo
1197 enddo
1198 enddo
1199c
1200 n = lx1*ly1*lz1*nelv
1201 call rzero(wx,n)
1202 call rzero(wy,n)
1203 call dsavg(wz)
1204c
1205 return
1206 end
1207c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.