c parep2.lbl.gov:$MDOCS/cl4_check/fread_cl4.sub 11/25/99
c from sy$seedis:[seedis.map_routines]fread_cl4.sub 8/31/96
c read a file in cl4 format

c	warning - maximum of 4000 polygons, 900 points per polygon
c
	subroutine read_cl4(iunitcl4,level,maxpoly,maxnpts,x0,y0,scalar,
     x	ixmin,ixmax,iymin,iymax,fconv,maxd,ngeoplus2,jarray,nptsarray,
     x	jxarray,jyarray,nbpt)

c	5 geocodes for each polygon
	integer jarray(5,1)
c	number of points in each polygon
	integer nptsarray(1)
c	point array for all polygons
 	integer jxarray(4000,1), jyarray(4000,1)
c	point array for one polygon
        integer jx(900),jy(900)

c	input variables
c	iunitcl4 = input CL4 file, already open
c	level = level indicator (11,21,31,41,43)
c	maxpoly = maximum number of polygons allowed
c	maxnpts = maximum number of points in any polygon

c	output variables
c	fconv = convergence RMS error
c	maxd = no. of polygons in CL3 input file
c	ngeo = no. of geocodes plus 2
c	x0,y0,scalar,ixmin,ixmax,iymin,iymax
c	nbpt = no. of non-boundary points

c	internal variables

	character*10 dummy

	character*80 line

      character*9  info1
      parameter (info1  = 'longitude') 
      character*8  info2
      parameter (info2  = 'latitude' )
      character*19 info3
      parameter (info3  = 'number of map units')

c	read header of CLOSE4 file

c 'DEMP input file in CLOSE4 format'
	read (iunitcl4,*)

c '  1             Average density rho0 (not used)'
	read (iunitcl4,*)

c '  1e-4          Convergence RMS error'
	read (iunitcl4,*) fconv

c maxd, ngeo+2, ' # poly, # codes (seq,st,cy,...poly)'
	read (iunitcl4,*) maxd, ngeoplus2

	if (maxd.gt.maxpoly) then
		write (6,601) maxd, maxpoly
601		format (1x,i5,' polygons found. maximum is ',i5)
cVMS		call exit
		call exit(1)
	endif

	x0=-999.
	y0=-99.
	scalar=-999
	ixmin=99999
	ixmax=-99999
	iymin=99999
	iymax=-99999

c definition of map units:
c x0 = longitude in degrees at (x=0)
c y0 = latitude  in degrees at (y=0)
c scalar = number of map units in one kilometer
c x = scalar * 111.195 * (longitude-x0) * cos (y0 * PI / 180)
c y = scalar * 111.195 * (latitude-y0)
c longitude = x0 + x/(111.195 * scalar * cos (y0 * PI / 180) )
c latitude  = y0 + y/(111.195 * scalar)
c PI = 3.14159

	j1=-999
	j2=-999
	j3=-999
	j4=-999
	j5=-999

	do 402 i=1,maxd

	jarray(1,i)=0
	jarray(2,i)=0
	jarray(3,i)=0
	jarray(4,i)=0
	jarray(5,i)=0
c			read input header for next polygon

			if (level.eq.11) read(iunitcl4,511)
     x			dummy,iseq,j1,j5
			if (level.eq.21) read(iunitcl4,521)
     x			dummy,iseq,j1,j2,j5
			if (level.eq.31) read(iunitcl4,531)
     x			dummy,iseq,j1,j2,j3,j5
			if (level.gt.40) read(iunitcl4,541)
     x			dummy,iseq,j1,j2,j3,j4,j5

c511			format(a11,3i6)
c521			format(a11,4i6)
c531			format(a11,5i6)
c541			format(a11,6i6)

511			format(a7,3i7)
521			format(a7,4i7)
531			format(a7,5i7)
541			format(a7,6i7)

	jarray(1,i)=j1
	jarray(2,i)=j2
	jarray(3,i)=j3
	jarray(4,i)=j4
	jarray(5,i)=j5

c ' targ area, pop, orig area'
	read(iunitcl4,*) orig
c ' present area, target magn'
	read(iunitcl4,*) target
c npts,  '        number of points'
	read(iunitcl4,*) npts
	nptsarray(i) = npts

	if (npts.gt.maxnpts) then
		write (6,602) npts, i, maxnpts
602		format (1x,i5,' points found in polygon ',i5,
     x		'. maximum is ',i5)
cVMS		call exit
		call exit(1)
	endif

	do 721 ii=1,npts
   		read(iunitcl4,*) jx(ii),jy(ii)
		jxarray(i,ii) = jx(ii)
		jyarray(i,ii) = jy(ii)
721	continue

402	continue

c "End of Boundary points  skipped line"
	read(iunitcl4,*)
c "			   skipped line"
	read(iunitcl4,*)
c " Non-Boundary points"
	read(iunitcl4,*)

	read(iunitcl4,*) nbpt
c "    401          Non-boundary points x,y, dx,dy"

	if (nbpt.gt.0) then
		do 403 i=1,nbpt
		read(iunitcl4,*)  ix, iy
		write(8,*) ix, iy
		ixmin=min(ixmin,ix)
		ixmax=max(ixmax,ix)
		iymin=min(iymin,iy)
		iymax=max(iymax,iy)
403		continue
	endif

c "End of Data"
	read(iunitcl4,*)

c     Find and save the longitude, latitude, and map units
 24   continue
c       read until the information is found or until and eof. 
        n = 0
        read(iunitcl4,'(a)',iostat=iostat,end= 26) line
        n = index(line,info1)
        if(n.ne.0) then
c          write(2,'(i,1x,a)') n,line
		read(line,*)  x0
	endif
        n = index(line,info2)
        if(n.ne.0) then
c          write(2,'(i,1x,a)') n,line
		read(line,*)  y0
	endif
        n = index(line,info3)
        if(n.ne.0) then
c          write(2,'(i,1x,a)') n,line
		read(line,*)  scalar
	endif

        goto 24

 26   continue

c "  -72.47850     longitude in degrees at (x=0)"
c	read(iunitcl4,*) x0
c "   43.87635     latitude  in degrees at (y=0)"
c	read(iunitcl4,*) y0
c "          10 number of map units in one kilometer'
c	read(iunitcl4,*) scalar

	end