% Creating met forcing output for ROMS: Create 3 files
create_files = 1;
load NCDCMeteorPage
tempC = (tempC - 32) ./9 .*5; % converting to celsius
greg_wind=datevec(jd);

%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];

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

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);
%    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);


grid_file='../FORCING/CALUS/calus_grid.nc';   % input: ROMS grid file
out_file1='../FORCING/CALUS/calus_cloudD.nc';   % output: ROMS forcing file
out_file2='../FORCING/CALUS/calus_bulkD.nc';
out_file3='../FORCING/CALUS/calus_windD.nc';
%out_file4='LG_swrad3.nc'

ang2 = nc_varget(grid_file,'angle');

Fname1 = out_file1;
Fname2 = out_file2;
Fname3 = out_file3;

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

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

    Fname3 = out_file3
    cmd=['! ncgen -o ' Fname3 ' wind_calusD.cdl '];
    eval(cmd);
end
% Read and write wind
 
disp('started wind')
m = length(jd_wind);
nc_varput(Fname3,'wind_time',jd_wind);

[nc,mc] = size(ang2); %mc = 150;  nc = 40;
for i=1:length(jd_wind);
wr=W(i).*exp(-sqrt(-1)*ang2);   % rotate wind to ROMS direction
WR(1,:,:) = wr;

nc_varput(Fname3,'Uwind',real(WR),[i-1 0 0],[1 nc mc]);
nc_varput(Fname3,'Vwind',imag(WR),[i-1 0 0],[1 nc mc]);
end
clear u v 
disp(' finished wind')
%---------------------------------------------------

if(1)
% write bulk vars

% pause
 nc_varput(Fname2,'qair_time',jd_wind);%-2440000;
 nc_varput(Fname2,'tair_time',jd_wind);%-2440000;
 nc_varput(Fname2,'pair_time',jd_wind);%-2440000;
 nc_varput(Fname1,'cloud_time',jd_wind);%-2440000;
tmp = ones(1,nc,mc);
tmpC = ones(1,nc,mc) .*0.1;
for i=1:length(jd_wind);
    tmp(:) = tempC(i);
 nc_varput(Fname2,'Tair',tmp,[i-1 0 0],[1 nc mc]);
    tmp(:) = relhum(i);
 nc_varput(Fname2,'Qair',tmp,[i-1 0 0],[1 nc mc]);
    tmp(:) = slp(i);
 nc_varput(Fname2,'Pair',tmp,[i-1 0 0],[1 nc mc]);
 nc_varput(Fname1,'cloud',tmpC,[i-1 0 0],[1 nc mc]);
end

disp(' finished Bulk vars')
%---------------------------------------------------
% swrad 


% for i=1:floor(length(jd_wind)/24);
%   f4{'srf_time'}(i)=jd_wind((i-1)*24+1);%-2440000;
%   f4{'swrad'}(i,:,:)=c_data(4,greg_wind((i-1)*24+1,2));
% end
% close(f4);
% disp(' done')
end