program site_to_txt

! note that in order for this program to work, you have to use files created with the modified version of calhornew, which outputs
! the header values to tell you the latitude and longitude range and spacing

! the calhornew program outputs starting from the northwest corner of the grid, and increments by longitude first

	implicit none
	
! xlatn,wlong,xlats,elong,dlat,dlong,nlat,nlon
	real :: north_lat, west_long, south_lat, east_long, lat_spacing, long_spacing, latitude, longitude, x, y
	integer :: number_times, ierr, counter, number_lat, number_long, long_counter, lat_counter, record_number, hundred_index
	character(len=180) :: filename, ice_model, earth_model, parameters_file, format_statement, outname, time_string, dummy
	real :: time
	integer :: time_index
	character(len=180), parameter :: outfile = "deformation.txt"

	real, allocatable, dimension(:) :: time_array
	real, allocatable, dimension(:,:,:) :: deform_values

	integer, parameter :: def_unit = 10,  out_unit = 30


	call getarg(1, filename)
	call getarg(2, dummy)
	read(dummy,*) time

	if(filename == "" .or. dummy == "") THEN
		write(6,*) "Usage: site_to_bin [filename] [ice model] [earth model]"
		stop
	endif

	open(unit=def_unit, file=filename, access="sequential", form="formatted", status="old")



	read(def_unit,*) number_times, north_lat, west_long, south_lat, east_long, lat_spacing, long_spacing, number_lat,&
		number_long


	allocate(time_array(number_times), stat=ierr)
	if(ierr /=0) THEN
		write(6,*) "problem allocating time_array", ierr
		stop
	endif

	allocate(deform_values(number_times, number_long, number_lat), stat=ierr)
	if(ierr /=0) THEN
		write(6,*) "problem allocating deform_values", ierr
		stop
	endif

	write(format_statement,*) "(10F6.2)"



	read(def_unit,format_statement) (time_array(counter), counter=1,number_times, 1)
	time_index = 0
	do counter = 1, number_times



		if(time_array(counter)*1000. == time) THEN
			time_index = counter
		endif
	end do

	if (time_index == 0) THEN
		write(6,*) "could not find time in def_mod file:", time
		stop
	endif


	do lat_counter = 1, number_lat
		do long_counter = 1, number_long

			read(def_unit, *) longitude, latitude, (deform_values(counter, long_counter, lat_counter), &
				counter=1, number_times, 1)

		end do



	end do




	open(unit=out_unit, file=outfile, access="sequential", form="formatted", status="replace")

	do lat_counter = 1, number_lat
		do long_counter = 1, number_long

			latitude = north_lat - real(lat_counter-1)*lat_spacing
			longitude =real(long_counter-1)*long_spacing -  west_long

			if(longitude >= 180) then
				longitude = longitude -360
			endif
			write(out_unit,*) longitude, latitude,deform_values(time_index, long_counter, lat_counter) - &
				deform_values(1, long_counter, lat_counter)


		end do

	end do


	close(out_unit)




	deallocate(time_array)
	deallocate(deform_values)

	close(unit=def_unit)

	write(6,*) "successfully extracted deformation"

contains



! calculate pi
real function pi()

	pi = 4.0 * atan(1.0)

end function pi

!************************************************************************************

! convert to radians
real function ctr()

	ctr = pi() / 180.0

end function ctr

end program site_to_txt
