load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

  begin

; print("==============> Type simulation (era/jra/nnrp1/nnrp2/snXX :")
; sim = systemfunc("read ifland; echo $ifland")


 
;-CREATE NECESSARY ARRAYS--------------------------------------------------------------------------
;--------------------------------------------------------------------------------------------------



names 	= (/ "Temperature",	"Dewpoint",	"WindSpeed",	 "WindDir",	"Sea-Level-Pres" /)
indices	= (/		6 ,		 7,     	  9,		10,		      11 /)





bias	= new( (/181,5/) , "float" )			;Shoot for 181 times and 5 variables
data	= new( (/181,11/) , "float" )			;Shoot for 181 times and 5 variables








;-LOOP THROUGH 205 STATIONS------------------------------------------------------------------------
;--------------------------------------------------------------------------------------------------
 stations = readAsciiTable("/import/wrkdir1/asemenov/process-model-output/ncdc/station-files/stations-name-205.txt", 1, "string", 0)
 coor 	  = readAsciiTable("/import/wrkdir1/asemenov/process-model-output/ncdc/station-files/stations-coords.txt", 2, "float", 0)
 


 do n=0,101-1

	lat 	= coor(n,0)
	lon 	= coor(n,1)
	station	= str_right_strip(stations(n,0))
	
	print("------ "+station+"------")

 	error = False
 

;-OPEN FILES---------------------------------------------------------------------------------------
;--------------------------------------------------------------------------------------------------

simfile        = "/import/wrkdir1/asemenov/process-model-output/ncdc/01-obtain-simulation-data/wrfout_d02_2011-03-15_00:00:00_"+station+".txt"
obsfile        = "/import/wrkdir1/asemenov/process-model-output/ncdc/station-data/"+station+".txt"


simdata	= readAsciiTable(simfile,12,"float",0)		;All have 61 times
obsdata	= readAsciiTable(obsfile,16,"float",0)		;Number of times vary but are between 8th at 12 and 11th at 00


;-GET DATA-----------------------------------------------------------------------------------------
;--------------------------------------------------------------------------------------------------
do v=0,4						;Loop through variables

	name 	= names(v)
	index	= indices(v)
	
	;print("Bias for "+name)

do nt=0,180



	dy_sim	= simdata(nt,2)
	hr_sim	= simdata(nt,3)
	var_sim	= simdata(nt,index)

	dy_obs	= obsdata(nt,2)
	hr_obs	= obsdata(nt,3)
	var_obs	= obsdata(nt,index+1)

;							print("days: "+dy_sim+", "+dy_obs)
;							print("hours: "+hr_sim+", "+hr_obs)
;
;if (v .eq. 3) then
;print("nt is "+nt)
;print("Obs: "+var_obs+" Sim: "+var_sim)
;end if


;--------------------CALCULATE BIAS----------------------------------------------------------------
;--------------------------------------------------------------------------------------------------

;     printVarSummary(dy_sim)
;     printVarSummary(dy_obs)
;     printVarSummary(hr_sim)
;     printVarSummary(hr_obs)

	if (dy_sim .eq. dy_obs .and. hr_sim .eq. hr_obs) then
	
		bias(nt,v)	= 	var_sim - var_obs
		
	else
	
		print("*****************************************")
		print("  Check the data.   Not paired in time.  ")
		error = True
	end if


delete(dy_sim)
delete(hr_sim)
delete(var_sim)
delete(dy_obs)
delete(hr_obs)
delete(var_obs)


end do							;End nt


delete(name)
delete(index)

							
end do							;End v



;## Correct for wind direction bias with magnitude greater than 180!
bias(:,3)	= where ( bias(:,3) .gt. 180,  360 - bias(:,3), bias(:,3) )
bias(:,3)  	= where ( bias(:,3) .lt. -180, 360 + bias(:,3), bias(:,3) )



do v=0,4
print(avg(bias(:,v)))
end do


;---------------------CREATE ARRAY TO BE PRINTED OUT TO FILE---------------------------------------
;--------------------------------------------------------------------------------------------------

data(:,0)	=	2011.
data(:,1)	=	03.
data(:,2)	=	obsdata(0:180,2)
data(:,3)	=	obsdata(0:180,3)
data(:,4)	=	lat
data(:,5)	=	lon
data(:,6:10)	=	bias

opt = True
;opt@fout = output_folder+"/bias/bias-"+station+"-"+sim+"-5vars.txt"
;opt@fout = output_folder+"/bias/bias-"+station+"-5vars.txt"
opt@fout = "/import/wrkdir1/asemenov/process-model-output/ncdc/03-calculate-stats/output/"+station+"-5vars.txt"


fmtx = "11f14.4"
write_matrix(data, fmtx , opt)

lat=0
lon=0
station=" "
data=0
bias=0
delete(obsdata)
delete(simdata)

end do							;Station loop

if (error .eq. True) then
	print("Go back and look at output. There is some error.")
end if

end