%%%=======================================================
%FileName: 	downscale_100percentiles.m
%Description: 	This script performs the statistical downscaling of RCM data
%		utilizing SDBC. 
%
%		SDBC is composed of the following steps: 
%		0) Grid adjustment: interpolation of curvlinear coordinates to rectangular grid
%		1) Epoch adjustment: removing the climatological differences in gcm data
%		2) Bias Correction:
%		3) Epoch adjustment: adding back in climatological differences in gcm data
%		
%Code Authors:
%		Original: Katherine Hegewisch(khegewisch@uidaho.edu)
%		Date:	  February, 2016
%		Version: 1 
%		Updates: 
%
function []=downscale_100percentiles(VAR_TO_RUN,MODEL_NUM,EXP_NUM,HIST_NUM,outputdir,writehistorical,writefuture,CMIP_NUM,sdMode);

	%sdMode

	%%%=======================================================
	%		PARAMETERS
 	%%%=======================================================
	 [outputdir,method_name,region_string,MODEL_NAME,EXP_NAME,HIST_NAME,VAR_NAME,VAR_MULT_OR_ADD,...
		OBS_COARSE_DIR,OBS_FINE_DIR,GCM_DIR,...
		COARSE_XX,COARSE_YY,...
                FINE_XX_TARGET,FINE_YY_TARGET,HISTORICAL_YEARS,FUTURE_YEARS,...
                EA_YR_WINDOW,EA_DAY_WINDOW,...
                BC_SEMI_DAY_WINDOW_MULT,BC_SEMI_DAY_WINDOW_ADD,...
                MACA_SEMI_DAY_WINDOW,MACA_NUM_BEST_MODELS,GRID_CHUNK_PTS,units,...
		RSDS_DIR,BC_MAX_PR_RATIO,EA_VARIABLES,DRY_DAY_THRES_COARSE,DRY_DAY_THRES_FINE,...
		VAR_OF_PREC,ERRORYES]=...
                parameters(MODEL_NUM,EXP_NUM,VAR_TO_RUN,outputdir,CMIP_NUM);
	DAYS = 365;

	f5=find(VAR_TO_RUN==VAR_OF_PREC);
	f10=find(VAR_TO_RUN==10);
	f1=find(VAR_TO_RUN==1);
	f2=find(VAR_TO_RUN==2);

	NUM_DOWNSCALE_VAR = length(VAR_TO_RUN);
	HISTORICAL_REF_YEAR = min(HISTORICAL_YEARS);
	FUTURE_REF_YEAR = min(FUTURE_YEARS);
	HISTORICAL_YEAR_INDICES=(HISTORICAL_YEARS-HISTORICAL_REF_YEAR+1)'; 
	FUTURE_YEAR_INDICES=(FUTURE_YEARS-FUTURE_REF_YEAR+1)' +max(HISTORICAL_YEAR_INDICES);
	YEARS_HISTORICAL=length(HISTORICAL_YEAR_INDICES);
	YEARS_FUTURE=length(FUTURE_YEAR_INDICES);
	YEARS_GCM=length(FUTURE_YEAR_INDICES)+length(HISTORICAL_YEAR_INDICES);


	%%%=======================================================
	fprintf('Downscaling: %s: %s\n',MODEL_NAME{MODEL_NUM},EXP_NAME{EXP_NUM} )
	%%%=======================================================
  	%BREAK REGION UP INTO SPATIAL CHUNKS
	%%%=======================================================
	   [FULL_FINE_XX,FULL_FINE_YY,X_ranges,Y_ranges]=...
                break_curvilinear_grid(OBS_FINE_DIR,FINE_XX_TARGET,FINE_YY_TARGET,...
                GRID_CHUNK_PTS,1);

	%%%=======================================================
	%      RUN OVER ALL VARIABLES  
	%%%=======================================================
        for var=1:NUM_DOWNSCALE_VAR
	%['variable=',num2str(var)]
           kk=VAR_TO_RUN(var);
	  %%%=======================================================
	   %     SAVE FILE 
	   %%%=======================================================
	   if(writehistorical==1);
	      exp=HIST_NUM; %historical
	      historical=save_file(FULL_FINE_XX,FULL_FINE_YY,...
                HISTORICAL_YEAR_INDICES,...
                method_name,MODEL_NAME,MODEL_NUM,HIST_NAME,exp,VAR_NAME,kk,outputdir,DAYS,...
		units{kk},HISTORICAL_YEARS);
	   end

	   if(writefuture==1);
		   exp=EXP_NUM; %future
		   future=save_file(FULL_FINE_XX,FULL_FINE_YY,...
			FUTURE_YEAR_INDICES,...
			method_name,MODEL_NAME,MODEL_NUM,EXP_NAME,exp,VAR_NAME,kk,outputdir,DAYS,...
			units{kk},FUTURE_YEARS);
	   end;

	   %%%=======================================================
	   %    RUN OVER GRID FILLING IN DATA 
	   %%%=======================================================
 	   for j=1:length(Y_ranges)-1;
		if(j~=length(Y_ranges)-1)
	      		 Y_range = [Y_ranges(j):1:Y_ranges(j+1)-1]; %lat index range
		else
	      		 Y_range = [Y_ranges(j):1:Y_ranges(j+1)]; %lat index range
		end
		for  i=1:length(X_ranges)-1;
[j i]
			if(i~=length(X_ranges)-1)
                        	X_range = [X_ranges(i):1:X_ranges(i+1)-1]; %lon index range
			else
                        	X_range = [X_ranges(i):1:X_ranges(i+1)]; %lon index range
			end
			FINE_YY=FULL_FINE_YY(Y_range);
			FINE_XX=FULL_FINE_XX(X_range);

			%arbitrary here for curvlinear
			NUM_FINE_XX=size(FINE_XX,2);        
			NUM_FINE_YY=size(FINE_YY,2);        

			%================================================
			%        GET OBS FINE DATA
			%================================================
			[obsdataFine,YEARS_OBS,FINE_XX,FINE_YY]=get_obsdataFine(OBS_FINE_DIR,...
			FINE_XX,FINE_YY,DAYS,kk,units,region_string);
			if(kk==VAR_OF_PREC);
				f = find(obsdataFine<DRY_DAY_THRES_FINE); obsdataFine(f) = 0.0;
			end
 		 %================================================
		%only do the following if there are non-NaNs in obs
		 fobs=find(~isnan(obsdataFine));
                if(~isempty(fobs) );
			  [gcmdataFine,YEARS_OBS]=get_curvilinear_rcmdata(...
				GCM_DIR,char(VAR_NAME(kk)),char(MODEL_NAME(MODEL_NUM)),...
				char(EXP_NAME(EXP_NUM)),region_string,FINE_XX,FINE_YY,units,...
				HISTORICAL_YEARS,FUTURE_YEARS,char(HIST_NAME(HIST_NUM)));
			
			 obsdataFine = permute(obsdataFine,[2 1 3 4]);
			 gcmdataFine = permute(gcmdataFine,[2 1 3 4]);
			%================================================
			%        BIAS CORRECT FINE GCM DATA W 100 PERCENTILES
			%================================================
			if(isempty(f10)); %don't bias correct wind speed.. use WRF instead
				 obsdataFine = permute(obsdataFine,[3 4 1 2]);
				 gcmdataFine = permute(gcmdataFine,[3 4 1 2]);
	
				 [gcmdataFine(:,:,:,:)]=get_bias_correction_100percentiles(obsdataFine,gcmdataFine, ...
                                               FINE_XX,FINE_YY,NUM_FINE_XX,NUM_FINE_YY,...
                                               HISTORICAL_YEAR_INDICES,FUTURE_YEAR_INDICES,...
                                                VAR_MULT_OR_ADD(kk),kk,...
                                                BC_SEMI_DAY_WINDOW_MULT,BC_SEMI_DAY_WINDOW_ADD,...
                                                BC_MAX_PR_RATIO,DRY_DAY_THRES_FINE);

			 	obsdataFine = permute(obsdataFine,[3 4 1 2]);
				 gcmdataFine = permute(gcmdataFine,[3 4 1 2]);

			end;
			obsdataFine = permute(obsdataFine,[2 1 3 4]);
			 gcmdataFine = permute(gcmdataFine,[2 1 3 4]);
			%================================================
		  	%      FINAL ADJUSTMENT 
		   	%==========================================
			  gcmdataFine=check_bounds_fine(gcmdataFine,kk,'fine',DRY_DAY_THRES_FINE,...
					RSDS_DIR,Y_range,X_range);

		        %=====================================
			%    CHECK on DTR 
			%=====================================
			%check on dtr<0
			if(kk==2);
				%have tasmax saved already... alter tasmin if tasmax<tasmin to tasmin =taxmax
			     if(writehistorical==1);
	   			exp=HIST_NUM;%historical
				historical_tmax=read_file(method_name,MODEL_NAME,MODEL_NUM,HIST_NAME,exp,...
				VAR_NAME,1,outputdir);
				%temp_mx = double(historical_tmax.data(:,:,Y_range,X_range)); %for curvilinear
				temp_mx = double(historical_tmax.data(:,:,X_range,Y_range));
				temp_mx =permute(temp_mx,[3 4 1 2]);
				temp_mn = gcmdataFine(:,:,:,HISTORICAL_YEAR_INDICES);
				f =find(temp_mx<temp_mn);
				if(~isempty(f));
					temp_mn(f) = temp_mx(f);
					gcmdataFine(:,:,:,HISTORICAL_YEAR_INDICES) = temp_mn;      
				end
			    end

			     if(writefuture==1);
	   			exp=EXP_NUM; %future
           			future_tmax=read_file(method_name,MODEL_NAME,MODEL_NUM,EXP_NAME,exp,...
				VAR_NAME,1,outputdir);
				%temp_mx = double(future_tmax.data(:,:,Y_range,X_range)); %for curvilinear
				temp_mx = double(future_tmax.data(:,:,X_range,Y_range));
				temp_mx =permute(temp_mx,[3 4 1 2]);
				temp_mn = gcmdataFine(:,:,:,FUTURE_YEAR_INDICES);
				f =find(temp_mx<temp_mn);
				if(~isempty(f));  
					temp_mn(f) = temp_mx(f);
					gcmdataFine(:,:,:,FUTURE_YEAR_INDICES) = temp_mn;      
				end
			   end;

			end;
			   %==========================================
			   %        SAVE DATA
			   %==========================================
			   gcmdataFine=permute(gcmdataFine,[3 4 1 2]);  %days,yrs,lat,lon
			   f = find(~isnan(gcmdataFine));
			   if(~isempty(f))
			       if(writehistorical==1);
				   %write historical gcm downscaled data 
				%historical.data(:,:,Y_range,X_range) =  ...   %curvilinear change
				historical.data(:,:,X_range,Y_range) =  ...
					single(gcmdataFine(:,HISTORICAL_YEAR_INDICES,:,:));
				   fprintf('Saved historical data, grid lat %i lon %i\n',j,i)
				end
				if(writefuture==1);
				   %write future gcm downscaled data
				  %future.data(:,:,Y_range,X_range) = ... %curvilinear change
				  future.data(:,:,X_range,Y_range) = ...
					 single(gcmdataFine(:,FUTURE_YEAR_INDICES,:,:));
				  fprintf('Saved future data too. \n')
				end;
	 		   end
			   %==========================================
		else;    
			'no obs points'
		end %fobs
		%==========================================
		%==========================================

		end %i
	     end %j

	end %var
end %of *.m
