%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Add N2O data in input CROCO files from Global Atlas (WOA or CARS)
%
%  N2O distribution from Nevison et al. (2003) formulation
%
%  Data input format (netcdf):
%     variable(T, Z, Y, X)
%     T : time [Months]
%     Z : Depth [m]
%     Y : Latitude [degree north]
%     X : Longitude [degree east]
%
%
%  Elodie Gutknecht, 2013
%  Gildas Cambon, 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
%
%  Title
%
crocotools_param
%%%%%%%%%%%%%%%%%%% END USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%
%
% options for write_time_attributes
insecond =  0 ;
add_cycle = 1 ;
Yorig='';
%
% Initialize Yorig if not provided
if ~exist('Yorig', 'var') ,  Yorig = []; , end
%
% Get time attributes
[time_unit_att,time_second_unit_att,calendar_att]=...
    get_time_attributes(Yorig);
%
%
%=========

if makeini
    nc=netcdf(grdname,'r');
    h = nc{'h'}(:);
    close(nc)
    % O2 initial conditions
    nc=netcdf(ininame,'r');
    theta_s = nc{'theta_s'}(:);
    theta_b = nc{'theta_b'}(:);
    Tcline = nc{'Tcline'}(:);
    N =  length(nc('s_rho'));
    vtransform=nc{'Vtransform'}(:);
    if  ~exist('vtransform')
        vtransform=1; %Old Vtransform
        disp([' NO VTRANSFORM parameter found'])
        disp([' USE TRANSFORM default value vtransform = 1'])
    end
    O2_ini   = nc{'O2'}(:);
    close(nc)
    type = 'initial conditions file' ;
    history = 'CROCO' ;
    [KK,LL,MM]=size(O2_ini);
    %
    zw_ini=zlevs(h,0.,theta_s,theta_b,Tcline,N,'w',vtransform);
    N2O_ini=zeros(KK,LL,MM);N2O_ini=NaN;
    for k=1:KK
        for j=1:LL
            for i=1:MM
                N2O_ini(k,j,i)=nevis_2003(squeeze(zw_ini(k,j,i)),squeeze(O2_ini(k,j,i)));
            end
        end
    end
    % Find NaN
    find(isnan(N2O_ini)==1);%
    %
    % open the ini file
    nc = netcdf(ininame,'write');
    % new variable
    %%redef(nc);
    %
    nc{'N2O'} = ncdouble('time','s_rho','eta_rho','xi_rho') ;
    nc{'N2O'}.long_name = 'Nitrous oxide';
    nc{'N2O'}.fields = 'NO2, scalar, series';
    %%endef(nc);
    % write new variable
    nc{'N2O'}(:,:,:) =  N2O_ini;
    % Synchronize on disk
    close(nc);
end

%==========
if makeclim
    nc=netcdf(grdname,'r');
    h = nc{'h'}(:);
    close(nc)
    % O2 climatological conditions
    nc=netcdf(clmname,'r');
    theta_s = nc{'theta_s'}(:);
    theta_b = nc{'theta_b'}(:);
    Tcline = nc{'Tcline'}(:);
    N =  length(nc('s_rho'));
    vtransform=nc{'Vtransform'}(:);
    if  ~exist('vtransform')
        vtransform=1; %Old Vtransform
        disp([' NO VTRANSFORM parameter found'])
        disp([' USE TRANSFORM default value vtransform = 1'])
    end
    O2_clm  = nc{'O2'}(:);
    O2_time = nc{'o2_time'}(:);
    O2_cycle = nc{'o2_time'}.cycle_length(:);
    zeta = nc{'zeta'}(:);
    close(nc)
    [TT,KK,LL,MM]=size(O2_clm);
    for itps=1:length(O2_time)
        zw_clm(itps,:,:,:)=zlevs(h,squeeze(zeta(itps,:,:)),theta_s,theta_b,Tcline,N,'w',vtransform);
    end
    N2O_clm=zeros(TT,KK,LL,MM); N2O_clm(:)=NaN;
    for t=1:TT
        for k=1:KK
            for j=1:LL
                for i=1:MM
                    zw= squeeze(zw_clm(t,k,j,i));
                    O2=squeeze(O2_clm(t,k,j,i));
                    N2O_clm(t,k,j,i)=nevis_2003(zw,02);
                end
            end
        end
    end
    % Find NaN
    find(isnan(N2O_clm)==1);
    %
    % add N20 in climatological file
    type = 'climatological conditions file' ;
    history = 'CROCO' ;
    %
    % open the clm file
    %
    nc = netcdf(clmname,'write');
    %
    % new variable
    %
    %%redef(nc);
    nc('n2o_time') = TT;
    nc{'n2o_time'} = ncdouble('n2o_time') ;
    nc{'n2o_time'}.long_name = ncchar('time climatological N2O');
    nc{'n2o_time'}.long_name = 'time climatological N2O';
    write_time_attributes(nc,'n2o_time',O2_cycle,time_unit_att,time_second_unit_att,...
        calendar_att,0,1);

    nc{'N2O'} = ncdouble('n2o_time','s_rho','eta_rho','xi_rho') ;
    nc{'N2O'}.long_name = ncchar('Nitrous oxide');
    nc{'N2O'}.long_name = 'Nitrous oxide';
    nc{'N2O'}.units = ncchar('mMol N2O m-3');
    nc{'N2O'}.units = 'mMol N2O m-3';
    nc{'N2O'}.fields = ncchar('NO2, scalar, series');
    nc{'N2O'}.fields = 'NO2, scalar, series';
    %%endef(nc);
    %
    % write new variable
    %
    nc{'n2o_time'}(:) =  O2_time(:);
    nc{'N2O'}(:,:,:,:) = N2O_clm(:);
    %
    % Synchronize on disk
    %
    close(nc);
end

%==========
if makebry
    % O2 boundary conditions
    nc=netcdf(grdname,'r');
    h = nc{'h'}(:);
    close(nc)

    nc=netcdf(bryname,'r');
    theta_s = nc{'theta_s'}(:);
    theta_b = nc{'theta_b'}(:);
    Tcline = nc{'Tcline'}(:);
    N =  length(nc('s_rho'));
    vtransform=nc{'Vtransform'}(:);
    if  ~exist('vtransform')
        vtransform=1; %Old Vtransform
        disp([' NO VTRANSFORM parameter found'])
        disp([' USE TRANSFORM default value vtransform = 1'])
    end
    O2_bry_west = nc{'O2_west'}(:);
    O2_bry_east = nc{'O2_east'}(:);
    O2_bry_south = nc{'O2_south'}(:);
    O2_bry_north = nc{'O2_north'}(:);
    O2_time = nc{'o2_time'}(:);
    O2_cycle = nc{'o2_time'}.cycle_length(:);
    zeta_west  = nc{'zeta_west'}(:);
    zeta_east  = nc{'zeta_east'}(:);
    zeta_south = nc{'zeta_south'}(:);
    zeta_north = nc{'zeta_north'}(:);
    close(nc)

    % BRY
    [TT,KK,MM]=size(O2_bry_south); [TT,KK,LL]=size(O2_bry_west);
    type='w';
    zw_bry_west  = zeros(TT,KK+1);
    zw_bry_east  = zeros(TT,KK+1);
    zw_bry_south = zeros(TT,KK+1);
    zw_bry_north = zeros(TT,KK+1);

    N2O_bry_east  = zeros(TT,KK,LL); N2O_bry_east(:)=NaN  ;
    N2O_bry_west  = zeros(TT,KK,LL); N2O_bry_west(:)=NaN  ;
    N2O_bry_south = zeros(TT,KK,MM); N2O_bry_south(:)=NaN ;
    N2O_bry_north = zeros(TT,KK,MM); N2O_bry_north(:)=NaN ;
    for t=1:TT
        for j=1:LL
            h_west  = squeeze(h(j,1));
            h_east  = squeeze(h(j,end));
            zw_bry_west(t,:,j) = zlevs_1d(h_west , squeeze(zeta_west(t,j)) , theta_s, theta_b, hc, N, type, vtransform);
            zw_bry_east(t,:,j)= zlevs_1d(h_east , squeeze(zeta_east(t,j)) , theta_s, theta_b, hc, N, type, vtransform);
        end
    end
    for t=1:TT
        for i=1:MM
            h_south  = squeeze(h(1,i));
            h_north  = squeeze(h(end,i));
            zw_bry_south(t,:,i)= zlevs_1d(h_south , squeeze(zeta_south(t,i)) , theta_s, theta_b, hc, N, type, vtransform);
            zw_bry_north(t,:,i)= zlevs_1d(h_north , squeeze(zeta_north(t,i)) , theta_s, theta_b, hc, N, type, vtransform);
        end
    end
    for t=1:TT
        for k=1:KK
            for j=1:LL
                zw_east=zw_bry_east(t,k,j);
                zw_west=zw_bry_west(t,k,j);
                O2_east=O2_bry_east(t,k,j);
                O2_west=O2_bry_west(t,k,j);
                N2O_bry_east(t,k,j)= nevis_2003(zw_east,O2_east);
                N2O_bry_west(t,k,j)= nevis_2003(zw_west,O2_west);
            end
        end
    end
    for t=1:TT
        for k=1:KK
            for i=1:MM
                zw_north=zw_bry_north(t,k,i);
                zw_south=zw_bry_south(t,k,i);
                O2_north=O2_bry_north(t,k,i);
                O2_south=O2_bry_south(t,k,i);
                N2O_bry_north(t,k,i)=nevis_2003(zw_north,O2_north);
                N2O_bry_south(t,k,i)=nevis_2003(zw_south,O2_south);
            end
        end
    end
    T= O2_time ;% time in days
    cycle=O2_cycle;

    %
    % add N20 in bryfile file
    %
    type = 'boundary conditions file' ;
    history = 'CROCO' ;
    %
    % open the bry file
    %
    nc = netcdf(bryname,'write');
    % new variable
    %
    %%redef(nc);
    nc('n2o_time') = length(T);
    nc{'n2o_time'} = ncdouble('n2o_time') ;
    nc{'n2o_time'}.long_name = ncchar('time climatological N2O');
    nc{'n2o_time'}.long_name = 'time climatological N2O';
    write_time_attributes(nc,'n2o_time',O2_cycle,time_unit_att,time_second_unit_att,...
                        calendar_att,0,1);

    nc{'N2O_east'} = ncdouble('n2o_time','s_rho','eta_rho') ;
    nc{'N2O_east'}.long_name = ncchar('Nitrous oxide');
    nc{'N2O_east'}.long_name = 'Nitrous oxide';
    nc{'N20_east'}.units = ncchar('mMol N2O m-3');
    nc{'N2O_east'}.units = 'mMol N2O m-3';
    nc{'N2O_east'}.fields = ncchar('NO2, scalar, series');
    nc{'N2O_east'}.fields = 'NO2, scalar, series';

    nc{'N2O_west'} = ncdouble('n2o_time','s_rho','eta_rho') ;
    nc{'N2O_west'}.long_name = ncchar('Nitrous oxide');
    nc{'N2O_west'}.long_name = 'Nitrous oxide';
    nc{'N2O_west'}.units = ncchar('mMol N2O m-3');
    nc{'N2O_west'}.units = 'mMol N2O m-3';
    nc{'N2O_west'}.fields = ncchar('NO2, scalar, series');
    nc{'N2O_west'}.fields = 'NO2, scalar, series';

    nc{'N2O_south'} = ncdouble('n2o_time','s_rho','xi_rho') ;
    nc{'N2O_south'}.long_name = ncchar('Nitrous oxide');
    nc{'N2O_south'}.long_name = 'Nitrous oxide';
    nc{'N20_south'}.units = ncchar('mMol N2O m-3');
    nc{'N2O_south'}.units = 'mMol N2O m-3';
    nc{'N2O_south'}.fields = ncchar('NO2, scalar, series');
    nc{'N2O_south'}.fields = 'NO2, scalar, series';

    nc{'N2O_north'} = ncdouble('n2o_time','s_rho','xi_rho') ;
    nc{'N2O_north'}.long_name = ncchar('Nitrous oxide');
    nc{'N2O_north'}.long_name = 'Nitrous oxide';
    nc{'N20_north'}.units = ncchar('mMol N2O m-3');
    nc{'N2O_north'}.units = 'mMol N2O m-3';
    nc{'N2O_north'}.fields = ncchar('NO2, scalar, series');
    nc{'N2O_north'}.fields = 'NO2, scalar, series';

    %
    % write new variable
    %
    nc{'n2o_time'}(:) =  T;
    nc{'N2O_east'}(:,:,:,:) =  N2O_bry_east;
    nc{'N2O_west'}(:,:,:,:) =  N2O_bry_west;
    nc{'N2O_south'}(:,:,:,:) = N2O_bry_south;
    nc{'N2O_north'}(:,:,:,:) = N2O_bry_north;
    %
    % Synchronize on disk
    %
    close(nc);
end






