%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Create and fill CROCO clim and bry files with OGCM
%  MERCATOR data.
%
%  On crocotools_param.m, available datasets:
%    SODA     -> Description :  1/4 deg SODA v2.2.4 monthly means
%                                                           (187101 -> 201012)
%    mercator -> 2 options for mercator_type parameter:
%
%                1 -->  1/12 deg Mercator global reanalysis (1993   -> )
%                                           (Monthly or Daily product)
%                2 -->  1/12 deg Mercator global analysis   (202206 -> )
%       
%            Note for 2:
%                    (7days from now -> now-1day : Near Real Time analysis)
%                    (more than 15days from now  : Best analysis)
%
%  Online reference to MERCATOR is at http://marine.copernicus.eu
%
%  Further Information:
%  http://www.croco-ocean.org
%
%  This file is part of CROCOTOOLS
%
%  CROCOTOOLS is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published
%  by the Free Software Foundation; either version 2 of the License,
%  or (at your option) any later version.
%
%  CROCOTOOLS is distributed in the hope that it will be useful, but
%  WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%  MA  02111-1307  USA
%
%  Copyright (c) 2005-2006 by Pierrick Penven
%  e-mail:Pierrick.Penven@ird.fr
%
%  Contributions of P. Marchesiello (IRD), J. Lefevre (IRD),
%                   and F. Colberg (UCT)
%
%  Updated    6-Sep-2006 by Pierrick Penven
%  Update     13-Sep-2009 by Gildas Cambon (IRD)
%  Update     14-March-2011 by Gildas Cambon & Serena Illig (IRD)
%  Update     24-November-2020 by Gildas Cambon(IRD)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%oct_start % to be used in batch mode %
clear all
close all
%%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%
%
% Common parameters
%
crocotools_param
%
disp('=============================================')
disp('Take care of OGCM defined in crocotools_param')
disp('=============================================')
%
%  MERCATOR : see oct_get_file_python_mercator.m for python-motu, url, and others options
%
itolap_tot=itolap_a + itolap_p;
disp(['Overlap before =',num2str(itolap_a)])
disp(['Overlap after =',num2str(itolap_p)])
disp(['Total overlap =',num2str(itolap_tot)])
disp(['...'])
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of user input  parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(['==================='])
disp([' Compute time axis '])
disp(['==================='])
if level==0
    nc_suffix='.nc';
else
    nc_suffix=['.nc.',num2str(level)];
    grdname=[grdname,'.',num2str(level)];
end
%
% Get the model grid
%
ncid = netcdf.open(grdname, 'NC_NOWRITE');
lon=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon_rho'));
lat=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat_rho'));
angle=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'angle'));
h=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'h'));
netcdf.close(ncid);
%
if ~strcmp(OGCM, 'mercator')
    error('Error: OGCM must be "mercator" in crocotools_param.m. Exiting script.');
  end
%
% Extract data over the internet
%
if Download_data==1
    %
    % Get the model limits
    %
    lonmin=min(min(lon));
    lonmax=max(max(lon));
    latmin=min(min(lat));
    latmax=max(max(lat));
    %
    % Download data with DODS (the download matlab routine depends on the OGCM)
    %
    disp('Download data...')
    % Loop on the years
    %
    for Y=Ymin:Ymax
        disp(['Processing year: ',num2str(Y)])
        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
            thedatemonth=['Y',num2str(Y),'M',num2str(M)];
            oct_download_mercator_python(pathCMC,user,password,mercator_type,...
                                    product_id, ...
                                    Y,M, ...
                                    lonmin,lonmax,latmin,latmax,hmax, ...
                                    OGCM_dir,OGCM_prefix,thedatemonth,Yorig)
        end
    end  %End loop over month and years for the downloading with python/motu client
end %End loop for Download_data
%
%------------------------------------------------------------------------------------
%
% Get the OGCM grid
%
ncid = netcdf.open([OGCM_dir,OGCM_prefix,'Y',num2str(Ymin),'M',num2str(Mmin),'.cdf'], 'NC_NOWRITE');
lonT=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lonT'));
latT=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latT'));
lonU=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lonU'));
latU=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latU'));
lonV=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lonV'));
latV=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'latV'));
Z=-netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'depth'));
NZ=length(Z);
NZ=NZ-rmdepth;
Z=Z(1:NZ);
netcdf.close(ncid);
%
% Initial file
% (the strategy is to oct_start at the begining of a month)
% it is possible to do some temporal interpolation...
% but I am too lazy. lets oct_start the first day of
% month Mmin of year Ymin... with the first data available.
%
if makeini==1
    if  ~exist('vtransform')
        vtransform=1; %Old Vtransform
        disp([' NO VTRANSFORM parameter found'])
        disp([' USE vtransform default value  Vtransfor = 1'])
    end
    ininame=[ini_prefix,'Y',num2str(Ymin),'M',num2str(sprintf(Mth_format,Mmin)),nc_suffix];
    %
    % Process the time in Yorig time (i.e days since Yorig-01-01)
    %
    tini=datenum(Ymin,Mmin,1)-datenum(Yorig,1,1);
    disp([' '])
    disp(['=> Create an initial file for ',datestr(tini+datenum(Yorig,1,1));])
    oct_create_inifile(ininame,grdname,CROCO_title,...
                    theta_s,theta_b,hc,N,...
                    tini,'clobber', vtransform,Yorig);
    ncid_ini = netcdf.open(ininame, 'NC_WRITE');
    oct_interp_OGCM(OGCM_dir,OGCM_prefix,Ymin,Mmin,Roa,interp_method,...
                lonU,latU,lonV,latV,lonT,latT,Z,1,...
                ncid_ini,[],lon,lat,angle,h,1,obc,vtransform)
    oct_nc_add_globatt(ininame,Yorig,Mmin,Dmin,Hmin,Min_min,Smin,OGCM);
    netcdf.close(ncid_ini);
end
%
% Clim and Bry files
%
if makeclim==1 | makebry==1
    if  ~exist('vtransform')
        vtransform=1; %Old Vtransform
        disp([' NO VTRANSFORM parameter found'])
        disp([' USE vtransform default value  Vtransfor = 1'])
    end
    %
    % 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
            %
            % Add 2 times step in the CROCO files: 1 at the beginning and 1 at the end
            %
            ncid = netcdf.open([OGCM_dir,OGCM_prefix,'Y',num2str(Y),'M',num2str(M),'.cdf'], 'NC_NOWRITE');
            OGCM_time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'time'));
            ntimes=length(OGCM_time);
            if ntimes==1
                dt = eomday(Y, M); % monthly files
                %disp(['Number of days of the given month: ',num2str(dt)]);
                if strcmp(OGCM, 'mercator')
                    dt = eomday(Y, M) + 1 ;  % add +1 days to avoid problems 
                                             % with mercator starting at first days month   
                end
                itolap_a=1; itolap_p=1;
                itolap_tot=itolap_a + itolap_p;
                disp(['Reduced overlap for monthly files'])
                disp(['...'])
            else
                dt=max(gradient(OGCM_time));
            end
            %
            %Fill the time axis
            %
            croco_time=0*(1:ntimes+itolap_tot);
            %Current month
            croco_time(itolap_a+1:end-itolap_p)=OGCM_time;
            %
            %Previous  month
            %
            disp(['==================================='])
            for aa= 1:itolap_a
                disp(['Compute beginning overlap, time index:',num2str(aa)])
                disp(['Add ',num2str(-(itolap_a + 1 - aa)), ' timestep dt'])
                disp(['--------'])
                croco_time(aa) = croco_time(itolap_a+1) - ((itolap_a + 1 - aa).* dt);
            end
            %
            %Next month
            %
            disp(['==================================='])
            for aa= 1:itolap_p
                disp(['Compute end overlap, time index:',num2str(ntimes+itolap_tot - itolap_p + aa)])
                disp(['Add ',num2str(aa), ' timestep dt'])
                disp(['--------'])
                croco_time(end - itolap_p +  aa   ) = croco_time(end - itolap_p) +  aa.* dt;
            end
            disp(['==================================='])
            netcdf.close(ncid);
            %-----------------------------------------------------
            %
            % Create and open the CROCO files
            %
            if makebry==1
                bryname=[bry_prefix,'Y',num2str(Y),...
                    'M',num2str(sprintf(Mth_format,M)),nc_suffix];
                disp([' '])
                disp(['=> Create bry file : ',  bryname ])
                oct_create_bryfile(bryname,grdname,CROCO_title,[1 1 1 1],...
                                theta_s,theta_b,hc,N,...
                                croco_time,0,'clobber',vtransform,Yorig);
                oct_nc_add_globatt(bryname,Yorig,Mmin,Dmin,Hmin,Min_min,Smin,OGCM);
                ncid_bry = netcdf.open(bryname, 'NC_WRITE');
            else
                ncid_bry=[];
            end
            if makeclim==1
                clmname=[clm_prefix,'Y',num2str(Y),...
                        'M',num2str(sprintf(Mth_format,M)),nc_suffix];
                oct_create_climfile(clmname,grdname,CROCO_title,...
                                theta_s,theta_b,hc,N,...
                                croco_time,0,'clobber',vtransform,Yorig);
                oct_nc_add_globatt(clmname,Yorig,Mmin,Dmin,Hmin,Min_min,Smin,OGCM);
                ncid_clm = netcdf.open(clmname, 'NC_WRITE');
            else
                ncid_clm=[];
            end
            %
            % Check if there are OGCM files for the previous Month
            %
            fname=[OGCM_dir,OGCM_prefix,'Y',num2str(Ym),'M',num2str(Mm),'.cdf'];
            if exist(fname)==0
                disp(['   No data for the previous month: using current month'])
                Mm=M;
                Ym=Y;
                tndx_OGCM=ones(itolap_a,1);
            else
                ncid = netcdf.open(fname, 'NC_NOWRITE');
                tndx_OGCM=[(netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 'time'))- (itolap_a -1) ):1: (netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 'time')))];
                netcdf.close(ncid);
            end
            %
            % Perform the interpolations for the previous month
            %
            disp(' Previous month :')
            disp('=================')
            for aa=1:itolap_a
                disp(['Beg overlap # ', num2str(aa),' -> tindex ',num2str(aa)])
                disp(['It. of prev month used for it= ',num2str(tndx_OGCM(aa))])
                oct_interp_OGCM(OGCM_dir,OGCM_prefix,Ym,Mm,Roa,interp_method,...
                            lonU,latU,lonV,latV,lonT,latT,Z,tndx_OGCM(aa),...
                            ncid_clm,ncid_bry,lon,lat,angle,h,aa,obc,vtransform)
            end
            %
            % Perform the interpolations for the current month
            %

            disp(' Current month :')
            disp('================')
            for tndx_OGCM=1:ntimes
                disp([' Time step : ',num2str(tndx_OGCM),' of ',num2str(ntimes),' :'])
                oct_interp_OGCM(OGCM_dir,OGCM_prefix,Y,M,Roa,interp_method,...
                            lonU,latU,lonV,latV,lonT,latT,Z,tndx_OGCM,...
                            ncid_clm,ncid_bry,lon,lat,angle,h,tndx_OGCM+itolap_a,obc,vtransform)
            end
            %
            % Read the OGCM file for the next month
            %
            fname=[OGCM_dir,OGCM_prefix,'Y',num2str(Yp),'M',num2str(Mp),'.cdf'];
            if exist(fname)==0
                disp(['   No data for the next month: using current month'])
                Mp=M;
                Yp=Y;
                for aa=1:itolap_p
                    tndx_OGCM(aa)=ntimes;
                end
            else
                for aa=1:itolap_p
                    tndx_OGCM(aa)=aa;
                end;
            end
            %
            % Perform the interpolations for the next month
            %
            disp(' Next month :')
            disp('=============')
            for aa=1:itolap_p
                disp(['End Overlap #',num2str(aa),' -> tindex ',num2str(ntimes+itolap_a+aa)])
                disp(['It. of next month used for it= ',num2str(tndx_OGCM(aa))])
                oct_interp_OGCM(OGCM_dir,OGCM_prefix,Yp,Mp,Roa,interp_method,...
                            lonU,latU,lonV,latV,lonT,latT,Z,tndx_OGCM(aa),...
                            ncid_clm,ncid_bry,lon,lat,angle,h,ntimes+itolap_a+aa,obc,vtransform)
            end
            %
            % Close the CROCO files
            %
            if ~isempty(ncid_clm)
                netcdf.close(ncid_clm);
            end
            if ~isempty(ncid_bry)
                netcdf.close(ncid_bry);
            end
            %
        end
    end
end
%
% Spin-up: (reproduce the first year 'SPIN_Long' times)
% just copy the files for the first year and change the time
%
if SPIN_Long>0
    %
    % Initial file
    %
    if makeini==1
        %
        % Copy the file
        %
        ininame=[ini_prefix,'Y',num2str(Ymin),'M',num2str(sprintf(Mth_format,Mmin)),nc_suffix];
        ininame2=[ini_prefix,'Y',num2str(Ymin-SPIN_Long),'M',num2str(sprintf(Mth_format,Mmin)),nc_suffix];
        disp(['Create ',ininame2])
        eval(['!cp ',ininame,' ',ininame2])
        %
        % Change the time
        %
        ncid = netcdf.open(ininame2, 'NC_WRITE');
        time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'scrum_time'))-365.*SPIN_Long*(24*3600);
        netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'scrum_time'), time);
        netcdf.close(ncid);
    end
    %
    M=Mmin-1;
    Y=Ymin-SPIN_Long;
    for month=1:12*SPIN_Long
        M=M+1;
        if M==13
            M=1;
            Y=Y+1;
        end
        %
        % Climatology files
        %
        if makeclim==1
            %
            % Copy the file
            %
            clmname=[clm_prefix,'Y',num2str(Ymin),'M',num2str(sprintf(Mth_format,M)),nc_suffix];
            clmname2=[clm_prefix,'Y',num2str(Y),'M',num2str(sprintf(Mth_format,M)),nc_suffix];
            disp(['Create ',clmname2])
            eval(['!cp ',clmname,' ',clmname2])
            %
            % Change the time
            %
            ncid = netcdf.open(clmname2, 'NC_WRITE');
            time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'tclm_time'))-365.*(Ymin-Y);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'tclm_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'temp_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'sclm_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'salt_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'uclm_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'vclm_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'v2d_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'v3d_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'ssh_time'), time);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'zeta_time'), time);
            netcdf.close(ncid);
        end
        %
        % Boundary files
        %
        if makebry==1
            %
            % Copy the file
            %
            bryname=[bry_prefix,'Y',num2str(Ymin),'M',num2str(sprintf(Mth_format,M)),nc_suffix];
            bryname2=[bry_prefix,'Y',num2str(Y),'M',num2str(sprintf(Mth_format,M)),nc_suffix];
            disp(['Create ',bryname2])
            eval(['!cp ',bryname,' ',bryname2])
            %
            % Change the time
            %
            ncid = netcdf.open(bryname2, 'NC_WRITE');
            time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bry_time'))-365.*(Ymin-Y);
            netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'bry_time'), time);
            netcdf.close(ncid);
        end
    end
end
%---------------------------------------------------------------
% Make a few plots
%---------------------------------------------------------------
if makeplot==1
    disp(' ')
    disp(' Make a few plots...')
    if makeini==1
        ininame=[ini_prefix,'Y',num2str(Ymin),'M',num2str(sprintf(Mth_format,Mmin)),nc_suffix];
        figure
        oct_test_clim(ininame,grdname,'temp',1,coastfileplot)
        figure
        oct_test_clim(ininame,grdname,'salt',1,coastfileplot)
    end
    if makeclim==1
        clmname=[clm_prefix,'Y',num2str(Y),'M',num2str(sprintf(Mth_format,M)),nc_suffix];
        figure
        oct_test_clim(clmname,grdname,'temp',1,coastfileplot)
        figure
        oct_test_clim(clmname,grdname,'salt',1,coastfileplot)
    end
    if makebry==1
        bryname=[bry_prefix,'Y',num2str(Y),'M',num2str(sprintf(Mth_format,M)),nc_suffix];
        figure
        oct_test_bry(bryname,grdname,'temp',1,obc)
        figure
        oct_test_bry(bryname,grdname,'salt',1,obc)
        figure
        oct_test_bry(bryname,grdname,'u',1,obc)
        figure
        oct_test_bry(bryname,grdname,'v',1,obc)
    end
end