%=========================================================================%
%% This routine takes soundings type data [lon] [lat] [data] which may be %
%% in the form of scatter data and then interpolates it to a ROMS grid    %
%% using an inverse-square weighting.  Thereafter, the interpolated fields%
%% are written into a *.mat file.                                         %
%%                                  Lyon Lanerolle NOAA/NOS/OCS/CSDL/MMAP %
%=========================================================================%

%=========================================================================%
%  Loop through each of the wet points in the grid file (after loading    %
%  in the lat, lon and mask values from the ROMS grid netcdf file) and    %
%  interpolate the data depths on to the points on the grid.  The inter-  %
%  polation method is one where the data points for the given radius of   %
%  influence are identified and then they are interpolated via a weighting%
%  method where the weights are the inverse of the square of the distance %
%  from the grid point of interest.                                       %
%=========================================================================%

m2deg=1/112000; % Conversion from meter to degrees 1 degree ~ 112000 m (112 km)

f=netcdf(grd_fl);

  lon_roms=f{'lon_rho'}(:);
  lat_roms=f{'lat_rho'}(:);
  msk_roms=f{'mask_rho'}(:);
  dx=1./f{'pm'}(:);dx=m2deg.*dx;
  dy=1./f{'pn'}(:);dy=m2deg.*dy;

close(f)

load MLW_correction.mat
MLW_corr=abs(q);
clear q

load MLLW_correction.mat
MLLW_corr=abs(q);
clear q

Iwet=find(msk_roms > 0.5);
Nwet=length(Iwet);

[Ny_roms,Nx_roms]=size(lon_roms);

h_roms=dpth_min.*ones(Ny_roms,Nx_roms);

disp('Now beginning the interpolation.......')
disp(' ')

icc=0;

for j=1:Ny_roms
  for i=1:Nx_roms
    
    if msk_roms(j,i) > 0.5

      icc=icc+1;

      lonr=lon_roms(j,i);
      latr=lat_roms(j,i);

      lon_min=lonr-0.5*dx(j,i);
      lon_max=lonr+0.5*dx(j,i);
      lat_min=latr-0.5*dy(j,i);
      lat_max=latr+0.5*dy(j,i);

      xv(1)=lon_min;xv(2)=lon_min;xv(3)=lon_max;xv(4)=lon_max;
      yv(1)=lat_max;yv(2)=lat_min;yv(3)=lat_min;yv(4)=lat_max;

      IN_mlw=inpolygon(lon_mlw,lat_mlw,xv,yv);Inn_mlw=find(IN_mlw > 0.5);
      IN_mllw=inpolygon(lon_mllw,lat_mllw,xv,yv);Inn_mllw=find(IN_mllw > 0.5);
      IN_enc=inpolygon(lon_enc,lat_enc,xv,yv);Inn_enc=find(IN_enc > 0.5);
      IN_paper=inpolygon(lon_paper,lat_paper,xv,yv);Inn_paper=find(IN_paper > 0.5);

      if length(Inn_mlw) > 0

        londat=lon_mlw(Inn_mlw);
        latdat=lat_mlw(Inn_mlw);
        dpthdat=dpth_mlw(Inn_mlw);
     
        dlon=londat-lonr;
        dlat=latdat-latr;

        ddis=sqrt(dlon.*dlon+dlat.*dlat);
        ddmin=mean(min(ddis));

        if ddmin <= 0.00001

          Idd=find(ddmin == ddis);
          if length(Idd) > 1
            h_roms(j,i)=mean(dpthdat(Idd));
          else
            h_roms(j,i)=dpthdat(Idd);
          end

        else

          ww=1./(ddis.*ddis);
          h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

        end

        h_roms(j,i)=h_roms(j,i)+MLW_corr(j,i);

      elseif length(Inn_mllw) > 0

        londat=lon_mllw(Inn_mllw);
        latdat=lat_mllw(Inn_mllw);
        dpthdat=dpth_mllw(Inn_mllw);
     
        dlon=londat-lonr;
        dlat=latdat-latr;

        ddis=sqrt(dlon.*dlon+dlat.*dlat);
        ddmin=mean(min(ddis));

        if ddmin <= 0.00001

          Idd=find(ddmin == ddis);
          if length(Idd) > 1
            h_roms(j,i)=mean(dpthdat(Idd));
          else
            h_roms(j,i)=dpthdat(Idd);
          end

        else

          ww=1./(ddis.*ddis);
          h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

        end

        h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

      elseif length(Inn_enc) > 0

        londat=lon_enc(Inn_enc);
        latdat=lat_enc(Inn_enc);
        dpthdat=dpth_enc(Inn_enc);
     
        dlon=londat-lonr;
        dlat=latdat-latr;

        ddis=sqrt(dlon.*dlon+dlat.*dlat);
        ddmin=mean(min(ddis));

        if ddmin <= 0.00001

          Idd=find(ddmin == ddis);
          if length(Idd) > 1
            h_roms(j,i)=mean(dpthdat(Idd));
          else
            h_roms(j,i)=dpthdat(Idd);
          end

        else

          ww=1./(ddis.*ddis);
          h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

        end

        h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

      elseif length(Inn_paper) > 0

        londat=lon_paper(Inn_paper);
        latdat=lat_paper(Inn_paper);
        dpthdat=dpth_paper(Inn_paper);
     
        dlon=londat-lonr;
        dlat=latdat-latr;

        ddis=sqrt(dlon.*dlon+dlat.*dlat);
        ddmin=mean(min(ddis));

        if ddmin <= 0.00001

          Idd=find(ddmin == ddis);
          if length(Idd) > 1
            h_roms(j,i)=mean(dpthdat(Idd));
          else
            h_roms(j,i)=dpthdat(Idd);
          end

        else

          ww=1./(ddis.*ddis);
          h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

        end

        h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

      else

        lon_min=lonr-dx(j,i);
        lon_max=lonr+dx(j,i);
        lat_min=latr-dy(j,i);
        lat_max=latr+dy(j,i);

        xv(1)=lon_min;xv(2)=lon_min;xv(3)=lon_max;xv(4)=lon_max;
        yv(1)=lat_max;yv(2)=lat_min;yv(3)=lat_min;yv(4)=lat_max;

        IN_mlw=inpolygon(lon_mlw,lat_mlw,xv,yv);Inn_mlw=find(IN_mlw > 0.5);
        IN_mllw=inpolygon(lon_mllw,lat_mllw,xv,yv);Inn_mllw=find(IN_mllw > 0.5);
        IN_enc=inpolygon(lon_enc,lat_enc,xv,yv);Inn_enc=find(IN_enc > 0.5);
        IN_paper=inpolygon(lon_paper,lat_paper,xv,yv);Inn_paper=find(IN_paper > 0.5);

        if length(Inn_mlw) > 0

          londat=lon_mlw(Inn_mlw);
          latdat=lat_mlw(Inn_mlw);
          dpthdat=dpth_mlw(Inn_mlw);
     
          dlon=londat-lonr;
          dlat=latdat-latr;

          ddis=sqrt(dlon.*dlon+dlat.*dlat);
          ddmin=mean(min(ddis));

          if ddmin <= 0.00001

            Idd=find(ddmin == ddis);
            if length(Idd) > 1
              h_roms(j,i)=mean(dpthdat(Idd));
            else
              h_roms(j,i)=dpthdat(Idd);
            end

          else

            ww=1./(ddis.*ddis);
            h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

          end

          h_roms(j,i)=h_roms(j,i)+MLW_corr(j,i);

        elseif length(Inn_mllw) > 0

          londat=lon_mllw(Inn_mllw);
          latdat=lat_mllw(Inn_mllw);
          dpthdat=dpth_mllw(Inn_mllw);
     
          dlon=londat-lonr;
          dlat=latdat-latr;

          ddis=sqrt(dlon.*dlon+dlat.*dlat);
          ddmin=mean(min(ddis));

          if ddmin <= 0.00001

            Idd=find(ddmin == ddis);
            if length(Idd) > 1
              h_roms(j,i)=mean(dpthdat(Idd));
            else
              h_roms(j,i)=dpthdat(Idd);
            end

          else

            ww=1./(ddis.*ddis);
            h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

          end

          h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

        elseif length(Inn_enc) > 0

          londat=lon_enc(Inn_enc);
          latdat=lat_enc(Inn_enc);
          dpthdat=dpth_enc(Inn_enc);
     
          dlon=londat-lonr;
          dlat=latdat-latr;

          ddis=sqrt(dlon.*dlon+dlat.*dlat);
          ddmin=mean(min(ddis));

          if ddmin <= 0.00001

            Idd=find(ddmin == ddis);
            if length(Idd) > 1
              h_roms(j,i)=mean(dpthdat(Idd));
            else
              h_roms(j,i)=dpthdat(Idd);
            end

          else

            ww=1./(ddis.*ddis);
            h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

          end

          h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

        elseif length(Inn_paper) > 0

          londat=lon_paper(Inn_paper);
          latdat=lat_paper(Inn_paper);
          dpthdat=dpth_paper(Inn_paper);
     
          dlon=londat-lonr;
          dlat=latdat-latr;

          ddis=sqrt(dlon.*dlon+dlat.*dlat);
          ddmin=mean(min(ddis));

          if ddmin <= 0.00001

            Idd=find(ddmin == ddis);
            if length(Idd) > 1
              h_roms(j,i)=mean(dpthdat(Idd));
            else
              h_roms(j,i)=dpthdat(Idd);
            end

          else

            ww=1./(ddis.*ddis);
            h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

          end

          h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

        else

          lon_min=lonr-1.5*dx(j,i);
          lon_max=lonr+1.5*dx(j,i);
          lat_min=latr-1.5*dy(j,i);
          lat_max=latr+1.5*dy(j,i);

          xv(1)=lon_min;xv(2)=lon_min;xv(3)=lon_max;xv(4)=lon_max;
          yv(1)=lat_max;yv(2)=lat_min;yv(3)=lat_min;yv(4)=lat_max;

          IN_mlw=inpolygon(lon_mlw,lat_mlw,xv,yv);Inn_mlw=find(IN_mlw > 0.5);
          IN_mllw=inpolygon(lon_mllw,lat_mllw,xv,yv);Inn_mllw=find(IN_mllw > 0.5);
          IN_enc=inpolygon(lon_enc,lat_enc,xv,yv);Inn_enc=find(IN_enc > 0.5);
          IN_paper=inpolygon(lon_paper,lat_paper,xv,yv);Inn_paper=find(IN_paper > 0.5);

          if length(Inn_mlw) > 0

            londat=lon_mlw(Inn_mlw);
            latdat=lat_mlw(Inn_mlw);
            dpthdat=dpth_mlw(Inn_mlw);
     
            dlon=londat-lonr;
            dlat=latdat-latr;

            ddis=sqrt(dlon.*dlon+dlat.*dlat);
            ddmin=mean(min(ddis));

            if ddmin <= 0.00001

              Idd=find(ddmin == ddis);
              if length(Idd) > 1
                h_roms(j,i)=mean(dpthdat(Idd));
              else
                h_roms(j,i)=dpthdat(Idd);
              end

            else

              ww=1./(ddis.*ddis);
              h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

            end

            h_roms(j,i)=h_roms(j,i)+MLW_corr(j,i);

          elseif length(Inn_mllw) > 0

            londat=lon_mllw(Inn_mllw);
            latdat=lat_mllw(Inn_mllw);
            dpthdat=dpth_mllw(Inn_mllw);
     
            dlon=londat-lonr;
            dlat=latdat-latr;

            ddis=sqrt(dlon.*dlon+dlat.*dlat);
            ddmin=mean(min(ddis));

            if ddmin <= 0.00001

              Idd=find(ddmin == ddis);
              if length(Idd) > 1
                h_roms(j,i)=mean(dpthdat(Idd));
              else
                h_roms(j,i)=dpthdat(Idd);
              end

            else

              ww=1./(ddis.*ddis);
              h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

            end

            h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

          elseif length(Inn_enc) > 0

            londat=lon_enc(Inn_enc);
            latdat=lat_enc(Inn_enc);
            dpthdat=dpth_enc(Inn_enc);
     
            dlon=londat-lonr;
            dlat=latdat-latr;

            ddis=sqrt(dlon.*dlon+dlat.*dlat);
            ddmin=mean(min(ddis));

            if ddmin <= 0.00001

              Idd=find(ddmin == ddis);
              if length(Idd) > 1
                h_roms(j,i)=mean(dpthdat(Idd));
              else
                h_roms(j,i)=dpthdat(Idd);
              end

            else

              ww=1./(ddis.*ddis);
              h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

            end

            h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

          elseif length(Inn_paper) > 0

            londat=lon_paper(Inn_paper);
            latdat=lat_paper(Inn_paper);
            dpthdat=dpth_paper(Inn_paper);
     
            dlon=londat-lonr;
            dlat=latdat-latr;

            ddis=sqrt(dlon.*dlon+dlat.*dlat);
            ddmin=mean(min(ddis));

            if ddmin <= 0.00001

              Idd=find(ddmin == ddis);
              if length(Idd) > 1
                h_roms(j,i)=mean(dpthdat(Idd));
              else
                h_roms(j,i)=dpthdat(Idd);
              end

            else

              ww=1./(ddis.*ddis);
              h_roms(j,i)=sum(ww.*dpthdat)./sum(ww);

            end

            h_roms(j,i)=h_roms(j,i)+MLLW_corr(j,i);

          else

            h_roms(j,i)=nan;

          end

        end

      end

      rat=icc/Nwet;
      if mod(icc-1,200) == 0
        disp(['Completed ',num2str(100.*rat),'% of the interpolation....'])
      end

    end

  end
end

bathy=h_roms;

save bathy_interp_output.mat bathy

disp('Depth data interpolation complete, now exiting.............')

