clc
clear all
close all

coordPath = ['/discover/nobackup/jkolassa/AMSRE_ASCAT_DA_experiment/AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02/output/SMAP_EASEv2_M09_CONUS/rc_out/'];


startYear = 2008;
stopYear = 2008;
startMonth = 1;
stopMonth = 12;
collectionID = 10;
nrOutVar = 46;
varIndex = 19;


%--------------------------------------------------------------------------
% COMPUTATIONS
%--------------------------------------------------------------------------

% account for data availability

startMonths = ones((stopYear-startYear+1),1);
startMonths(1) = startMonth;
stopMonths = 12*ones((stopYear-startYear+1),1);
stopMonths(end) = stopMonth;
totalMonths = sum(stopMonths-startMonths+1);

% account for effect of leap year

leapYears = 1992:4:2013;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];

% get variable names and units according to collection ID and number of
% outputs variables

[data_tag, data_units] = get_data_tag(collectionID, nrOutVar, nrOutVar);

% read no DA coordinate data

tileCoordFile = [coordPath 'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02.ldas_tilecoord.bin'];
[tile_coord ] = read_tilecoord( tileCoordFile, 0 );

tileGridFile = [coordPath 'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02.ldas_tilegrids.bin'];
[ tile_grid_g, tile_grid_d ] = read_tilegrids( tileGridFile );

obsPoints = length(tile_coord.com_lat);

plot_lats = tile_coord.com_lat(1:obsPoints);
plot_lons = tile_coord.com_lon(1:obsPoints);


yearIndex = 0;
figureIndex = 0;



% read data for the entire period

for y = startYear%:stopYear
    y
    
    yearStr = num2str(y);
    yearIndex = yearIndex + 1;
    
    if any(y==leapYears)
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end
    
    for m = startMonths(yearIndex):stopMonths(yearIndex)
        m
        
        if m < 10
            monthStr = ['0' num2str(m)];
        else
            monthStr = num2str(m);
        end
        
        for d = 20%:monthLength(m)
            
            figureIndex = figureIndex + 1;
            
            if d < 10
                dayStr = ['0' num2str(d)];
            else
                dayStr = num2str(d);
            end
            
            xIndex = 0;
            
            for i = 2:5
                
                xIndex = xIndex+1;
                
                xcorr = ['0' num2str(i)];
                
                
                catPath = ['/discover/nobackup/jkolassa/AMSRE_ASCAT_DA_experiment/AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr' xcorr '/output/SMAP_EASEv2_M09_CONUS/cat/ens_avg/'];
                anaPath = ['/discover/nobackup/jkolassa/AMSRE_ASCAT_DA_experiment/AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr' xcorr '/output/SMAP_EASEv2_M09_CONUS/ana/ens_avg/'];
                
                
                obsFcstAna_File_6 = [anaPath 'Y' yearStr '/M' monthStr '/' ...
                    'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr' xcorr '.ens_avg.ldas_ObsFcstAna.'...
                    yearStr ...
                    monthStr dayStr '_0600z.bin'];
                
                obsFcstAna_File_9 = [anaPath 'Y' yearStr '/M' monthStr '/' ...
                    'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr' xcorr '.ens_avg.ldas_ObsFcstAna.'...
                    yearStr ...
                    monthStr dayStr '_0900z.bin'];
                
                [date_time_6,              ...
                    obs_assim_6,              ...
                    obs_species_6,            ...
                    obs_tilenum_6,            ...
                    obs_lon_6,                ...
                    obs_lat_6,                ...
                    obs_obs_6,                ...
                    obs_obsvar_6,             ...
                    obs_fcst_6,               ...
                    obs_fcstvar_6,            ...
                    obs_ana_6,                ...
                    obs_anavar_6              ...
                    ] =                      ...
                    read_ObsFcstAna(obsFcstAna_File_6);
                
                [date_time_9,              ...
                    obs_assim_9,              ...
                    obs_species_9,            ...
                    obs_tilenum_9,            ...
                    obs_lon_9,                ...
                    obs_lat_9,                ...
                    obs_obs_9,                ...
                    obs_obsvar_9,             ...
                    obs_fcst_9,               ...
                    obs_fcstvar_9,            ...
                    obs_ana_9,                ...
                    obs_anavar_9              ...
                    ] =                      ...
                    read_ObsFcstAna(obsFcstAna_File_9);
                
                lonVec = [obs_lon_6; obs_lon_9];
                latVec = [obs_lat_6; obs_lat_9];
                obsVec = [obs_obs_6; obs_obs_9];
                fcstVec = [obs_fcst_6; obs_fcst_9];
                obsStd = [obs_obsvar_6; obs_obsvar_9];
                obsStd = sqrt(obsStd);
                fcstStd = [obs_fcstvar_6; obs_fcstvar_9];
                fcstStd = sqrt(fcstStd);
                
                obsIndex = xIndex;
                fcstIndex = 4+xIndex;
                
                figure(figureIndex)
                subplot(2,4,obsIndex)
                hist(obsStd,20)
                axis([0 0.05 0 40])
                title(['obs std (xcorr = 0.' num2str(i) ')'])
%                 map_plot_filled_sub(latVec,lonVec,obsStd,['obs std (xcorr = 0.' num2str(i) ')'],...
%                     'soil moisture [m^3 m^{-3}]',[-0.1 0.1])
                subplot(2,4,fcstIndex)
                hist(fcstStd,20)
                axis([0 0.05 0 40])
                title(['fcst std (xcorr = 0.' num2str(i) ')'])
%                 map_plot_filled_sub(latVec,lonVec,fcstStd,['fcst std (xcorr = 0.' num2str(i) ')'],...
%                     'soil moisture [m^3 m^{-3}]',[-0.1 0.1])
                
                
                
            end
            
        end
    end
end



% fprintf(ofp, '%10d%10d%13.5e%13.5e%13.5e', ...
% % 	  N_tile, N_poly, edge_min, edge_max, edge_dx );
% 
% format long
%   
% scalePath = '/discover/nobackup/projects/gmao/ssd/land/l_data/AMSR_ASCAT_NN/AMSR_ASCAT/scaling/';
% scaleCoeffFile = [scalePath 'AMSRE_ASCAT.cdf.poly.hscale_1.25_2007-2011.all_monthsSM_1.dat'];
% coeff_data = dlmread(scaleCoeffFile);
% % fid = fopen(scaleCoeffFile);
% % coeff_data = fscanf(fid,'%10d %10d %13.5e %13.5e %13.5e','delimiter','\n');
% % fclose(fid);
% % 
% % statFile = [scalePath 'AMSRE_ASCAT_DA_noDA_AMSRE_ASCAT.stats.hscale_1.25_2007-2011.all_monthsSM_1.dat'];
% % stat_data = dlmread(statFile);
% 
% % fid = fopen(statFile);
% % stat_data = fscanf(fid,'%10d%7d%13.5e%13.5e%13.5e%13.5e%13.5e%8.3f%8.3f');
% % fclose(fid);
% amsre_ascat_coord_path = '/discover/nobackup/jkolassa/AMSRE_ASCAT_DA_experiment/';
% amsre_ascat_coord_path2 = '/discover/nobackup/jkolassa/AMSRE_ASCAT_DA_experiment/AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr0/output/SMAP_EASEv2_M09_CONUS/rc_out/';
% 
% amsre_ascat_tileCoordFile = [amsre_ascat_coord_path 'AMSRE_ASCAT_DA_noDA_AMSRE_ASCAT.ldas_tilecoord.bin'];
% [amsre_ascat_tile_coord ] = read_tilecoord( amsre_ascat_tileCoordFile, 0 );
% 
% amsre_ascat_tileCoordFile2 = [amsre_ascat_coord_path2 'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr0.ldas_tilecoord.bin'];
% [amsre_ascat_tile_coord2 ] = read_tilecoord( amsre_ascat_tileCoordFile2, 0 );
% 
% amsre_ascat_tileGridFile = [amsre_ascat_coord_path 'AMSRE_ASCAT_DA_noDA_AMSRE_ASCAT.ldas_tilegrids.bin'];
% [ amsre_ascat_tile_grid_g, amsre_ascat_tile_grid_d ] = read_tilegrids( amsre_ascat_tileGridFile );
% 
% amsre_ascat_Points = length(amsre_ascat_tile_coord.com_lat);
% % 
% amsre_ascat_plot_lats = amsre_ascat_tile_coord.com_lat(1:amsre_ascat_Points);
% amsre_ascat_plot_lons = amsre_ascat_tile_coord.com_lon(1:amsre_ascat_Points);
% 
% tileNum = amsre_ascat_tile_coord2.tile_id;
% tileIndex = ismember(coeff_data(:,1),tileNum);
% 
% new_coeffData = [coeff_data(tileIndex,:)];
% 
% testFile = 'test.dat';
% dlmwrite(testFile,new_coeffData,'%6d%13.5f%13.5f%13.5f%13.5f%13.5f%13.5f%13.5f%13.5f');
% testData = load(testFile);
% 
% 
% saveFile = [scalePath 'AMSRE_ASCAT_small.cdf.poly.hscale_1.25_2007-2011.all_monthsSM_1.dat'];
% fid = fopen(saveFile);
% fprintf(fid,'%10d%10d%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e\n',coeff_data(1,:));
% fprintf(fid,'%10f%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e',new_coeffData);
% fclose(fid);
% 
% 
% data = load(saveFile);
% 
% % fid = fopen(saveFile);
% header = fscanf(fid,'%10d%10d%13.5e%13.5e%13.5e\n',1);
% test_coeff_data = fscanf(fid,'%10d%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e%13.5e\n');
% fclose(fid);
% % 
% readPath = '/discover/nobackup/jkolassa/AMSRE_ASCAT_DA_experiment/AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02/output/SMAP_EASEv2_M09_CONUS/ana/ens_avg/Y2007/M09/';
% 
% readFile1 = [readPath 'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02.ens_avg.ldas_ObsFcstAna.20070929_0300z.bin'];
% readFile2 = [readPath 'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02.ens_avg.ldas_ObsFcstAna.20070929_0600z.bin'];
% readFile3 = [readPath 'AMSRE_ASCAT_DA_NN_AMSRE_ASCAT_xcorr02.ens_avg.ldas_ObsFcstAna.20070929_0900z.bin'];
% 
% [date_time1,              ...
% 	  obs_assim1,              ...
% 	  obs_species1,            ...
% 	  obs_tilenum1,            ...
% 	  obs_lon1,                ...
% 	  obs_lat1,                ...
% 	  obs_obs1,                ...
%  	  obs_obsvar1,             ...
% 	  obs_fcst1,               ...
%  	  obs_fcstvar1,            ...
% 	  obs_ana1,                ...
%  	  obs_anavar1              ...
%  	 ] =                      ...
%     read_ObsFcstAna(readFile1);
% 
% [date_time2,              ...
% 	  obs_assim2,              ...
% 	  obs_species2,            ...
% 	  obs_tilenum2,            ...
% 	  obs_lon2,                ...
% 	  obs_lat2,                ...
% 	  obs_obs2,                ...
%  	  obs_obsvar2,             ...
% 	  obs_fcst2,               ...
%  	  obs_fcstvar2,            ...
% 	  obs_ana2,                ...
%  	  obs_anavar2              ...
%  	 ] =                      ...
%     read_ObsFcstAna(readFile2);
% 
% [date_time3,              ...
% 	  obs_assim3,              ...
% 	  obs_species3,            ...
% 	  obs_tilenum3,            ...
% 	  obs_lon3,                ...
% 	  obs_lat3,                ...
% 	  obs_obs3,                ...
%  	  obs_obsvar3,             ...
% 	  obs_fcst3,               ...
%  	  obs_fcstvar3,            ...
% 	  obs_ana3,                ...
%  	  obs_anavar3              ...
%  	 ] =                      ...
%     read_ObsFcstAna(readFile3);


