% ----------------------------- % 
% Download CNAPS
% Note: Create variable names and then "eval" to execute bash commands 
% Written by: Taylor Shropshire 
% Date 10/07/22
% ------------------------------ % 
clear; clc; close all;

down_flag=0;
if down_flag==1
% ############################################### % 
% ################## Download CNAPS ############# % 
% ############################################### % 

% Download info 
out_path='/gpfs_backup/he_data/he/tashrops/projects/ROFFS/CNAPS/';
url='http://oceanus.meas.ncsu.edu:8080/thredds/dodsC/fmrc/useast2_coawst_roms/';
file_name='COAWST-ROMS_SWAN_Forecast_Model_Run_Collection_best.ncd';
dir_out='/gpfs_backup/he_data/he/tashrops/projects/ROFFS/CNAPS_opendap/'
file_form='.ascii';

xind=[0:481];   xdim_ori=length(xind);
yind=[0:401];   ydim_ori=length(yind);
zind=[0:35];    zdim_ori=length(zind);
tind=[0:1791];  tdim_ori=length(tind);

get_grid=0;
if get_grid==1
% ------------------------------ % 
% ----- Download Grid Info------ % 
% ------------------------------ % 
% Get longitude
xrng=['%5B',num2str(xind(1)),':1:',num2str(xind(end)),'%5D'];
yrng=['%5B',num2str(yind(1)),':1:',num2str(yind(end)),'%5D'];
var_rng=[xrng,yrng];
file_down=[url,file_name,file_form,'?','lon_rho',var_rng];
file_out=[dir_out,'grid/lon_rho.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Get latitude 
xrng=['%5B',num2str(xind(1)),':1:',num2str(xind(end)),'%5D'];
yrng=['%5B',num2str(yind(1)),':1:',num2str(yind(end)),'%5D'];
var_rng=[xrng,yrng];
file_down=[url,file_name,file_form,'?','lat_rho',var_rng];
file_out=[dir_out,'grid/lat_rho.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Get time 
trng=['%5B',num2str(tind(1)),':1:',num2str(tind(end)),'%5D'];
file_down=[url,file_name,file_form,'?','time_run',trng];
file_out=[dir_out,'grid/time_run.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Get bathy 
xrng=['%5B',num2str(xind(1)),':1:',num2str(xind(end)),'%5D'];
yrng=['%5B',num2str(yind(1)),':1:',num2str(yind(end)),'%5D'];
var_rng=[xrng,yrng];
file_down=[url,file_name,file_form,'?','h',var_rng];
file_out=[dir_out,'grid/bathy_rho.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Get mask
xrng=['%5B',num2str(xind(1)),':1:',num2str(xind(end)),'%5D'];
yrng=['%5B',num2str(yind(1)),':1:',num2str(yind(end)),'%5D'];
var_rng=[xrng,yrng];
file_down=[url,file_name,file_form,'?','mask_rho',var_rng];
file_out=[dir_out,'grid/mask_rho.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Get depth_t 
trng=['%5B',num2str(zind(1)),':1:',num2str(zind(end)),'%5D'];
file_down=[url,file_name,file_form,'?','s_rho',trng];
file_out=[dir_out,'grid/s_rho.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Get depth_w
trng=['%5B',num2str(zind(1)),':1:',num2str(zind(end)),'%5D'];
file_down=[url,file_name,file_form,'?','s_w',trng];
file_out=[dir_out,'grid/s_w.bin'];
eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);

% Example 
% http://oceanus.meas.ncsu.edu:8080/thredds/dodsC/fmrc/useast2_coawst_roms/COAWST-ROMS_SWAN_Forecast_Model_Run_Collection_best.ncd.ascii?temp%5B0:1:1%5D%5B0:1:35%5D%5B0:1:10%5D%5B0:1:10%5D
% ---------------------------------- % 
% ---------------------------------- % 
% ---------------------------------- % 
end % if get_grid==1

% Read files
file_out=[dir_out,'grid/lon_rho.bin'];
data_all=importdata(file_out);
xc=data_all.data;
file_out=[dir_out,'grid/lat_rho.bin'];
data_all=importdata(file_out);
yc=data_all.data;
file_out=[dir_out,'grid/bathy_rho.bin'];
data_all=importdata(file_out);
bathy=data_all.data;
file_out=[dir_out,'grid/mask_rho.bin'];
data_all=importdata(file_out);
mask=data_all.data;

% Set NaN for mask
mask(find(mask==0))=NaN;
mask_ori=mask;

% Simulation Date
file_out=[dir_out,'grid/time.bin'];
data_all=importdata(file_out);
model_date=data_all.data;
date_ref=datenum(2022,01,24);
model_date=date_ref+(model_date/24);

% Depth 
file_out=[dir_out,'grid/s_rho.bin'];
data_all=importdata(file_out);
dfrac_t=data_all.data;
file_out=[dir_out,'grid/s_w.bin'];
data_all=importdata(file_out);
dfrac_w=data_all.data;

% Change sign
dfrac_t=dfrac_t*-1;
dfrac_w=dfrac_w*-1;

% Extract
[indy,indx]=find(xc>=-82 & xc<=-71 & yc>=28 & yc<=39);
x0=min(indx); 		x1=max(indx);		dx=1;
y0=min(indy);		y1=max(indy);		dy=1;
z0=1;			z1=zdim_ori;		dz=1;
xc=xc(y0:dy:y1,x0:dx:x1);
yc=yc(y0:dy:y1,x0:dx:x1);
mask=mask(y0:dy:y1,x0:dx:x1);
bathy=bathy(y0:dy:y1,x0:dx:x1);
bathy2=bathy.*mask;

% Set dims
xdim=size(xc,2);
ydim=size(yc,1);
tdim=length(model_date);
zdim=zdim_ori;

% Plot for reference
figure
pcolor(xc,yc,mask);shading flat; colorbar; hold on;
plot(cc(:,1),cc(:,2),'k');
plot(lonR,latR,'.r','markersize',20);
axis equal; axis([min(lonR)-4 max(lonR)+4 min(latR)-4 max(latR)+4]);
xlabel('Longtidue'); ylabel('Latitude'); set(gca,'FontSize',14);

% ------------------------------------------------ % 
% --------------- Down 3D Files ------------------ % 
% ------------------------------------------------ % 
var_name={'temp','salt'};
var_name_MOD={'TMP','SAL'};
var_name_long={'Temperature','Salinity'};
var_units={'deg C','PSU (g/kg)'};

% Switched for some reason.
xx0=y0;
yy0=x0;
xx1=y1;
yy1=x1;

% Set download domain limits
xrng=['%5B',num2str(xx0-1),':',num2str(dx),':',num2str(xx1-1),'%5D'];
yrng=['%5B',num2str(yy0-1),':',num2str(dy),':',num2str(yy1-1),'%5D'];
zrng=['%5B',num2str(z0-1),':',num2str(dz),':',num2str(z1-1),'%5D'];

exist_flag=[];
% =============================== % 
% ========= Main Loop =========== % 
% =============================== % 
%field_all=[];
for tt=1:length(uni_date)

display(' *********************** ');
display(' *********************** ');
display(' *********************** ');
display(['Reading ' datestr(uni_date(tt))]); tic; %pause;
display(['Day File = ', num2str(tt), '/',num2str(length(uni_date))]);

field_all=[];
for vv=1:length(var_name)
display(' ');
display(['Var = ', cell2mat(var_name(vv))]);

% ---- Get File Name ------ % 
ind_tt=find(model_date==uni_date(tt));
t0=ind_tt;	
t1=ind_tt+1;
trng=['%5B',num2str(t0-1),':1:',num2str(t0-1),'%5D'];
var_rng=[trng,zrng,xrng,yrng];

var_str=cell2mat(var_name(vv));
file_down=[url,file_name,file_form,'?',var_str,var_rng];

yy=num2str(uni_year(tt),'%02.f');
mm=num2str(uni_mon(tt),'%02.f');
dd=num2str(uni_day(tt),'%02.f');
fn=[var_str,'_',yy,mm,dd,'.bin'];
file_out=[dir_out,'bin/',var_str,'/',fn];
% ------------------------- % 

% ***** Download ******* % 
display('Downloading');
%eval(['! wget -O ' file_out ' ' ['"',file_down,'"']]);
eval(['! wget -O ' file_out ' ' ['"',file_down,'"'] ' -q']);
%eval(['! wget -O ' file_out ' ' ['"',file_down,'"'] ' -q --show-progress']);
%eval(['! wget -q --show-progress -O ' file_out ' ' ['"',file_down,'"']]);
%eval(['! wget --progress=bar -O ' file_out ' ' ['"',file_down,'"']]);
% *********************** %

% ----- Read and modify array to 3D ---- % 
data_all=importdata(file_out);
file_flag=0;

% !!!!!!!!! Check if CNAPS field exist !!!!!!!!! % 
if ~isempty(data_all)
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! % 

file_flag=1;
field=data_all.data;
field(find(field==1e37))=NaN;

del_flag=0;
if del_flag==1
% Option to remove file
eval(['!rm -f ' file_out]);
end 

field2=[];
ind1=1; ind2=ydim;
for zz=size(field,1)/ydim:-1:1;
var=field(ind1:ind2,:);        field2(:,:,zz)=var;
ind1=ind2+1;                    ind2=ind1+(ydim-1);
end
field=field2; clear field2; 
% ------------------------------------- % 

% *************************************** % 
% ***** Interpolate to z-level ********** % 
% *************************************** % 
% note: optional... could just use first sigma layer. Better to be 
% consistant with z-level

% -----Get depth interface ------- % 
var_rng=[trng,xrng,yrng];
var_str='zeta';
file_down=[url,file_name,file_form,'?',var_str,var_rng];
fn_zeta=[var_str,'_',yy,mm,dd,'.bin'];
file_zeta=[dir_out,'bin/',var_str,'/',fn_zeta];
display('Downloading SSH');
eval(['! wget -q -O ' file_out ' ' ['"',file_down,'"']]);
data_all=importdata(file_out);
zeta=data_all.data;
zeta(find(zeta==1e37))=NaN;

tot_H=bathy+zeta;
depth_t=[];	depth_w=[]; 
for zz=1:zdim
depth_t(:,:,zdim-(zz-1))=tot_H*dfrac_t(zz);
depth_w(:,:,zdim-(zz-1))=tot_H*dfrac_w(zz);
end
% ----------------------------------- % 

% ------- Interpolate ------------ % 
display('Interpolate to z-level');
dz1=10;			zlim1=100;
dz2=50;			zlim2=300;
dz3=200;		zlim3=1500;
dz4=500;		zlim4=5000;
int1=[0:dz1:zlim1];
int2=[zlim1+dz2:dz2:zlim2];
int3=[zlim2+dz3:dz3:zlim3];
int4=[zlim3+dz4:dz4:zlim4-dz4];
depth_zl=[int1,int2,int3,int4,zlim4-1];

zdim_int=length(depth_zl);

field_int=zeros(ydim,xdim,length(depth_zl))*NaN;
for yy=1:ydim
for xx=1:xdim

if mask(yy,xx)==1
prof_v=squeeze(field(yy,xx,:));
prof_d=squeeze(depth_t(yy,xx,:));
pred=interp1(prof_d,prof_v,depth_zl,'linear','extrap');
ind=find(depth_zl>=bathy(yy,xx) | pred<0);
pred(ind)=NaN;
field_int(yy,xx,:)=pred;
end 

end 
end 
% ------------------------------------- % 

% Record for a single netCDF
field2=[];
for zz=1:zdim_int
field2(:,:,zz)=field_int(:,:,zz)';
end
field_all(:,:,:,vv)=field2;
%field_all(:,:,:,vv,tt)=field2;

else 	%if ~isempty(data_all)
display(' ');
display('!!!!!!!!!!!!!!!!!!!!!!!!!!');
display('!!!!!!!!!!!!!!!!!!!!!!!!!!');
display('!!!!!!!!!!!!!!!!!!!!!!!!!!');
display('CNAPS Field Does Not Exist'); 
display('!!!!!!!!!!!!!!!!!!!!!!!!!!');
display('!!!!!!!!!!!!!!!!!!!!!!!!!!');
display('!!!!!!!!!!!!!!!!!!!!!!!!!!');
display(' ');
end % if ~isempty(data_all)

exist_flag(1,tt)=file_flag;
% *************************************** % 
% *************************************** % 
% *************************************** % 
end % for vv=1:length(var_name)

%display('Stop Before Writting NetCDF - Pause'); pause;

if del_flag==0 & file_flag==1
% ******************************************** % 
% ******* Write netCDF File ****************** % 
% ******************************************** %
display(' '); display('Writting NetCDF');  display(' ');
% note: optional... could just interpolate to the ROFFS points
% don't run if deleting ascii file after reading.

% ----- Grid File ----- % 
if tt==1
nc_path='/gpfs_backup/he_data/he/tashrops/projects/ROFFS/CNAPS_opendap/netCDF/';
yy=num2str(uni_year(tt),'%02.f');
mm=num2str(uni_mon(tt),'%02.f');
dd=num2str(uni_day(tt),'%02.f');
nc_name=['grid/CNAPS_grid.nc'];
nc_file=[nc_path,nc_name];
eval(['! rm -f ' nc_file]);

nccreate(nc_file,'XC','Dimensions',{'X', xdim, 'Y', ydim});
ncwriteatt(nc_file, 'XC', 'standard_name', 'XC');
ncwriteatt(nc_file, 'XC', 'long_name', 'Longitude');
ncwriteatt(nc_file, 'XC', 'units', 'degrees west');
ncwriteatt(nc_file, 'XC', 'CoordinateAxisType', 'XC');
ncwrite(nc_file,'XC',xc');
nccreate(nc_file,'YC','Dimensions',{'X', xdim, 'Y', ydim});
ncwriteatt(nc_file, 'YC', 'standard_name', 'YC');
ncwriteatt(nc_file, 'YC', 'long_name', 'Latitude');
ncwriteatt(nc_file, 'YC', 'units', 'degrees north');
ncwriteatt(nc_file, 'YC', 'CoordinateAxisType', 'YC');
ncwrite(nc_file,'YC',yc');
nccreate(nc_file,'ZC','Dimensions',{'Z', zdim_int});
ncwriteatt(nc_file, 'ZC', 'standard_name', 'ZC');
ncwriteatt(nc_file, 'ZC', 'long_name', 'rho depth');
ncwriteatt(nc_file, 'ZC', 'units', 'm');
ncwriteatt(nc_file, 'ZC', 'CoordinateAxisType', 'ZC');
ncwrite(nc_file,'ZC',depth_zl);
var_str='BATHY';
nccreate(nc_file,var_str,'Dimensions',{'X', xdim,'Y',ydim});
ncwriteatt(nc_file, var_str, 'standard_name', var_str);
ncwriteatt(nc_file, var_str, 'long_name', 'Bathymetry');
ncwriteatt(nc_file, var_str, 'units', 'm');
ncwriteatt(nc_file, var_str, 'missing_value', 'NaN');
ncwrite(nc_file,var_str,bathy2');
end 
% ----------------------- % 

% ----- Dynamic Fields ----- % 
nc_path='/gpfs_backup/he_data/he/tashrops/projects/ROFFS/CNAPS_opendap/netCDF/';
yy=num2str(uni_year(tt),'%02.f');
mm=num2str(uni_mon(tt),'%02.f');
dd=num2str(uni_day(tt),'%02.f');
nc_name=['fields/CNAPS_',yy,mm,dd,'.nc'];
nc_file=[nc_path,nc_name];
eval(['! rm -f ' nc_file]);

var_str='SSH';
nccreate(nc_file,var_str,'Dimensions',{'X', xdim,'Y',ydim});
ncwriteatt(nc_file, var_str, 'standard_name', var_str);
ncwriteatt(nc_file, var_str, 'long_name', 'Sea Surface Height');
ncwriteatt(nc_file, var_str, 'units', 'm');
ncwriteatt(nc_file, var_str, 'missing_value', 'NaN');
ncwrite(nc_file,var_str,zeta');

for vv=1:length(var_name_MOD);
var_str=cell2mat(var_name_MOD(vv));
nccreate(nc_file,var_str,'Dimensions',{'X', xdim,'Y',ydim,'Z',zdim_int});
ncwriteatt(nc_file, var_str, 'standard_name', var_str);
ncwriteatt(nc_file, var_str, 'long_name', cell2mat(var_name_long(vv)));
ncwriteatt(nc_file, var_str, 'units', cell2mat(var_units(vv)));
ncwriteatt(nc_file, var_str, 'missing_value', 'NaN');
ncwrite(nc_file,var_str,field_all(:,:,:,vv));
end 
% ------------------------- % 

% ******************************************** % 
% ******************************************** % 
% ******************************************** % 
end % if del_flag==0
 
toc;
display('Done Reading One Day'); %pause;
display(' *********************** ');
display(' *********************** ');
display(' *********************** ');
display(' '); 

end % for tt=1:length(uni_date)

display('Done Downloading NetCDF'); pause; 
% ======================================== % 
% ======================================== % 
% ======================================== % 
% ---------------------------------------------- % 
% ---------------------------------------------- % 
% ---------------------------------------------- % 
% ####################################################### % 
% ####################################################### % 
% ####################################################### % 
end % if down_flag==1



















