clear;clc 

% Read in files 

path='/JRA-55/model_resolution/';
for i=1:2 
    disp(i) 
    if i==1
        var_abbr='054_pwat'; 
        var_name='PWAT_GDS4_EATM';
    end 
    
    if i==2
        var_abbr='227_cw'; 
        var_name='CW_GDS4_EATM';
    end 
    
    if i==3
        var_abbr='058_cice'; 
        var_name='CICE_GDS4_EATM';
    end 
    
    % Weird things with this data: skip first timestep (all missing) and use 3 hrs offset from initial time instead of 6 hrs
    data_in1=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2003010100_2003033121.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in2=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2003040100_2003063021.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in3=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2003070100_2003093021.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in4=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2003100100_2003123121.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in5=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2004010100_2004033121.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in6=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2004040100_2004063021.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in7=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2004070100_2004093021.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in8=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2004100100_2004123121.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in9=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2005010100_2005033121.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in10=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2005040100_2005063021.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in11=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2005070100_2005093021.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);
    data_in12=ncread([path 'fcst_column.' var_abbr '.reg_tl319.2005100100_2005123121.nc'],var_name,[1,1,1,2],[Inf,Inf,1,Inf]);

    data=cat(4,data_in1,data_in2,data_in3,data_in4,data_in5,data_in6,data_in7,data_in8,data_in9,data_in10,data_in11,data_in12);
    data=squeeze(data);
    size(data)
    
    data(:,:,1697:1700)=[]; % remove leap day 2004 
    data_dims=size(data);     
    skip=1:4:data_dims(3); dims=size(skip); nDays=dims(2);
    
    data_daily_mean=zeros(data_dims(1),data_dims(2),nDays);
    
    % Take daily means
    for j=1:nDays
        day_start=skip(j);
        data_daily_mean(:,:,j)=nanmean(data(:,:,day_start:(day_start+3)),3);
    end 
    
    if i==1
        pwat=data_daily_mean;
    end 
    
    if i==2
        cw=data_daily_mean;
    end 
    
    if i==3
        cice=data_daily_mean; 
    end 
    
end 

total_cloud=cw; % Can change this to add in cloud ice (cice) if needed

% Save 
final_dims=size(total_cloud);

% filename='fcst_column.pwat_cwpluscice.reg_tl319.2003-2005.nc';
% nccreate(filename,'pwat','Dimensions',{'i',final_dims(1),'j',final_dims(2),'time',final_dims(3)},'Format','netcdf4_classic');
% nccreate(filename,'cwpluscice','Dimensions',{'i',final_dims(1),'j',final_dims(2),'time',final_dims(3)},'Format','netcdf4_classic');
% ncwriteatt(filename,'/','creation_time',datestr(now));
% ncwriteatt(filename,'/','_FillValue',9.96921e36);
% ncwrite(filename,'pwat',pwat);
% ncwrite(filename,'cwpluscice',total_cloud);

filename='fcst_column.pwat_cw.reg_tl319.2003-2005.nc';
nccreate(filename,'pwat','Dimensions',{'i',final_dims(1),'j',final_dims(2),'time',final_dims(3)},'Format','netcdf4_classic');
nccreate(filename,'cw','Dimensions',{'i',final_dims(1),'j',final_dims(2),'time',final_dims(3)},'Format','netcdf4_classic');
ncwriteatt(filename,'/','creation_time',datestr(now));
ncwriteatt(filename,'/','_FillValue',9.96921e36);
ncwrite(filename,'pwat',pwat);
ncwrite(filename,'cw',cw);
