%-----------------------------------------------------------------------
%  Set out river name, location, and runoff direction.
%-----------------------------------------------------------------------
%
%  This script CREATES a river runoff forcing file
%
%  Notice that for now the river identification number is that of the
%  river count.  This variable is needed to fool ROMS generic IO
%  interphase; it can be any real number since will never used.
%
%  The River flag can have any of the following values:
%
%      River.flag(:) = 0,  All Tracer source/sink are off.
%      River.flag(:) = 1,  Only temperature is on.
%      River.flag(:) = 2,  Only salinity is on.
%      River.flag(:) = 3,  Both temperature and salinity are on. 
%      River.flag(:) = 4,  Temperature, salinity and sediment are on.
%
%_________________________________________________________________________

Fname = '../FORCING/swf_riverD.nc'

N=20;       % number of vertical levels
f=3;        % flag for with sediment or not (3 for temp+ salinity, 4 for + sediment.)

% settting up non-time-stepping values specific to the river
Name='River1'; lstr=length(Name);
River.Name(1,1:lstr)=Name;
River.Ypos(1)=79;  
River.Xpos(1)=197;
River.dir(1)=1;
River.num(1)=1;
River.flag(1)=f;
River.qfac(1)=-1;    % to modify the direction of flux.

Name='River2'; lstr=length(Name);
River.Name(2,1:lstr)=Name;
River.Ypos(2)=119;
River.Xpos(2)=197;
River.dir(2)=0;
River.num(2)=2;
River.flag(2)=f;
River.qfac(2)=1;    % to modify the direction of flux.

Nrivers=length(River.dir);

fname = '/work/users/abever/FUGATE/FORCING/swf_windD.nc';
wind_time = nc_varget(fname,'wind_time');
wind_time = [wind_time(1)-1;wind_time];

Nrec=length(wind_time);
% River.mud1=zeros(Nrivers,N,Nrec);
% River.mud2=zeros(Nrivers,N,Nrec);
% River.mud3=zeros(Nrivers,N,Nrec);
% River.mud4=zeros(Nrivers,N,Nrec);
% River.mud5=zeros(Nrivers,N,Nrec);
% for i=1:Nrivers
%     for j=1:N
%         River.mud4(i,j,:)=dis_gpL4*.9;
%         River.mud5(i,j,:)=dis_gpL4*0.1;
%         River.mud6(i,j,:)=dis_gpL4*0.01;
%     end
% end

River.trans = ones(Nrivers,Nrec);
River.trans = River.trans .* 50;

for i=1:Nrivers
River.trans(i,:)=River.qfac(i)*River.trans(i,:);
    
River.vshape(i,N  )=0.10;
River.vshape(i,N-1)=0.10;
River.vshape(i,N-2)=0.10;
River.vshape(i,N-3)=0.10;
River.vshape(i,N-4)=0.10;
River.vshape(i,N-5)=0.10;
River.vshape(i,N-6)=0.10;
River.vshape(i,N-7)=0.075;
River.vshape(i,N-8)=0.075;
River.vshape(i,N-9)=0.05;
River.vshape(i,N-10)=0.05;
River.vshape(i,N-11)=0.01;
River.vshape(i,N-12)=0.01;
River.vshape(i,N-13)=0.01;
River.vshape(i,N-14)=0.01;
River.vshape(i,N-15)=0.01;
end


River.temp = ones([Nrivers,N, Nrec]).*6;
River.salt = ones([Nrivers,N, Nrec]).*0;
River.time = wind_time;

cmd=['! ncgen -o ' Fname ' ./river.cdl'];
eval(cmd);


nc = netcdf(Fname,'write')

nc{'river'}(:)=[River.num];
nc{'river_Xposition'}(:)=River.Ypos(:);
nc{'river_Eposition'}(:)=River.Xpos(:);
nc{'river_direction'}(:)=River.dir(:);
nc{'river_flag'}(:)=River.flag(:);
nc{'river_time'}(1:Nrec)=River.time(:);
nc{'river_transport'}(:)=permute(River.trans,[2 1]);
nc{'river_Vshape'}(:)=permute(River.vshape,[2 1]);
nc{'river_temp'}(:)=permute(River.temp,[3 2 1]);
nc{'river_salt'}(:)=permute(River.salt,[3 2 1]);
% nc{'river_mud_01'}(:)=permute(River.mud1,[3 2 1]);
% nc{'river_mud_02'}(:)=permute(River.mud2,[3 2 1]);
% nc{'river_mud_03'}(:)=permute(River.mud3,[3 2 1]);
% nc{'river_mud_04'}(:)=permute(River.mud4,[3 2 1]); 
% nc{'river_mud_05'}(:)=permute(River.mud5,[3 2 1]); 
% nc{'river_mud_06'}(:)=permute(River.mud6,[3 2 1]); 
% nc{'river_mud7'}(:)=permute(River.mud7,[3 2 1]);
% nc{'river_mud8'}(:)=permute(River.mud8,[3 2 1]); 
nc.date(:)=datestr(now);
close(nc);





