% edited by jdm 20 Mar 2014


clear all;
close all;

%--------------------------------------------------------------------------
% Common Variables
%--------------------------------------------------------------------------

% Date Time Data
startyear=2008;
startmonth=03;
startday=16;
starthour=09;
startminute=00;

endyear=2008;
endmonth=03;
endday=16;
endhour=18;
endminute=00;

% Lat and Lon bounds
ymax=42.1;
ymin=40.1;
xmax=-111.3;
xmin=-113.7;

% Radar Directory Selection
%--------------------------------------------------------------------------

% Radar data in netcdf format from Trevor and Kristens climos are saved in:
radardir = '/uufs/chpc.utah.edu/common/home/steenburgh-group1/KMTX/level3/events/netcdf/';

% Radar data for the 27 Oct 2010 case
% radardir='/uufs/chpc.utah.edu/common/home/steenburgh-group4/trevor/solpex/iop1/radar/level3/netcdf/';

% Radar data in netcdf format from 2010-2011 my pure banded non-banded
% cases: now includes 27 Oct 2010
%radardir = '/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/matlab/radarprecip/cases/';

% WRF Directory Selection
%--------------------------------------------------------------------------

wrfdir = '/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.5/';

% Case Type
% case_type = 'banded';
case_type = 'nonbanded';

wrfdir = [wrfdir case_type];

% Ensure month day and hour strings are two digits i.e. 02 not 2
  startmonthstr=num2str(startmonth);
  if startmonth<10
    startmonthstr=['0' num2str(startmonth)];
  end
  
  startdaystr=num2str(startday);
  if startday<10
    startdaystr=['0' num2str(startday)];
  end
  
  starthourstr=num2str(starthour);
  if starthour<10
    starthourstr=['0' num2str(starthour)];
  end
  
  startminutestr=num2str(startminute);
  if startminute<10
    startminutestr=['0' num2str(startminute)];
  end

  endmonthstr=num2str(endmonth);
  if endmonth<10
    endmonthstr=['0' num2str(endmonth)];
  end
  
  enddaystr=num2str(endday);
  if endday<10
    enddaystr=['0' num2str(endday)];
  end
  
  endhourstr=num2str(endhour);
  if endhour<10
    endhourstr=['0' num2str(endhour)];
  end
  
  endminutestr=num2str(endminute);
  if endminute<10
    endminutestr=['0' num2str(endminute)];
  end
  
% Create path and filename strings.
wrfpath1 = [num2str(startyear) '_' startmonthstr '_' startdaystr];

wrfoutdtg1 = [num2str(startyear) '-' startmonthstr '-' startdaystr '_' starthourstr ':' startminutestr ':00' ];

wrfpath2 = [num2str(endyear) '_' endmonthstr '_' enddaystr];

wrfoutdtg2 = [num2str(endyear) '-' endmonthstr '-' enddaystr '_' endhourstr ':' endminutestr ':00' ];

%--------------------------------------------------------------------------
% Compute NEXRAD SWE Estimate
%--------------------------------------------------------------------------

startnum=datenum(startyear,startmonth,startday,starthour,startminute,0);
endnum=datenum(endyear,endmonth,endday,endhour,endminute,0);

% Extract data and compute SWE total.

prevtime=0;
thistime=startnum;
while thistime<=endnum
  thisvec=datevec(thistime);
  thisyear=thisvec(1);
  thismonth=thisvec(2);
  thisday=thisvec(3);
  thishour=thisvec(4);
  thisminute=thisvec(5);
  monthstr=num2str(thismonth);
  
  if thismonth<10
    monthstr=['0' num2str(thismonth)];
  end
  
  daystr=num2str(thisday);
  if thisday<10
    daystr=['0' num2str(thisday)];
  end
  
  hourstr=num2str(thishour);
  if thishour<10
    hourstr=['0' num2str(thishour)];
  end
  
  minutestr=num2str(thisminute);
  if thisminute<10
    minutestr=['0' num2str(thisminute)];
  end
  
  radarfile=[radardir num2str(thisyear) '/' 'KMTX_N0R_' num2str(thisyear) monthstr daystr '_' hourstr minutestr '00.nc']; % All other cases. Now includes 27 Oct.
%   radarfile=[radardir '/' 'KMTX_N0R_' num2str(thisyear) monthstr daystr '_' hourstr minutestr '00.nc']; % 27 Oct 2010 case
  disp(['Searching for ' radarfile]);  
  disp(['Searching for ' radarfile]);
  if exist(radarfile)
    disp(['Found ' radarfile])
    if prevtime==0;
      firsttime=thistime;
      prevtime=thistime;
      dbz1 = (ncread(radarfile,'value'))'; %Trevors data
%      dbz1 = (ncread(radarfile,'bref'))';  %My data and 27 Oct 2010
      longitude = ncread(radarfile,'lon');
      latitude = ncread(radarfile,'lat');
      dbz1(isnan(dbz1))=0;
      dbz1(dbz1<0)=0;
      xind=find(longitude>=xmin & longitude<=xmax);
      yind=find(latitude>=ymin & latitude<=ymax);
      dbz1=dbz1(yind,xind);
      ptotal=dbz1*0;
    end
    if thistime>firsttime;
      timeelapsed=(thistime-prevtime)*24.0; %time elapsed in hours
      prevtime=thistime;
      dbz2 = (ncread(radarfile,'value'))'; %Trevors data
%       dbz2 = (ncread(radarfile,'bref'))'; %My data and 27 Oct 2010
      longitude = ncread(radarfile,'lon');
      latitude = ncread(radarfile,'lat');
      dbz2(isnan(dbz2))=0;
      dbz2(dbz2<0)=0;
      xind=find(longitude>=xmin & longitude<=xmax);
      yind=find(latitude>=ymin & latitude<=ymax);
      dbz2=dbz2(yind,xind);
      longitude=longitude(xind);
      latitude=latitude(yind);
      meandbz=0.5*(dbz1+dbz2);

      prate=((10.0.^(meandbz/10.0))/75.0).^0.5;
      prate(find(meandbz==0))=0;
      pstep=prate*timeelapsed;
      ptotal=ptotal+pstep;
      dbz1=dbz2;
    end
  end
  thistime=thistime+1/1440;
end

%--------------------------------------------------------------------------
% Compute WRF Accumulated Precip
%--------------------------------------------------------------------------

wrffile1=[wrfdir '/' wrfpath1 '/control/wrfout/wrfout_d03_' wrfoutdtg1];

wrffile2=[wrfdir '/' wrfpath2 '/control/wrfout/wrfout_d03_' wrfoutdtg2];

% extract variables from wrfout NetCDFs
xlat = double(ncread(wrffile1,'XLAT'));

xlong = double(ncread(wrffile1,'XLONG'));
 
ptota2 = double(ncread(wrffile2,'RAINNC'));
% ptota2 = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.5/banded/2007_12_27/control/wrfout/wrfout_d03_2007-12-28_00:00:00','RAINNC'));
% ptota2 = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.5/banded/2008_01_15/control/wrfout/wrfout_d03_2008-01-16_20:00:00','RAINNC'));
% ptota2 = double(ncread('/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/WRF_3.5/banded/2010_12_31/control/wrfout/wrfout_d03_2011-01-01_04:00:00','RAINNC'));
ptota1 = double(ncread(wrffile1,'RAINNC'));



% Difference total precipitation fields.
ptota = ptota2-ptota1;


%--------------------------------------------------------------------------
% Subset and interpolate WRF data to NEXRAD grid
%--------------------------------------------------------------------------

% Define the NEXRAD data grid
[longitude1,latitude1]  = meshgrid(longitude,latitude);

% Make the WRF lat, lon, and accumulation into vectors for triscatteredinterp
y1=reshape(xlat,numel(xlat),1);
x1=reshape(xlong,numel(xlong),1);
var1=reshape(ptota,numel(ptota),1);

% Interpolating function
F1 = TriScatteredInterp(x1,y1,var1);

% Map data from the interpolating function to NEXRAD grid defined above. 
intwrf = F1(longitude1,latitude1);

%--------------------------------------------------------------------------
% Create NetCDF files.  One for observations (NEXRAD) and one for forecast (WRF) 
%--------------------------------------------------------------------------

% Determine the number of grid points of the common grid
lat = size(latitude,1);
lon = size(longitude,1);
lat_ll = min(latitude);
lon_ll = min(longitude);

% Create NetCDF file name
%--------------------------------------------------------------------------
NEXRADname = ['NEXRAD' num2str(thisyear) monthstr daystr '.nc'];

% Create NEXRAD netcdf
%--------------------------------------------------------------------------
nccreate(NEXRADname,'lat','Format','classic','Dimensions',{'lat',lat,'lon',lon});
nccreate(NEXRADname,'lon','Format','classic','Dimensions',{'lat','lon'});
nccreate(NEXRADname,'APCP','Format','classic','Dimensions',{'lat','lon'});
ncwrite(NEXRADname,'lat',latitude1);
ncwrite(NEXRADname,'lon',longitude1);
ncwrite(NEXRADname,'APCP',ptotal);

% Write net cdf variable attributes to spoof MET
%--------------------------------------------------------------------------

% Define variables for integer attributes
%--------------------------------------------------------------------------

% UNIX Epoch Time same for obs.  I use date from the
% tutorial data to spoof the code.
X = int32(1123459200);

% Accumulation time.  I use date from the
% tutorial data to spoof the code.
Z = int32(43200);

% Delta lat/lon
dll = latitude(2) - latitude(1);

% Define variable attributes
%--------------------------------------------------------------------------
ncwriteatt(NEXRADname,'APCP','name','APCP');
ncwriteatt(NEXRADname,'APCP','long_name','Total precipitation');
ncwriteatt(NEXRADname,'APCP','level','A12');
ncwriteatt(NEXRADname,'APCP','units','mm');
ncwriteatt(NEXRADname,'APCP','FillValue',-9999.);
ncwriteatt(NEXRADname,'APCP','init_time','DTG matches tutorial to spoof MET');
ncwriteatt(NEXRADname,'APCP','init_time_ut',X);
ncwriteatt(NEXRADname,'APCP','valid_time','DTG matches tutorial to spoof MET');
ncwriteatt(NEXRADname,'APCP','valid_time_ut',X);
ncwriteatt(NEXRADname,'APCP','accum_time','Time matches tutorial to spoof MET');
ncwriteatt(NEXRADname,'APCP','accum_time_sec',Z);

% Define global attributes
%--------------------------------------------------------------------------
ncwriteatt(NEXRADname,'/','FileOrigins','pcp_combine format spoofed with MATLAB. precipnetcdf.m');
ncwriteatt(NEXRADname,'/','MET_version','V4.1');
ncwriteatt(NEXRADname,'/','MET_tool','NONE');
ncwriteatt(NEXRADname,'/','RunCommand','Accumulation periods vary by case');
ncwriteatt(NEXRADname,'/','Projection','LatLon');
ncwriteatt(NEXRADname,'/','lat_ll',lat_ll);
ncwriteatt(NEXRADname,'/','lon_ll',lon_ll);
ncwriteatt(NEXRADname,'/','delta_lat',dll);
ncwriteatt(NEXRADname,'/','delta_lon',dll);
ncwriteatt(NEXRADname,'/','Nlat',lat);
ncwriteatt(NEXRADname,'/','Nlon',lon);

% Create NetCDF file name
%--------------------------------------------------------------------------
WRFname = ['WRF' num2str(thisyear) monthstr daystr '.nc'];

% Create WRF netcdf.
%--------------------------------------------------------------------------
nccreate(WRFname,'lat','Format','classic','Dimensions',{'lat',lat,'lon',lon});
nccreate(WRFname,'lon','Format','classic','Dimensions',{'lat','lon'});
nccreate(WRFname,'APCP','Format','classic','Dimensions',{'lat','lon'});

ncwrite(WRFname,'lat',latitude1);
ncwrite(WRFname,'lon',longitude1);
ncwrite(WRFname,'APCP',intwrf);

% Write variable attributes
%--------------------------------------------------------------------------

% Set variables for integer attributes
%--------------------------------------------------------------------------

% UNIX Epoch Time different for forecast.  I use date from the
% tutorial data to spoof the code.
X = int32(1123372800);
Y = int32(1123459200);

% Define variable attributes
%--------------------------------------------------------------------------
ncwriteatt(WRFname,'APCP','name','APCP');
ncwriteatt(WRFname,'APCP','long_name','Total precipitation');
ncwriteatt(WRFname,'APCP','level','A12');
ncwriteatt(WRFname,'APCP','units','mm');
ncwriteatt(WRFname,'APCP','FillValue',-9999.);
ncwriteatt(WRFname,'APCP','init_time','DTG matches tutorial to spoof MET');
ncwriteatt(WRFname,'APCP','init_time_ut',X);
ncwriteatt(WRFname,'APCP','valid_time','DTG matches tutorial to spoof MET');
ncwriteatt(WRFname,'APCP','valid_time_ut',Y);
ncwriteatt(WRFname,'APCP','accum_time','Time matches tutorial to spoof MET');
ncwriteatt(WRFname,'APCP','accum_time_sec',Z);

% Define global attributes
%--------------------------------------------------------------------------
ncwriteatt(WRFname,'/','FileOrigins','pcp_combine format spoofed with MATLAB. precipnetcdf.m');
ncwriteatt(WRFname,'/','MET_version','V4.1');
ncwriteatt(WRFname,'/','MET_tool','NONE');
ncwriteatt(WRFname,'/','RunCommand','Accumulation periods vary by case');
ncwriteatt(WRFname,'/','Projection','LatLon');
ncwriteatt(WRFname,'/','lat_ll',lat_ll);
ncwriteatt(WRFname,'/','lon_ll',lon_ll);
ncwriteatt(WRFname,'/','delta_lat',dll);
ncwriteatt(WRFname,'/','delta_lon',dll);
ncwriteatt(WRFname,'/','Nlat',lat);
ncwriteatt(WRFname,'/','Nlon',lon);

%--------------------------------------------------------------------------
% Plot accumulations from NEXRAD and regridded WRF to verify the code worked. 
%--------------------------------------------------------------------------

% Set up figure and plot data
V = [1,5];
figsize=12;
figure(1)
set(gcf,'Renderer','Zbuffer'); %OpenGL causes a seg fault.  This avoids OpenGL
axesm('MapProjection','tranmerc','Grid','on','MapLatLimit',[ymin ymax],...
    'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on');
xfigsize=figsize;
yfigsize=figsize*(((ymax-ymin)/(xmax-xmin))/cos(((ymax+ymin)/2)/57.3));
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [xfigsize yfigsize]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 xfigsize yfigsize]);
contourm(latitude1,longitude1,intwrf,V, 'b','linewidth',2);
hold on;
contourm(latitude1,longitude1,ptotal,V,'k','linewidth',2);

% Read in and plot elevation data
elev = (ncread('/uufs/chpc.utah.edu/common/home/horel-group/archive/terrain/KMTX_terrain.nc','elevation'))';
longitude2 = double(ncread('/uufs/chpc.utah.edu/common/home/horel-group/archive/terrain/KMTX_terrain.nc','longitude'));
latitude2 = double(ncread('/uufs/chpc.utah.edu/common/home/horel-group/archive/terrain/KMTX_terrain.nc','latitude'));
xind=find(longitude2>=xmin & longitude2<=xmax);
yind=find(latitude2>=ymin & latitude2<=ymax);
elev=elev(yind,xind);
longitude2=longitude2(xind);
latitude2=latitude2(yind);
contourm(latitude2,longitude2,elev,1500:300:3500,'color',[0.5 0.25 0],'linewidth',1);

% Read in and plot Lakes and state boundaries.
load conus;
S2 = shaperead('hydrogp020','UseGeoCoords',true); %Lakes and other Waterbodies
levels = [S2.F_CODE];
L1 = S2(levels==1);
geoshow([L1.Lat],[L1.Lon],'linewidth',1.5,'Color','black');
geoshow(statelat,statelon,'DisplayType','line','Color','k','linewidth',1.5,'linestyle','-');

% Save Figure.

figname = [num2str(thisyear) monthstr daystr hourstr '.png'];

print('-dpng','-r200',figname);

