clc;clear all;

%addpath(genpath('/home/hipub/code/COAWST_WRF4.0.3/Tools'))
addpath(genpath('/data/brian/Tools'))
%typhoon
% hisnc = '/data/hipub/temp/COAWST_WRF4.0.3/Projects/f3owf_test/20220519.typ/ocean_his_coupled.nc'
% %winter
% hisnc = '/data/hipub/temp/COAWST_WRF4.0.3/Projects/f3owf_test/20220517.winter/ocean_his_coupled.nc'
% %summer
% hisnc = '/data/hipub/temp/COAWST_WRF4.0.3/Projects/f3owf_test/20220507.summer/ocean_his_coupled.nc';
%model = '/data/brian/grid_rom/rom_20m/domin1/backfill_grid/backfill_testmud/coawst_ini.nc';

% hisnc='N_roms_grid.nc'
modelgrid = 'N_roms_grid.nc'
modelwrite = 'N_coawst_ini.nc'

wdr = '/data3/yiwen/matlab_coawst/';

eval(['cd ',wdr])

tic
%get his_ocean grid
numX=ncread(hisnc,'lon_rho');
numY=ncread(hisnc,'lat_rho');
his_lon=ncread(hisnc,'lon_rho',[1 1],[length(numX)-1 1]);
his_lat=ncread(hisnc,'lat_rho',[1 1],[1 length(numY)-1]);
his_time=ncread(hisnc,'ocean_time');
hisu_lon=ncread(hisnc,'lon_u',[1 1],[length(numX)-1 1]);
hisu_lat=ncread(hisnc,'lat_u',[1 1],[1 length(numY)-1]);
hisv_lon=ncread(hisnc,'lon_v',[1 1],[length(numX)-1 1]);
hisv_lat=ncread(hisnc,'lat_v',[1 1],[1 length(numY)-1]);
%lim of 
lon = ncread(modelgrid,'lon_rho');
lat = ncread(modelgrid,'lat_rho');
ulon = lon(1:length(lon(:,1))-1,:);ulat=lat(1:length(lon(:,1))-1,:);
vlon = lon(:,1:length(lat(1,:))-1);vlat=lat(:,1:length(lat(1,:))-1);
%hisu_lat=ncread(hisnc,'lat_u')
%hisu_lon=ncread(hisnc,'lon_u')

%xl=min(min(lon));xr=max(max(lon));
%yb=min(min(lat));yt=max(max(lat));

%xg=his_lon;
%yg=his_lat;
%xg1=hisv_lon;
%yg1=hisv_lat;
%[ii] = find(xg>=xl & xg<=xr);
%[jj] = find(yg>=yb & yg<=yt);

%[iiu] = find(xg1>=xl & xg1<=xr);
%[jju] = find(yg1>=yb & yg1<=yt);

%ocean_time=ncread(hisnc,'ocean_time');
%read data ocean_his_coupled
temp_his = ncread(hisnc,'temp');
salt_his = ncread(hisnc,'salt');
u_his = ncread(hisnc,'u');
v_his = ncread(hisnc,'v');
%coawst_ini
temp_ini = ncread(model,'temp');
salt_ini = ncread(model,'salt');
u_ini = ncread(model,'u');
v_ini = ncread(model,'v');


%[slon,slat] = meshgrid(his_lon(ii),his_lat(jj));
%[uslon,uslat] = meshgrid(hisu_lon(iiu),hisu_lat(jju));
%[vslon,vslat] = meshgrid(his_lon(ii(1:length(ii)-1)),his_lat(jj));


%temp = griddata(his_lon,his_lat,temp_his(1:799,1:799,1,2),lon,lat)
%u_ini = griddata(hisu_lon,hisu_lat,u_his(1:799,1:799,1,2),ulon,ulat)
%v_ini = griddata(hisv_lon,hisv_lat,v_his(1:799,1:799,1,2),vlon,vlat)


%figure;pcolor(his_lon(ii,1),his_lat(1,jj),temp_his(jj,ii,1,2));shading interp;colorbar

%figure;pcolor(ulon,ulat,u_ini);shading interp;colorbar;
%figure;pcolor(hisu_lon(iiu,1),hisu_lat(1,jju),u_his(jju,iiu,1,2));shading interp;colorbar;

%figure;pcolor(vlon,vlat,v_ini);shading interp;colorbar;
%figure;pcolor(hisu_lon(iiu,1),hisv_lat(1,jju),v_his(jju,iiu,1,2));shading interp;colorbar;

for i = (1:15)
    temp_ini(:,:,i) = griddata(his_lon,his_lat,temp_his(1:799,1:799,i,73),lon,lat);
    salt_ini(:,:,i) = griddata(his_lon,his_lat,salt_his(1:799,1:799,i,73),lon,lat);
    u_ini(:,:,i) = griddata(hisu_lon,hisu_lat,u_his(1:799,1:799,i,73),ulon,ulat);
    v_ini(:,:,i) = griddata(hisv_lon,hisv_lat,v_his(1:799,1:799,i,73),vlon,vlat);
end

for i = (1:15)
    temp_ini(:,:,i) = griddata(his_lon,his_lat,temp_his(1:799,1:799,i,5),lon,lat);
    salt_ini(:,:,i) = griddata(his_lon,his_lat,salt_his(1:799,1:799,i,5),lon,lat);
    u_ini(:,:,i) = griddata(hisu_lon,hisu_lat,u_his(1:799,1:799,i,5),ulon,ulat);
    v_ini(:,:,i) = griddata(hisv_lon,hisv_lat,v_his(1:799,1:799,i,5),vlon,vlat);
end


    
ncwrite(modelwrite,'salt',salt_ini);
ncwrite(modelwrite,'temp',temp_ini);
ncwrite(modelwrite,'u',u_ini);
ncwrite(modelwrite,'v',v_ini);

ncwrite(modelwrite,'ocean_time',17)

d50_ini = ncread(modelwrite,'grain_diameter');
g_d = ncread(modelwrite,'grain_density');
s_vel = ncread(modelwrite,'settling_vel');
e_stre = ncread(modelwrite,'erosion_stress');
mf_01 = ncread(modelwrite,'mudfrac_01');


[1.5]/1000
[1.0]/1000; 
2850
ncwrite(modelwrite,'grain_diameter',d50_ini*5);
ncwrite(modelwrite,'grain_density',(g_d/2350)*2850);
ncwrite(modelwrite,'settling_vel',s_vel*0.8)
ncwrite(modelwrite,'erosion_stress',e_stre*0.001);
ncwrite(modelwrite,'mudfrac_01',mf_01+1);
ncwrite(modelwrite,'mudmass_01',mf_01+13250);
ncwrite(modelwrite,'bed_porosity',mf_01+0.5);
ncwrite(modelwrite,'bed_thickness',mf_01+10);


a = ncread(modelwrite,'settling_vel')
sustr1 = ncread(hisnc,'sustr');
svstr1 = ncread(hisnc,'svstr');
sustr = sustr1(:,1:799,69:121);
svstr = svstr1(1:799,:,69:121);
his_time = linspace(0,13,53);
his_time = linspace(0,19,77);
his_time = linspace(0,29,117);

lon_u = ncread(hisnc,'lon_u');
lat_u = ncread(hisnc,'lat_u');
lon = lon_u(1:799,1:799);
lat = lat_u(1:799,1:799);

lon_v = ncread(hisnc,'lon_v');
lat_v = ncread(hisnc,'lat_v');

create_roms_forcings(lon,lat,his_time,'/data3/yiwen/matlab_coawst/climatology_output/N_forcing.nc','sustr','svstr')

