Actual source code: user.F

  1: !---------------------------------------------------------------
  2: ! The following subroutines are from node2t.f in the original
  3: ! code - D. K. Kaushik (1/17/97)
  4: !---------------------------------------------------------------
  5: !
  6: !
  7: ! 2-D Navier Stokes on Unstructured Grids

  9: !============================== FORLINK ==============================72
 10: !
 11: !  FORLINK establishes links between FORTRAN common blocks and C
 12: !
 13: !=====================================================================72
 14:       subroutine FORLINK()

 16:       implicit none
 17: #include <petsc/finclude/petscsys.h>
 18:       PetscScalar  title(20),beta,alpha
 19:       PetscScalar  Re,dt,tot,res0,resc
 20:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
 21:       PetscScalar  cfl1,cfl2
 22:       integer nsmoth,iflim,itran,nbtran,jupdate,                        &
 23:      &        nstage,ncyct,iramp,nitfo
 24:       PetscScalar  gtol
 25:       integer icycle,nsrch,ilu0,ifcn
 26:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                 &
 27:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
 28:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,         &
 29:      &             nstage,ncyct,iramp,nitfo
 30:       common/gmcom/gtol,icycle,nsrch,ilu0,ifcn

 32:       call CLINK(title,cfl1,gtol)
 33: !
 34: ! End of subroutine FORLINK
 35: !
 36:       return
 37:       end

 39: !============================== Block_Initialization  ======== =======72
 40: !
 41: ! Initializes the common blocks members for turbulence model
 42: !
 43: !=====================================================================72
 44:       block data Block_Initialization
 45:       implicit none
 46: #include <petsc/finclude/petscsys.h>
 47:       PetscScalar vkar,cmu,ce1,ce2
 48:       PetscScalar aplus1,aplus2,turbinf
 49:       PetscScalar cb1,sig,cb2,cw1,cw2
 50:       PetscScalar cw3,cv1,ct1,ct2,ct3,ct4
 51:       common/turb/vkar,cmu,ce1,ce2,aplus1,aplus2,turbinf
 52:       common/spalrt/cb1,sig,cb2,cw1,cw2,cw3,cv1,ct1,ct2,ct3,ct4
 53:       data vkar,cmu,ce1,ce2/0.41,0.09,1.2,2.0 /
 54:       data aplus1,aplus2,turbinf/26.0,10.0,0.1 /
 55:       data cb1,sig,cb2,cw2,cw3/0.1355,0.66667,0.622,0.3,2.0/
 56: !
 57: ! Comment out old coefficients
 58: !      data cv1,ct1,ct2,ct3,ct4/7.1,1.0,2.0,1.1,2.0/
 59: !
 60:        data cv1,ct1,ct2,ct3,ct4/7.1d0,1.0d0,2.0d0,1.2d0,0.5d0/

 62:       end

 64: !================================== INIT =============================72
 65: !
 66: ! Initializes the flow field
 67: !
 68: !=====================================================================72
 69:       subroutine INIT(nnodes, qvec, turbre, amut,nvnode, ivnode,irank)
 70:       implicit none
 71: #include <petsc/finclude/petscsys.h>
 72:       PetscScalar  turbre(1),amut(1)
 73:       integer ivnode(1),nnodes,nvnode,irank
 74:       PetscScalar  title(20),beta,alpha
 75:       PetscScalar  Re,dt,tot,res0,resc
 76:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
 77:       PetscScalar  gamma,gm1,gp1,gm1g,gp1g,ggm1
 78:       PetscScalar  p0,rho0,c0,u0,v0,w0,et0,h0,pt0
 79:       PetscScalar  vkar,cmu,ce1,ce2
 80:       PetscScalar  aplus1,aplus2,turbinf
 81:       PetscScalar  cb1,sig,cb2,cw1,cw2,cw3
 82:       PetscScalar  cv1,ct1,ct2,ct3,ct4
 83:       PetscScalar  cfl1,cfl2
 84:       integer nsmoth,iflim,itran,nbtran,jupdate,                          &
 85:      &        nstage,ncyct,iramp,nitfo
 86:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                   &
 87:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
 88:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
 89:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
 90:       common/turb/vkar,cmu,ce1,ce2,aplus1,aplus2,turbinf
 91:       common/spalrt/cb1,sig,cb2,cw1,cw2,cw3,cv1,ct1,ct2,ct3,ct4
 92:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,           &
 93:      &             nstage,ncyct,iramp,nitfo
 94:       PetscScalar  pi,conv,rmu,chi,fv1
 95:       Integer i,k,n
 96: #if defined(INTERLACING)
 97:        PetscScalar qvec(4,nnodes)
 98: #define qnod(i,j) qvec(i,j)
 99: #else
100:        PetscScalar qvec(nnodes,4)
101: #define qnod(i,j) qvec(j,i)
102: #endif
103: !
104:       cw1 = cb1/vkar/vkar + (1.0d0 + cb2)/sig
105:       if (ivisc.gt.4d0)turbinf = 1.341946d0
106: !      if (ivisc.gt.4)turbinf = 0.5
107:       if (ivisc.gt.4.and.itran.eq.1)turbinf = 0.1d0
108: !
109: ! Note that for Spalarts model, I use turbinf as the freestream value of
110: ! the dependent variable just as in the Baldwin-Barth model.
111: ! The constant is set so that in the freestream, nu_t=0.009 (1.341946)
112: !
113: !     print *, "I am in INIT"
114:       res0 = 1.0d0
115:       resc = 1.0d0

117:       gamma = 1.0d0
118:       gm1   = 1.0d0
119:       gp1   = 1.0d0
120:       gm1g  = 1.0d0
121:       gp1g  = 1.0d0
122:       ggm1  = 1.0d0

124:       pi = 4.0d0*datan(1.0d0)
125:       conv = 180.0d0/pi
126:       p0    = 1.0d0
127: #if defined(CFL3D_AXIS)
128:       u0    = cos(alpha/conv)
129:       v0    = 0.0d0
130:       w0    = sin(alpha/conv)
131: #else
132:       u0    = cos(alpha/conv)
133:       v0    = sin(alpha/conv)
134:       w0    = 0.0d0
135: #endif
136:       do n = 1,nnodes
137:          qnod(1,n) = p0
138:          qnod(2,n) = u0
139:          qnod(3,n) = v0
140:          qnod(4,n) = w0
141:       enddo
142:       if (ivisc.eq.3) then
143:        do n = 1,nnodes
144:          turbre(n)  = 0.0d0
145:          amut(n)    = 0.0d0
146:        enddo
147:       endif
148: !
149: ! If viscous, zero out the velocity on the surface
150: !
151: !     print *, "Just Before Viscous"

153:       do 9010 i = 1,nvnode
154: !
155: ! Compute the velocity normal to the surface
156: !
157:         k       = ivnode(i)
158:         qnod(2,k) = 0.0d0
159:         qnod(3,k) = 0.0d0
160:         qnod(4,k) = 0.0d0
161:  9010 continue
162: !
163: ! If turbulent, initialize turbre
164: !
165:       if (ivisc.eq.3.or.ivisc.eq.4) then
166:          if (irank .eq. 0) write(10,110)
167:          if (irank .eq. 0) write(10,120)vkar,cmu,ce1,ce2
168:          if (irank .eq. 0) write(10,130)aplus1,aplus2,turbinf
169:          do 1010 n = 1,nnodes
170:             turbre(n) = turbinf
171:             amut(n)   = cmu*turbinf
172:  1010    continue
173:       end if

175:       if (ivisc.eq.5.or.ivisc.eq.6) then
176:          if (irank .eq. 0) write(10,110)
177:          if (irank .eq. 0) write(10,140)vkar,cb1,sig,cb2
178:          if (irank .eq. 0) write(10,150)cw1,cw2,cw3,cv1
179:          if (irank .eq. 0) write(10,160)ct1,ct2,ct3,ct4
180:          do 1020 n = 1,nnodes
181:             turbre(n) = turbinf
182:             rmu = 1.0d0
183:             chi = turbre(n)/rmu
184:             fv1 = chi**3/(chi**3 + cv1**3)
185:             amut(n)   = fv1*turbinf
186:  1020    continue
187:       end if
188: !     print *, "I am out of INIT"

190:       return
191:   110 format(1h ,'Parameters for turbulence model')
192:   120 format(1h ,'k=',f10.5,' cmu=',f10.5,' ce1=',f10.5,'ce2=',f10.5)
193:   130 format(1h ,'aplus1',f10.5,' aplus2=',f10.5,' turbinf=',f10.5)
194:   140 format(1h ,'k=',f10.5,' cb1=',f10.5,' sig=',f10.5,'cb2=',f10.5)
195:   150 format(1h ,'cw1=',f10.5,' cw2=',f10.5,' cw3=',f10.5,' cv1=',f10.5)
196:   160 format(1h ,'ct1=',f10.5,' ct2=',f10.5,' ct3=',f10.5,' ct4=',f10.5)
197: !
198: ! End of subroutine INIT
199: !

201:       end

203: !================================ READR1 ====================================
204: !
205: !  Reads input parameters
206: !
207: !============================================================================
208:       subroutine READR1(ileast, irank)

210:       implicit none
211: #include <petsc/finclude/petscsys.h>
212:       PetscScalar  title(20),beta,alpha
213:       PetscScalar  Re,dt,tot,res0,resc
214:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
215:       PetscScalar  cfl1,cfl2
216:       integer nsmoth,iflim,itran,nbtran,jupdate,                                &
217:      &        nstage,ncyct,iramp,nitfo
218:       PetscScalar  gtol
219:       integer icycle,nsrch,ilu0,ifcn
220:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                         &
221:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
222:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,                 &
223:      &             nstage,ncyct,iramp,nitfo
224:       common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
225:       integer ileast, irank
226:       integer i

228:       read(7,10)(title(i),i=1,20)
229:       if (irank .eq. 0) write(10,11)(title(i),i=1,20)

231:       read(7,10)

233:       read(7,24)mseq,ihane,ivisc,ileast,iflim,jupdate,ifcn
234:       if (irank .eq. 0) write(10,25)mseq,ihane,ivisc
235:       if (irank .eq. 0) write(10,28)ileast,iflim,jupdate

237:       read(7,10)

239:       read(7,12)beta,alpha,Re
240:       if (irank .eq. 0) write(10,13)beta,alpha,Re

242:       read(7,10)

244:       read(7,26)cfl2,dt,irest,itran,nbtran
245:       if (irank .eq. 0) write(10,27)cfl2,dt,irest

247:       if (irank .eq. 0) then
248:        if (ivisc.eq.5.or.ivisc.eq.6) write(10,123)itran,nbtran
249:       endif

251:    10 format(20a4)
252:    11 format(1h ,20a4)
253:    12 format(2f10.5,e14.7)
254:    13 format(1h ,'beta = ',f10.5,' Alpha = ',f10.5,' Re = ',e14.7)
255:    24 format(i10,i10,i10,i10,i10,i10,i10)
256:    25 format(1h ,'mseq = ',i3,' ihane = ',i3,' ivisc=',i3)
257:    26 format(2f10.5,3i10)
258:    27 format(1h ,' cfl2= ',e14.7,' dt= ',f10.5,'irest= ',i5)
259:    28 format(1h ,'ileast= ',i5,' iflim= ',i5,' jupdate= ',i5)
260:   123 format(1h ,'itran = ',i5,' nbtran= ',i5)

262:       return
263:       end

265: !================================ TECFLO =============================72
266: !
267: !  Writes a formatted file that contains an input file for TECPLOT
268: !
269: !=====================================================================72
270:       subroutine TECFLO(nnodes,                                            &
271:      &                  nnbound,nvbound,nfbound,                           &
272:      &                  nnfacet,nvfacet,nffacet,                           &
273:      &                  nsnode,nvnode,nfnode,                              &
274:      &                  title,x,y,z,q,                                     &
275:      &                  nnpts,nntet,nvpts,nvtet,nfpts,nftet,               &
276:      &                  f2ntn,f2ntv,f2ntf,isnode,ivnode,ifnode,            &
277:      &                  irank)
278: !
279:       implicit none
280: #include <petsc/finclude/petscsys.h>
281:       PetscScalar pinf
282:       PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
283:       PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
284:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
285:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
286:       integer nnodes,nnbound,nvbound,nfbound,                              &
287:      &        nnfacet,nvfacet,nffacet,                                     &
288:      &        nsnode,nvnode,nfnode,irank
289:       integer nntet(1),nnpts(1)
290:       integer nvtet(1),nvpts(1)
291:       integer nftet(1),nfpts(1)
292:       integer f2ntn(nnfacet,4),f2ntv(nvfacet,4),f2ntf(nffacet,4)
293:       integer isnode(1),ivnode(1),ifnode(1)
294:       PetscScalar  x(1),y(1),z(1),q(4,nnodes)
295:       integer i,j,i1,i2,i3,ic,isp,ist

297:       character c4*4,title*20

299:       pinf = p0
300: !
301: !   + write tecplot header
302: !
303:       write(13,'(a)') 'TITLE="' // title // '"'
304:       write(13,'(a)') 'VARIABLES="X     ","Y     ","Z     ","RHO   ",',     &
305:      &     '"U     ","V     ","W     ","P/Pinf"'
306: !
307: !   + do a zone-title, so we can keep track of things
308: !   + start with solid-wall boundary surfaces
309: !
310:       isp   = 1
311:       ist   = 1
312:       do 10 i=1,nnbound
313:         write(c4,"(i4)") i
314:         if (i .ge.    0 .and. i .le.    9) ic = 4
315:         if (i .ge.   10 .and. i .le.   99) ic = 3
316:         if (i .ge.  100 .and. i .le.  999) ic = 2
317:         if (i .ge. 1000 .and. i .le. 9999) ic = 1

319:         write(13,1000) "nn." // c4(ic:4),nnpts(i),nntet(i)

321:         write(13,1010) (x(isnode(j))      ,j=isp,isp+nnpts(i)-1)
322:         write(13,1010) (y(isnode(j))      ,j=isp,isp+nnpts(i)-1)
323:         write(13,1010) (z(isnode(j))      ,j=isp,isp+nnpts(i)-1)
324:         write(13,1010) (1.0                 ,j=isp,isp+nnpts(i)-1)
325:         write(13,1010) (q(2,isnode(j))      ,j=isp,isp+nnpts(i)-1)
326:         write(13,1010) (q(3,isnode(j))      ,j=isp,isp+nnpts(i)-1)
327:         write(13,1010) (q(4,isnode(j))      ,j=isp,isp+nnpts(i)-1)
328:         write(13,1010) (q(1,isnode(j))/pinf ,j=isp,isp+nnpts(i)-1)

330:         do 30 j=ist,ist+nntet(i)-1
331:           i1 = f2ntn(j,1) - isp + 1
332:           i2 = f2ntn(j,2) - isp + 1
333:           i3 = f2ntn(j,3) - isp + 1

335:           write(13,1020) i1,i2,i3,i3
336:    30   continue

338:         isp = isp + nnpts(i)
339:         ist = ist + nntet(i)

341:    10 continue
342: !
343: !   + start with viscous boundary surfaces
344: !
345:       isp   = 1
346:       ist   = 1
347:       do 70 i=1,nvbound
348:         write(c4,"(i4)") i
349:         if (i .ge.    0 .and. i .le.    9) ic = 4
350:         if (i .ge.   10 .and. i .le.   99) ic = 3
351:         if (i .ge.  100 .and. i .le.  999) ic = 2
352:         if (i .ge. 1000 .and. i .le. 9999) ic = 1

354:         write(13,1000) "nv." // c4(ic:4),nvpts(i),nvtet(i)

356:         write(13,1010) (x(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
357:         write(13,1010) (y(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
358:         write(13,1010) (z(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
359:         write(13,1010) (1.0                 ,j=isp,isp+nvpts(i)-1)
360:         write(13,1010) (q(2,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
361:         write(13,1010) (q(3,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
362:         write(13,1010) (q(4,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
363:         write(13,1010) (q(1,ivnode(j))/pinf ,j=isp,isp+nvpts(i)-1)

365:         do 90 j=ist,ist+nvtet(i)-1
366:            i1 = f2ntv(j,1) - isp + 1
367:            i2 = f2ntv(j,2) - isp + 1
368:            i3 = f2ntv(j,3) - isp + 1

370:            write(13,1020) i1,i2,i3,i3
371:    90   continue

373:         isp = isp + nvpts(i)
374:         ist = ist + nvtet(i)

376:    70 continue
377: !
378: !   + do the far-field boundary surfaces
379: !
380:       isp   = 1
381:       ist   = 1
382:       do 40 i=1,nfbound
383:          write(c4,"(i4)") i
384:          if (i .ge.    0 .and. i .le.    9) ic = 4
385:          if (i .ge.   10 .and. i .le.   99) ic = 3
386:          if (i .ge.  100 .and. i .le.  999) ic = 2
387:          if (i .ge. 1000 .and. i .le. 9999) ic = 1

389:          write(13,1000) "ff." // c4(ic:4),nfpts(i),nftet(i)

391:          write(13,1010) (x(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
392:          write(13,1010) (y(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
393:          write(13,1010) (z(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
394:          write(13,1010) (1.0                 ,j=isp,isp+nfpts(i)-1)
395:          write(13,1010) (q(2,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
396:          write(13,1010) (q(3,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
397:          write(13,1010) (q(4,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
398:          write(13,1010) (q(1,ifnode(j))/pinf ,j=isp,isp+nfpts(i)-1)

400:          do 60 j=ist,ist+nftet(i)-1
401:             i1 = f2ntf(j,1) - isp + 1
402:             i2 = f2ntf(j,2) - isp + 1
403:             i3 = f2ntf(j,3) - isp + 1

405:             write(13,1020) i1,i2,i3,i3
406:    60    continue

408:          isp = isp + nfpts(i)
409:          ist = ist + nftet(i)

411:    40  continue

413: !
414: !      End of subroutine TECFLO
415: !
416:  1000  format('ZONE T="',a,'", Z=0, I=',i6,', J=',i6,', F=FEBLOCK')
417:  1010  format(1P10E13.5)
418:  1020  format(4I10)

420:        return
421:        end

423: !================================ FORCE  =============================72
424: !
425: !  Calculates the forces
426: !  Modified - D. K. Kaushik (1/15/97)
427: !  Added new parameters - clift, cdrag, cmom, irank, nvertices
428: !
429: !=====================================================================72
430:       subroutine FORCE(nnodes,nedge,                                    &
431:      &                 isnode,ivnode,                                   &
432:      &                 nnfacet,f2ntn,nnbound,                           &
433:      &                 nvfacet,f2ntv,nvbound,                           &
434:      &                 evec,qvec,xyz,sface_bit,vface_bit,               &
435:      &                 clift, cdrag, cmom,irank,nvertices)
436:       implicit none
437: #include <petsc/finclude/petscsys.h>
438:       PetscScalar  title(20),beta,alpha
439:       PetscScalar  Re,dt,tot,res0,resc
440:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
441:       PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
442:       PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
443:       PetscScalar rms(1001),clw(1001),cdw(1001),                             &
444:      &       cmw(1001),xres(1001)
445:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                 &
446:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
447:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
448:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
449:       common/history/rms,clw,cdw,cmw,xres

451:       integer nnodes,nedge,nnfacet,nvfacet,                             &
452:      &        nnbound,nvbound,                                          &
453:      &        irank,nvertices
454:       integer isnode(1),ivnode(1)
455:       integer f2ntn(nnfacet,4)
456:       integer f2ntv(nvfacet,4)
457:       integer sface_bit(nnfacet), vface_bit(nvfacet)

459:       PetscScalar clift, cdrag, cmom
460:       PetscScalar clloc,cdloc,cmloc
461:       PetscScalar clglo,cdglo,cmglo
462:       PetscScalar cl_ghost, cd_ghost, cm_ghost
463:       PetscScalar coef_loc(3), coef_glo(3)
464:       PetscScalar pi,conv,csa,sna,xr,yr,zr
465:       PetscScalar x1,y1,z1,x2,y2,z2,x3,y3
466:       PetscScalar z3,xmid,ymid,zmid
467:       PetscScalar p1,p2,p3,u1,u2,u3,v1,v2,v3
468:       PetscScalar w1,w2,w3,ax,ay,az,bx,by,bz
469:       PetscScalar xnorm,ynorm,znorm
470:       PetscScalar dcx,dcy,dcz,cp,press
471:       integer icount,n,node1,node2,node3,ierr

473: #if defined(INTERLACING)
474:       PetscScalar qvec(4,nvertices)
475:       PetscScalar xyz(3,nvertices)
476:       integer evec(2,nedge)
477: #define qnod(i,j) qvec(i,j)
478: #define x(i) xyz(1,i)
479: #define y(i) xyz(2,i)
480: #define z(i) xyz(3,i)
481: #define eptr(j,i) evec(i,j)
482: #else
483:       PetscScalar qvec(nvertices,4)
484:       PetscScalar xyz(nvertices,3)
485:       integer evec(nedge,2)
486: #define qnod(i,j) qvec(j,i)
487: #define x(i) xyz(i,1)
488: #define y(i) xyz(i,2)
489: #define z(i) xyz(i,3)
490: #define eptr(i,j) evec(i,j)
491: #endif
492: !
493:       pi = 4.d0*atan(1.d0)
494:       conv = 180.d0/pi
495:       csa=cos(alpha/conv)
496:       sna=sin(alpha/conv)
497: !
498: !  initialize forces to zero
499: !
500:       clw(ntt) = 0.0d0
501:       cdw(ntt) = 0.0d0
502:       cmw(ntt) = 0.0d0
503:       cl_ghost = 0.0d0
504:       cd_ghost = 0.0d0
505:       cm_ghost = 0.0d0
506:       clloc = 0.0d0
507:       cdloc = 0.0d0
508:       cmloc = 0.0d0
509:       clglo = 0.0d0
510:       cdglo = 0.0d0
511:       cmglo = 0.0d0
512:       coef_glo(1) = 0.0d0
513:       coef_glo(2) = 0.0d0
514:       coef_glo(3) = 0.0d0

516:       xr = 0.25d0
517:       yr = 0.0d0
518:       zr = 0.0d0

520:       do 30 n = 1, nnfacet
521:                node1 = isnode(f2ntn(n,1))
522:                node2 = isnode(f2ntn(n,2))
523:                node3 = isnode(f2ntn(n,3))

525:                x1    = x(node1)
526:                y1    = y(node1)
527:                z1    = z(node1)

529:                x2    = x(node2)
530:                y2    = y(node2)
531:                z2    = z(node2)

533:                x3 = x(node3)
534:                y3 = y(node3)
535:                z3 = z(node3)

537:                ax = x2 - x1
538:                ay = y2 - y1
539:                az = z2 - z1

541:                bx = x3 - x1
542:                by = y3 - y1
543:                bz = z3 - z1
544: !
545: !  norm points outward, away from grid interior.
546: !  norm magnitude is area of surface triangle.
547: !
548:                xnorm =-0.5d0*(ay*bz - az*by)
549:                ynorm = 0.5d0*(ax*bz - az*bx)
550:                znorm =-0.5d0*(ax*by - ay*bx)

552:                p1 = qnod(1,node1)
553:                u1 = qnod(2,node1)
554:                v1 = qnod(3,node1)
555:                w1 = qnod(4,node1)

557:                p2 = qnod(1,node2)
558:                u2 = qnod(2,node2)
559:                v2 = qnod(3,node2)
560:                w2 = qnod(4,node2)

562:                p3 = qnod(1,node3)
563:                u3 = qnod(2,node3)
564:                v3 = qnod(3,node3)
565:                w3 = qnod(4,node3)

567:                press = (p1 + p2 + p3)/3.0d0
568:                cp    = 2.d0*(press - 1.0d0)
569: !
570:                dcx = cp*xnorm
571:                dcy = cp*ynorm
572:                dcz = cp*znorm

574:                xmid = (x1 + x2 + x3)
575:                ymid = (y1 + y2 + y3)
576:                zmid = (z1 + z2 + z3)

578:                if (sface_bit(n) .eq. 1) then
579: #if defined(CFL3D_AXIS)

581:                 clloc = clloc - dcx*sna     + dcz*csa
582:                 cdloc = cdloc + dcx*csa + dcz*sna
583:                 cmloc = cmloc                                           &
584:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx

586: #else
587:                 clloc = clloc - dcx*sna + dcy*csa
588:                 cdloc = cdloc + dcx*csa + dcy*sna
589:                 cmloc = cmloc                                           &
590:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx
591: #endif
592:                endif

594:  30    continue
595: !
596: ! Viscous boundary
597: !
598:        do 60 n = 1, nvfacet
599:                node1 = ivnode(f2ntv(n,1))
600:                node2 = ivnode(f2ntv(n,2))
601:                node3 = ivnode(f2ntv(n,3))

603:                x1    = x(node1)
604:                y1    = y(node1)
605:                z1    = z(node1)

607:                x2    = x(node2)
608:                y2    = y(node2)
609:                z2    = z(node2)

611:                x3 = x(node3)
612:                y3 = y(node3)
613:                z3 = z(node3)

615:                ax = x2 - x1
616:                ay = y2 - y1
617:                az = z2 - z1

619:                bx = x3 - x1
620:                by = y3 - y1
621:                bz = z3 - z1
622: !
623: !  norm points outward, away from grid interior.
624: !  norm magnitude is area of surface triangle.
625: !
626:                xnorm =-0.5d0*(ay*bz - az*by)
627:                ynorm = 0.5d0*(ax*bz - az*bx)
628:                znorm =-0.5d0*(ax*by - ay*bx)

630:                p1    = qnod(1,node1)
631:                u1    = qnod(2,node1)
632:                v1    = qnod(3,node1)
633:                w1    = qnod(4,node1)

635:                p2    = qnod(1,node2)
636:                u2    = qnod(2,node2)
637:                v2    = qnod(3,node2)
638:                w2    = qnod(4,node2)

640:                p3    = qnod(1,node3)
641:                u3    = qnod(2,node3)
642:                v3    = qnod(3,node3)
643:                w3    = qnod(4,node3)

645:                press = (p1 + p2 + p3)/3.0d0
646:                cp    = 2.d0*(press -1.0d0)

648:                dcx = cp*xnorm
649:                dcy = cp*ynorm
650:                dcz = cp*znorm

652:                xmid = (x1 + x2 + x3)
653:                ymid = (y1 + y2 + y3)
654:                zmid = (z1 + z2 + z3)

656:                if (vface_bit(n) .eq. 1) then
657:                 clloc = clloc - dcx*sna + dcy*csa
658:                 cdloc = cdloc + dcx*csa + dcy*sna
659:                 cmloc = cmloc                                            &
660:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx
661:                endif

663:  60   continue
664: !
665:       icount = 3
666:       coef_loc(1) = clloc
667:       coef_loc(2) = cdloc
668:       coef_loc(3) = cmloc
669:       call MPI_ALLREDUCE(coef_loc,coef_glo,icount,                       &
670:      &                   MPIU_SCALAR, MPIU_SUM,                           &
671:      &                   MPI_COMM_WORLD,ierr)
672: !     call MPI_ALLREDUCE(cdloc,cdglo,icount,
673: !    &                   MPIU_REAL_PRECISION, MPI_SUM,
674: !    &                   MPI_COMM_WORLD,ierr)
675: !     call MPI_ALLREDUCE(cmloc,cmglo,icount,
676: !    &                   MPIU_REAL_PRECISION, MPI_SUM,
677: !    &                   MPI_COMM_WORLD,ierr)
678: !
679:        clw(ntt) = coef_glo(1)
680:        cdw(ntt) = coef_glo(2)
681:        cmw(ntt) = coef_glo(3)
682:        clift = clw(ntt)
683:        cdrag = cdw(ntt)
684:        cmom = cmw(ntt)
685: !
686: ! Update the timing information - Added by D. K. Kaushik (1/17/97)
687:       xres(ntt) = tot
688: !
689: ! End of subroutine FORCE
690: !
691:       return
692:       end

694: !================================ DELTAT2 ============================72
695: !
696: ! Calculate a time step for each cell
697: ! Note that this routine assumes conservative variables
698: !
699: !=====================================================================72
700:       subroutine DELTAT2(nnodes,nedge,qvec,cdt,                          &
701:      &                  xyz,vol,xyzn,evec,                               &
702:      &                  sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,             &
703:      &                  nsnode,nvnode,nfnode,isnode,ivnode,ifnode,       &
704:      &                  LocalTS,irank,nvertices)
705: !
706:       implicit none
707: #include <petsc/finclude/petscsys.h>
708:       PetscScalar  title(20),beta,alpha
709:       PetscScalar  Re,dt,tot,res0,resc
710:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
711:       PetscScalar  gamma,gm1,gp1,gm1g,gp1g,ggm1
712:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                  &
713:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
714:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
715:       integer nnodes,nedge,nsnode,nvnode,nfnode,LocalTS,irank,nvertices
716:       PetscScalar  vol(1)
717:       PetscScalar  sxn(1),syn(1),szn(1)
718:       PetscScalar  vxn(1),vyn(1),vzn(1)
719:       PetscScalar  fxn(1),fyn(1),fzn(1)
720:       PetscScalar  cdt(1)
721:       integer isnode(1),ivnode(1),ifnode(1)
722:       integer n,i,node1,node2,inode
723:       integer ierr
724:       PetscLogDouble flops
725:       PetscScalar xnorm,ynorm,znorm,area
726:       PetscScalar u1,v1,w1,u2,v2,w2,u,v,c,w
727:       PetscScalar ubar,Vn,term
728: #if defined(INTERLACING)
729:       PetscScalar qvec(4,nvertices)
730:       PetscScalar xyz(3,nvertices)
731:       PetscScalar xyzn(4,nedge)
732:       integer evec(2,nedge)
733: #define qnod(i,j) qvec(i,j)
734: #define x(i) xyz(1,i)
735: #define y(i) xyz(2,i)
736: #define z(i) xyz(3,i)
737: #define xn(i) xyzn(1,i)
738: #define yn(i) xyzn(2,i)
739: #define zn(i) xyzn(3,i)
740: #define rl(i) xyzn(4,i)
741: #define eptr(j,i) evec(i,j)
742: #else
743:       PetscScalar qvec(nvertices,4)
744:       PetscScalar xyz(nvertices,3)
745:       PetscScalar xyzn(nedge,4)
746:       integer evec(nedge,2)
747: #define qnod(i,j) qvec(j,i)
748: #define x(i) xyz(i,1)
749: #define y(i) xyz(i,2)
750: #define z(i) xyz(i,3)
751: #define xn(i) xyzn(i,1)
752: #define yn(i) xyzn(i,2)
753: #define zn(i) xyzn(i,3)
754: #define rl(i) xyzn(i,4)
755: #define eptr(i,j) evec(i,j)
756: #endif
757: !
758: !  If local time steping, loop over faces
759: !  and calculate time step as cdt = V/(sum(|u.n| +c.area)
760: !  This is time step for cfl=1. We will multiply by cfl number later
761: !
762: ! First loop over nodes and zero out cdt
763: !
764:       flops = 0.0
765:       if (LocalTS .gt. 0) then
766:          do 1000 i = 1,nvertices
767:             cdt(i) = 0.0d0
768:  1000    continue
769: !
770:          do 1020 n = 1, nedge
771:                node1 = eptr(n,1)
772:                node2 = eptr(n,2)
773: !
774: ! Get normal to face
775: !
776:                xnorm = xn(n)
777:                ynorm = yn(n)
778:                znorm = zn(n)
779:                area  = rl(n)
780: !
781:                xnorm = xnorm*area
782:                ynorm = ynorm*area
783:                znorm = znorm*area
784: !
785: !/*
786: ! xnorm = xnormal * area of face
787: ! ynorm = ynormal * area of face
788: ! znorm = znormal * area of face
789: !*/
790:                u1   = qnod(2,node1)
791:                v1   = qnod(3,node1)
792:                w1   = qnod(4,node1)

794:                u2   = qnod(2,node2)
795:                v2   = qnod(3,node2)
796:                w2   = qnod(4,node2)
797: !
798: ! Get average values on face
799: !
800:                u    = 0.5d0*(u1 + u2)
801:                v    = 0.5d0*(v1 + v2)
802:                w    = 0.5d0*(w1 + w2)
803:                ubar = xn(n)*u + yn(n)*v + zn(n)*w
804:                c    = sqrt(ubar*ubar + beta)

806:                term = abs(u*xnorm + v*ynorm + w*znorm) + c*area
807:                cdt(node1) = cdt(node1) + term
808:                cdt(node2) = cdt(node2) + term
809: !
810:  1020    continue
811:          flops = flops + 27.0*nedge
812: !
813: ! Now loop over boundaries and close the contours
814: !
815:          do 1030 i = 1,nsnode
816:             inode = isnode(i)
817: !
818: ! Get the normal
819: !
820:             xnorm = sxn(i)
821:             ynorm = syn(i)
822:             znorm = szn(i)
823:             area  = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)

825:             u   = qnod(2,inode)
826:             v   = qnod(3,inode)
827:             w   = qnod(4,inode)
828:             ubar= (u*xnorm + v*ynorm + w*znorm)/area
829:             c   = sqrt(ubar*ubar + beta)
830: !
831:             Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
832:             cdt(inode) = cdt(inode) + Vn
833: !
834:  1030    continue
835:          flops = flops + 24.0*nsnode

837: !
838: ! Now viscous faces
839: !
840: !
841:          do 1040 i = 1,nvnode
842:             inode = ivnode(i)
843: !
844: ! Get the normal
845: !
846:             xnorm = vxn(i)
847:             ynorm = vyn(i)
848:             znorm = vzn(i)
849:             area  = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
850: !
851: ! Because velocities are zero, c is just sqrt(beta)
852: !
853:             c  = sqrt(beta)
854:             Vn = c*area

856:             cdt(inode) = cdt(inode) + Vn
857: !
858:  1040    continue
859:          flops = flops + 9.0*nvnode
860: !
861: ! Now far field
862: !
863:          do 1050 i = 1,nfnode
864:             inode = ifnode(i)
865: !
866: ! Get the normal
867: !
868:             xnorm = fxn(i)
869:             ynorm = fyn(i)
870:             znorm = fzn(i)
871:             area  = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)

873:             u   = qnod(2,inode)
874:             v   = qnod(3,inode)
875:             w   = qnod(4,inode)
876:             ubar= (u*xnorm + v*ynorm + w*znorm)/area
877:             c = sqrt(ubar*ubar + beta)

879:             Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
880:             cdt(inode) = cdt(inode) + Vn
881:  1050    continue
882:          flops = flops + 24.0*nfnode

884:          do 1060 n = 1,nvertices
885:             cdt(n) = vol(n)/cdt(n)
886:  1060    continue
887:          flops = flops + nvertices
888:       else
889: !
890: ! If not doing local time stepping just set cdt=1
891: !
892:          do 1070 n = 1,nvertices
893:             cdt(n) = 1.0d0
894:  1070    continue
895:       end if
896:       call PetscLogFlops(flops,ierr)
897: !
898: ! End of subroutine DELTAT2
899: !
900:       return
901:       end

903: !================================= FLUX ======================================
904: !
905: ! Calculates the fluxes on the face and performs the flux balance
906: !
907: !=============================================================================
908: #if defined(_OPENMP)
909: #if defined(HAVE_EDGE_COLORING)
910:       subroutine FLUX(nnodes, ncell, nedge,                                     &
911:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
912:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
913:      &                nnfacet,f2ntn,nnbound,                                    &
914:      &                nvfacet,f2ntv,nvbound,                                    &
915:      &                nffacet,f2ntf,nfbound,                                    &
916:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
917:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
918:      &                max_threads,                                              &
919:      &                ncolor,ncount,                                            &
920:      &                irank,nvertices)
921: #elif defined(HAVE_REDUNDANT_WORK)
922:       subroutine FLUX(nnodes, ncell, nedge,                                     &
923:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
924:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
925:      &                nnfacet,f2ntn,nnbound,                                    &
926:      &                nvfacet,f2ntv,nvbound,                                    &
927:      &                nffacet,f2ntf,nfbound,                                    &
928:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
929:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
930:      &                max_threads,                                              &
931:      &                resd,                                                     &
932:      &                irank,nvertices)
933: #else
934:       subroutine FLUX(nnodes, ncell, nedge,                                     &
935:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
936:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
937:      &                nnfacet,f2ntn,nnbound,                                    &
938:      &                nvfacet,f2ntv,nvbound,                                    &
939:      &                nffacet,f2ntf,nfbound,                                    &
940:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
941:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
942:      &                max_threads,                                              &
943:      &                nedgeAllThr,part_thr,nedge_thr,edge_thr,xyzn_thr,         &
944:      &                irank,nvertices)
945: #endif
946: #else
947:       subroutine FLUX(nnodes, ncell, nedge,                                     &
948:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
949:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
950:      &                nnfacet,f2ntn,nnbound,                                    &
951:      &                nvfacet,f2ntv,nvbound,                                    &
952:      &                nffacet,f2ntf,nfbound,                                    &
953:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
954:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
955:      &                irank,nvertices)
956: #endif
957:       implicit none
958: #include <petsc/finclude/petscsys.h>
959:       PetscScalar  title(20),beta,alpha
960:       PetscScalar  Re,dt,tot,res0,resc
961:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
962:       PetscScalar  gamma,gm1,gp1,gm1g,gp1g,ggm1
963:       PetscScalar  p0,rho0,c0,u0,v0,w0,et0,h0,pt0
964:       PetscScalar  rms(1001),clw(1001),cdw(1001)
965:       PetscScalar  cmw(1001),xres(1001)
966:       PetscScalar  cfl1,cfl2
967:       integer nsmoth,iflim,itran,nbtran,jupdate,                                &
968:      &        nstage,ncyct,iramp,nitfo
969:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                         &
970:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
971:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
972:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
973:       common/history/rms,clw,cdw,cmw,xres
974:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,                 &
975:      &             nstage,ncyct,iramp,nitfo
976:       integer nnodes, ncell, nedge,                                             &
977:      &        nsface, nvface, nfface,                                           &
978:      &        nsnode, nvnode, nfnode,                                           &
979:      &        nnfacet,nvfacet,nffacet,                                          &
980:      &        nnbound,nvbound,nfbound,                                          &
981:      &        irank,nvertices
982:       PetscScalar  sxn(1),syn(1),szn(1)
983:       PetscScalar  vxn(1),vyn(1),vzn(1)
984:       PetscScalar  fxn(1),fyn(1),fzn(1)
985:       PetscScalar  phi(nvertices,4)
986:       integer isface(1),ivface(1),ifface(1)
987:       integer isnode(1),ivnode(1),ifnode(1)
988:       integer f2ntn(nnfacet,4)
989:       integer f2ntv(nvfacet,4)
990:       integer f2ntf(nffacet,4)
991:       integer i,n,node1,node2,node3,inode
992:       integer ierr,first_edge,last_edge
993:       PetscLogDouble flops
994:       PetscScalar  xmean,ymean,zmean,xnorm,ynorm
995:       PetscScalar  znorm,rlen,area,dot,size
996:       PetscScalar  X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3
997:       PetscScalar  p1,p2,p3
998:       PetscScalar  rx,ry,rz,pL,uL,wL,vL
999:       PetscScalar  pR,uR,vR,wR
1000:       PetscScalar  ubarL,ubarR,ubar,p,u,v,w,c
1001:       PetscScalar  c2,c20
1002:       PetscScalar  ubar0,pi,ui,vi,wi,unorm
1003:       PetscScalar  phi1,phi2,phi3,phi4,phi5
1004:       PetscScalar  phi6,phi7,phi8,phi9
1005:       PetscScalar  eig1,eig2,eig3,eig4
1006:       PetscScalar  dp,du,dv,dw
1007:       PetscScalar  ti11,ti21,ti31,ti41
1008:       PetscScalar  dv1,dv2,dv3,dv4
1009:       PetscScalar  r11,r21,r31,r41
1010:       PetscScalar  r12,r22,r32,r42
1011:       PetscScalar  r13,r23,r33,r43
1012:       PetscScalar  r14,r24,r34,r44
1013:       PetscScalar  t1,t2,t3,t4
1014:       PetscScalar  t11,t21,t31,t41
1015:       PetscScalar  t12,t22,t32,t42
1016:       PetscScalar  t13,t23,t33,t43
1017:       PetscScalar  t14,t24,t34,t44
1018:       PetscScalar  fluxp1,fluxp2,fluxp3,fluxp4
1019:       PetscScalar  fluxm1,fluxm2,fluxm3,fluxm4
1020:       PetscScalar  res1,res2,res3,res4
1021:       PetscScalar  rhs1,rhs2,rhs3,rhs4
1022:       PetscScalar  c68,c18,second
1023:       PetscScalar  ax,ay,az,bx,by,bz
1024:       PetscScalar  pa,pb,pc,ub,vb,wb
1025: #if defined(INTERLACING)
1026:       PetscScalar qvec(4,nvertices)
1027:       PetscScalar grad(3,4,nvertices)
1028:       PetscScalar resvec(4,nnodes)
1029:       PetscScalar xyz(3,nvertices)
1030:       PetscScalar xyzn(4,nedge)
1031:       integer evec(2,nedge)
1032: #define qnod(i,j) qvec(i,j)
1033: #define res(i,j) resvec(i,j)
1034: #define gradx(i,j) grad(1,i,j)
1035: #define grady(i,j) grad(2,i,j)
1036: #define gradz(i,j) grad(3,i,j)
1037: #define x(i) xyz(1,i)
1038: #define y(i) xyz(2,i)
1039: #define z(i) xyz(3,i)
1040: #define xn(i) xyzn(1,i)
1041: #define yn(i) xyzn(2,i)
1042: #define zn(i) xyzn(3,i)
1043: #define ra(i) xyzn(4,i)
1044: #define eptr(j,i) evec(i,j)
1045: #else
1046:       PetscScalar qvec(nvertices,4)
1047:       PetscScalar grad(nvertices,4,3)
1048:       PetscScalar resvec(nnodes,4)
1049:       PetscScalar xyz(nvertices,3)
1050:       PetscScalar xyzn(nedge,4)
1051:       integer evec(nedge,2)
1052: #define qnod(i,j) qvec(j,i)
1053: #define res(i,j) resvec(j,i)
1054: #define gradx(i,j) grad(j,i,1)
1055: #define grady(i,j) grad(j,i,2)
1056: #define gradz(i,j) grad(j,i,3)
1057: #define x(i) xyz(i,1)
1058: #define y(i) xyz(i,2)
1059: #define z(i) xyz(i,3)
1060: #define xn(i) xyzn(i,1)
1061: #define yn(i) xyzn(i,2)
1062: #define zn(i) xyzn(i,3)
1063: #define ra(i) xyzn(i,4)
1064: #define eptr(i,j) evec(i,j)
1065: #endif
1066: #if defined(_OPENMP)
1067:        integer my_thread_id,omp_get_thread_num,                             &
1068:      &         max_threads,omp_get_num_threads,                             &
1069:      &         first_node,last_node,chunk
1070: #if defined(HAVE_EDGE_COLORING)
1071:        integer ncolor,ncount(1)
1072: #elif defined(HAVE_REDUNDANT_WORK)
1073:        PetscScalar resd(4,nnodes)
1074: #undef res
1075: #define res(a,b) resd(a,b)
1076: #else
1077:        integer nedgeAllThr,part_thr(nnodes),nedge_thr(0:max_threads),         &
1078:      &     edge_thr(2,nedgeAllThr)
1079:        PetscScalar xyzn_thr(4,nedgeAllThr)
1080: #undef eptr
1081: #define eptr(a,b) edge_thr(b,a)
1082: #undef xn
1083: #undef yn
1084: #undef zn
1085: #undef ra
1086: #define xn(i) xyzn_thr(1,i)
1087: #define yn(i) xyzn_thr(2,i)
1088: #define zn(i) xyzn_thr(3,i)
1089: #define ra(i) xyzn_thr(4,i)
1090: #endif
1091: #endif
1092: !      logging variables
1093:        integer flux_event, fun3d_class, flag
1094:        character * 16 flux_label
1095:        data flag/-1/,flux_label/'FluxEvaluation  '/
1096:        save flux_event, fun3d_class, flag, flux_label
1097:        if (flag .eq. -1) then
1098:           call PetscClassIdRegister('PetscFun3d',fun3d_class,ierr)
1099:           call PetscLogEventRegister(flux_label,                        &
1100:      &         fun3d_class,flux_event,ierr)
1101:           flag = 1
1102:        endif
1103:        call PetscLogEventBegin(flux_event,ierr)
1104:        flops = 0.0
1105:        second = 1.0d0
1106: !      if (ntt.le.nitfo)second = 0.0
1107: !      print *, 'Second is' , second
1108: #if defined(_OPENMP)
1109: #if defined(HAVE_EDGE_COLORING)
1110:        first_edge = 1
1111:        do i = 1,ncolor
1112:              last_edge = first_edge + ncount(i) - 1
1113: #elif defined(HAVE_REDUNDANT_WORK)
1114:        do i = 1,nnodes
1115:           do j = 1,4
1116:              resd(j,i) = 0.0
1117:           enddo
1118:        enddo
1119:        first_edge = 1
1120:        last_edge = nedge
1121: #endif
1122: #else
1123:        first_edge = 1
1124:        last_edge = nedge
1125: #endif
1126: !
1127: ! Loop over all the faces and calculate flux
1128: !
1129: !$omp  parallel default(shared) private(n,node1,node2,xmean,ymean,           &
1130: !$omp& zmean,xnorm,ynorm,znorm,rlen,dot,x1,y1,z1,size,x2,y2,z2,rx,           &
1131: !$omp& ry,rz,pL,uL,vL,wL,ubarL,pR,uR,vR,wR,ubarR,p,u,v,w,ubar,phi1,          &
1132: !$omp& phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,c2,c,eig1,eig2,eig3,          &
1133: !$omp& eig4,dp,du,dv,dw,ti11,ti21,ti31,ti41,dv1,dv2,dv3,dv4,r11,r21,         &
1134: !$omp& r31,r41,r12,r22,r32,r42,r13,r23,r33,r43,r14,r24,r34,r44,t1,t2,        &
1135: !$omp& t3,t4,fluxp1,fluxp2,fluxp3,fluxp4,fluxm1,fluxm2,fluxm3,fluxm4,        &
1136: !$omp& res1,res2,res3,res4,c68,c18,inode)                                    &
1137: !$omp& private(x3,y3,z3,p1,p2,p3,ax,ay,az,bx,by,bz,pa,pb,pc,area,c20,        &
1138: !$omp& c0,ubar0,t11,t21,t31,t41,t12,t22,t32,t42,t13,t23,t33,t43,t14,         &
1139: !$omp& t24,t34,t44,pi,ui,vi,wi,unorm,rhs1,rhs2,rhs3,rhs4,ub,vb,wb,node3)     &
1140: !$omp& private(my_thread_id,chunk)                                           &
1141: !$omp& reduction(+:flops)                                                    &
1142: #if defined(_OPENMP)
1143: #if defined(HAVE_EDGE_COLORING)
1144: #elif defined(HAVE_REDUNDANT_WORK)
1145: !$omp& firstprivate(resd)
1146: #else
1147: !$omp& private(first_edge,last_edge)
1148: #endif
1149: #endif
1150: !!$omp1 shared(evec,qvec,resvec,xyz,xyzn,grad,nedge,nnodes,second,flops,
1151: !!$omp2 title,beta,alpha,Re,dt,tot,res0,resc,ntt,mseq,ivisc,irest,icyc,
1152: !!$omp3 ihane,ntturb,gamma,gm1,gp1,gm1g,gp1g,ggm1,p0,rho0,c0,u0,v0,w0,
1153: !!$omp4 et0,h0,pt0,rms,clw,cdw,cmw,xres,cfl1,cfl2,nsmoth,iflim,
1154: !!$omp5 itran,nbtran,jupdate,nstage,ncyct,iramp,nitfo)
1155: !!$      my_thread_id = omp_get_thread_num()
1156: !!$omp   critical
1157: !!$      print *, 'My thread id on processor ',irank,' is ',my_thread_id
1158: !!$omp   end critical
1159: #if defined(_OPENMP)
1160: #if defined(HAVE_EDGE_COLORING)
1161: !$omp do
1162: #elif defined(HAVE_REDUNDANT_WORK)
1163: !$omp do
1164: #else
1165: !     tot_threads = omp_get_num_threads()
1166:       my_thread_id = omp_get_thread_num()
1167: !     chunk = nnodes/tot_threads
1168: !     first_node = my_thread_id*chunk+1
1169: !     if (my_thread_id .eq. (size-1)) then
1170: !        last_node = nnodes
1171: !     else
1172: !        last_node = (my_thread_id+1)*chunk
1173: !     endif
1174:       first_edge = nedge_thr(my_thread_id)
1175:       last_edge = nedge_thr(my_thread_id+1)-1

1177: !      print *, 'On thread ', my_thread_id,', first_edge = ',first_edge,        &
1178: !     &         ', last_edge = ',last_edge

1180: #endif
1181: #else
1182: #endif
1183:       do n = first_edge, last_edge
1184:         node1 = eptr(n,1)
1185:         node2 = eptr(n,2)
1186: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
1187:         if ((part_thr(node1) .eq. my_thread_id)                                 &
1188:      &       .or.(part_thr(node2) .eq. my_thread_id)) then
1189: #endif
1190: !
1191: ! Calculate unit normal to face and length of face
1192: !
1193:           xmean = .5d0*(x(node1) + x(node2))
1194:           ymean = .5d0*(y(node1) + y(node2))
1195:           zmean = .5d0*(z(node1) + z(node2))
1196:           xnorm  = xn(n)
1197:           ynorm  = yn(n)
1198:           znorm  = zn(n)
1199:           rlen   = ra(n)
1200: !
1201: ! Now lets get our other 2 vectors
1202: ! For first vector, use {1,0,0} and subtract off the component
1203: ! in the direction of the face normal. If the inner product of
1204: ! {1,0,0} is close to unity, use {0,1,0}
1205: !
1206:          dot = xnorm
1207:          if (abs(dot).lt.0.95d0)then
1208:           X1 = 1. - dot*xnorm
1209:           Y1 =    - dot*ynorm
1210:           Z1 =    - dot*znorm
1211:          else
1212:           dot = ynorm
1213:           X1 =    - dot*xnorm
1214:           Y1 = 1.0d0 - dot*ynorm
1215:           Z1 =    - dot*znorm
1216:          end if
1217: !
1218: ! Normalize the first vector
1219: !
1220:          size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
1221:          X1 = X1/size
1222:          Y1 = Y1/size
1223:          Z1 = Z1/size
1224: !
1225: ! Take cross-product of normal and V1 to get V2
1226: !
1227:          X2 = ynorm*Z1 - znorm*Y1
1228:          Y2 = znorm*X1 - xnorm*Z1
1229:          Z2 = xnorm*Y1 - ynorm*X1
1230: !
1231: ! Get variables on "left" and "right" side of face
1232: !
1233:          rx     = second*(xmean - x(node1))
1234:          ry     = second*(ymean - y(node1))
1235:          rz     = second*(zmean - z(node1))
1236:          pL = qnod(1,node1) + gradx(1,node1)*rx                           &
1237:      &                       + grady(1,node1)*ry                           &
1238:      &                       + gradz(1,node1)*rz
1239:          uL = qnod(2,node1) + gradx(2,node1)*rx                           &
1240:      &                       + grady(2,node1)*ry                           &
1241:      &                       + gradz(2,node1)*rz
1242:          vL = qnod(3,node1) + gradx(3,node1)*rx                           &
1243:      &                       + grady(3,node1)*ry                           &
1244:      &                       + gradz(3,node1)*rz
1245:          wL = qnod(4,node1) + gradx(4,node1)*rx                           &
1246:      &                       + grady(4,node1)*ry                           &
1247:      &                       + gradz(4,node1)*rz

1249:          ubarL  = xnorm*uL + ynorm*vL + znorm*wL
1250: !        c2L = ubarl*ubarL + beta
1251: !        cL  = sqrt(c2L)
1252: !
1253:          rx     = second*(xmean - x(node2))
1254:          ry     = second*(ymean - y(node2))
1255:          rz     = second*(zmean - z(node2))
1256:          pR   = qnod(1,node2) + gradx(1,node2)*rx                         &
1257:      &                         + grady(1,node2)*ry                         &
1258:      &                         + gradz(1,node2)*rz
1259:          uR   = qnod(2,node2) + gradx(2,node2)*rx                         &
1260:      &                         + grady(2,node2)*ry                         &
1261:      &                         + gradz(2,node2)*rz
1262:          vR   = qnod(3,node2) + gradx(3,node2)*rx                         &
1263:      &                         + grady(3,node2)*ry                         &
1264:      &                         + gradz(3,node2)*rz
1265:          wR   = qnod(4,node2) + gradx(4,node2)*rx                         &
1266:      &                         + grady(4,node2)*ry                         &
1267:      &                         + gradz(4,node2)*rz
1268:          ubarR  = xnorm*uR + ynorm*vR + znorm*wR
1269: !        c2R = ubarR*ubarR + beta
1270: !        cR  = sqrt(c2R)
1271: !
1272: ! Compute averages
1273: !
1274:           p = .5d0*(pL + pR)
1275:           u = .5d0*(uL + uR)
1276:           v = .5d0*(vL + vR)
1277:           w = .5d0*(wL + wR)

1279:           ubar  = xnorm*u + ynorm*v + znorm*w
1280:           phi1  = xnorm*beta + u*ubar
1281:           phi2  = ynorm*beta + v*ubar
1282:           phi3  = znorm*beta + w*ubar
1283:           phi4  = Y2*phi3 - Z2*phi2
1284:           phi5  = Z2*phi1 - X2*phi3
1285:           phi6  = X2*phi2 - Y2*phi1
1286:           phi7  = Z1*phi2 - Y1*phi3
1287:           phi8  = X1*phi3 - Z1*phi1
1288:           phi9  = Y1*phi1 - X1*phi2
1289:           c2    = ubar*ubar + beta
1290:           c     = sqrt(c2)
1291: !
1292: ! Now compute eigenvalues, eigenvectors, and strengths
1293: !
1294:           eig1 = abs(ubar)
1295:           eig2 = abs(ubar)
1296:           eig3 = abs(ubar + c)
1297:           eig4 = abs(ubar - c)
1298: !
1299:           dp = pr - pl
1300:           du = ur - ul
1301:           dv = vr - vl
1302:           dw = wr - wl
1303: !
1304: ! Components of T(inverse) (I will divide by c2 later)
1305: !
1306:           ti11 = -(u*phi4 + v*phi5 + w*phi6)/beta
1307:           ti21 = -(u*phi7 + v*phi8 + w*phi9)/beta
1308:           ti31 =  .5d0*(c-ubar)/beta
1309:           ti41 = -.5d0*(c+ubar)/beta
1310: !
1311: ! jumps (T(inverse)*dq)
1312:           dv1 = (ti11*dp + phi4*du + phi5*dv + phi6*dw)/c2
1313:           dv2 = (ti21*dp + phi7*du + phi8*dv + phi9*dw)/c2
1314:           dv3 = .5d0*(2.d0*ti31*dp + xnorm*du + ynorm*dv + znorm*dw)/c2
1315:           dv4 = .5d0*(2.d0*ti41*dp + xnorm*du + ynorm*dv + znorm*dw)/c2
1316: !
1317: ! Now get elements of T (call it r for now)
1318: !
1319:           r11 = 0.0d0
1320:           r21 = X1
1321:           r31 = Y1
1322:           r41 = Z1

1324:           r12 = 0.0d0
1325:           r22 = X2
1326:           r32 = Y2
1327:           r42 = Z2

1329:           r13 = c*beta
1330:           r23 = xnorm*beta + u*(ubar + c)
1331:           r33 = ynorm*beta + v*(ubar + c)
1332:           r43 = znorm*beta + w*(ubar + c)
1333: !
1334:           r14 = -c*beta
1335:           r24 = xnorm*beta + u*(ubar - c)
1336:           r34 = ynorm*beta + v*(ubar - c)
1337:           r44 = znorm*beta + w*(ubar - c)
1338: !
1339: ! Calculate T* |lambda| *T(inverse)
1340: !
1341:           t1 = eig1*r11*dv1 + eig2*r12*dv2 + eig3*r13*dv3 + eig4*r14*dv4
1342:           t2 = eig1*r21*dv1 + eig2*r22*dv2 + eig3*r23*dv3 + eig4*r24*dv4
1343:           t3 = eig1*r31*dv1 + eig2*r32*dv2 + eig3*r33*dv3 + eig4*r34*dv4
1344:           t4 = eig1*r41*dv1 + eig2*r42*dv2 + eig3*r43*dv3 + eig4*r44*dv4
1345: !
1346: ! Modify to calculate .5(fl +fr) from nodes
1347: ! instead of extrapolated ones
1348: !
1349: !         pL    = qnod(1,node1)
1350: !         uL    = qnod(2,node1)
1351: !         vL    = qnod(3,node1)
1352: !         wL    = qnod(4,node1)
1353: !         ubarL = xnorm*uL + ynorm*vL + znorm*wL
1354: !
1355:           fluxp1 = rlen*beta*ubarL
1356:           fluxp2 = rlen*(uL*ubarL + xnorm*pL)
1357:           fluxp3 = rlen*(vL*ubarL + ynorm*pL)
1358:           fluxp4 = rlen*(wL*ubarL + znorm*pL)
1359: !
1360: ! Now the right side
1361: !
1362: !         pR    = qnod(1,node2)
1363: !         uR    = qnod(2,node2)
1364: !         vR    = qnod(3,node2)
1365: !         wR    = qnod(4,node2)
1366: !         ubarR = xnorm*uR + ynorm*vR + znorm*wR
1367: !
1368:           fluxm1 = rlen*beta*ubarR
1369:           fluxm2 = rlen*(uR*ubarR + xnorm*pR)
1370:           fluxm3 = rlen*(vR*ubarR + ynorm*pR)
1371:           fluxm4 = rlen*(wR*ubarR + znorm*pR)
1372: !
1373:           res1 = 0.5d0*(fluxp1 + fluxm1 - rlen*t1)
1374:           res2 = 0.5d0*(fluxp2 + fluxm2 - rlen*t2)
1375:           res3 = 0.5d0*(fluxp3 + fluxm3 - rlen*t3)
1376:           res4 = 0.5d0*(fluxp4 + fluxm4 - rlen*t4)
1377: !
1378:           flops = flops + 318.0

1380: !
1381: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
1382:           if ((part_thr(node1).eq.my_thread_id)                                 &
1383:      &        .and.(node1.le.nnodes)) then
1384: #else
1385:           if (node1.le.nnodes) then
1386: #endif
1387:            res(1,node1) = res(1,node1) + res1
1388:            res(2,node1) = res(2,node1) + res2
1389:            res(3,node1) = res(3,node1) + res3
1390:            res(4,node1) = res(4,node1) + res4
1391:            flops = flops + 4.0
1392:           endif
1393: !
1394: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
1395:           if ((part_thr(node2).eq.my_thread_id)                                 &
1396:      &        .and.(node2.le.nnodes)) then
1397: #else
1398:           if (node2.le.nnodes) then
1399: #endif
1400:            res(1,node2) = res(1,node2) - res1
1401:            res(2,node2) = res(2,node2) - res2
1402:            res(3,node2) = res(3,node2) - res3
1403:            res(4,node2) = res(4,node2) - res4
1404:            flops = flops + 4.0
1405:           endif
1406: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK)&& !defined(HAVE_EDGE_COLORING)
1407:          endif
1408: #endif
1409: !
1410:       enddo
1411: #if defined(_OPENMP)
1412: #if defined(HAVE_EDGE_COLORING)
1413: !$omp enddo
1414: !$omp end parallel
1415:        first_edge = first_edge + ncount(i)
1416:       enddo
1417: #elif defined(HAVE_REDUNDANT_WORK)
1418: !$omp enddo
1419: !$omp critical (residual_update)
1420:        do n = 1,nnodes
1421:           resvec(1,n)=resvec(1,n)+res(1,n)
1422:           resvec(2,n)=resvec(2,n)+res(2,n)
1423:           resvec(3,n)=resvec(3,n)+res(3,n)
1424:           resvec(4,n)=resvec(4,n)+res(4,n)
1425:        enddo
1426:        flops = flops + 4.0*nnodes
1427: !$omp end critical (residual_update)
1428: !$omp end parallel
1429: #undef res
1430: #define res(a,b) resvec(a,b)
1431: #else
1432: !$omp end parallel
1433: #endif
1434: #else
1435: #endif
1436:       call PetscLogFlops(flops,ierr)
1437:       call PetscLogEventEnd(flux_event,ierr)
1438: !     c68 = 6./8.
1439: !     c18 = 1./8.
1440:       c68 = 0.75d0
1441:       c18 = 0.125d0
1442: !
1443: ! Close contour over the boundaries
1444: ! First do inviscid faces
1445: !
1446:       do n = 1, nnfacet
1447:                node1 = isnode(f2ntn(n,1))
1448:                node2 = isnode(f2ntn(n,2))
1449:                node3 = isnode(f2ntn(n,3))

1451:                x1 = x(node1)
1452:                y1 = y(node1)
1453:                z1 = z(node1)
1454:                p1 = qnod(1,node1)

1456:                x2 = x(node2)
1457:                y2 = y(node2)
1458:                z2 = z(node2)
1459:                p2 = qnod(1,node2)

1461:                x3 = x(node3)
1462:                y3 = y(node3)
1463:                z3 = z(node3)
1464:                p3 = qnod(1,node3)

1466:                ax = x2 - x1
1467:                ay = y2 - y1
1468:                az = z2 - z1

1470:                bx = x3 - x1
1471:                by = y3 - y1
1472:                bz = z3 - z1
1473: !
1474: ! Normal points away from grid interior.
1475: ! Magnitude is 1/3 area of surface triangle.
1476: !
1477:                xnorm =-0.5d0*(ay*bz - az*by)/3.d0
1478:                ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
1479:                znorm =-0.5d0*(ax*by - ay*bx)/3.d0

1481:                pa = c68*p1 + c18*(p2 + p3)
1482:                pb = c68*p2 + c18*(p3 + p1)
1483:                pc = c68*p3 + c18*(p1 + p2)
1484: !
1485:                flops = flops + 35.0
1486:           if (node1.le.nnodes) then
1487:                res(2,node1) = res(2,node1) + xnorm*pa
1488:                res(3,node1) = res(3,node1) + ynorm*pa
1489:                res(4,node1) = res(4,node1) + znorm*pa
1490:                flops = flops + 6.0
1491:           endif

1493:           if (node2.le.nnodes) then
1494:                res(2,node2) = res(2,node2) + xnorm*pb
1495:                res(3,node2) = res(3,node2) + ynorm*pb
1496:                res(4,node2) = res(4,node2) + znorm*pb
1497:                flops = flops + 6.0
1498:           endif

1500:           if (node3.le.nnodes) then
1501:                res(2,node3) = res(2,node3) + xnorm*pc
1502:                res(3,node3) = res(3,node3) + ynorm*pc
1503:                res(4,node3) = res(4,node3) + znorm*pc
1504:                flops = flops + 6.0
1505:           endif

1507:       enddo
1508: !
1509: ! Now viscous faces
1510: !
1511:       do n = 1, nvfacet
1512:                node1 = ivnode(f2ntv(n,1))
1513:                node2 = ivnode(f2ntv(n,2))
1514:                node3 = ivnode(f2ntv(n,3))

1516:                x1 = x(node1)
1517:                y1 = y(node1)
1518:                z1 = z(node1)
1519:                p1 = qnod(1,node1)

1521:                x2 = x(node2)
1522:                y2 = y(node2)
1523:                z2 = z(node2)
1524:                p2 = qnod(1,node2)

1526:                x3 = x(node3)
1527:                y3 = y(node3)
1528:                z3 = z(node3)
1529:                p3 = qnod(1,node3)

1531:                ax = x2 - x1
1532:                ay = y2 - y1
1533:                az = z2 - z1

1535:                bx = x3 - x1
1536:                by = y3 - y1
1537:                bz = z3 - z1
1538: !
1539: ! norm point away from grid interior.
1540: ! norm magnitude is 1/3 area of surface triangle.
1541: !
1542:                xnorm =-0.5d0*(ay*bz - az*by)/3.d0
1543:                ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
1544:                znorm =-0.5d0*(ax*by - ay*bx)/3.d0

1546:                pa = c68*p1 + c18*(p2 + p3)
1547:                pb = c68*p2 + c18*(p3 + p1)
1548:                pc = c68*p3 + c18*(p1 + p2)
1549: !
1550:                flops = flops + 35.0
1551: !
1552:           if (node1.le.nnodes) then
1553:                 res(2,node1) = res(2,node1) + xnorm*pa
1554:                 res(3,node1) = res(3,node1) + ynorm*pa
1555:                 res(4,node1) = res(4,node1) + znorm*pa
1556:                 flops = flops + 6.0

1558:           endif

1560:           if (node2.le.nnodes) then
1561:                 res(2,node2) = res(2,node2) + xnorm*pb
1562:                 res(3,node2) = res(3,node2) + ynorm*pb
1563:                 res(4,node2) = res(4,node2) + znorm*pb
1564:                 flops = flops + 6.0
1565:           endif

1567:           if (node3.le.nnodes) then
1568:                 res(2,node3) = res(2,node3) + xnorm*pc
1569:                 res(3,node3) = res(3,node3) + ynorm*pc
1570:                 res(4,node3) = res(4,node3) + znorm*pc
1571:                 flops = flops + 6.0
1572:           endif

1574:       enddo
1575: !
1576: ! The next section of code is for when you dont care about
1577: ! preserving linear data on boundary. Also, doing this
1578: ! matches the left hand side when not doing Newton-Krylov
1579: ! Usually just go around unless you are just experimenting
1580: !
1581:       goto 1025
1582: !
1583: ! Loop over the boundaries
1584: ! First do inviscid nodes
1585: !
1586: !     do 1020 i = 1,nsnode
1587: !       inode = isnode(i)
1588: !
1589: !       xnorm = sxn(i)
1590: !       ynorm = syn(i)
1591: !       znorm = szn(i)
1592: !
1593: !       p = qnod(1,inode)
1594: !
1595: !       if (inode .le. nnodes) then
1596: !        res(2,inode) = res(2,inode) + xnorm*p
1597: !        res(3,inode) = res(3,inode) + ynorm*p
1598: !        res(4,inode) = res(4,inode) + znorm*p
1599: !       endif
1600: !
1601: !1020 continue
1602: !
1603: ! Now viscous nodes
1604: !
1605: !     do 1030 i = 1,nvnode
1606: !       inode = ivnode(i)
1607: !
1608: !       xnorm = vxn(i)
1609: !       ynorm = vyn(i)
1610: !       znorm = vzn(i)
1611: !
1612: !       p = qnod(1,inode)
1613: !
1614: !       if (inode .le. nnodes) then
1615: !        res(2,inode) = res(2,inode) + xnorm*p
1616: !        res(3,inode) = res(3,inode) + ynorm*p
1617: !        res(4,inode) = res(4,inode) + znorm*p
1618: !       endif
1619: !
1620: !1030 continue
1621: !
1622:  1025 continue
1623: !
1624: ! Now do far-field
1625: !
1626:        do i = 1,nfnode
1627:          inode   = ifnode(i)
1628: !
1629: !/*
1630: ! Get normal and "other" 2 vectors. Remember that fxn,fyn and fzn
1631: ! has the magnitude of the face contained in it.
1632: !*/
1633:          xnorm   = fxn(i)
1634:          ynorm   = fyn(i)
1635:          znorm   = fzn(i)
1636:          area    = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
1637:          xnorm   = xnorm/area
1638:          ynorm   = ynorm/area
1639:          znorm   = znorm/area
1640: !
1641: ! Now lets get our other 2 vectors
1642: ! For first vector, use {1,0,0} and subtract off the component
1643: ! in the direction of the face normal. If the inner product of
1644: ! {1,0,0} is close to unity, use {0,1,0}
1645: !
1646:          dot = xnorm
1647:          if (abs(dot).lt.0.95d0)then
1648:           X1 = 1.d0 - dot*xnorm
1649:           Y1 =    - dot*ynorm
1650:           Z1 =    - dot*znorm
1651:          else
1652:           dot = ynorm
1653:           X1 =    - dot*xnorm
1654:           Y1 = 1.d0 - dot*ynorm
1655:           Z1 =    - dot*znorm
1656:          end if
1657: !
1658: ! Normalize the first vector (V1)
1659: !
1660:          size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
1661:          X1 = X1/size
1662:          Y1 = Y1/size
1663:          Z1 = Z1/size
1664: !
1665: ! Take cross-product of normal with V1 to get V2
1666: !
1667:          X2 = ynorm*Z1 - znorm*Y1
1668:          Y2 = znorm*X1 - xnorm*Z1
1669:          Z2 = xnorm*Y1 - ynorm*X1

1671: !
1672: ! Calculate elements of T and T(inverse) evaluated at freestream
1673: !
1674:          ubar0 = xnorm*u0 + ynorm*v0 + znorm*w0
1675:          c20   = ubar0*ubar0 + beta
1676:          c0    = sqrt(c20)
1677:          phi1  = xnorm*beta + u0*ubar0
1678:          phi2  = ynorm*beta + v0*ubar0
1679:          phi3  = znorm*beta + w0*ubar0
1680:          phi4  = Y2*phi3 - Z2*phi2
1681:          phi5  = Z2*phi1 - X2*phi3
1682:          phi6  = X2*phi2 - Y2*phi1
1683:          phi7  = Z1*phi2 - Y1*phi3
1684:          phi8  = X1*phi3 - Z1*phi1
1685:          phi9  = Y1*phi1 - X1*phi2

1687:          t11 = 0.0d0
1688:          t21 = X1
1689:          t31 = Y1
1690:          t41 = Z1

1692:          t12 = 0.0d0
1693:          t22 = X2
1694:          t32 = Y2
1695:          t42 = Z2

1697:          t13 =  c0*beta
1698:          t23 = xnorm*beta + u0*(ubar0 + c0)
1699:          t33 = ynorm*beta + v0*(ubar0 + c0)
1700:          t43 = znorm*beta + w0*(ubar0 + c0)

1702:          t14 = -c0*beta
1703:          t24 = xnorm*beta + u0*(ubar0 - c0)
1704:          t34 = ynorm*beta + v0*(ubar0 - c0)
1705:          t44 = znorm*beta + w0*(ubar0 - c0)

1707:          ti11 = -(u0*phi4 + v0*phi5 + w0*phi6)/beta
1708:          ti21 = -(u0*phi7 + v0*phi8 + w0*phi9)/beta
1709:          ti31 =  .5d0*(c0 - ubar0)/beta
1710:          ti41 = -.5d0*(c0 + ubar0)/beta
1711: !
1712: ! Now, get the variables on the "inside"
1713: !
1714:          pi      = qnod(1,inode)
1715:          ui      = qnod(2,inode)
1716:          vi      = qnod(3,inode)
1717:          wi      = qnod(4,inode)
1718:          unorm   = xnorm*ui + ynorm*vi + znorm*wi
1719: !
1720: ! If ubar is negative, take the reference condition from outside
1721: !
1722: !
1723:          if (unorm.gt.0.0d0)then
1724:           pr = pi
1725:           ur = ui
1726:           vr = vi
1727:           wr = wi
1728:          else
1729:           pr = p0
1730:           ur = u0
1731:           vr = v0
1732:           wr = w0
1733:          end if
1734: !
1735: ! Set rhs
1736: !
1737:          rhs1 = (ti11*pr + phi4*ur + phi5*vr + phi6*wr)/c20
1738:          rhs2 = (ti21*pr + phi7*ur + phi8*vr + phi9*wr)/c20
1739:          rhs3 = .5d0*(2.d0*ti31*pi + xnorm*ui + ynorm*vi + znorm*wi)/c20
1740:          rhs4 = .5d0*(2.d0*ti41*p0 + xnorm*u0 + ynorm*v0 + znorm*w0)/c20
1741: !
1742: ! Now do matrix multiplication to get values on boundary
1743: !
1744:          pb =                       t13*rhs3 + t14*rhs4
1745:          ub = t21*rhs1 + t22*rhs2 + t23*rhs3 + t24*rhs4
1746:          vb = t31*rhs1 + t32*rhs2 + t33*rhs3 + t34*rhs4
1747:          wb = t41*rhs1 + t42*rhs2 + t43*rhs3 + t44*rhs4

1749:          ubar = xnorm*ub + ynorm*vb + znorm*wb
1750:          flops = flops + 180.0
1751: !
1752:         if (inode.le.nnodes) then
1753:           res(1,inode) = res(1,inode)+area*beta*ubar
1754:           res(2,inode) = res(2,inode)+area*(ub*ubar + xnorm*pb)
1755:           res(3,inode) = res(3,inode)+area*(vb*ubar + ynorm*pb)
1756:           res(4,inode) = res(4,inode)+area*(wb*ubar + znorm*pb)
1757:           flops = flops + 18.0
1758:         endif
1759: !
1760:        enddo
1761:        call PetscLogFlops(flops,ierr)
1762: !
1763: ! End of subroutine FLUX
1764: !
1765:       return
1766:       end

1768: !---------------------------------------------------------------
1769: ! The following subroutines are from node3t.f in the original
1770: ! code - D. K. Kaushik (1/17/97)
1771: !---------------------------------------------------------------
1772: !
1773: !=============================== SUMGS ===============================72
1774: !
1775: ! Gets the weights for calculating gradients using least squares
1776: !
1777: !=====================================================================72
1778:       subroutine SUMGS(nnodes,nedge,evec,xyz,rxy,irank,nvertices)

1780:       implicit none
1781: #include <petsc/finclude/petscsys.h>
1782:       integer nnodes,nedge,irank,nvertices
1783: #if defined(INTERLACING)
1784:       PetscScalar xyz(3,nvertices)
1785:       PetscScalar rxy(7,nnodes)
1786:       integer evec(2,nedge)
1787: #define x(i) xyz(1,i)
1788: #define y(i) xyz(2,i)
1789: #define z(i) xyz(3,i)
1790: #define r11(i) rxy(1,i)
1791: #define r12(i) rxy(2,i)
1792: #define r13(i) rxy(3,i)
1793: #define r22(i) rxy(4,i)
1794: #define r23(i) rxy(5,i)
1795: #define r33(i) rxy(6,i)
1796: #define r44(i) rxy(7,i)
1797: #undef eptr
1798: #define eptr(j,i) evec(i,j)
1799: #else
1800:       PetscScalar xyz(nvertices,3)
1801:       PetscScalar rxy(nnodes,7)
1802:       integer evec(nedge,2)
1803: #define x(i) xyz(i,1)
1804: #define y(i) xyz(i,2)
1805: #define z(i) xyz(i,3)
1806: #define r11(i) rxy(i,1)
1807: #define r12(i) rxy(i,2)
1808: #define r13(i) rxy(i,3)
1809: #define r22(i) rxy(i,4)
1810: #define r23(i) rxy(i,5)
1811: #define r33(i) rxy(i,6)
1812: #define r44(i) rxy(i,7)
1813: #define eptr(i,j) evec(i,j)
1814: #endif
1815:       integer i,n,node1,node2,ierr
1816:       PetscScalar x1,y1,z1,x2,y2,z2,dx,dy,dz
1817:       PetscScalar dist,weight,w2,w11,w22,w33
1818:       PetscScalar r12r11,r13r11,r23r22,rmess

1820: !
1821: ! Initialize all the rij to 0.0
1822: !
1823:       do 1000 i = 1,nnodes
1824:          r11(i) = 0.0
1825:          r12(i) = 0.0
1826:          r13(i) = 0.0
1827:          r22(i) = 0.0
1828:          r23(i) = 0.0
1829:          r33(i) = 0.0
1830:          r44(i) = 0.0
1831:  1000 continue
1832: !
1833: ! Now loop over the edges and accumulate the r
1834: !
1835:       do 1020 n = 1, nedge
1836: !
1837:           node1 = eptr(n,1)
1838:           node2 = eptr(n,2)
1839: !
1840:           x1 = x(node1)
1841:           y1 = y(node1)
1842:           z1 = z(node1)
1843:           x2 = x(node2)
1844:           y2 = y(node2)
1845:           z2 = z(node2)
1846: !
1847:           dx = x2 - x1
1848:           dy = y2 - y1
1849:           dz = z2 - z1
1850:           dist = sqrt(dx*dx + dy*dy + dz*dz)
1851: !         weight = 1.0d0/dist
1852:           weight = 1.0d0
1853:           w2 = weight*weight
1854: !
1855:           if (node1 .le. nnodes) then
1856:            r11(node1) = r11(node1) + (x2 - x1)*(x2 - x1)*w2
1857:            r12(node1) = r12(node1) + (x2 - x1)*(y2 - y1)*w2
1858:            r13(node1) = r13(node1) + (x2 - x1)*(z2 - z1)*w2
1859:           endif
1860: !
1861:           if (node2 .le. nnodes) then
1862:            r11(node2) = r11(node2) + (x1 - x2)*(x1 - x2)*w2
1863:            r12(node2) = r12(node2) + (x1 - x2)*(y1 - y2)*w2
1864:            r13(node2) = r13(node2) + (x1 - x2)*(z1 - z2)*w2
1865:           endif
1866:  1020 continue
1867: !
1868: !/*
1869: ! Now calculate ||x(:)|| = r11 by taking the square root
1870: ! Also divide r12 and r13 by ||x(:)||
1871: !*/
1872:       do 1030 i = 1,nnodes
1873:         r11(i) = sqrt(r11(i))
1874:         r12(i) = r12(i)/r11(i)
1875:         r13(i) = r13(i)/r11(i)
1876:  1030 continue
1877: !
1878: !/*
1879: ! Now calculate r22 and r23
1880: !*/
1881:       do 1050 n = 1, nedge
1882: !
1883:           node1 = eptr(n,1)
1884:           node2 = eptr(n,2)
1885: !
1886:           x1 = x(node1)
1887:           y1 = y(node1)
1888:           z1 = z(node1)
1889:           x2 = x(node2)
1890:           y2 = y(node2)
1891:           z2 = z(node2)
1892: !
1893:           dx = x2 - x1
1894:           dy = y2 - y1
1895:           dz = z2 - z1
1896:           dist = sqrt(dx*dx + dy*dy + dz*dz)
1897: !         weight = 1.0/dist
1898:           weight = 1.0d0
1899:           dx = weight*dx
1900:           dy = weight*dy
1901:           dz = weight*dz
1902:           w2 = weight*weight
1903: !
1904:           if (node1 .le. nnodes) then
1905:            r22(node1) = r22(node1) +                                         &
1906:      &                  (dy-dx*r12(node1)/r11(node1))**2
1907:            r23(node1) = r23(node1) + dz*                                     &
1908:      &                  (dy-dx*r12(node1)/r11(node1))
1909:           endif
1910: !
1911:           if (node2 .le. nnodes) then
1912:            r22(node2) = r22(node2) +                                         &
1913:      &                  (-dy + dx*r12(node2)/r11(node2))**2
1914:            r23(node2) = r23(node2) - dz*                                     &
1915:      &                  (-dy+dx*r12(node2)/r11(node2))
1916:           endif
1917:  1050 continue
1918: !
1919: !/*
1920: ! Now finish getting r22 and r23
1921: !*/
1922:       do 1060 i = 1,nnodes
1923:         r22(i) = sqrt(r22(i))
1924:         r23(i) = r23(i)/r22(i)
1925:  1060 continue
1926: !
1927: !/*
1928: ! Now all we have to do is get r33
1929: !*/
1930:       do 1080 n = 1, nedge
1931: !
1932:           node1 = eptr(n,1)
1933:           node2 = eptr(n,2)
1934: !
1935:           x1 = x(node1)
1936:           y1 = y(node1)
1937:           z1 = z(node1)
1938:           x2 = x(node2)
1939:           y2 = y(node2)
1940:           z2 = z(node2)
1941: !
1942:           dx = x2 - x1
1943:           dy = y2 - y1
1944:           dz = z2 - z1
1945:           dist = sqrt(dx*dx + dy*dy + dz*dz)
1946: !         weight = 1.0/dist
1947:           weight = 1.0d0
1948:           dx = weight*dx
1949:           dy = weight*dy
1950:           dz = weight*dz
1951:           w2 = weight*weight
1952: !
1953:           if (node1 .le. nnodes) then
1954:            r33(node1) = r33(node1) +                                        &
1955:      &                  (dz-dx*r13(node1)/r11(node1)-                       &
1956:      &                  r23(node1)/r22(node1)*                              &
1957:      &                  (dy - dx*r12(node1)/                                &
1958:      &                   r11(node1)))**2
1959:           endif
1960: !
1961:           if (node2 .le. nnodes) then
1962:            r33(node2) = r33(node2) +                                        &
1963:      &                  (-dz+dx*r13(node2)/r11(node2)-                      &
1964:      &                  r23(node2)/r22(node2)*                              &
1965:      &                  (-dy + dx*r12(node2)/                               &
1966:      &                  r11(node2)))**2
1967:           endif
1968: !
1969:  1080 continue
1970: !
1971: !/*
1972: ! Now just get the magnitude of r33
1973: !*/
1974:       do 1090 i = 1,nnodes
1975:         r33(i) = sqrt(r33(i))
1976:  1090 continue
1977: !
1978: !/*
1979: ! Added by Dinesh Kaushik (6/27/97)
1980: ! The following addition changes the meaning of r11 .. r33. r44
1981: ! is the new parameter introduced by me. The new definitions
1982: ! are taken from LSTGS (where these parameters
1983: ! are used finally).
1984: ! r11->w11
1985: ! r22->w22
1986: ! r33->w33
1987: ! r12->r12r11
1988: ! r13->r13r11
1989: ! r23->r23r22
1990: ! r44->rmess */

1992:       do i = 1,nnodes
1993:            w11 = 1.d0/(r11(i)*r11(i))
1994:            w22 = 1.d0/(r22(i)*r22(i))
1995:            w33 = 1.d0/(r33(i)*r33(i))
1996:            r12r11 = r12(i)/r11(i)
1997:            r13r11 = r13(i)/r11(i)
1998:            r23r22 = r23(i)/r22(i)
1999:            rmess=(r12(i)*r23(i)-                                             &
2000:      &            r13(i)*r22(i))/                                            &
2001:      &           (r11(i)*r22(i)*                                             &
2002:      &            r33(i)*r33(i))

2004:            r11(i) = w11
2005:            r22(i) = w22
2006:            r33(i) = w33
2007:            r12(i) = r12r11
2008:            r13(i) = r13r11
2009:            r23(i) = r23r22
2010:            r44(i) = rmess
2011:       enddo
2012: !
2013: ! Finished with SUMGS
2014: !
2015:       return
2016:       end

2018: !================================= LSTGS =============================72
2019: !
2020: ! Calculates the Gradients at the nodes using weighted least squares
2021: ! This subroutine solves using Gram-Schmidt
2022: !
2023: !=====================================================================72
2024:       subroutine LSTGS(nnodes,nedge,evec,                                &
2025:      &                 qvec,grad,xyz,                                    &
2026:      &                 rxy,irank,nvertices)
2027: !
2028:       implicit none
2029: #include <petsc/finclude/petscsys.h>
2030:       integer nnodes,nedge,irank,nvertices
2031: !
2032: #if defined(INTERLACING)
2033:       PetscScalar qvec(4,nvertices)
2034:       PetscScalar xyz(3,nvertices)
2035:       PetscScalar rxy(7,nnodes)
2036:       PetscScalar grad(3,4,nnodes)
2037:       integer evec(2,nedge)
2038: #define qnod(i,j) qvec(i,j)
2039: #define gradx(i,j) grad(1,i,j)
2040: #define grady(i,j) grad(2,i,j)
2041: #define gradz(i,j) grad(3,i,j)
2042: #define x(i) xyz(1,i)
2043: #define y(i) xyz(2,i)
2044: #define z(i) xyz(3,i)
2045: #define r11(i) rxy(1,i)
2046: #define r12(i) rxy(2,i)
2047: #define r13(i) rxy(3,i)
2048: #define r22(i) rxy(4,i)
2049: #define r23(i) rxy(5,i)
2050: #define r33(i) rxy(6,i)
2051: #define r44(i) rxy(7,i)
2052: #define eptr(j,i) evec(i,j)
2053: #else
2054:       PetscScalar qvec(nvertices,4)
2055:       PetscScalar xyz(nvertices,3)
2056:       PetscScalar rxy(nnodes,7)
2057:       PetscScalar grad(nnodes,4,3)
2058:       integer evec(nedge,2)
2059: #define qnod(i,j) qvec(j,i)
2060: #define gradx(i,j) grad(j,i,1)
2061: #define grady(i,j) grad(j,i,2)
2062: #define gradz(i,j) grad(j,i,3)
2063: #define x(i) xyz(i,1)
2064: #define y(i) xyz(i,2)
2065: #define z(i) xyz(i,3)
2066: #define r11(i) rxy(i,1)
2067: #define r12(i) rxy(i,2)
2068: #define r13(i) rxy(i,3)
2069: #define r22(i) rxy(i,4)
2070: #define r23(i) rxy(i,5)
2071: #define r33(i) rxy(i,6)
2072: #define r44(i) rxy(i,7)
2073: #define eptr(i,j) evec(i,j)
2074: #endif
2075:       PetscScalar  title(20),beta,alpha
2076:       PetscScalar  Re,dt,tot,res0,resc
2077:       integer i,ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2078:       PetscScalar  cfl1,cfl2
2079:       integer nsmoth,iflim,itran,nbtran,jupdate,                          &
2080:      &        nstage,ncyct,iramp,nitfo
2081:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                   &
2082:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2083:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,           &
2084:      &             nstage,ncyct,iramp,nitfo
2085:       integer n,node1,node2,ierr
2086:       PetscLogDouble flops
2087:       PetscScalar  dx1,dy1,dz1,dx2,dy2,dz2,                                    &
2088:      &        dq1,dq2,dq3,dq4,                                            &
2089:      &        weight,w11,r12r11,r13r11,w22,r23r22,w33,rmess,              &
2090:      &        coef1,coef2,termx,termy,termz
2091: !     logging variables
2092: !     integer flag
2093: !     integer grad_event,node1_event,node2_event
2094: !     character * 16 grad_label, node1_label, node2_label
2095: !     data flag/-1/,grad_label/'GRAD            '/
2096: !     data node1_label/'NODE1           '/
2097: !     data node2_label/'NODE2           '/
2098: !     save grad_event, grad_label, flag
2099: !     save node1_event,node2_event,node1_label, node2_label

2101: !     if (flag .eq. -1) then
2102: !        call PetscLogEventRegister(grad_event,grad_label,ierr)
2103: !        call PetscLogEventRegister(node1_event,node1_label,ierr)
2104: !        call PetscLogEventRegister(node2_event,node2_label,ierr)
2105: !        flag = 1
2106: !     endif
2107: !     call PetscLogEventBegin(grad_event,0,0,0,0,ierr)

2109:       flops = 0.0
2110: ! For checking out the code input a linear distribution
2111: !
2112: !     write(6,700)nnodes,nedge
2113: ! 700 format(1h ,'nnodes=',i5,' nedge=',i5)
2114: !     do 1001 i = 1,nnodes
2115: !     write(6,800)i,x(i),y(i),z(i)
2116: ! 800 format(1h ,'i x y z=',i5,3(f10.5,1x))
2117: !       qnod(1,i) = 1.0*x(i) +  2.0*y(i) + 3.0*z(i)
2118: !       qnod(2,i) = 3.0*x(i) +  4.0*y(i) + 6.0*z(i)
2119: !       qnod(3,i) = 5.0*x(i) +  6.0*y(i) + 9.0*z(i)
2120: !       qnod(4,i) = 7.0*x(i) +  8.0*y(i) + 12.0*z(i)
2121: !       qnod(i,5) = 9.0*x(i) + 10.0*y(i) + 15.0*z(i)
2122: !1001 continue
2123: !
2124: ! Zero out the gradients
2125: !
2126:       do 1000 i = 1,nnodes
2127:          gradx(1,i) = 0.0
2128:          grady(1,i) = 0.0
2129:          gradz(1,i) = 0.0

2131:          gradx(2,i) = 0.0
2132:          grady(2,i) = 0.0
2133:          gradz(2,i) = 0.0

2135:          gradx(3,i) = 0.0
2136:          grady(3,i) = 0.0
2137:          gradz(3,i) = 0.0

2139:          gradx(4,i) = 0.0
2140:          grady(4,i) = 0.0
2141:          gradz(4,i) = 0.0
2142:  1000 continue

2144: ! If second order, loop over all the faces accumulate sums
2145: !
2146: !     nitfo = 0
2147: !     if (ntt.gt.nitfo.or.ivisc.gt.0)then
2148: !     if (1 .gt. 0) then
2149:        do 1020 n = 1, nedge
2150:          node1 = eptr(n,1)
2151:          node2 = eptr(n,2)
2152: !        if ((node1 .le. nnodes).or.(node2 .le. nnodes)) then
2153:            dx1 = x(node2) - x(node1)
2154:            dy1 = y(node2) - y(node1)
2155:            dz1 = z(node2) - z(node1)
2156: !
2157: !          dist = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1)
2158: !          weight = 1.0/dist
2159:            weight = 1.0d0
2160: !          w2 = weight*weight
2161: !
2162:            dx1 = weight*dx1
2163:            dy1 = weight*dy1
2164:            dz1 = weight*dz1

2166:            flops = flops + 6.0
2167: !
2168: !         call PetscLogEventBegin(node1_event,0,0,0,0,ierr)
2169:           if (node1 .le. nnodes) then
2170:            dq1 = weight*(qnod(1,node2) - qnod(1,node1))
2171:            dq2 = weight*(qnod(2,node2) - qnod(2,node1))
2172:            dq3 = weight*(qnod(3,node2) - qnod(3,node1))
2173:            dq4 = weight*(qnod(4,node2) - qnod(4,node1))
2174: !
2175: !          w11 = 1./(r11(node1)*r11(node1))
2176: !          w22 = 1./(r22(node1)*r22(node1))
2177: !          w33 = 1./(r33(node1)*r33(node1))
2178: !          r12r11 = r12(node1)/r11(node1)
2179: !          r13r11 = r13(node1)/r11(node1)
2180: !          r23r22 = r23(node1)/r22(node1)
2181: !          rmess  = (r12(node1)*r23(node1) - r13(node1)*r22(node1))/
2182: !    1              (r11(node1)*r22(node1)*r33(node1)*r33(node1))
2183: !
2184:            w11 = r11(node1)
2185:            r12r11 = r12(node1)
2186:            r13r11 = r13(node1)
2187:            w22 = r22(node1)
2188:            r23r22 = r23(node1)
2189:            w33 = r33(node1)
2190:            rmess  = r44(node1)
2191: !
2192:            coef1  = dy1 - dx1*r12r11
2193:            coef2  = dz1 - dx1*r13r11 - r23r22*coef1
2194:            termx = dx1*w11 - w22*r12r11*coef1 + rmess*coef2
2195:            termy = w22*coef1 - r23r22*w33*coef2
2196:            termz = w33*coef2
2197: !
2198:            gradx(1,node1) = gradx(1,node1) + termx*dq1
2199:            grady(1,node1) = grady(1,node1) + termy*dq1
2200:            gradz(1,node1) = gradz(1,node1) + termz*dq1
2201: !
2202:            gradx(2,node1) = gradx(2,node1) + termx*dq2
2203:            grady(2,node1) = grady(2,node1) + termy*dq2
2204:            gradz(2,node1) = gradz(2,node1) + termz*dq2
2205: !
2206:            gradx(3,node1) = gradx(3,node1) + termx*dq3
2207:            grady(3,node1) = grady(3,node1) + termy*dq3
2208:            gradz(3,node1) = gradz(3,node1) + termz*dq3
2209: !
2210:            gradx(4,node1) = gradx(4,node1) + termx*dq4
2211:            grady(4,node1) = grady(4,node1) + termy*dq4
2212:            gradz(4,node1) = gradz(4,node1) + termz*dq4
2213: !
2214:            flops = flops + 49.0
2215:           endif
2216: !         call PetscLogEventEnd(node1_event,0,0,0,0,ierr)
2217: !
2218: ! Now do the other node
2219: !
2220: !         call PetscLogEventBegin(node2_event,0,0,0,0,ierr)
2221:           if (node2 .le. nnodes) then
2222:            dx2 = -dx1
2223:            dy2 = -dy1
2224:            dz2 = -dz1
2225: !
2226:            dq1 = weight*(qnod(1,node1) - qnod(1,node2))
2227:            dq2 = weight*(qnod(2,node1) - qnod(2,node2))
2228:            dq3 = weight*(qnod(3,node1) - qnod(3,node2))
2229:            dq4 = weight*(qnod(4,node1) - qnod(4,node2))
2230: !
2231: !          w11 = 1./(r11(node2)*r11(node2))
2232: !          w22 = 1./(r22(node2)*r22(node2))
2233: !          w33 = 1./(r33(node2)*r33(node2))
2234: !          r12r11 = r12(node2)/r11(node2)
2235: !          r13r11 = r13(node2)/r11(node2)
2236: !          r23r22 = r23(node2)/r22(node2)
2237: !          rmess  = (r12(node2)*r23(node2) - r13(node2)*r22(node2))/
2238: !    &              (r11(node2)*r22(node2)*r33(node2)*r33(node2))

2240:            w11 = r11(node2)
2241:            r12r11 = r12(node2)
2242:            r13r11 = r13(node2)
2243:            w22 = r22(node2)
2244:            r23r22 = r23(node2)
2245:            w33 = r33(node2)
2246:            rmess  = r44(node2)
2247: !
2248:            coef1  = dy2 - dx2*r12r11
2249:            coef2  = dz2 - dx2*r13r11 - r23r22*coef1
2250:            termx = dx2*w11 - w22*r12r11*coef1 + rmess*coef2
2251:            termy = w22*coef1 - r23r22*w33*coef2
2252:            termz = w33*coef2
2253: !
2254:            gradx(1,node2) = gradx(1,node2) + termx*dq1
2255:            grady(1,node2) = grady(1,node2) + termy*dq1
2256:            gradz(1,node2) = gradz(1,node2) + termz*dq1
2257: !
2258:            gradx(2,node2) = gradx(2,node2) + termx*dq2
2259:            grady(2,node2) = grady(2,node2) + termy*dq2
2260:            gradz(2,node2) = gradz(2,node2) + termz*dq2
2261: !
2262:            gradx(3,node2) = gradx(3,node2) + termx*dq3
2263:            grady(3,node2) = grady(3,node2) + termy*dq3
2264:            gradz(3,node2) = gradz(3,node2) + termz*dq3
2265: !
2266:            gradx(4,node2) = gradx(4,node2) + termx*dq4
2267:            grady(4,node2) = grady(4,node2) + termy*dq4
2268:            gradz(4,node2) = gradz(4,node2) + termz*dq4
2269: !
2270:            flops = flops + 52.0
2271:           endif
2272: !         call PetscLogEventEnd(node2_event,0,0,0,0,ierr)
2273: !       endif
2274: !
2275:  1020  continue
2276: !     end if
2277:       call PetscLogFlops(flops,ierr)
2278: !     call PetscLogEventEnd(grad_event,0,0,0,0,ierr)
2279: !
2280: ! End of LSTGS
2281: !
2282:       return
2283:       end

2285: !=================================== GETRES ==========================72
2286: !
2287: ! Calculates the residual
2288: ! Last Modified - D. K. Kaushik 1/23/97
2289: ! I have eliminated the input variables which were not needed -
2290: ! dq, A, B, iupdate
2291: !
2292: !=====================================================================72
2293: #if defined(_OPENMP)
2294: #if defined(HAVE_EDGE_COLORING)
2295:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
2296:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
2297:      &                  ileast,isnode,ivnode,ifnode,                        &
2298:      &                  nnfacet,f2ntn,nnbound,                              &
2299:      &                  nvfacet,f2ntv,nvbound,                              &
2300:      &                  nffacet,f2ntf,nfbound,                              &
2301:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
2302:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
2303:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
2304:      &                  max_threads,                                        &
2305:      &                  ncolor,ncount,                                      &
2306:      &                  LocalTS,irank, nvertices)
2307: #elif defined(HAVE_REDUNDANT_WORK)
2308:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
2309:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
2310:      &                  ileast,isnode,ivnode,ifnode,                        &
2311:      &                  nnfacet,f2ntn,nnbound,                              &
2312:      &                  nvfacet,f2ntv,nvbound,                              &
2313:      &                  nffacet,f2ntf,nfbound,                              &
2314:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
2315:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
2316:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
2317:      &                  max_threads,                                        &
2318:      &                  resd,                                               &
2319:      &                  LocalTS,irank, nvertices)
2320: #else
2321:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
2322:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
2323:      &                  ileast,isnode,ivnode,ifnode,                        &
2324:      &                  nnfacet,f2ntn,nnbound,                              &
2325:      &                  nvfacet,f2ntv,nvbound,                              &
2326:      &                  nffacet,f2ntf,nfbound,                              &
2327:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
2328:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
2329:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
2330:      &                  max_threads,                                        &
2331:      &                  nedgeAllThr,part_thr,nedge_thr,                     &
2332:      &                  edge_thr,xyzn_thr,                                  &
2333:      &                  LocalTS,irank, nvertices)
2334: #endif
2335: #else
2336:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
2337:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
2338:      &                  ileast,isnode,ivnode,ifnode,                        &
2339:      &                  nnfacet,f2ntn,nnbound,                              &
2340:      &                  nvfacet,f2ntv,nvbound,                              &
2341:      &                  nffacet,f2ntf,nfbound,                              &
2342:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
2343:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
2344:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
2345:      &                  LocalTS,irank, nvertices)
2346: #endif
2347: !
2348:       implicit none
2349: #include <petsc/finclude/petscsys.h>
2350: !
2351:       integer nnodes, ncell, nedge,                                         &
2352:      &        nbface,ileast,ires,                                           &
2353:      &        nsface, nvface, nfface,                                       &
2354:      &        nsnode, nvnode, nfnode,                                       &
2355:      &        nnfacet,nvfacet,nffacet,                                      &
2356:      &        nnbound,nvbound,nfbound,                                      &
2357:      &        LocalTS,irank,nvertices
2358:       integer evec(nedge,2)
2359:       integer isface(1),ivface(1),ifface(1)
2360:       integer isnode(1),ivnode(1),ifnode(1)
2361:       integer c2n(ncell,4),c2e(ncell,6)
2362:       integer f2ntn(nnfacet,4)
2363:       integer f2ntv(nvfacet,4)
2364:       integer f2ntf(nffacet,4)
2365: !
2366:       PetscScalar us(nbface,3,4),vs(nbface,3,4)
2367:       PetscScalar as(nbface,3,4)
2368:       PetscScalar sxn(1),syn(1),szn(1)
2369:       PetscScalar vxn(1),vyn(1),vzn(1)
2370:       PetscScalar fxn(1),fyn(1),fzn(1)
2371:       PetscScalar xyz(nvertices,3),area(nvertices)
2372:       PetscScalar xyzn(nedge,4)
2373:       PetscScalar turbre(1),slen(1)
2374:       PetscScalar qvec(nvertices,4)
2375:       PetscScalar phi(nvertices,4)
2376:       PetscScalar cdt(nvertices)
2377:       PetscScalar grad(3,4,nvertices)
2378: !     real dq(nnodes,4)
2379: !     real r11(nvertices),r12(nvertices),r13(nvertices)
2380: !     real r22(nvertices),r23(nvertices),r33(nvertices)
2381:       PetscScalar amut(nnodes)
2382: !
2383:       integer i
2384:       PetscScalar  title(20),beta,alpha
2385:       PetscScalar  Re,dt,tot,res0,resc
2386:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2387:       PetscScalar  cfl1,cfl2
2388:       integer nsmoth,iflim,itran,nbtran,jupdate,                             &
2389:      &        nstage,ncyct,iramp,nitfo
2390:       PetscScalar  gtol
2391:       integer icycle,nsrch,ilu0,ifcn
2392:       PetscScalar  rms(1001),clw(1001),cdw(1001),                                 &
2393:      &        cmw(1001),xres(1001)
2394:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                      &
2395:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2396:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,              &
2397:      &             nstage,ncyct,iramp,nitfo
2398:       common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
2399:       common/history/rms,clw,cdw,cmw,xres
2400: !     logging variables
2401: !     integer flux_event, flag, ierr
2402: !     integer delta2_event
2403: !     character * 16 flux_label
2404: !     character * 16 delta2_label
2405: !     data flag/-1/,flux_label/'FluxEvaluation  '/
2406: !     data delta2_label/'DELTA2          '/
2407: !     save flux_event, flag, flux_label
2408: !     save delta2_event,delta2_label
2409: #if defined(INTERLACING)
2410:        PetscScalar resvec(4,nnodes)
2411: #undef res
2412: #define res(i,j) resvec(i,j)
2413: #else
2414:        PetscScalar resvec(nnodes,4)
2415: #define res(i,j) resvec(j,i)
2416: #endif
2417: #if defined(_OPENMP)
2418:        integer max_threads
2419: #if defined(HAVE_EDGE_COLORING)
2420:        integer ncolor,ncount(1)
2421: #elif defined(HAVE_REDUNDANT_WORK)
2422:        PetscScalar resd(4,nnodes)
2423: #else
2424:        integer nedgeAllThr,part_thr(nnodes),nedge_thr(0:max_threads),         &
2425:      &         edge_thr(2,nedgeAllThr)
2426:        PetscScalar xyzn_thr(4,nedgeAllThr)
2427: #endif
2428: #endif
2429: !

2431: !     if (flag .eq. -1) then
2432: !        call PetscLogEventRegister(delta2_event,delta2_label,ierr)
2433: !        call PetscLogEventRegister(flux_event,flux_label,ierr)
2434: !        flag = 1
2435: !     endif
2436: !
2437: ! Calculate the time step
2438: !
2439:       if (ires.eq.1) goto 888
2440: !     call PetscLogEventBegin(delta2_event,0,0,0,0,ierr)
2441:       call DELTAT2(nnodes,nedge,qvec,cdt,                                   &
2442:      &             xyz,area,xyzn,evec,                                      &
2443:      &             sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,                     &
2444:      &             nsnode,nvnode,nfnode,isnode,ivnode,ifnode,               &
2445:      &             LocalTS,irank,nvertices)
2446: !     call PetscLogEventEnd(delta2_event,0,0,0,0,ierr)
2447:   888 continue
2448: !
2449: !/*
2450: !   Calculate the gradients
2451: !   ----Kyle seems to recommend only LSTGS for gradients,
2452: !   so I have commented the GETGRAD call - DKK (1/17/97)
2453: !
2454: !     if (ileast.eq.0) then
2455: !        call GETGRAD(nnodes,ncell,nedge,nsface,nvface,nfface,
2456: !    &                isface,ivface,ifface,eptr,ncolor,ncount,
2457: !    &                qnod,gradx,grady,x,y,
2458: !    &                area,wx,wy,xn,yn,rl)
2459: !
2460: !     else if (ileast.eq.4) then
2461: !     if (ileast.eq.4) then
2462: !        call LSTGS(nnodes,nedge,eptr,
2463: !    &              qnod,gradx,grady,gradz,xyz,
2464: !    &              r11,r12,r13,r22,r23,r33,irank,nvertices)
2465: !     end if
2466: !
2467: ! zero out residuals (viscous residuals are zeroed in vfluxnew)
2468: !*/
2469:       do 1002 i = 1,nnodes
2470:           res(1,i)=0.0d0
2471:           res(2,i)=0.0d0
2472:           res(3,i)=0.0d0
2473:           res(4,i)=0.0d0

2475:           phi(i,1)=1.0d0
2476:           phi(i,2)=1.0d0
2477:           phi(i,3)=1.0d0
2478:           phi(i,4)=1.0d0
2479:  1002   continue
2480: !
2481: !/*
2482: !
2483: ! If not doing Newton-Krylov and iflim=1 call the Flux Limiter
2484: !
2485: !     if (iflim.eq.1.and.ifcn.ne.1)then
2486: !      call TIMLIM(nnodes,nedge,qnod,res,dq,phi,ncolor,ncount,
2487: !    &             gradx,grady,gradz,x,y,z,eptr)
2488: !
2489: ! If we used the limiter we need to zero out the residual again
2490: ! since we used it for scratch space
2491: !
2492: !     do 1003 i = 1,nnodes
2493: !         res(1,i)=0.0
2494: !         res(2,i)=0.0
2495: !         res(3,i)=0.0
2496: !         res(4,i)=0.0
2497: !1003   continue
2498: !
2499: !     end if
2500: !*/
2501: !   Split the fluxes and perform the flux balance
2502: !
2503: !       call PetscLogEventBegin(flux_event,0,0,0,0,ierr)

2505: #if defined(_OPENMP)
2506: #if defined(HAVE_EDGE_COLORING)
2507:       call FLUX(nnodes, ncell, nedge,                                     &
2508:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
2509:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
2510:      &          nnfacet,f2ntn,nnbound,                                    &
2511:      &          nvfacet,f2ntv,nvbound,                                    &
2512:      &          nffacet,f2ntf,nfbound,                                    &
2513:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
2514:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
2515:      &          max_threads,                                              &
2516:      &          ncolor,ncount,                                            &
2517:      &          irank,nvertices)
2518: #elif defined(HAVE_REDUNDANT_WORK)
2519:       call FLUX(nnodes, ncell, nedge,                                     &
2520:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
2521:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
2522:      &          nnfacet,f2ntn,nnbound,                                    &
2523:      &          nvfacet,f2ntv,nvbound,                                    &
2524:      &          nffacet,f2ntf,nfbound,                                    &
2525:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
2526:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
2527:      &          max_threads,                                              &
2528:      &          resd,                                                     &
2529:      &          irank,nvertices)
2530: #else
2531:       call FLUX(nnodes, ncell, nedge,                                     &
2532:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
2533:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
2534:      &          nnfacet,f2ntn,nnbound,                                    &
2535:      &          nvfacet,f2ntv,nvbound,                                    &
2536:      &          nffacet,f2ntf,nfbound,                                    &
2537:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
2538:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
2539:      &          max_threads,                                              &
2540:      &          nedgeAllThr,part_thr,nedge_thr,edge_thr,xyzn_thr,         &
2541:      &          irank,nvertices)
2542: #endif
2543: #else
2544:       call FLUX(nnodes, ncell, nedge,                                     &
2545:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
2546:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
2547:      &          nnfacet,f2ntn,nnbound,                                    &
2548:      &          nvfacet,f2ntv,nvbound,                                    &
2549:      &          nffacet,f2ntf,nfbound,                                    &
2550:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
2551:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
2552:      &          irank,nvertices)
2553: #endif
2554: !          call PetscLogEventEnd(flux_event,0,0,0,0,ierr)
2555: !/*
2556: ! calculate viscous fluxes
2557: !
2558: !     if (ivisc.gt.0) then
2559: !           call VISRHS (nnodes,ncell,nedge,
2560: !           call EDGEVIS(nnodes,ncell,nedge,
2561: !    &                  nsnode,nvnode,nfnode,isnode,ivnode,ifnode,
2562: !    &                  nsface,nvface,nfface,isface,ivface,ifface,
2563: !    &                  nnfacet,f2ntn,nnbound,ncolorn,countn,
2564: !    &                  nvfacet,f2ntv,nvbound,ncolorv,countv,
2565: !    &                  nffacet,f2ntf,nfbound,ncolorf,countf,
2566: !    &                  nccolor,nccount,
2567: !    &                  eptr,c2n,c2e,
2568: !    &                  sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,
2569: !    &                  x,y,z,gradx,grady,gradz,
2570: !    &                  qnod,amut,res,phi)
2571: !     end if
2572: !*/
2573: ! End of subroutine GETRES
2574: !
2575:       return
2576:       end

2578: !---------------------------------------------------------------
2579: ! The following subroutine is from node4t.f in the original
2580: ! code - D. K. Kaushik (1/17/97)
2581: !---------------------------------------------------------------
2582: !
2583: !=============================================================================
2584: !
2585: !  Opens files for I/O
2586: !
2587: !=============================================================================
2588:       SUBROUTINE OPENM(irank)
2589: !
2590: !  TAPE7  -- input:  mach number, angle of attack etc..
2591: !  TAPE9  -- input:  reads restart file
2592: !  TAPE10 -- output: residual history
2593: !  TAPE11 -- output: writes restart file
2594: !  TAPE12 -- output: writes residual and lift for plotting
2595: !  TAPE13 -- output: writes flowfield for contour plotting
2596: !
2597:       implicit none
2598:       integer  irank
2599: !     OPEN(UNIT= 7,FILE='home/petsc/datafiles/fun3dgrid/ginput.faces',            &
2600: !    &form='formatted',STATUS='OLD')

2602: !     OPEN(UNIT=9,FILE='framer.bin',
2603: !    &form='unformatted',STATUS='old')

2605:       if (irank .eq. 0) OPEN(UNIT= 10,FILE='frame.out3',                          &
2606:      &form='formatted',STATUS='unknown')

2608: !     OPEN(UNIT= 11,FILE='frame.bin',
2609: !    +form='unformatted',STATUS='unknown')

2611: !     OPEN(UNIT= 12,FILE='frame.plt',
2612: !    +form='formatted',STATUS='unknown')

2614: !     OPEN(UNIT= 13,FILE='frame.tec',
2615: !    +form='formatted',STATUS='unknown')

2617: !     OPEN(UNIT= 14,FILE='frame.fast.g',
2618: !    +form='unformatted',STATUS='unknown')

2620: !     OPEN(UNIT= 15,FILE='frame.fast.q',
2621: !    +form='unformatted',STATUS='unknown')

2623:       return
2624:       end
2625: !
2626: !
2627: !===================================================================
2628: !
2629: ! Get the IA, JA, and IAU arrays
2630: !
2631: !===================================================================
2632:       subroutine GETIA(nnodes,nedge,evec,ia,ideg,irank)
2633:       implicit none
2634: #include <petsc/finclude/petscsys.h>
2635:       integer nnodes,nedge,irank
2636:       integer ia(1),ideg(1)
2637: #if defined(INTERLACING)
2638:        integer evec(2,nedge)
2639: #define eptr(j,i) evec(i,j)
2640: #else
2641:        integer evec(nedge,2)
2642: #define eptr(i,j) evec(i,j)
2643: #endif
2644:        integer i,node1,node2
2645: !
2646: ! First get the degree of each node using ideg as a dummy array
2647: !
2648:       do 1000 i = 1,nnodes
2649:         ideg(i) = 0
2650:  1000 continue
2651: !
2652:       do 1010 i = 1,nedge
2653:         node1 = eptr(i,1)
2654:         node2 = eptr(i,2)
2655: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
2656:         ideg(node1) = ideg(node1) + 1
2657:         ideg(node2) = ideg(node2) + 1
2658: #else
2659:         if (node1 .le. nnodes) ideg(node1) = ideg(node1) + 1
2660:         if (node2 .le. nnodes) ideg(node2) = ideg(node2) + 1
2661: #endif
2662:  1010 continue
2663: !
2664: ! Now we can fill the ia array fairly easily
2665: !
2666:       ia(1) = 1
2667:       do 1020 i = 1,nnodes
2668:         ia(i+1) = ia(i) + ideg(i) + 1
2669: !       write(9,100)i,ideg(i)
2670: ! 100   format(1h ,'deg(',i6,')=',i6)
2671:  1020 continue
2672: !
2673:       return
2674:       end
2675: !
2676: !===================================================================
2677: !
2678: ! Get the IA, JA, and IAU arrays
2679: !
2680: !===================================================================
2681:       subroutine GETJA(nnodes,nedge,evec,ia,ja,iwork,irank)
2682:       implicit none
2683: #include <petsc/finclude/petscsys.h>
2684:       integer nnodes,nedge,irank
2685:       integer ia(1),ja(1),iwork(1)
2686: #if defined(INTERLACING)
2687:        integer evec(2,nedge)
2688: #define eptr(j,i) evec(i,j)
2689: #else
2690:        integer evec(nedge,2)
2691: #define eptr(i,j) evec(i,j)
2692: #endif
2693:        integer i,index,node1,node2,index1,index2,                          &
2694:      &         istart,iend
2695: !     open(unit=90,file='map.out',status='UNKNOWN')
2696: !
2697: ! Now we need to get the JA array
2698: ! First fill the diagonal places
2699: !
2700:       do 1040 i = 1,nnodes
2701:         index = ia(i)
2702:         ja(index) = i
2703:         iwork(i) = 1
2704:  1040 continue
2705: !
2706:       do 1030 i = 1,nedge
2707:         node1 = eptr(i,1)
2708:         node2 = eptr(i,2)
2709: !
2710: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
2711:         index1 = ia(node1) + iwork(node1)
2712:         iwork(node1) = iwork(node1) + 1
2713:         ja(index1) = node2
2714:         index2 = ia(node2) + iwork(node2)
2715:         iwork(node2) = iwork(node2) + 1
2716:         ja(index2) = node1
2717: #else
2718:         if (node1 .le. nnodes) then
2719:           index1 = ia(node1) + iwork(node1)
2720:           iwork(node1) = iwork(node1) + 1
2721:           ja(index1) = node2
2722:         endif
2723:         if (node2 .le. nnodes) then
2724:           index2 = ia(node2) + iwork(node2)
2725:           iwork(node2) = iwork(node2) + 1
2726:           ja(index2) = node1
2727:         endif
2728: #endif
2729:  1030 continue
2730: !
2731: ! Now lets sort all our "bins" and get the correct one on the diagonal
2732: !
2733:       do 1050 i = 1,nnodes
2734:         istart = ia(i)
2735:         iend   = ia(i+1) - 1
2736: !       write(9,200)i,istart,iend
2737: ! 200   format(1h ,'Sorting ',i6,' istart iend = ',i6,1x,i6)
2738:         call SORTER(istart,iend,ja,i)
2739:  1050 continue
2740: !
2741: ! Now get the "fhelp" array which will assist in assembling
2742: ! the flux Jacobians into the correct location in the alu array
2743: !
2744: !     write(90,*) 'fhelp array'
2745: !     do 1060 i = 1,nedge
2746: !       node1 = eptr(i,1)
2747: !       node2 = eptr(i,2)
2748: !
2749: ! First take care of node1
2750: !
2751: !       idiag = iau(node1)
2752: !
2753: ! If the offdiagonal term is ordered later in the ja array
2754: !
2755: !       if (node2.gt.node1)then
2756: !        jstart = idiag + 1
2757: !        jend   = ia(node1+1) - 1
2758: !       else
2759: !        jstart = ia(node1)
2760: !        jend   = idiag -1
2761: !       end if
2762: !
2763: !        do 1070 j = jstart,jend
2764: !          if (ja(j).eq.node2)fhelp(i,1) = j
2765: !1070    continue
2766: !
2767: !
2768: ! Now take care of node2
2769: !
2770: !       idiag = iau(node2)
2771: !
2772: ! If the offdiagonal term is ordered later in the ja array
2773: !
2774: !       if (node1.gt.node2)then
2775: !        jstart = idiag + 1
2776: !        jend   = ia(node2+1) - 1
2777: !       else
2778: !        jstart = ia(node2)
2779: !        jend   = idiag -1
2780: !       end if
2781: !
2782: !       do 1080 j = jstart,jend
2783: !         if (ja(j).eq.node1)fhelp(i,2) = j
2784: !1080   continue
2785: !       write(90,*) i,fhelp(i,1),fhelp(i,2)
2786: !1060 continue
2787: !      close(90)
2788: !
2789:       return
2790:       end
2791: !
2792: !
2793: !===================================================================
2794: !
2795: ! Sort each of our bins
2796: !
2797: !===================================================================
2798:       subroutine SORTER(istart,iend,ja,inode)
2799:       implicit none
2800: #include <petsc/finclude/petscsys.h>
2801:       integer istart,iend,inode
2802:       integer ja(1)
2803: !
2804:       integer min,minsave,jsave,i,j
2805:       do 1000 i = istart,iend
2806:         min = ja(i)
2807:         minsave = ja(i)
2808:         jsave = i
2809:         do 1010 j = i+1,iend
2810:           if (ja(j).lt.min)then
2811:             min = ja(j)
2812:             jsave = j
2813:           end if
2814:  1010   continue
2815:         ja(i) = min
2816:         ja(jsave) = minsave
2817: !       if (ja(i).eq.inode)iau(inode) = i
2818:  1000 continue
2819: !
2820:       return
2821:       end