program tridemp

c parep2.lbl.gov:$MDOCS/tridemp/tridemp.f 10/10/96
c convert a CL4 DEMP map to TRI format
c for each polygon:

c	warning - maximum of 2200 polygons, 900 points per polygon
c
	character argv1*79
	character argv2*79
	character argv3*79
	character argv4*79
	character argv9*79
	character argv8*79
	character argv7*79

c input file 1: pre-DEMP CL4 (one line per polygon)
c input file 2: post-DEMP CL4 (one line per polygon)
c input file 3: pre-DEMP TRI (one line per triangle)   (7)+(11)
c output file 4: post-DEMP TRI (one line per triangle) (7)+(11)
c output file 9: population and areas for whole map  (.mappop)
c	line 1: pre-demp  x0,y0,scalar,xmin,xmax,ymin,ymax (7)
c	line 2: post-demp x0,y0,scalar,xmin,xmax,ymin,ymax (7)
c	line 3: pop, 2x predemp area, 2x postdemp area, npoly, ngeo (5)
c output file 8: population and areas (one line per geocode) (.geopop)
c	seq,  pop, 2x predemp area, 2x postdemp area, 4 geocodes, npoly (8)
c output file 7: population and areas (one line per polygon) (.polypop)
c	polyseq, pop, 2x predemp area, 2x postdemp area, 4 geocodes (8)

c tridemp.exe input.1 input.2 input.3 output.4 output.9 output.8 output.7

c	number of points in current polygon
	integer npts

c	from input CL4 files: oldx, oldy, newx, newy
	integer kxy(4,10000)
c	max number of points including duplicates
	data kmax /10000/
c	max number of polygons
	data maxpoly /2200/
c	max number of points in any one polygon
      	data maxnpts /900/
c	point array for all polygons - file 1
	integer jxarray(2200,900), jyarray(2200,900)
c	point array for all polygons - file 2
	integer kxarray(2200,900), kyarray(2200,900)

c	5 geocodes for each polygon - maximum of 2200 polygons
	integer jarray(5,2200)
c	number of points in each polygon
	integer nptsarray(2200)
c	pop of each polygon - max of 2200 polygons
	real poparray(2200)
c	2 x area of each polygon - integer
	integer jareaarray(2200)
	integer kareaarray(2200)


	iunit1=1
	iunit2=2
	iunit3=3
	iunit4=4
	iunit7=7
	iunit8=8
	iunit9=9


c	debug output
	iotmsg=6

c	implemented only for 4-geocode levels eg clsc,nmcdtr80,etc
	level=41

c ---------------------------------------------------------
c input file #1 - pre-DEMP Cl4
c file 1: jxarray,jyarray,jxmin,jxmax,jymin,jymax,jareaarray
	call getarg(1,argv1)
	write(iotmsg,*) argv1
	open (iunit1,file=argv1,type='old')
	call read_cl4a(iunit1,level,maxpoly,maxnpts,x0,y0,scalar,
     x	jxmin,jxmax,jymin,jymax,fconv,maxd,ngeoplus2,jarray,nptsarray,
     x	jxarray,jyarray,nbpt,poparray,jareaarray)

c nickel header as stored in pre-DEMP file
	write(iotmsg,710)x0,y0,scalar,jxmin,jxmax,jymin,jymax
c 2 x area of first polygon
	jarea = -jareaarray(1)
	write(iotmsg,1000) jarea
	write(iotmsg,1000) maxd
c ---------------------------------------------------------
c all variables from file #1 are overwritten except the coordinate values
c input file #2 - post-DEMP Cl4
c file 2: kxarray,kyarray,kxmin,kxmax,kymin,kymax,kareaarray

	call getarg(2,argv2)
	write(iotmsg,*) argv2
	open (iunit2,file=argv2,type='old')
	call read_cl4a(iunit2,level,maxpoly,maxnpts,x0,y0,scalar,
     x	kxmin,kxmax,kymin,kymax,fconv,maxd,ngeoplus2,jarray,nptsarray,
     x	kxarray,kyarray,nbpt,poparray,kareaarray)

c nickel header as stored in post-DEMP file
	write(iotmsg,710)x0,y0,scalar,kxmin,kxmax,kymin,kymax
c 2 x area of first polygon
	karea = -kareaarray(1)
	write(iotmsg,1000) karea
	write(iotmsg,1000) maxd

c -------------------------------------------------------------------
c       store old and new coordinates from files 1 and 2

c	counter for array of stored coordinates
	k=0

	do 3402 i=1,maxd
c	number of points in polygon i
	npts = nptsarray(i)
	do 3403 j=1,npts
	if (k.ge.kmax) call exit(1)
	k=k+1
c	oldx, oldy, newx, newy
	kxy(1,k)=jxarray(i,j)
	kxy(2,k)=jyarray(i,j)
	kxy(3,k)=kxarray(i,j)
	kxy(4,k)=kyarray(i,j)
3403	continue
3402	continue
	ktot=k
	write(iotmsg,1000) ktot

c -------------------------------------------------------------------
	call getarg(5,argv9)
	write(iotmsg,*) argv9
	open (iunit9,file=argv9)

	call getarg(6,argv8)
	write(iotmsg,*) argv8
	open (iunit8,file=argv8)

	call getarg(7,argv7)
	write(iotmsg,*) argv7
	open (iunit7,file=argv7)

c	for unit 9 (one line for whole map)
	mappoly=0
	geopopmap=0.
	jareamap=0
	kareamap=0
	
c	for unit 8 (one line per geocode)
	igeo=0
	npoly=0
	geopop=0.
	jareageo=0
	kareageo=0
	j1old=0
	j2old=0
	j3old=0
	j4old=0

c	whenever a new geocode is encountered, write to unit 8
c	geopop,jareageo,kareageo and old geocodes

c	skip first polygon (external boundary)
	do 2000 i=2,maxd

c	population, orig area, demp area for each polygon
	write (7,2001) i,poparray(i),jareaarray(i),kareaarray(i),
     x	(jarray(j,i),j=1,4)
2001	format (i6,f15.5,2i10,5i7)

	j1=jarray(1,i)
	j2=jarray(2,i)
	j3=jarray(3,i)
	j4=jarray(4,i)
	if (j1.ne.j1old.or.j2.ne.j2old
     x	.or.j3.ne.j3old.or.j4.ne.j4old) then
	    if (npoly.gt.0) write (8,2001) igeo,geopop,jareageo,kareageo,
     x		j1old,j2old,j3old,j4old,npoly
		npoly=0
		geopop=0
		jareageo=0
		kareageo=0
		igeo=igeo+1
	endif

	mappoly=mappoly+1
	geopopmap=geopopmap+poparray(i)
	jareamap=jareamap+jareaarray(i)
	kareamap=kareamap+kareaarray(i)

	npoly=npoly+1
	geopop=geopop+poparray(i)
	jareageo=jareageo+jareaarray(i)
	kareageo=kareageo+kareaarray(i)

	j1old=j1
	j2old=j2
	j3old=j3
	j4old=j4
2000	continue
	if (npoly.gt.0) write (8,2001) igeo,geopop,jareageo,kareageo,
     x		j1old,j2old,j3old,j4old,npoly

c	write headers to file 9
	write(iunit9,710) x0,y0,scalar,jxmin,jxmax,jymin,jymax
	write(iotmsg,710) x0,y0,scalar,jxmin,jxmax,jymin,jymax
	write(iunit9,710) x0,y0,scalar,kxmin,kxmax,kymin,kymax
	write(iotmsg,710) x0,y0,scalar,kxmin,kxmax,kymin,kymax
c	write one line to file 9
	write(iunit9,2002) geopopmap,jareamap,kareamap,mappoly,igeo
	write(iotmsg,2002) geopopmap,jareamap,kareamap,mappoly,igeo
2002	format (f18.5,2i20,2i7)

c -------------------------------------------------------------------
c	read file 3 - input TRI file
	
c	read and write header line
	call getarg(3,argv3)
	write(iotmsg,*) argv3
	open (iunit3,file=argv3,type='old')
	read(iunit3,710) x0,y0,scalar,ixmin,ixmax,iymin,iymax
	write(iotmsg,710)x0,y0,scalar,ixmin,ixmax,iymin,iymax

c	area of input and output map
	iareatot = 0
	jareatot = 0

	call getarg(4,argv4)
	write(iotmsg,*) argv4
	open (iunit4,file=argv4)
c	header line of file 4 same as file 2
	write(iunit4,710) x0,y0,scalar,kxmin,kxmax,kymin,kymax

5000	continue
	read(iunit3,3004,end=999) ixa,iya,ixb,iyb,ixc,iyc,iarea,
     x		i1,i2,i3,i4
	iareatot = iareatot+iarea

c	translate old to new coords

	ierra=1
	jxa=0
	jya=0
	do 6001 k=1,ktot
	if (ixa.eq.kxy(1,k).and.iya.eq.kxy(2,k)) then
		ierra=0
		jxa=kxy(3,k)
		jya=kxy(4,k)
		goto 60011
	endif
6001	continue
60011	continue

	ierrb=1
	jxb=0
	jyb=0
	do 6002 k=1,ktot
	if (ixb.eq.kxy(1,k).and.iyb.eq.kxy(2,k)) then
		ierrb=0
		jxb=kxy(3,k)
		jyb=kxy(4,k)
		goto 60021
	endif
6002	continue
60021	continue

	ierrc=1
	jxc=0
	jyc=0
	do 6003 k=1,ktot
	if (ixc.eq.kxy(1,k).and.iyc.eq.kxy(2,k)) then
		ierrc=0
		jxc=kxy(3,k)
		jyc=kxy(4,k)
		goto 60031
	endif
6003	continue
60031	continue

	jarea=0
	if (ierra+ierrb+ierrc .eq. 0) then
	jarea=-(jxa*jyb-jxb*jya + 
     x		jxb*jyc-jxc*jyb + jxc*jya-jxa*jyc)
	jareatot=jareatot+jarea
	endif

	write (iunit4,3004) jxa,jya,jxb,jyb,jxc,jyc,jarea,i1,i2,i3,i4
3004    format(6i7,i14, 1x,i4,1x,i3,1x,i3,1x,i6)

	goto 5000

999	continue
c	area from file 3
	write(iotmsg,1000) iareatot
c	header line of file 4 same as file 2
	write(iotmsg,710) x0,y0,scalar,kxmin,kxmax,kymin,kymax
c	area from file 4
	write(iotmsg,1000) jareatot
1000	format (i20)

	call exit(0)

	include 'format.inc'

	end

	include 'fread_cl4a.sub'