;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	takes one input file, and returns the name of 3 input files that
;;	it is split into (though it will actually create as many input 
;;	files as are in the binary file originally)
;;
;;	9/12/00 also creates a seperate file that stores all lat, lon, and
;;	date data with the station name in ascii format
;;
;;	1/27/01	if it finds a prec, tmin, and tmax file with the same name it 
;;	calls centclimate with those files, to create a .wth file no longer
;;	returns file names
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro splitfile, fullinput


type=strarr(6)
type(0)='prec'
type(1)='Error1'
type(2)='Error2'
type(4)='tmax'
type(5)='tmin'

filelist=strarr(6)
filelist(0)='1' & filelist(4)='2'

curpoint=0l
count=0
name='                       '
lat='    '
lon='     '
mo=0b
yr=0b
year = 0

openw, keyun, /get, string(fullinput + '.key')
openr, un, /get, fullinput

while not eof(un) do begin

	yr=0b
	name='                       '
	lat='    '
	lon='     '
	point_lun, un, curpoint+45
	readu, un, lat
	readu, un, lon
	point_lun, un, curpoint+55
	readu, un, mo
	readu, un, yr

	if mo ne 1 then yr=yr+1
	year=fix(yr)+1800

	c=0b
	point_lun, un, curpoint+54
	readu, un, c

	point_lun, un, curpoint+2
	readu, un, name
	print, name, lat, '  ', lon, year
	name=strcompress(name, /remove_all)

	out=name+'.wth'

	if c gt 5 then c=2
	filelist(c) = name

	name=string(name+type(c)+'.bin')

	point_lun, un, curpoint+65
	i=0
	readu, un, i &	byteorder, i

	if i le 0 or i gt 200 then begin 
		print, 'ERROR', i, name
		curpoint=curpoint+819l
;		break
	endif else begin
		data=bytarr(819l*i)
		point_lun, un, curpoint

		readu, un, data
		openw, oun, /get, name
		writeu, oun, data
		close, oun &	free_lun, oun

		curpoint=long(curpoint+(819l*i))

		if filelist(0) eq filelist(4) and filelist(4) eq filelist(5) then begin
			;wait, 1 ;;ensure we finished writing?
			print, 'Creating .wth file for ', filelist(0), ' with ', i, ' years'
			printf, keyun, out,' ',lat,'  ',lon,' ',year

			centclimate, string(filelist(0),'prec.bin'), $
					string(filelist(4),'tmax.bin'), $
					string(filelist(5), 'tmin.bin'), out
		endif
	endelse
endwhile

close, keyun, un
free_lun, keyun, un

end;;splitfiles


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Sets the initial file pointer positition and gets available year
;;	info.  
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function initfile, un, startyr, totyrs


point_lun, un, 55

;;read in the first year available
startmo=0b
startyr=0b
endmo=0b
endyr=0b
years=0

readu, un, startmo
readu, un, startyr
readu, un, endmo
readu, un, endyr

point_lun, un, 65
readu, un, years
byteorder, years

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;set the startpointer
;;	if the first month is not january than skip the first year (go to the 
;;	first complete year
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

if (endyr-startyr+1 gt years) then begin
	print, endyr, startyr, years, endmo, startmo
	free_lun, un
	return, -1
endif

case startmo of
	1 : begin 
		startyr=fix(startyr)+1800
		point_lun, un, 75
	    end
	else : begin
		startyr=fix(startyr)+1801
		point_lun, un, 819+75
	    end
endcase



;; need the end year to calculate total years so do a similar trick to find the
;; last year that has all twelve months? (though we can take a year that finished in Oct?)
if endmo gt 10 then begin
	endyr=fix(endyr)+1800
endif else endyr=fix(endyr)+1799

;;calculate the number of years we have to read in
	totyrs=endyr-startyr+1
print, startyr, endyr, totyrs, startmo, endmo

return, 0

end;;initfile



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;sets all the file pointers to the first year they all have info for
;;	probably these should all be the same, but just in case. 
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro setpointers, un1, un2, un3, st1, st2, st3

;;let us know if and why we actually have to do this
	print, st1, st2, st3

	stlast=max([st1, st2, st3])

;;calculateshow many years to offset each unit
	offset1=stlast-st1
	offset2=stlast-st2
	offset3=stlast-st3
	if offset1 ne offset2 or offset2 ne offset3 then print, 'ERROR',offset1,offset2,offset3

;;finds the current position to offset from
	cur1=(fstat(un1)).cur_ptr
	cur2=(fstat(un2)).cur_ptr
	cur3=(fstat(un3)).cur_ptr

;;calculates the new position
	new1=cur1+(819l*offset1)
	new2=cur2+(819l*offset2)
	new3=cur3+(819l*offset3)

;;sets the file pointer positions
	point_lun, un1, new1
	point_lun, un2, new2
	point_lun, un3, new3
end;setpointers



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	average all days in each month and return one years monthly ave
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getmonths, curyr, averageit

	yearbymonth=fltarr(12)

	for mo=0, 11 do begin
		tmpmo=curyr(*,mo)
		index=where(tmpmo ge -100)

		if n_elements(index) le 20 then begin yearbymonth(mo)=-99.99
		endif else if averageit then begin yearbymonth(mo)=mean(tmpmo(index))
		endif else yearbymonth(mo)=total(tmpmo(index))

	endfor

	return, yearbymonth
end;;getmonths



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	reads in ncdc binary data and outputs a data file that century
;;	can read
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro centclimate, inprec, intmax, intmin, outFile

if n_elements(inprec) eq 0 then begin
	print, ''
	Print, '* * * * * * * * * * * * * * * * * * *'
	print, ''
	print, 'Calling sequence is :'
	Print, 'centclimate, precip file, tmax file, tmin file, output file'
	print, ''
	print, 'or if all of your climate data is in one file then :'
	print, 'centclimate, inputfile'
	print, ''
	Print, '* * * * * * * * * * * * * * * * * * *'
	print, ''
	return
endif


;; This is probably the standard first procedure now.  takes a large binary
;;  file and splits it into three binary files, then calls centclimate on those
if n_elements(intmax) eq 0 then begin 
	splitfile, inprec
	return
endif

;;open files
	openr, prcun, /get, inprec
	openr, tmxun, /get, intmax
	openr, tmnun, /get, intmin

;;set initial filepointer position and get year info
if initfile(prcun, precstart, precyears) eq -1 then begin 
	free_lun, prcun, tmxun, tmnun
	return
endif
if initfile(tmxun, tmaxstart, tmaxyears) eq -1 then begin 
	free_lun, prcun, tmxun, tmnun
	return
endif
if initfile(tmnun, tminstart, tminyears) eq -1 then begin 
	free_lun, prcun, tmxun, tmnun
	return
endif


;;if necessary, reset the file pointers so they all point to the same year, assumes
;; overlapping coverage
if (precstart ne tmaxstart) or (precstart ne tminstart) then $
	setpointers, prcun, tmxun, tmnun, precstart, tmaxstart, tminstart


totyears=min([tmaxyears, tminyears, precyears])
styear=max([precstart, tmaxstart, tminstart])

;;set up arrays to read number of years *number of months
	tmax=fltarr(12, totyears, /nozero)
	tmin=fltarr(12, totyears, /nozero)
	prec=fltarr(12, totyears, /nozero)


tmpyr=intarr(31,12)
tmpmo=intarr(31)

tmaxbase=(fstat(tmxun)).cur_ptr
tminbase=(fstat(tmnun)).cur_ptr
precbase=(fstat(prcun)).cur_ptr
print, tmaxbase, tminbase, precbase

;; read data from files
for i=0, totyears-1 do begin
	readu, tmxun, tmpyr
	byteorder, tmpyr
	tmax(*,i)=getmonths(tmpyr, 1)
	point_lun, tmxun, tmaxbase+(819l*(i+1))

	readu, tmnun, tmpyr
	byteorder, tmpyr
	tmin(*,i)=getmonths(tmpyr, 1)
	point_lun, tmnun, tminbase+(819l*(i+1))

	readu, prcun, tmpyr
	byteorder, tmpyr
	prec(*,i)=getmonths(tmpyr, 0)
	point_lun, prcun, precbase+(819l*(i+1))
endfor

index=where(tmax eq -99.99)
tmax=(tmax-32)*(5./9.)	;; convert Fahrenheit to Celsius
if index(0) ne -1 then tmax(index)=-99.99

index=where(tmin eq -99.99)
tmin=(tmin-32)*(5./9.)	;; convert Fahrenheit to Celsius
if index(0) ne -1 then tmin(index)=-99.99

index=where(prec eq -99.99)
prec=prec*2.54/100		;; convert inches to cm
if index(0) ne -1 then prec(index)=-99.99

;; close files
	close, tmxun, tmnun, prcun
	free_lun, tmxun, tmnun, prcun


;;open our output file
	openw, oun, /get, outFile


;;for each year, create a string for each variable with all the months
for i=0, totyears-1 do begin

	year=string(format='(I4)', (i+styear))

	precdata=string(format='(F7.2)', prec(*,i))
	curline=string('prec'+'  '+year+precdata(0)+precdata(1)+precdata(2) $
	+precdata(3)+precdata(4)+precdata(5)+precdata(6)+precdata(7) $
	+precdata(8)+precdata(9)+precdata(10)+precdata(11))
	printf, oun, curline

	mindata=string(format='(F7.2)', tmin(*,i))
	curline=string('tmin'+'  '+year+mindata(0)+mindata(1)+mindata(2) $
	+mindata(3)+mindata(4)+mindata(5)+mindata(6)+mindata(7) $
	+mindata(8)+mindata(9)+mindata(10)+mindata(11))
	printf, oun, curline


	maxdata=string(format='(F7.2)', tmax(*,i))
	curline=string('tmax'+'  '+year+maxdata(0)+maxdata(1)+maxdata(2) $
	+maxdata(3)+maxdata(4)+maxdata(5)+maxdata(6)+maxdata(7) $
	+maxdata(8)+maxdata(9)+maxdata(10)+maxdata(11))
	printf, oun, curline

endfor
close, oun
free_lun, oun

end