%------------------------------------------------------------------------------------------------
%   a program to vertically extrapolate/interpolate the 
%    repeated measurements at same time (different vertical location) to get
%   11 vertical levels of data using Chesapeake Bay Program measurements
%   i.e. get 11 values for one meausring time at one site, and intend to use the results 
%   (monthly sampled) to interpolate to the whole grid system in 3D 
%-------------------------------------------------------------------------------------------------

clear all;
startup;

indir='./2sepvar';
outdir='./4vertextrap';
if(~exist(outdir,'dir'))
  eval(['mkdir ' outdir]);
end

%read the grid file and find the depth of the station 
grd_file=['../model-results/ChesROMS/grid/chesroms_grd_smoother_msl.nc'];

ncgrid=netcdf(grd_file);
  lat_rho=ncgrid{'lat_rho'}(:);
  lon_rho=ncgrid{'lon_rho'}(:);
  h=ncgrid{'h'}(:);
  mask=ncgrid{'mask_rho'}(:);  %land/ocean mask, 1 for ocean, 0 for land;
  Inan=find(mask<0.5);   %find all the land locations
  h(Inan)=nan;   %set the depth to be nan for all land locations
close(ncgrid);

%--figdir----
figdir='./graph/';
if(~exist(figdir,'dir'))
  eval(['mkdir ' figdir]);
end

%try only doing DO for now%

parname=input('please input the parameter name (upper case): (''DO'', ''WTEMP'',''SALINITY'',''PO4'' ...) ');

iyear_plot=input('please input year to plot:(e.g. 1999, input 0 if no plots wanted)');

infilenames=dir([indir '/' char(parname) '_*.txt']); 

nf=length(infilenames);   %find number files to process

if ~isdir(outdir) 
     eval(['mkdir ' outdir]);   %make the directory
end

for i=1:nf
 
   a=load([indir '/' infilenames(i).name]);   %load values

   outfile=[outdir '/' infilenames(i).name];  %output averaged data

   display(['processing file number' num2str(i)  ':' infilenames(i).name]);


   [nr,nc]=size(a);    %number of measurements(records) in the file and number of columns

   n=1;

   k=1;

   sum=0;

   fl=fopen(outfile,'w');   %open output file 

   nstart=1;      %start index of repeated measurment
   nend=nstart;   %end index of repeated measurement
  
   while nend <= nr 

      mm=a(nstart,2);  %month
      dd=a(nstart,3);  %day
      yyyy=a(nstart,1); %year
      hh=a(nstart,4);   %hour
      mnt=a(nstart,5);  %minute

      k=nend+1;   

      while k <= nr 
            mm1=a(k,2);
            dd1=a(k,3);
            yyyy1=a(k,1);
            hh1=a(k,4);
            mnt1=a(k,5);
%
% measured at same time yet different depth of the same site
            if mm1==mm && dd1==dd && yyyy1==yyyy  && hh1 == hh && mnt1==mnt 

%
% measured at the same day, different time, different depth of the same site
%
%            if mm1==mm && dd1==dd && yyyy1==yyyy   

               nend=k;  %find this record is measured at the same time as before
               k=k+1;   %go to next line 

            else
               break;   %break the loop, go out of the loop immediately because
            end         %kend is already the last of the group of measurements
      end

%     nstart
%     nend

     average_nutr=mean(a(nstart:nend,7));   %calculate the averaged nutrient
                                              %of this group of meaurements

%   find depth of that station using grid info
     lat_sta=a(nstart,8);
     lon_sta=a(nstart,9);
     depth_sta=griddata(lat_rho,lon_rho,h,lat_sta, lon_sta,'nearest');
         
     depth_data=a(nstart:nend,6);   %the depths of the measurements (postive downward from surface)
     nutr_data=a(nstart:nend,7);    %the nutrient concentration of the measurements
 
                                    %  make sure thhe depth_data elements are distinct
     [depth_uniq Iuniq Juniq]=unique(depth_data);
     nutr_uniq=nutr_data(Iuniq);
     depth_data=[]; nutr_data=[];  %clear and set back
     depth_data=depth_uniq; nutr_data=nutr_uniq;
     clear('depth_uniq','nutr_uniq','Iuniq','Juniq');      %clear garbages 

     if(isnan(depth_sta))   %if found the location is actually land 
        depth_sta=max(depth_data);
     else
        depth_sta=max(max(depth_data),depth_sta); 
     end

     depth_min=0;
     depth_max=depth_sta;   %read this from the grid file      

     dz=(depth_max-depth_min)/10.0;
     depth_inter=(depth_max:-dz:depth_min); %depth levels to be interpolated/extrapolated
					    %depth_inter(1) for bottom, depth_inter(Nz) for surface

     Nz=length(depth_inter);               %number of vertical levels, should be 11

     if (length(depth_data)>=2) %more than 2 data points vertically, do the interpolation
       vert_nutr=interp1(depth_data,nutr_data,depth_inter,'nearest','extrap');  %do the interpolation/extrapolation
                                                                             %using the nearest vertical layers

%----------Make plots of them here------
       if (yyyy==iyear_plot) 
           hfig=figure('position',[1, 1, 800, 600],'visible','off');clf; 
           hold on;
           hdata=plot(nutr_data, -depth_data ,'ro','MarkerSize',10,'MarkerFaceColor','m'); 
           hvert=plot(vert_nutr, -depth_inter,'b^','MarkerSize',10,'MarkerFaceColor','g');
           legend([hdata,hvert],'measurments','measurments interpolated vertically');
           data_range=max( max(nutr_data),max(vert_nutr'));
           if data_range< 0.0001
              data_range = 0.0001 ;
           end
           xlim([-data_range*1.1,data_range*1.1]) ;
           ylim([-max(max(depth_data),max(depth_inter'))-1, 1]);  %from 1m below bottom to 1m above water)
           xlabel(parname);
           ylabel('depth below surface (m)');
           title( [parname ' of day ' datestr(datenum([yyyy mm dd hh mnt 0])) ]);
           saveas(gcf,[figdir  infilenames(i).name '_' num2str(yyyy,'%04d') ... 
                                       '_' num2str(mm  ,'%02d') ...
                                       '_' num2str(dd  ,'%02d') ... 
                                       '_' num2str(hh  ,'%02d') ...
                                       '_' num2str(mnt ,'%02d') '.eps'],'eps');
           saveas(gcf,[figdir  infilenames(i).name '_' num2str(yyyy,'%04d') ... 
                                       '_' num2str(mm  ,'%02d') ... 
                                       '_' num2str(dd  ,'%02d') ...
                                       '_' num2str(hh  ,'%02d') ... 
                                       '_' num2str(mnt ,'%02d') '.jpg'],'jpg');
           saveas(gcf,[figdir  infilenames(i).name '_' num2str(yyyy,'%04d') ... 
                                       '_' num2str(mm  ,'%02d') ... 
                                       '_' num2str(dd  ,'%02d') ...
                                       '_' num2str(hh  ,'%02d') ... 
                                       '_' num2str(mnt ,'%02d') 'fig'],'fig');
           delete(hfig);
        end
%----------

     else 
       vert_nutr=ones([Nz,1])*nutr_data(1,1);  %otherwise, just copy them to Nz levels
     end

     b(1:8)=0;

     b(1)=yyyy;
     b(2)=mm;
     b(3)=dd;
     b(4)=hh;
     b(5)=mnt;
     b(6)=average_nutr;
     b(7)=a(nstart,8);    %lat
     b(8)=a(nstart,9);    %lon
%
% print the exact time of the measurement and the depth averaged data
%     fprintf(fl,'%4d %2d %2d %2d %2d %2d  %6.2f %6.4f %6.4f \n', b(1),b(2),b(3),b(4),b(5),0,b(6),b(7),b(8));
%

%
% print the day of the measurement, the time is of measurement is assumed to be
% (hh+hh1)/2 : (mnt+mnt1)/2, i.e. the middle of the time span of the measurments of that day
% finally print the averaged value of that site on that day (differnt time, different depth)

 
     hmid= (hh+hh1)/2.0;  %hh:mnt is the first measurement time of that site on that day
     minmid= (mnt+mnt1)/2.0;%hh1:mnt1 is the last measurement time of that site on that day
     tmid=hmid*60+minmid;

     b(4)=floor(tmid/60);  %hour of the day
     b(5)=floor(rem(tmid,60));    %minutes of that hour

     for kz=1:Nz
       fprintf(fl,'%4d %2d %2d %2d %2d %2d  %6.2f %6.2f %6.4f %6.4f \n', b(1),b(2),b(3),b(4),b(5),0,depth_inter(kz),vert_nutr(kz),b(7),b(8));
     end

     nstart=nend+1;    %start line of the next group of measurements
     nend=nstart; 

  end
  fclose(fl);

end

