clear all
close all
%clc
warning off all 
addpath(genpath('/Users/bouch/Documents/MATLAB/'))
addpath(genpath('/Users/bouch/Desktop/'))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
BLK_dir 	   = '/Users/bouch/Desktop/PROJECTS/MODELING/CROCO/Config_test_EPAC/INPUT_FILES/BLK/';
FRC_dir 	   = '/Users/bouch/Desktop/PROJECTS/MODELING/CROCO/Config_test_EPAC/INPUT_FILES/FRC2/';
grdname 	   = '/Users/bouch/Desktop/PROJECTS/MODELING/CROCO/Config_test_EPAC/INPUT_FILES/Grid_refined_UCLAtools/croco_epac12s3_grd.nc';

BLK_prefix    = 'croco_epac12s3_blk_CFSRv2'; %% files pre-processed by process_SODA_outputs.csh (==> DODS format)
FRC_prefix    = 'croco_epac12s3_frc_CFSRv2'; %% files names with bartropic U and V used to make CROCRO bry filers
%
Ymin          = 2015;          % first forcing year
Ymax          = 2017;          % last  forcing year
Mmin          = 1;             % first forcing month
Mmax          = 12;            % last  forcing month
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get Grid info:
% =============
%
nc=netcdf(grdname,'r');
L1=length(nc('xi_psi'));
M1=length(nc('eta_psi'));
close(nc);
Lp1=L1+1;
Mp1=M1+1;


% Loop on the years and the months
for Y=Ymin:Ymax
  if Y==Ymin 
    mo_min=Mmin;
  else
    mo_min=1;
  end
  if Y==Ymax
    mo_max=Mmax;
  else
    mo_max=12;
  end
  for M=mo_min:mo_max
    disp(' ')
    disp(['Processing  year ',num2str(Y),...
    ' - month ',num2str(M)])
    disp(' ')
    %
    Mm=M-1;Ym=Y;
    if Mm==0
Mm=12;
Ym=Y-1;
    end
    Mp=M+1;Yp=Y;
    if Mp==13
Mp=1;
Yp=Y+1;
    end
    %
    blkname=([BLK_dir,BLK_prefix,'_Y',num2str(Y),'M',num2str(M),'.nc']);
    frcname=([FRC_dir,FRC_prefix,'_Y',num2str(Y),'M',num2str(M),'.nc']);


% Get wind stress and time from BLK files:
% ========================================
%

disp('    Get Wind Stress & Time from BLK file')
nc2=netcdf(blkname);
time=nc2{'bulk_time'}(:);
tauy=nc2{'svstr'}(:);
taux=nc2{'sustr'}(:);
close(nc2)

disp('    Create FRC file')

nc = netcdf(frcname, 'clobber');
redef(nc);


%
%  Create dimensions
%

nc('xi_u') = L1;
nc('eta_u') = Mp1;
nc('xi_v') = Lp1;
nc('eta_v') = M1;
nc('xi_rho') = Lp1;
nc('eta_rho') = Mp1;
nc('xi_psi') = L1;
nc('eta_psi') = M1;
nc('sms_time') = length(time);
%
%  Create variables and attributes
%



nc{'sms_time'} = ncdouble('sms_time');
nc{'sms_time'}.long_name = ncchar('surface momentum stress time');
nc{'sms_time'}.long_name = 'surface momentum stress time';
nc{'sms_time'}.units = ncchar('days');
nc{'sms_time'}.units = 'days';

nc{'sustr'} = ncdouble('sms_time', 'eta_u', 'xi_u');
nc{'sustr'}.long_name = ncchar('surface u-momentum stress');
nc{'sustr'}.long_name = 'surface u-momentum stress';
nc{'sustr'}.units = ncchar('Newton meter-2');
nc{'sustr'}.units = 'Newton meter-2';

nc{'svstr'} = ncdouble('sms_time', 'eta_v', 'xi_v');
nc{'svstr'}.long_name = ncchar('surface v-momentum stress');
nc{'svstr'}.long_name = 'surface v-momentum stress';
nc{'svstr'}.units = ncchar('Newton meter-2');
nc{'svstr'}.units = 'Newton meter-2';

endef(nc);
%
% Write time variables
%

% for iii=1:length(time)
% nc{'sms_time'}(iii) = time(iii);
% nc{'shf_time'}(:) = 0;
% nc{'swf_time'}(:) = 0;
% nc{'sst_time'}(:) = 0;
% nc{'srf_time'}(:) = 0;
% nc{'sss_time'}(:) = 0;
% nc{'wwv_time'}(:) = 0;
%
% nc{'sustr'}(iii,:,:) = squeeze(taux(iii,:,:));
% nc{'svstr'}(iii,:,:) = squeeze(tauy(iii,:,:));
% end

for iii=1:length(time)
nc{'sms_time'}(:) = time;
nc{'shf_time'}(:) = 0;
nc{'swf_time'}(:) = 0;
nc{'sst_time'}(:) = 0;
nc{'srf_time'}(:) = 0;
nc{'sss_time'}(:) = 0;
nc{'wwv_time'}(:) = 0;

nc{'sustr'}(:) = taux;
nc{'svstr'}(:) = tauy;

close(nc)

end
end