%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  oct_make_bulk: Build a CROCO bulk file using octave-netcdf procedural API
%
clear all
close all
%%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%
%crocotools_param should set filenames and variables used below
crocotools_param
as_consts

sat_file     =[coads_dir,'sat.cdf']; sat_name='sat';
sst_file     =[coads_dir,'sst.cdf']; sst_name='sst';
airdens_file =[coads_dir,'airdens.cdf']; airdens_name='airdens';
_u3_file     =[coads_dir,'u3.cdf']; u3_name='u3';

% use other file names as in original

% Read in the grid (procedural)
ng = netcdf.open(grdname,'NC_NOWRITE');
Lp = netcdf.inqDim(ng, netcdf.inqDimID(ng,'xi_rho'));
Mp = netcdf.inqDim(ng, netcdf.inqDimID(ng,'eta_rho'));
lon = netcdf.getVar(ng, netcdf.inqVarID(ng,'lon_rho'));
lat = netcdf.getVar(ng, netcdf.inqVarID(ng,'lat_rho'));
lonu = netcdf.getVar(ng, netcdf.inqVarID(ng,'lon_u'));
latu = netcdf.getVar(ng, netcdf.inqVarID(ng,'lat_u'));
lonv = netcdf.getVar(ng, netcdf.inqVarID(ng,'lon_v'));
latv = netcdf.getVar(ng, netcdf.inqVarID(ng,'lat_v'));
angle = netcdf.getVar(ng, netcdf.inqVarID(ng,'angle'));
netcdf.close(ng);

cosa = cos(angle); sina = sin(angle);

% Create the forcing file (oct_create_bulk handles procedural create)
oct_create_bulk(blkname,grdname,CROCO_title,coads_time,coads_cycle,Yorig);

% Open bulk file for writing
nc = netcdf.open(blkname,'NC_WRITE');
var_tair = netcdf.inqVarID(nc,'tair');
var_rhum = netcdf.inqVarID(nc,'rhum');
var_prate = netcdf.inqVarID(nc,'prate');
var_radlw = netcdf.inqVarID(nc,'radlw');
var_radlw_in = netcdf.inqVarID(nc,'radlw_in');
var_radsw = netcdf.inqVarID(nc,'radsw');
var_wspd = netcdf.inqVarID(nc,'wspd');
var_uwnd = netcdf.inqVarID(nc,'uwnd');
var_vwnd = netcdf.inqVarID(nc,'vwnd');
var_sustr = netcdf.inqVarID(nc,'sustr');
var_svstr = netcdf.inqVarID(nc,'svstr');

% Loop on time and fill variables using oct_ext_data
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  data = oct_ext_data(sat_file,sat_name,tindex,lon,lat,time,Roa,1);
  netcdf.putVar(nc,var_tair,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],data);
end
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  data = 0.01*oct_ext_data(rh_file,rh_name,tindex,lon,lat,time,Roa,1);
  netcdf.putVar(nc,var_rhum,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],data);
end
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  data = 0.8*oct_ext_data(precip_file,precip_name,tindex,lon,lat,time,Roa,1);
  netcdf.putVar(nc,var_prate,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],data);
end
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  radlw = oct_ext_data(lrf_file,lrf_name,tindex,lon,lat,time,Roa,1);
  netcdf.putVar(nc,var_radlw,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],radlw);
  sst = oct_ext_data(sst_file,sst_name,tindex,lon,lat,time,Roa,1);
  lwup = emiss_lw.*sigmaSB.*((sst+CtoK).^4);
  radlw_in = -(radlw - lwup);
  netcdf.putVar(nc,var_radlw_in,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],radlw_in);
end
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  data = oct_ext_data(srf_file,srf_name,tindex,lon,lat,time,Roa,1);
  netcdf.putVar(nc,var_radsw,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],data);
end

% wind and uwnd/vwnd
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  wspd = oct_ext_data(w3_file,w3_name,tindex,lon,lat,time,Roa,1);
  netcdf.putVar(nc,var_wspd,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],wspd);
end
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  uwnd = oct_ext_data(u3_file,u3_name,tindex,lon,lat,time,Roa,1);
  vwnd = oct_ext_data(v3_file,v3_name,tindex,lon,lat,time,Roa,1);
  u10 = rho2u_2d(uwnd.*cosa + vwnd.*sina);
  v10 = rho2v_2d(vwnd.*cosa - uwnd.*sina);
  netcdf.putVar(nc,var_uwnd,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],u10);
  netcdf.putVar(nc,var_vwnd,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],v10);
end

% stress
for tindex=1:length(coads_time)
  time = netcdf.getVar(nc, netcdf.inqVarID(nc,'bulk_time'), tindex-1, 1);
  tx = oct_ext_data(taux_file,taux_name,tindex,lon,lat,time,Roa,2);
  ty = oct_ext_data(tauy_file,tauy_name,tindex,lon,lat,time,Roa,2);
  sustr = rho2u_2d(tx.*cosa + ty.*sina);
  svstr = rho2v_2d(ty.*cosa - tx.*sina);
  netcdf.putVar(nc,var_sustr,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],sustr);
  netcdf.putVar(nc,var_svstr,[tindex-1 0 0],[1 size(lon,1) size(lon,2)],svstr);
end

netcdf.close(nc);

% plotting calls unchanged (call original test_forcing)
if makeplot==1
  test_forcing(blkname,grdname,'tair',[1 4 7 10],3,coastfileplot)
end

return
