program merge_ice_models

! this program read the models created using "icesheet", and combines them to create a file that can be read into calsea.
! note that the most recent time that represents "0" must be listed first in the file "dates". Also, the dates must be in decending order.

	implicit none
	
	integer :: istat, date, number_dates, number_lat, number_long, file_counter, lat_index, long_index, format_int_value
	integer :: lat_counter, long_counter, date_counter, counter, num_long_ice, time_counter, thickness_int
	real :: min_lat, min_long, max_lat, max_long, longitude, latitude, thickness, colat, colong, read_date
	integer, dimension(:,:,:), allocatable :: model_values
	integer, dimension(:), allocatable :: dates_array
	character(len=180) :: filename, ice_format, date_format, input_format

	! this is the spacing I used for making the ice models
	real, parameter :: lat_spacing = 0.25, long_spacing=0.5 

	! this value is the oldest date in the model, assumed to be similar to present

	integer, parameter :: old_date = 200000
	integer, parameter :: out_num = 30


	character(len=180) :: outfile, folder

	call getarg(1,outfile)
	call getarg(2,folder)

	if(outfile == "" .or. folder == "") THEN
		write(6,*) "Usage: merge_ice_models [output file] [folder containing ice files]"
		stop
	endif

	! initialize the minimum and maximum coordinates. These values are in the center of Hudson Bay, so they should be far away from the extreme values

	min_lat = 61.
	min_long = -88.
	max_lat = 61.
	max_long = -88.
	number_dates = 0

	open(unit=10, file="times_to_calculate", access="sequential", status="old", form="formatted")

	initial_read: do

		read(10,*, iostat=istat) read_date

		if(istat /= 0) THEN
			exit initial_read
		endif

		date = int(read_date*1000.)

		filename = ""
		if(date == 0) THEN
			format_int_value = 1
		else
			format_int_value = int(log10(real(date))) + 1
		endif


		write(input_format,*)  "(A",len_trim(folder),",A1,I", format_int_value, ",A17)"


		write(filename,input_format) adjustl(trim(folder)),"/", date, "/model_values.txt"

		write(6,*) filename
		number_dates = number_dates+1

		! open the model file for the particular date

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

		initial_model_read: do

			read(20,*, iostat=istat) longitude, latitude, thickness

			if(istat < 0 ) THEN
				exit initial_model_read
			elseif (istat > 0) THEN
				cycle initial_model_read
			endif

			if(longitude < min_long) THEN
				min_long = longitude
			endif

			if(longitude > max_long) THEN
				max_long = longitude
			endif

			if(latitude < min_lat) THEN
				min_lat = latitude
			endif

			if(latitude > max_lat) THEN
				max_lat = latitude
			endif


		end do initial_model_read

		close(unit=20)

	end do initial_read

	! with the extreme values of the model known, read in the ice values

	number_lat = nint( (max_lat-min_lat) / lat_spacing) + 1
	number_long = nint( (max_long-min_long) / long_spacing) + 1
	number_dates = number_dates + 1 ! add in an extra date to represent the oldest time when things were zero


	allocate(model_values(number_dates, number_long, number_lat))

	allocate(dates_array(number_dates))
	model_values = 0
	dates_array = 0
	rewind(unit=10)

	! restart reading the model files, putting everything into the array

	second_read: do file_counter = number_dates - 1, 1, -1

		read(10,*, iostat=istat) read_date

		if(istat /= 0) THEN
			exit second_read
		endif
		date = int(read_date*1000.)
		dates_array(file_counter) = date

		filename = ""
		if(date == 0) THEN
			format_int_value = 1
		else
			format_int_value = int(log10(real(date))) + 1
		endif


!		write(input_format,*)  "(A",len_trim(folder),",I", format_int_value, ",A17)"
		write(input_format,*)  "(A",len_trim(folder),",A1,I", format_int_value, ",A17)"

!		write(filename,input_format) adjustl(trim(folder)), date, "/model_values.txt"
		write(filename,input_format) adjustl(trim(folder)),"/", date, "/model_values.txt"

		! open the model file for the particular date

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

		second_model_read: do

			read(20,*, iostat=istat) longitude, latitude, thickness

			if(istat < 0) THEN
				exit second_model_read
			elseif(istat > 0) THEN
				cycle second_model_read
			endif

			lat_index = nint((latitude - min_lat)/lat_spacing) + 1
			long_index = nint((longitude - min_long)/long_spacing) + 1

			thickness_int = nint(thickness)
!			if(file_counter == 1) THEN

!				model_values(:,long_index, lat_index) =thickness_int


!			else
				if (model_values(file_counter,long_index, lat_index) < thickness_int) THEN
					model_values(file_counter,long_index, lat_index) = thickness_int 
				endif
	!			write(6,*) file_counter,long_index, lat_index, model_values(file_counter,long_index, lat_index)
!			endif

		end do second_model_read
		
		close(unit=20)
	end do second_read

	close(unit=10)

	dates_array(number_dates) = old_date

!	do date_counter = 1, number_dates
!		do lat_counter = 1, number_lat, 1
!			do long_counter = 1, number_long, 1
!				if(model_values(1,long_counter,lat_counter) /=0) THEN
!				write(656,*) dates_array(date_counter), real(long_counter-1)*long_spacing+min_long, 90. - &
!					(real(lat_counter-1)*lat_spacing+min_lat), &
!					model_values(date_counter,long_counter,lat_counter), model_values(1,long_counter,lat_counter),&
!				 model_values(date_counter,long_counter,lat_counter) - model_values(1,long_counter,lat_counter)
!				endif
!			end do
!		end do
!	end do

	! make the array relative to present

	do date_counter = 2, number_dates

		model_values(date_counter,:,:) = model_values(date_counter,:,:) - model_values(1,:,:)
		
!		do lat_counter = 1, number_lat, 1
!			do long_counter = 1, number_long, 1
!				if(model_values(date_counter,long_counter, lat_counter) < -100) THEN
!					model_values(date_counter,long_counter, lat_counter) = 0
!				endif
!			end do
!		end do
	
	end do




	! as per convention in calsea, thickness values are negative

	model_values = -model_values

	! construct the ice model

	write(ice_format,*)  "(f8.3,", number_dates-1, "I6)"
	
	write(6,*) ice_format

	write(date_format,*) "(", number_dates-1, "f8.3)"
	open(unit=out_num, file=outfile, access="sequential", form="formatted", status="replace")

	write(out_num,'(A80)') adjustl(outfile)
	write(out_num,'(f8.3,1X,f8.3)') 90.-max_lat, 90.-min_lat
	write(out_num,'(f8.3,1X,f8.3)') long_spacing, lat_spacing
	write(out_num,'(I3)') number_dates - 1
	write(out_num,date_format) (real(dates_array(counter))/1000., counter=number_dates,2,-1)


	lat_loop: do lat_counter = 1, number_lat, 1

		if (sum(model_values(2:number_dates,:,lat_counter)) == 0) THEN ! don't need these values in the file
			cycle lat_loop
		else
			num_long_ice = 0
			do long_counter = 1, number_long, 1

				if (sum(model_values(2:number_dates,long_counter,lat_counter)) /= 0) THEN 
					num_long_ice = num_long_ice + 1
				endif
			end do

			colat = 90. - (real(lat_counter-1)*lat_spacing+min_lat)

			write(out_num,'(f8.3,1X,I6)') colat, num_long_ice

			do long_counter = 1, number_long, 1

				if (sum(model_values(2:number_dates,long_counter,lat_counter)) /= 0) THEN

					colong = real(long_counter-1)*long_spacing+min_long

					if(colong < 0) THEN
						colong = 360. + colong
					endif

					write(out_num,ice_format) colong, &
					     (model_values(time_counter,long_counter, lat_counter), &
!						time_counter = zero_index, num_times,1)
						time_counter = number_dates, 2,-1)
!					write(6,*) (model_values(time_counter,long_counter, lat_counter), &
!						time_counter = zero_index, num_times,1)
!						time_counter = number_dates, 2,-1)
				endif

			end do

		endif

	end do lat_loop
	
	close(out_num)
	deallocate(dates_array, model_values)

end program merge_ice_models
