program ext
implicit none
c
c mnel = maximum number of elements
c mncon = maximum number of connections
c mv = maximum number of values (p,t,sg,...,fx,fy,...)
c nhsh = size of hash table
c
integer mnel
integer mncon
integer mcat
integer mv
integer nhsh
parameter (mnel = 32000)
parameter (mncon = 120000)
parameter (mcat = 140000)
parameter (mv = 80)
parameter (nhsh = 256)
common/id/id
common/xy/ x
common/ec/ ec
common/e/ e
common/cc/ cc
common/hsh/ hsh
common/ind/ ind
common/v/ v
common/zn/ zn
common/zo/ zo
common/cat/ cat
character*5 id(mnel)
double precision x(3,mnel)
integer ec(mnel)
integer e(2,mncon)
integer cc(2,mncon)
integer hsh(nhsh)
integer ind(mnel)
double precision a(mncon)
double precision v(mv,mnel)
double precision pv
integer c1(mv)
integer cl(mv)
integer uv(20)
integer cat(4,mcat)
integer zn(mnel)
integer zo(mnel)
character*4 cprnt
character*140 file,line,lc,lc2
character*400 v1,want
double precision denom
double precision dx,dy,dz
integer i
integer i1
integer i11
integer i1l
integer i2
integer i21
integer i2l
integer iargc
integer ic
integer ihash
integer ii1
integer iil
integer il
integer iprnt
integer itype
integer iwant
integer ixyz
integer izone
integer j
integer k
integer l
integer l0
integer l1
integer l9
integer ll
integer lprnt
integer lv
integer nc
integer ncat
integer nd
integer ne
integer ne1
integer ne2
integer nn
integer nu1
integer nu2
integer nv
integer nxyz
c
c read arguments
c
ixyz = 0
izone = 0
l = iargc()
i = 0
iwant = 0
10 i = i + 1
if(i .le. l) then
call getarg(i, file)
if(file(1:1) .eq. '-') then
if(file .eq. '-yz') then
ixyz = 1
else if(file .eq. '-xz') then
ixyz = 2
else if(file .eq. '-xy') then
ixyz = 3
else if(file .eq. '-xyz') then
ixyz = 4
else if(file .eq. '-z' .and. i .lt. l) then
izone = 1
i = i + 1
call getarg(i, file)
open(unit = 3,file=file,status='old')
else if(file .eq. '-w' .and. i .lt. l) then
i = i + 1
call getarg(i, file)
open(unit = 1,file=file,status='old')
want(1:1) = ' '
read(1,'(a)') want(2:)
close(unit=1)
call inspect(want,want,l0,iwant,k)
else
stop 'unknown argument'
endif
goto 10
else if(i .lt. l) then
stop 'too many arguments'
endif
else
write(6,'(a,$)') 'Input file name: '
read(5,'(a)') file
endif
c
c read mesh file
c
call rmesh(ne,id,x,ec,e,a,cc,hsh,ind,mnel,mncon,nhsh,ixyz)
do i = 1, ne
zo(i) = i
zn(i) = i
enddo
ne1 = ne
c
c compute catalog
c
if(ixyz .lt. 4 .and. izone .eq. 0) then
call catalog(ncat,cat,x,ec,e,cc,ne,mcat)
endif
c
c loop through each printed time step (count in iprnt)
c locate beginning of printout and read time
c
iprnt = 0
open (unit=1,file=file,status='old')
i11 = 1
20 if(lc(i11:i11+9) .ne. 'total time') then
read(1,'(a)',end=240) line
line(1:1) = ' '
call inspect(line,lc,i11,l,i1l)
goto 20
endif
read(1,'(a)',end=240) line
iprnt = iprnt + 1
c
c get printout count (iprnt) as character for use in file names
c (cprnt) -- first nonzero character is in lprnt
c
write(cprnt,'(i4)') iprnt
call inspect(cprnt,cprnt,lprnt,l1,k)
write(6,30) cprnt(lprnt:4),line(3:20)
30 format(' Found printout(',a,'): ',a)
c
c open output file
c
line = file
call inspect(line,line,l1,l,k)
l = l + 1
line(l:l+5-lprnt) = '.'//cprnt(lprnt:4)
l = l + 5 - lprnt
open (unit=2,file=line(1:l),status='unknown')
c
c look through the printout
c
nv = 0
40 read(1,'(a)',end=100) line
line(1:1) = ' '
call inspect(line,lc,i11,l,i1l)
if(lc(i11:i11+9) .eq. 'total time') goto 110
if(lc(i11:i11+3) .ne. 'elem') goto 40
c
c choose first variable names
c
if(nv .eq. 0) then
ic = ichar(line(i11:i11))
nxyz = 2
v1 = 'Variables = x y z'
lv = 11
v1(lv+2:lv+2) = char(ic+19)
if(ixyz .ne. 1) lv = lv + 2
v1(lv+2:lv+2) = char(ic+20)
if(ixyz .ne. 2) lv = lv + 2
v1(lv+2:lv+2) = char(ic+21)
if(ixyz .ne. 3) lv = lv + 2
if(lv .eq. 17) nxyz = 3
endif
itype = 6
i2l = i1l
50 call inspect(lc(i2l+1:l),lc(i2l+1:l),ii1,l9,iil)
ii1 = ii1 + i2l
iil = iil + i2l
if(lc(ii1:iil) .ne. 'index') then
if(itype .ne. 6) goto 40
itype = 5
i21 = ii1
i2l = iil
if(lc(ii1:iil) .eq. 'elem2') goto 50
goto 40
endif
il = iil
nn = 0
nd = 0
60 if(il .lt. l) then
call inspect(lc(il+1:l),lc(il+1:l),i1,l9,ll)
i1 = i1 + il
il = ll + il
nn = nn + 1
c1(nn) = i1
cl(nn) = il
uv(nn) = 1
c
c include a check as to whether this variable lc(i1:il)
c is wanted here
c
if(iwant .eq. 0) goto 70
if(index(want,lc(i1-1:il+1)) .ne. 0) goto 70
uv(nn) = 0
nd = nd + 1
goto 60
70 if(lc(i1:i1+2) .eq. 'vel') uv(nn) = -1
l9 = il + 2 - i1
v1(lv+1:lv+l9) = line(i1-1:il)
i2 = lv + 2
if(itype .eq. 5 .and. v1(i2+3:i2+3) .eq. '(') then
v1(i2+1:lv+l9) = v1(i2+4:lv+l9)
l9 = l9 - 3
if(l9 .gt. 6 .and. v1(i2+4:i2+4) .eq. ')') then
v1(i2+2:lv+l9) = v1(i2+5:lv+l9)
l9 = l9 - 3
endif
else if(itype .eq. 6 .and. v1(i2+1:i2+1) .eq. '(') then
v1(i2+1:lv+l9) = v1(i2+2:lv+l9)
l9 = l9 - 2
endif
if(l9 .gt. itype) l9 = itype
lv = lv + l9
if(itype .eq. 6) goto 60
i2 = i2 - 1
lv = lv + 1
v1(lv:lv) = 'x'
if(ixyz .ne. 1) then
v1(lv+1:lv+l9) = v1(i2:lv-1)
lv = lv + l9 + 1
endif
v1(lv:lv) = 'y'
if(ixyz .eq. 2) then
v1(lv:lv) = 'z'
else if(ixyz .ne. 3) then
v1(lv+1:lv+l9) = v1(i2:i2+l9-1)
lv = lv + l9 + 1
v1(lv:lv) = 'z'
endif
goto 60
endif
nd = nn - nd
if(itype .eq. 5) nd = nd*nxyz
do i = 1, nd
do nu1 = 1, ne
v(i + nv,nu1) = 0.
enddo
enddo
i = nv
if(nv + nd .gt. mv) goto 110
nv = nv + nd
nd = i
c
c a line of @'s ends this section
c blank lines are ignored as are the page break headings
c
l0 = 0
nu2 = -1
80 read(1,'(a)') line
if(line(2:21) .eq. '@@@@@@@@@@@@@@@@@@@@') goto 40
nu1 = ihash(line(i11:i1l),id,ind,hsh,ne,nhsh)
if(itype .eq. 5) then
nu2 = ihash(line(i21:i1l),id,ind,hsh,ne,nhsh)
if(l0 .eq. 0 .and. nu1 .eq. 0 .and. nu2 .eq. 0 .and.
1 line(i11:i11+1) .eq. ' ') then
nu1 = ihash(line(i11+2:i1l+2),id,ind,hsh,ne,nhsh)
nu2 = ihash(line(i21+2:i1l+2),id,ind,hsh,ne,nhsh)
if(nu1 .eq. 0 .or. nu2 .eq. 0) goto 80
i11 = i11 + 2
i1l = i1l + 2
i21 = i21 + 2
i2l = i2l + 2
endif
endif
if(nu1 .eq. 0 .or. nu2 .eq. 0) goto 80
if(l0 .eq. 0) then
c
c determine which columns contain the numbers
c
call inspect(line(ii1:),lc,l0,l,k)
l = l + ii1 - 1
if(cl(nn).lt.l) cl(nn) = l
c1(nn+1) = l + 1
cl(nn+1) = l + 1
iil = ii1 + k - 1
il = iil
do i = 1, nn
i1 = il
il = cl(i+1)
call inspect(line(i1+1:il),lc,l0,l,k)
il = i1 + k
if(k .gt. 40) k = 40
i1 = il - k
c1(i) = i1 + 1
cl(i) = il
enddo
endif
if(itype .eq. 5) then
c
c extra flow printout computation
c
dx = x(1,nu1) - x(1,nu2)
dy = x(2,nu1) - x(2,nu2)
dz = x(3,nu1) - x(3,nu2)
denom = dx*dx + dy*dy
if(nxyz .eq. 3) denom = denom + dz*dz
if(denom .le. 0.) goto 80
nc = ec(nu1)
do while (nc .ne. 0 .and. nu2 .ne. e(2,nc))
i = (nu1 - e(1,nc))/(e(2,nc) - e(1,nc))
nc = cc(i+1,nc)
enddo
if(nc .eq. 0) then
stop 'no connection'
c gchen: next 2 lines
c print *, 'Missing connection'
c go to 80
endif
c
c compute 1/2 unit vector (dx,dy) from second element to first
c (this is the positive flow direction).
c 1/2 is because the flow is added both entering and leaving
c the element.
c
denom = sqrt(denom)
denom = denom + denom
dx = dx/denom
dy = dy/denom
dz = dz/denom
endif
i = nd
do j = 1, nn
if(uv(j) .eq. 0) goto 90
i1 = c1(j)
il = cl(j)
k = i1 + 40 - il
lc(1:k-1) = ' '
lc(k:40) = line(i1:il)
read(lc,'(f40.0)') pv
if(itype .eq. 6) then
i = i + 1
v(i,nu1) = pv
goto 90
endif
denom = 1.
if(uv(j) .gt. 0) denom = a(nc)
c --- added by gchen
c denom=1.0
c
v(i+1,nu1) = v(i+1,nu1) + dx*pv/denom
v(i+1,nu2) = v(i+1,nu2) + dx*pv/denom
v(i+2,nu1) = v(i+2,nu1) + dy*pv/denom
v(i+2,nu2) = v(i+2,nu2) + dy*pv/denom
if(nxyz .eq. 3) then
v(i+3,nu1) = v(i+3,nu1) + dz*pv/denom
v(i+3,nu2) = v(i+3,nu2) + dz*pv/denom
endif
i = i + nxyz
90 enddo
goto 80
c
c data has all been collected
c
100 iprnt = -1
110 write(2,'(a)') v1(1:lv)
l = 0
if(izone .eq. 0) goto 190
v1 = 'Zone'
lv = 4
rewind 3
120 ncat = 0
ne1 = 0
do i = 1, ne
zo(i) = 0
zn(i) = 0
enddo
130 l = 0
read(3,'(a)',end=150) line
call inspect(line,lc2,l0,l,k)
if(lc2(l0:l0+4) .ne. 'zone') then
ncat = ncat + 1
ne2 = ne1
l0 = 2
do i = 1, 4
nu1 = ihash(line(l0:l0+4),id,ind,hsh,ne,nhsh)
if(nu1 .eq. 0) then
write(6,140) line(l0:l0+4)
140 format(' Unknown element "',A5,'" in zone file')
ncat = ncat - 1
goto 130
endif
l0 = l0 + 6
j = zn(nu1)
if(j .eq. 0) then
ne2 = ne2 + 1
j = ne2
zn(nu1) = j
zo(j) = nu1
endif
cat(i,ncat) = j
enddo
ne1 = ne2
goto 130
endif
150 if(ne1 .ne. 0) goto 170
160 if(l .eq. 0) goto 230
v1 = line(l0:l)
lv = l - l0 + 1
goto 120
170 write(2,180) v1(1:lv),ne1,ncat
180 format(a,' F=FEPOINT I=',i5,' J=',i5)
goto 220
190 if(ncat .ne. 0) then
write(2,200) ne,ncat
200 format('Zone F=FEPOINT I=',i5,' J=',i5)
else
write(2,210) ne
210 format('Zone F=POINT I=',i5)
endif
220 do i = 1, ne1
j = zo(i)
write(2,'(1p,50e14.6)') (x(k,j),k=1,nxyz),(v(k,j),k=1,nv)
enddo
do i = 1, ncat
write(2,'(4i6)') (cat(j,i),j=1,4)
enddo
if(l .ne. 0) goto 160
c
c close the output file
c
230 close (unit=2)
if(iprnt .ge. 0) goto 20
c
c close input and end
c
240 close (unit=1)
stop
end