%-----------------------------------------------------------------------
%  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/CALUS/calus_river.nc'
create = 0;

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

% Set if making the river forcing file for the Caloosahatchee River or the
% whole alrger grid.
calus = 1;

% settting up non-time-stepping values specific to the river
if(calus)
    Name='River1'; lstr=length(Name);
    River.Name(1,1:lstr)=Name;
    River.Ypos(1)=149;
    River.Xpos(1)=13;
    River.dir(1)=0;
    River.num(1)=1;
    River.flag(1)=f;
    River.qfac(1)=-1;    % to modify the direction of flux.
else
    % 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.
end

Nrivers=length(River.dir);

%Loading in the data and setting the date
data = load('S79_2008.DAT');
yr = data(:,1);
mn = data(:,2);
day = data(:,3);
fw = data(:,4);

tmp = ones(length(yr),1).*0;
Rtime = julian([yr mn day tmp tmp tmp]) - 2440000;
River.time = Rtime;

Nrec=length(Rtime);

% converting to m3/s
fw = fw .* 12^3;   %ft3 to in3
fw = fw .* 2.54^3; % in3 to mm3
fw = fw ./ 100^3;  % cm3 to m3

% Settting freshwater discharge, no longer using a constant values, instead
% using the acual discharge
%River.trans = ones(Nrivers,Nrec);
%River.trans = River.trans .* 5;
River.trans = fw;

% Set accordingly if using river sediment. Set as suspended sediment
% concentration in kg/m^3 (g/L).
% 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

% Setting the vertical freshwater distribution, here setting it so 1/N*100
% percent of the total discharge comes from each verticle layer. This is
% the same as setting a uniform with depth river discharge.
River.vshape = ones(N,Nrivers) .* 1/N
for i=1:Nrivers
    
% Setting the river direction. River dir and the sign of the freshwater
% discahrge dictate the direction of river discharge within the model grid.
River.trans(i,:)=River.qfac(i)*River.trans(i,:);

% An alternate way to specify the verticle discharge profile
% 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

% Setting temperature and Salinity.
River.temp = ones(Nrec,Nrivers,N).*6;
River.salt = ones(Nrec,Nrivers,N).*0;

% Creating the forcing file.
if(create)
    cmd=['! ncgen -o ' Fname ' ./river_calus.cdl'];
    eval(cmd);
end

% Setting the river time.
for i=1:length(Rtime);
    nc_varput(Fname,'river_time',River.time(i),i-1,1);
end

% Writing the values to the netcdf forcing file.
nc_varput(Fname,'river',[River.num]);
nc_varput(Fname,'river_Xposition',River.Ypos(:));
nc_varput(Fname,'river_Eposition',River.Xpos(:));
nc_varput(Fname,'river_direction',River.dir(:));
nc_varput(Fname,'river_flag',River.flag(:));
nc_varput(Fname,'river_transport',River.trans);
nc_varput(Fname,'river_Vshape',River.vshape);
tmp = ones(
nc_varput(Fname,'river_temp',permute(River.temp,[1 3 2]));
nc_varput(Fname,'river_salt',permute(River.salt,[1 3 2]));
% 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);