%Author: Joseph Zhang
%Generate ampl. and phases for elev, u,v from FES2014 using linear interp directly
%from nc files.
%Requires inputs: (1) list of open bnd lon/lat (generated by gen_fg.f90)
%                     longitude in [-180,180) or [0,360]
%                 (2) *.nc (in this dir)
%                 (3) iflag: 1: generate elev only; 2: elev, u,v
%  Also need to first download FES2014 elev and velocity data at AVISO site
%Outputs: ap.out.[1,2,...] (for bctides.in; see below for freq order)
clear all; close all;

iflag=2;
const={'s2','m2','n2','k2','k1','p1','o1','q1'};

%List of open bnd seg's that contain lon/lat of open bnd nodes
open_ll={'fg.bp.1','fg.bp.2','fg.bp.3'};

for iseg=1:length(open_ll)
%-----------------------------------------------------------
disp(['doing segment # ' num2str(iseg)]);
open=load(open_ll{iseg}); %ID,lon,lat of open bnd nodes
npt=size(open,1);

%miss_value=1.e10; %junk value
for ifl=1:2*iflag-1 %loop over elev, u,v
  for i=1:length(const)
    if(ifl==1)
      ncid = netcdf.open(['fes2014b_elevations_extrapolated/ocean_tide_extrapolated/' const{i} '.nc'],'NC_NOWRITE');
    elseif(ifl==2)
      ncid = netcdf.open(['eastward_velocity/' const{i} '.nc'],'NC_NOWRITE');
    elseif(ifl==3)
      ncid = netcdf.open(['northward_velocity/' const{i} '.nc'],'NC_NOWRITE');
    else
      error('Unknown iflag')
    end
    vid=netcdf.inqVarID(ncid,'lat'); %1D
    lat = netcdf.getVar(ncid, vid); %ascending order (-90,90)
    vid=netcdf.inqVarID(ncid,'lon');
    lon = netcdf.getVar(ncid, vid); %ascending order [0,360)
    nx=length(lon); ny=length(lat);
  
    if(ifl==1)
      vid=netcdf.inqVarID(ncid,'amplitude'); %(nx,ny)
    elseif(ifl==2)
      vid=netcdf.inqVarID(ncid,'Ua'); %(nx,ny)
    else
      vid=netcdf.inqVarID(ncid,'Va'); %(nx,ny)
    end
    amp0= netcdf.getVar(ncid, vid);
    amp0=amp0/100; %to meters or m/s
  
    if(ifl==1)
      vid=netcdf.inqVarID(ncid,'phase');
    elseif(ifl==2)
      vid=netcdf.inqVarID(ncid,'Ug');
    else
      vid=netcdf.inqVarID(ncid,'Vg');
    end
    pha0= netcdf.getVar(ncid, vid); %degr
  
    %Put phases into [0,360)
    indx=find(pha0<0);
    pha0(indx)=pha0(indx)+360;
    indx=find((pha0<0 | pha0>360) & pha0<1.4e10); %not junk but out of range
    if(sum(indx) ~=0) 
      pha0(indx)
      error('Phase out of bound');
    end;
  
  %  ap=zeros(npt,2);
    for j=1:npt
      lon2=open(j,2); if(lon2<0); lon2=lon2+360; end;
      lat2=open(j,3);
      I=find(lon>=lon2,1); %find 1st entry
      J=find(lat>=lat2,1);
      if(I<=1 || J<=1) 
        disp('Failed to find an interval'); [j lon2 lat2 I J]
        error('Bomb out'); 
      end
  
      ratx=(lon(I)-lon2)/(lon(I)-lon(I-1));
      raty=(lat(J)-lat2)/(lat(J)-lat(J-1));
      %Check junks
      amp_max=max([amp0(I-1,J-1);amp0(I,J-1);amp0(I-1,J);amp0(I,J)]);
      amp_min=min([amp0(I-1,J-1);amp0(I,J-1);amp0(I-1,J);amp0(I,J)]);
      if(amp_max>100)
        if(amp_min>100 || amp_min<0)
          [i j I J]
          error('All junks for amp:'); 
        else %use min
          amp(j,i,ifl)=amp_min;
        end
      else
        tmp1=ratx*amp0(I-1,J-1)+(1-ratx)*amp0(I,J-1);
        tmp2=ratx*amp0(I-1,J)+(1-ratx)*amp0(I,J);
        amp(j,i,ifl)=tmp1*raty+tmp2*(1-raty);
      end %amp_max
  
  %    if(ap(i,1)>2)
  %      disp('Amp>2m');
  %      [i j amp_max amp_min]
  %    end
  
      %Make sure there is no jump in phase range
      wild(1:4)=[pha0(I-1,J-1); pha0(I,J-1); pha0(I,J); pha0(I-1,J)]; %counter-clockwise
      pha_max=max(wild);
      pha_min=min(wild);
  
      if(pha_max>370) %junk
        if(pha_min<0 ||  pha_min>370)
          [i j I J]
          error('All junk phases:'); 
        end
        phase(j,i,ifl)=pha_min;
      else
        for k=1:4
          if(abs(wild(k)-pha_max)>180)
            wild(k)=wild(k)+360;
            if(abs(wild(k)-pha_max)>180)
              [i j k I J]
              wild(1:4)
              error('Phases still jump');
            end
          end
        end %for k
  
        tmp1=ratx*wild(1)+(1-ratx)*wild(2);
        tmp2=ratx*wild(4)+(1-ratx)*wild(3);
        phase(j,i,ifl)=tmp1*raty+tmp2*(1-raty);
      end %pha_max
  
    end %for j=1:npt
  end %for i - freq's
end %for ifl

%Output
fid=fopen(['ap.out.' num2str(iseg)],'w');
for i=1:length(const)
  fprintf(fid,'%s\n',const{i});
  fprintf(fid,'%f %f\n',[amp(:,i,1) phase(:,i,1)]');
end %for
  
if(iflag==2)
  for i=1:length(const)
    fprintf(fid,'%s\n',const{i});
    fprintf(fid,'%f %f %f %f\n',[amp(:,i,2) phase(:,i,2) amp(:,i,3) phase(:,i,3)]');
  end %for
end
fclose(fid);

clear open lon* lat* amp* pha*;
%-----------------------------------------------------------
end %for iseg
