%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script creates three different files that are used to force the
% meteorlogical variables within the ROMS model. Any variables not set here
% are set using the analytical files (ex. ana_rain.h) when the code is
% compiled. The netcdf files are created using the unix/linux ncgen
% command with .cdl files that are specific to each of the three forcing
% files.
%
%           Aaron Bever, November 2009
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

create_files = 1;  % To create the netcdf files (1 for yes, 0 to write into
                   % already created files.

% Load in the Data and convert the variables that need converted                   
load NCDCMeteorPage
tempC = (tempC - 32) ./9 .*5; % converting to celsius
greg_wind=datevec(jd);


% Set the starting and ending time for the files. This is currently set to
% make forcing files that encompass the entire wet and dry runs so there
% are not two sets of files. ROMS should search through the times and get
% the correct ones to match with the modeling time.
%DRY SEASON
%start = [2008 6 1 1 0 0];
%stop = [2008 8 1 0 0 0];
%WET SEASON
%start = [2008 8 1 1 0 0];
%stop = [2008 10 1 0 0 0];
%ALL
start = [2008 3 1 1 0 0];
stop = [2008 11 1 0 0 0];

is = near(jd,datenum(start));
ie = near(jd,datenum(stop));

% The data had bad values and nan's in it. This reinterpolates the data
% back onto it's own time grid, thus replacing the nan's with
% interpolations between the good values.
jd_wind=julian(greg_wind(is:ie,:))-2440000;
wdir=wdir/360*2*pi;
W=ws(is:ie).*exp(sqrt(-1).*wdir(is:ie));
W(1) = 1+i;
    p=find(isnan(W)==0);
    W = interp1(jd_wind(p),W(p),jd_wind);
relhum = relhum(is:ie);
    p=find(isnan(relhum)==0);
    relhum = interp1(jd_wind(p),relhum(p),jd_wind);
%    relhum(:) = mean(relhum);
slp = slp(is:ie);
    p=find(isnan(slp)==0);
    slp = interp1(jd_wind(p),slp(p),jd_wind);
% Precipitation didn't look good, or I wasn't sure how to use it.
%    slp(:) = mean(slp);
% pcp01 = pcp01(is:ie);
%     p=find(isnan(pcp01)==0);
%     pcp01 = interp1(jd_wind(p),pcp01(p),jd_wind);
tempC = tempC(is:ie);
    p=find(isnan(tempC)==0);
    tempC = interp1(jd_wind(p),tempC(p),jd_wind);
%    tempC(:) = mean(tempC);


% Set the file names.
grid_file='../FORCING/CALUS/calus_grid.nc';   % input: ROMS grid file
Fname1='calus_cloud.nc'   % output: ROMS forcing file
Fname2='calus_bulk.nc'
Fname3='calus_wind.nc'
%Fname4='LG_swrad3.nc'  % no radiation in data

% the grid angle is needed to rotate the winds onto ROMS grid coordinates
ang=nc_varget(grid_file,'angle');

if(create_files)
%%%% generate the netcdf file using ncgen
cmd=['! ncgen -o ' Fname1 ' cloud_calus.cdl '];
eval(cmd);

cmd=['! ncgen -o ' Fname2 ' bulk_calus.cdl '];
eval(cmd);

cmd=['! ncgen -o ' Fname3 ' wind_calus.cdl '];
eval(cmd);
end
 
% Add time variable to the netcdf file to initialize it to it's total size
m = length(jd_wind);
nc_varput(Fname3,'wind_time',jd_wind,0,m)

% Adding the wind to the forcing file
disp('Writing to wind file.')
[M,N] = size(ang);
for i=1:length(jd_wind);
 wr(1,:,:)=W(i)*exp(-sqrt(-1)*ang);   % rotate wind to ROMS direction
 nc_varput(Fname3,'Uwind',real(wr),[i-1 0 0],[1 M N]);
 nc_varput(Fname3,'Vwind',imag(wr),[i-1 0 0],[1 M N]);
end
disp(' finished wind')

%---------------------------------------------------

% Setting time variables
nc_varput(Fname2,'qair_time',jd_wind,0,m)
nc_varput(Fname2,'tair_time',jd_wind,0,m)
nc_varput(Fname2,'pair_time',jd_wind,0,m)
nc_varput(Fname1,'cloud_time',jd_wind,0,m)
% Settting Cloud cover at a low fraction (0<C>1);
nc_varput(Fname1,'cloud',ones(m,M,N0.1,[0 0 0],[m M N])

for i=1:length(jd_wind);
 nc_varput(Fname2,'Tair',ones(1,M,N).*tempC(i),[i-1 0 0],[1 M N]);
 nc_varput(Fname2,'Qair',ones(1,M,N).*relhum(i),[i-1 0 0],[1 M N]);
 nc_varput(Fname2,'Pair',ones(1,M,N).*slp(i),[i-1 0 0],[1 M N]);
end

disp(' finished Bulk vars')