	PROGRAM EXTRACT_DEF

! This program reads an deformation file created by calhornew and outputs 
! each time slice into an XYZ file that can be used to plot the ice sheet
! at various times

	IMPLICIT NONE

	CHARACTER(len=255) :: infoline, argname, filename, outname, outnamecomb, dummy, outlog, format1

	CHARACTER(len=1)	:: ext1
	CHARACTER(len=2)	:: ext2
	CHARACTER(len=3)	:: ext3, extfin
	integer, parameter :: max_times = 300, max_number_coords = 720
	real			:: ncolat, scolat, latsp, longsp, latitude, longitude, colat, colong
	integer		:: numtimes, outnamelength, numlong, k, l, m, n, j, numcolat, counter
	logical		:: tester

	real, dimension(max_times)	:: times

	real, dimension(max_number_coords, max_number_coords, max_times)	:: ivalues

	integer, dimension(max_number_coords)	:: longvalue

	


	call getarg(1, filename)
!	read(argname,*) filename

	call getarg(2,argname)
	read(argname,*) outname

	write(6,*) "extracting", filename
	write(6,*) "time", outname

! find the length of the output name

	tester=.true.
	outnamelength=0

	DO WHILE (tester)
	  outnamelength=outnamelength+1
	
	  if (outname(outnamelength:outnamelength) .eq. ' ') THEN
	   tester=.false.
	  ENDIF



	END DO

 

	call read_def(filename, times, numcolat, numlong, numtimes, latsp, longsp, colat, colong, ivalues,&
		max_times,max_number_coords)

! write out results


	outlog = outname(1:outnamelength-1) // "/" // outname(1:outnamelength-1) // "." // "log"
	open (unit=12, file=outlog, form="formatted", access="sequential", status="new")

	write(12,*) latsp, longsp
	counter = 1
	DO n =numtimes, 1, -1

		write(format1,'(A2,I1,A1)') "(i", int(log10(real(counter)))+1, ")"

		write(dummy,format1) counter


		read (dummy,*) extfin

		outnamecomb = outname(1:outnamelength-1) // "/" // outname(1:outnamelength-1) // "." // extfin

!		write(6,*) n, times(n)

		write(12,*) counter, times(n)

		open (unit=11, file=outnamecomb, access="sequential", form="formatted", status="new")

		DO m = 1, numcolat, 1

			latitude = real(m-1) * latsp + colat

			DO j=1, numlong, 1

				longitude = real(j-1) * longsp + colong

				if(longitude < 180.) THEn
					write(11, '(F8.3,1X,F7.3,1X,F6.0)') longitude, 90.-latitude,  ivalues(m, j, n)

				else

					write(11, '(F8.3,1X,F7.3,1X,F6.0)') longitude-360, 90.-latitude,  ivalues(m, j, n)

				endif


			enddo

		enddo

		close (unit=11)

		counter = counter + 1

	enddo



	close (unit=12)

	STOP
	END
		
