clc
clear all
close all

read_path = '/discover/nobackup/projects/geoscm/fzeng/ALEXI/';
save_path = '/discover/nobackup/projects/geoscm/fzeng/ALEXI/monthly/';
start_year = 2015;
stop_year = 2020;

leapYears = 1980:4:2022;
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];

mat_index = 0;

im_alexi_m36=214
jm_alexi_m36=98

latcenter1=[-59.5718,-59.0241,-58.4849,-57.9538,-57.4304,-56.9143,-56.4052,-55.9027,-55.4067,-54.9167,-54.4326,-53.954014,-53.4809,-53.0129,-52.5499,-52.0917,-51.6381,-51.1889,-50.7441,-50.3033,-49.8666,-49.4338,-49.0047,-48.5792,-48.1572,-47.7386,-47.3233,-46.9113,-46.5023,-46.0963,-45.6933,-45.2931,-44.8957,-44.5010,-44.1089,-43.7193,-43.3323,-42.9476,-42.5653,-42.1853,-41.8075,-41.4319,-41.0585,-40.6871,-40.3177,-39.9504,-39.5849,-39.2214,-38.8596,-38.4997,-38.1416,-37.7851,-37.4304,-37.0773,-36.7258,-36.3759,-36.0275,-35.6806,-35.3352,-34.9912,-34.6487,-34.3075,-33.9677,-33.6292,-33.2921,-32.9562,-32.6215,-32.2881,-31.9558,-31.6248,-31.2949,-30.9661,-30.6384,-30.3118,-29.9863,-29.6618,-29.3384,-29.0159,-28.6944,-28.3739,-28.0543,-27.7357,-27.4180,-27.1011,-26.7852,-26.4701,-26.1558,-25.8424,-25.5297,-25.2179,-24.9068,-24.5966,-24.2870,-23.9782,-23.6701,-23.3628,-23.0561,-22.7501,-22.4448,-22.1401,-21.8361,-21.5327,-21.2299,-20.9277,-20.6261,-20.3251,-20.0247,-19.7248,-19.4255,-19.1267,-18.8285,-18.5308,-18.2335,-17.9368,-17.6405,-17.3447,-17.0494,-16.7546,-16.4602,-16.1662,-15.8726,-15.5795,-15.2867,-14.9944,-14.7025,-14.4109,-14.1197,-13.8289,-13.5384,-13.2483,-12.9585,-12.6690,-12.3799,-12.0910,-11.8025,-11.5143,-11.2263,-10.9386,-10.6513,-10.3641,-10.0772,-9.79061,-9.50423,-9.21807,-8.93214,-8.64643,-8.36093,-8.07564,-7.79054,-7.50563,-7.22090,-6.93634,-6.65195,-6.36773,-6.08365,-5.79972,-5.51593,-5.23228,-4.94874,-4.66533,-4.38203,-4.09883,-3.81573,-3.53272,-3.24979,-2.96694,-2.68417,-2.40145,-2.11879,-1.83619,-1.55362,-1.27110,-0.988599,-0.706126,-0.423669,-0.141222];
  latcenter2=[0.141222,0.423669,0.706126,0.988599,1.27110,1.55362,1.83619,2.11879,2.40145,2.68417,2.96694,3.24979,3.53272,3.81573,4.09883,4.38203,4.66533,4.94874,5.23228,5.51593,5.79972,6.08365,6.36773,6.65195,6.93634,7.22090,7.50563,7.79054,8.07564,8.36093,8.64643,8.93214,9.21807,9.50423,9.79061,10.0772,10.3641,10.6513,10.9386,11.2263,11.5143,11.8025,12.0910,12.3799,12.6690,12.9585,13.2483,13.5384,13.8289,14.1197,14.4109,14.7025,14.9944,15.2867,15.5795,15.8726,16.1662,16.4602,16.7546,17.0494,17.3447,17.6405,17.9368,18.2335,18.5308,18.8285,19.1267,19.4255,19.7248,20.0247,20.3251,20.6261,20.9277,21.2299,21.5327,21.8361,22.1401,22.4448,22.7501,23.0561,23.3628,23.6701,23.9782,24.2870,24.5966,24.9068,25.2179,25.5297,25.8424,26.1558,26.4701,26.7852,27.1011,27.4180,27.7357,28.0543,28.3739,28.6944,29.0159,29.3384,29.6618,29.9863,30.3118,30.6384,30.9661,31.2949,31.6248,31.9558,32.2881,32.6215,32.9562,33.2921,33.6292,33.9677,34.3075,34.6487,34.9912,35.3352,35.6806,36.0275,36.3759,36.7258,37.0773,37.4304,37.7851,38.1416,38.4997,38.8596,39.2214,39.5849,39.9504,40.3177,40.6871,41.0585,41.4319,41.8075,42.1853,42.5653,42.9476,43.3323,43.7193,44.1089,44.5010,44.8957,45.2931,45.6933,46.0963,46.5023,46.9113,47.3233,47.7386,48.1572,48.5792,49.0047,49.4338,49.8666,50.3033,50.7441,51.1889,51.6381,52.0917,52.5499,53.0129,53.4809,53.9540,54.4326,54.9167,55.4067,55.9027,56.4052,56.9143,57.4304,57.9538,58.4849,59.0241,59.5718,60.1286,60.6948,61.2711,61.8582,62.4566,63.0672,63.6908,64.3284,64.9810,65.6498,66.3363,67.0421,67.7688,68.5188,69.2945,70.0989,70.9357,71.8094,72.7255,73.6913,74.7160,75.8126,76.9991,78.3030,79.7690,81.4803];

latcenter = [latcenter1 latcenter2];
loncenter = -160.394:0.373444:(-160.394+900*0.373444);


for y = start_year:stop_year

    y_str = num2str(y,'%04i');

    % define month length based on leap year status
    if sum(y==leapYears)>0
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

        m_str = num2str(m,'%02i');

        alexi_mat = NaN*ones(length(latcenter),length(loncenter),31);

        for d = 1:monthLength(m)
            dt = datetime([y m d]);
            doy = day(dt,'dayofyear');
            doy_str = num2str(doy,'%03i');

            read_file = [read_path y_str '_M36_conus/EDAY_ALEXI_V10E_' y_str doy_str '_M36_conus.dat'];
           
            fid = fopen(read_file,'r');
            data = fread(fid,[jm_alexi_m36 im_alexi_m36],'real*4');
            fclose(fid);
            data(data==-9999) = NaN;

            alexi_mat(245:(245+jm_alexi_m36-1),55:(55+im_alexi_m36-1),d) = data;

            datenum_vec(mat_index) = datenum(y,m,1);
            year_vec(mat_index) = y;
            month_vec(mat_index) = m;

        end 

        alexi_monthly_avg = nanmean(alexi_mat,3);

        save_file = [save_path 'ALEXI_ET_' y_str m_str '.nc'];

        nccreate(save_file,"ET","Dimensions",{"lat",size(alexi_monthly_avg,1),"lon",size(alexi_monthly_avg,2)}, ...
                   "FillValue","disable");
        nccreate(save_file,"latitude","Dimensions",{"lat",size(alexi_monthly_avg,1)}, ...
                   "FillValue","disable");
        nccreate(save_file,"longitude","Dimensions",{"lon",size(alexi_monthly_avg,2)}, ...
                   "FillValue","disable");


        ncwrite(save_file,'ET',alexi_monthly_avg);
        ncwriteatt(save_file,'ET',"long_name","monthly average evapotranspiration");
        ncwriteatt(save_file,'ET',"units","mm day^-1");

        ncwrite(save_file,'latitude',latcenter);
        ncwriteatt(save_file,'latitude',"long_name","latitude");
        ncwriteatt(save_file,'latitude',"units","degrees North");

        ncwrite(save_file,'longitude',loncenter);
        ncwriteatt(save_file,'longitude',"long_name","longitude");
        ncwriteatt(save_file,'longitude',"units","degrees East");
 
    end
end

evap = zeros(length(latcenter),length(loncenter));
evap(245:(245+jm_alexi_m36-1),55:(55+im_alexi_m36-1)) = data;

[lon_mat,lat_mat] = meshgrid(loncenter,latcenter);

%load coastlines;
%[latcells, loncells] = polysplit(coastlat, coastlon);
%
%latlim = [-60 90];
%lonlim = [-180 180];
%
%cmap_mat = cmocean('algae',12);
%%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
%cmap = cmap_mat;
%s1 = quantile(evap(:),0.01);
%s2 = quantile(evap(:),0.99);
%caxis_min = s1;
%caxis_max = s2;
%
%figure
%ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                                             
%    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
%framem off;
%gridm on;
%setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                                              
%    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                                                
%    'MLabelParallel','south','FontSize',12)
%geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
%pcolorm(lat_mat,lon_mat,evap(:,:))
%caxis([caxis_min caxis_max]);
%set(gcf,'Colormap',cmap);
%plotm(coastlat, coastlon,'k')
%hold on
%%geoshow(states,'facealpha',0);                                                                                     
%set(gca,'XLim',lonlim);
%set(gca,'YLim',latlim);
%tightmap
%title(['ALEXI ET [mm/day]'])
%h = colorbar;
%set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
%set(gcf,'paperunits','centimeters')
%set(gcf,'paperposition',[1 1 17 15])
%saveas(gcf,'test.png')



% convert data to monthly means


