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


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

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

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

232:       read(7,10)

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

238:       read(7,10)

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

243:       read(7,10)

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

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

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

263:       return
264:       end



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


301:       character c4*4,title*20

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

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

325:         write(13,1010) (x(isnode(j))      ,j=isp,isp+nnpts(i)-1)
326:         write(13,1010) (y(isnode(j))      ,j=isp,isp+nnpts(i)-1)
327:         write(13,1010) (z(isnode(j))      ,j=isp,isp+nnpts(i)-1)
328:         write(13,1010) (1.0                 ,j=isp,isp+nnpts(i)-1)
329:         write(13,1010) (q(2,isnode(j))      ,j=isp,isp+nnpts(i)-1)
330:         write(13,1010) (q(3,isnode(j))      ,j=isp,isp+nnpts(i)-1)
331:         write(13,1010) (q(4,isnode(j))      ,j=isp,isp+nnpts(i)-1)
332:         write(13,1010) (q(1,isnode(j))/pinf ,j=isp,isp+nnpts(i)-1)

334:         do 30 j=ist,ist+nntet(i)-1
335:           i1 = f2ntn(j,1) - isp + 1
336:           i2 = f2ntn(j,2) - isp + 1
337:           i3 = f2ntn(j,3) - isp + 1

339:           write(13,1020) i1,i2,i3,i3
340:    30   continue

342:         isp = isp + nnpts(i)
343:         ist = ist + nntet(i)

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

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

360:         write(13,1010) (x(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
361:         write(13,1010) (y(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
362:         write(13,1010) (z(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
363:         write(13,1010) (1.0                 ,j=isp,isp+nvpts(i)-1)
364:         write(13,1010) (q(2,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
365:         write(13,1010) (q(3,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
366:         write(13,1010) (q(4,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
367:         write(13,1010) (q(1,ivnode(j))/pinf ,j=isp,isp+nvpts(i)-1)

369:         do 90 j=ist,ist+nvtet(i)-1
370:            i1 = f2ntv(j,1) - isp + 1
371:            i2 = f2ntv(j,2) - isp + 1
372:            i3 = f2ntv(j,3) - isp + 1

374:            write(13,1020) i1,i2,i3,i3
375:    90   continue

377:         isp = isp + nvpts(i)
378:         ist = ist + nvtet(i)

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

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

395:          write(13,1010) (x(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
396:          write(13,1010) (y(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
397:          write(13,1010) (z(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
398:          write(13,1010) (1.0                 ,j=isp,isp+nfpts(i)-1)
399:          write(13,1010) (q(2,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
400:          write(13,1010) (q(3,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
401:          write(13,1010) (q(4,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
402:          write(13,1010) (q(1,ifnode(j))/pinf ,j=isp,isp+nfpts(i)-1)

404:          do 60 j=ist,ist+nftet(i)-1
405:             i1 = f2ntf(j,1) - isp + 1
406:             i2 = f2ntf(j,2) - isp + 1
407:             i3 = f2ntf(j,3) - isp + 1

409:             write(13,1020) i1,i2,i3,i3
410:    60    continue

412:          isp = isp + nfpts(i)
413:          ist = ist + nftet(i)

415:    40  continue

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

424:        return
425:        end

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

455:       integer nnodes,nedge,nnfacet,nvfacet,                             &
456:      &        nnbound,nvbound,                                          &
457:      &        irank,nvertices
458:       integer isnode(1),ivnode(1)
459:       integer f2ntn(nnfacet,4)
460:       integer f2ntv(nvfacet,4)
461:       integer sface_bit(nnfacet), vface_bit(nvfacet)

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

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

520:       xr = 0.25d0
521:       yr = 0.0d0
522:       zr = 0.0d0

524:       do 30 n = 1, nnfacet
525:                node1 = isnode(f2ntn(n,1))
526:                node2 = isnode(f2ntn(n,2))
527:                node3 = isnode(f2ntn(n,3))

529:                x1    = x(node1)
530:                y1    = y(node1)
531:                z1    = z(node1)

533:                x2    = x(node2)
534:                y2    = y(node2)
535:                z2    = z(node2)

537:                x3 = x(node3)
538:                y3 = y(node3)
539:                z3 = z(node3)

541:                ax = x2 - x1
542:                ay = y2 - y1
543:                az = z2 - z1

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

556:                p1 = qnod(1,node1)
557:                u1 = qnod(2,node1)
558:                v1 = qnod(3,node1)
559:                w1 = qnod(4,node1)

561:                p2 = qnod(1,node2)
562:                u2 = qnod(2,node2)
563:                v2 = qnod(3,node2)
564:                w2 = qnod(4,node2)

566:                p3 = qnod(1,node3)
567:                u3 = qnod(2,node3)
568:                v3 = qnod(3,node3)
569:                w3 = qnod(4,node3)

571:                press = (p1 + p2 + p3)/3.0d0
572:                cp    = 2.d0*(press - 1.0d0)
573: !
574:                dcx = cp*xnorm
575:                dcy = cp*ynorm
576:                dcz = cp*znorm

578:                xmid = (x1 + x2 + x3)
579:                ymid = (y1 + y2 + y3)
580:                zmid = (z1 + z2 + z3)

582:                if (sface_bit(n) .eq. 1) then
583: #if defined(CFL3D_AXIS)

585:                 clloc = clloc - dcx*sna     + dcz*csa
586:                 cdloc = cdloc + dcx*csa + dcz*sna
587:                 cmloc = cmloc                                           &
588:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx

590: #else
591:                 clloc = clloc - dcx*sna + dcy*csa
592:                 cdloc = cdloc + dcx*csa + dcy*sna
593:                 cmloc = cmloc                                           &
594:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx
595: #endif
596:                endif

598:  30    continue
599: !
600: ! Viscous boundary
601: !
602:        do 60 n = 1, nvfacet
603:                node1 = ivnode(f2ntv(n,1))
604:                node2 = ivnode(f2ntv(n,2))
605:                node3 = ivnode(f2ntv(n,3))

607:                x1    = x(node1)
608:                y1    = y(node1)
609:                z1    = z(node1)

611:                x2    = x(node2)
612:                y2    = y(node2)
613:                z2    = z(node2)

615:                x3 = x(node3)
616:                y3 = y(node3)
617:                z3 = z(node3)

619:                ax = x2 - x1
620:                ay = y2 - y1
621:                az = z2 - z1

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

634:                p1    = qnod(1,node1)
635:                u1    = qnod(2,node1)
636:                v1    = qnod(3,node1)
637:                w1    = qnod(4,node1)

639:                p2    = qnod(1,node2)
640:                u2    = qnod(2,node2)
641:                v2    = qnod(3,node2)
642:                w2    = qnod(4,node2)

644:                p3    = qnod(1,node3)
645:                u3    = qnod(2,node3)
646:                v3    = qnod(3,node3)
647:                w3    = qnod(4,node3)

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

652:                dcx = cp*xnorm
653:                dcy = cp*ynorm
654:                dcz = cp*znorm

656:                xmid = (x1 + x2 + x3)
657:                ymid = (y1 + y2 + y3)
658:                zmid = (z1 + z2 + z3)

660:                if (vface_bit(n) .eq. 1) then
661:                 clloc = clloc - dcx*sna + dcy*csa
662:                 cdloc = cdloc + dcx*csa + dcy*sna
663:                 cmloc = cmloc                                            &
664:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx
665:                endif

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


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

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

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

830:             u   = qnod(2,inode)
831:             v   = qnod(3,inode)
832:             w   = qnod(4,inode)
833:             ubar= (u*xnorm + v*ynorm + w*znorm)/area
834:             c   = sqrt(ubar*ubar + beta)
835: !
836:             Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
837:             cdt(inode) = cdt(inode) + Vn
838: !
839:  1030    continue
840:          flops = flops + 24.0*nsnode

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

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

878:             u   = qnod(2,inode)
879:             v   = qnod(3,inode)
880:             w   = qnod(4,inode)
881:             ubar= (u*xnorm + v*ynorm + w*znorm)/area
882:             c = sqrt(ubar*ubar + beta)

884:             Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
885:             cdt(inode) = cdt(inode) + Vn
886:  1050    continue
887:          flops = flops + 24.0*nfnode

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


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

1183: !      print *, 'On thread ', my_thread_id,', first_edge = ',first_edge,        &
1184: !     &         ', last_edge = ',last_edge

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

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

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

1330:           r12 = 0.0d0
1331:           r22 = X2
1332:           r32 = Y2
1333:           r42 = Z2

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

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

1457:                x1 = x(node1)
1458:                y1 = y(node1)
1459:                z1 = z(node1)
1460:                p1 = qnod(1,node1)

1462:                x2 = x(node2)
1463:                y2 = y(node2)
1464:                z2 = z(node2)
1465:                p2 = qnod(1,node2)

1467:                x3 = x(node3)
1468:                y3 = y(node3)
1469:                z3 = z(node3)
1470:                p3 = qnod(1,node3)

1472:                ax = x2 - x1
1473:                ay = y2 - y1
1474:                az = z2 - z1

1476:                bx = x3 - x1
1477:                by = y3 - y1
1478:                bz = z3 - z1
1479: !
1480: ! Normal points away from grid interior.
1481: ! Magnitude is 1/3 area of surface triangle.
1482: !
1483:                xnorm =-0.5d0*(ay*bz - az*by)/3.d0
1484:                ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
1485:                znorm =-0.5d0*(ax*by - ay*bx)/3.d0

1487:                pa = c68*p1 + c18*(p2 + p3)
1488:                pb = c68*p2 + c18*(p3 + p1)
1489:                pc = c68*p3 + c18*(p1 + p2)
1490: !
1491:                flops = flops + 35.0
1492:           if (node1.le.nnodes) then
1493:                res(2,node1) = res(2,node1) + xnorm*pa
1494:                res(3,node1) = res(3,node1) + ynorm*pa
1495:                res(4,node1) = res(4,node1) + znorm*pa
1496:                flops = flops + 6.0
1497:           endif

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

1506:           if (node3.le.nnodes) then
1507:                res(2,node3) = res(2,node3) + xnorm*pc
1508:                res(3,node3) = res(3,node3) + ynorm*pc
1509:                res(4,node3) = res(4,node3) + znorm*pc
1510:                flops = flops + 6.0
1511:           endif

1513:       enddo
1514: !
1515: ! Now viscous faces
1516: !
1517:       do n = 1, nvfacet
1518:                node1 = ivnode(f2ntv(n,1))
1519:                node2 = ivnode(f2ntv(n,2))
1520:                node3 = ivnode(f2ntv(n,3))

1522:                x1 = x(node1)
1523:                y1 = y(node1)
1524:                z1 = z(node1)
1525:                p1 = qnod(1,node1)

1527:                x2 = x(node2)
1528:                y2 = y(node2)
1529:                z2 = z(node2)
1530:                p2 = qnod(1,node2)

1532:                x3 = x(node3)
1533:                y3 = y(node3)
1534:                z3 = z(node3)
1535:                p3 = qnod(1,node3)

1537:                ax = x2 - x1
1538:                ay = y2 - y1
1539:                az = z2 - z1

1541:                bx = x3 - x1
1542:                by = y3 - y1
1543:                bz = z3 - z1
1544: !
1545: ! norm point away from grid interior.
1546: ! norm magnitude is 1/3 area of surface triangle.
1547: !
1548:                xnorm =-0.5d0*(ay*bz - az*by)/3.d0
1549:                ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
1550:                znorm =-0.5d0*(ax*by - ay*bx)/3.d0

1552:                pa = c68*p1 + c18*(p2 + p3)
1553:                pb = c68*p2 + c18*(p3 + p1)
1554:                pc = c68*p3 + c18*(p1 + p2)
1555: !
1556:                flops = flops + 35.0
1557: !
1558:           if (node1.le.nnodes) then
1559:                 res(2,node1) = res(2,node1) + xnorm*pa
1560:                 res(3,node1) = res(3,node1) + ynorm*pa
1561:                 res(4,node1) = res(4,node1) + znorm*pa
1562:                 flops = flops + 6.0

1564:           endif

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

1573:           if (node3.le.nnodes) then
1574:                 res(2,node3) = res(2,node3) + xnorm*pc
1575:                 res(3,node3) = res(3,node3) + ynorm*pc
1576:                 res(4,node3) = res(4,node3) + znorm*pc
1577:                 flops = flops + 6.0
1578:           endif

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

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

1693:          t11 = 0.0d0
1694:          t21 = X1
1695:          t31 = Y1
1696:          t41 = Z1

1698:          t12 = 0.0d0
1699:          t22 = X2
1700:          t32 = Y2
1701:          t42 = Z2

1703:          t13 =  c0*beta
1704:          t23 = xnorm*beta + u0*(ubar0 + c0)
1705:          t33 = ynorm*beta + v0*(ubar0 + c0)
1706:          t43 = znorm*beta + w0*(ubar0 + c0)

1708:          t14 = -c0*beta
1709:          t24 = xnorm*beta + u0*(ubar0 - c0)
1710:          t34 = ynorm*beta + v0*(ubar0 - c0)
1711:          t44 = znorm*beta + w0*(ubar0 - c0)

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

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

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

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

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

1998:       do i = 1,nnodes
1999:            w11 = 1.d0/(r11(i)*r11(i))
2000:            w22 = 1.d0/(r22(i)*r22(i))
2001:            w33 = 1.d0/(r33(i)*r33(i))
2002:            r12r11 = r12(i)/r11(i)
2003:            r13r11 = r13(i)/r11(i)
2004:            r23r22 = r23(i)/r22(i)
2005:            rmess=(r12(i)*r23(i)-                                             &
2006:      &            r13(i)*r22(i))/                                            &
2007:      &           (r11(i)*r22(i)*                                             &
2008:      &            r33(i)*r33(i))

2010:            r11(i) = w11
2011:            r22(i) = w22
2012:            r33(i) = w33
2013:            r12(i) = r12r11
2014:            r13(i) = r13r11
2015:            r23(i) = r23r22
2016:            r44(i) = rmess
2017:       enddo
2018: !
2019: ! Finished with SUMGS
2020: !
2021:       return
2022:       end


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

2108: !     if (flag .eq. -1) then
2109: !        call PetscLogEventRegister(grad_event,grad_label,ierr)
2110: !        call PetscLogEventRegister(node1_event,node1_label,ierr)
2111: !        call PetscLogEventRegister(node2_event,node2_label,ierr)
2112: !        flag = 1
2113: !     endif
2114: !     call PetscLogEventBegin(grad_event,0,0,0,0,ierr)

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

2138:          gradx(2,i) = 0.0
2139:          grady(2,i) = 0.0
2140:          gradz(2,i) = 0.0

2142:          gradx(3,i) = 0.0
2143:          grady(3,i) = 0.0
2144:          gradz(3,i) = 0.0

2146:          gradx(4,i) = 0.0
2147:          grady(4,i) = 0.0
2148:          gradz(4,i) = 0.0
2149:  1000 continue

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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