function oct_reformat_ECMWF(Ymin,Ymax,Mmin,Mmax,lonmin,lonmax,latmin,latmax,...
    ECMWF_dir,Yorig,My_ECMWF_dir)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extract data from global oct_netcdf format saved on my computer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Illig, 2010, from oct_download_NCEP
%  Updated    January 2016 (E. Cadier and S. Illig)
%

if nargin < 1
    Ymin=2007;
    Ymax=2008;
    Yorig=1900;
    Mmin=12;
    Mmax=1;
    lonmin=10;
    lonmax=30;
    latmin=-40;
    latmax=-20;
    ECMWF_dir='DATA/NCEP_Peru/';
    My_ECMWF_dir='../NCEP_REA1/';
end
%
% Definitions of names and directories
%
disp(['===================='])
disp(['Direct CONVERT Procedure'])
disp(['===================='])

disp(['ECMWF directory:'])
ecmwf_url=My_ECMWF_dir
catalog={'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...    
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ... 
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    'EI_ecmwf_' ...
    ''};
vnames={'SSTK'...  % surface land-sea mask [1=land; 0=sea]
    'T2M' ...      % 2 m temp. [k]
    'SSTK' ...     % surface temp. [k]
    'U10M' ...     % 10 m u wind [m/s]
    'V10M' ...     % 10 m v wind [m/s]
    'Q' ...        % 2 m specific humidity [kg/kg]
    'SWH' ...      % significant wave height [m]
    'MWD' ...      % mean wave direction [degree]
    'PP1D' ...     % peak wave period [s]
    'STR' ...      % surface thermal radiation [w/m^2.s] -> [w/m^2]
    'STRD' ...     % surface downward thermal radiation [w/m^2.s] -> [w/m^2]
    'SSR' ...      % surface solar radiation [w/m^2.s] -> [w/m^2]
    'TP' ...       % surface precipitation rate [m]->[kg/m^2/s]
    'EWSS' ...     % east-west surface stress [N m-2 s] -> [N m-2]
    'NSSS' ...     % north-south surface stress [N m-2 s] -> [N m-2]
    };
fnames={'sst'...  % surface land-sea mask [1=land; 0=sea]
    't2m' ...      % 2 m temp. [k]
    'sst' ...     % surface temp. [k]
    'u10' ...     % 10 m u wind [m/s]
    'v10' ...     % 10 m v wind [m/s]
    'q' ...        % 2 m specific humidity [kg/kg]
    'swh' ...      % significant wave height [m]
    'mwd' ...      % mean wave direction [degree]
    'pp1d' ...     % peak wave period [s]
    'str' ...      % surface thermal radiation [w/m^2.s] -> [w/m^2]
    'strd' ...     % surface downward thermal radiation [w/m^2.s] -> [w/m^2]
    'ssr' ...      % surface solar radiation [w/m^2.s] -> [w/m^2]
    'tp' ...       % surface precipitation rate [m]->[kg/m^2/s]
    'ewss' ...     % east-west surface stress [N m-2 s] -> [N m-2]
    'nsss' ...     % north-south surface stress [N m-2 s] -> [N m-2]
    };

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Common OpenDAP FTP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp([' '])
disp(['Get ECMWF data from ',num2str(Ymin),' to ',num2str(Ymax)])
disp(['From ',ecmwf_url]);
disp([' '])
disp(['Minimum Longitude: ',num2str(lonmin)])
disp(['Maximum Longitude: ',num2str(lonmax)])
disp(['Minimum Latitude: ',num2str(latmin)])
disp(['Maximum Latitude: ',num2str(latmax)])
disp([' '])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the directory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['Making output data directory ',ECMWF_dir])
eval(['!mkdir ',ECMWF_dir])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find a subset of the ECMWF grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ifile=[char(ecmwf_url),char(catalog(1)),char(vnames(1)),'_',sprintf('%04d',Ymin),sprintf('%02d',Mmin),'.nc'];
[i1min,i1max,i2min,i2max,i3min,i3max,jmin,jmax,lon,lat]=...
 oct_get_ECMWF_subgrid(ifile,lonmin,lonmax,latmin,latmax);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global loop on variable names
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:length(vnames)
    disp(['=========================='])
    disp(['VNAME IS ',char(vnames(k))]);
    disp(['=========================='])
    %
    % 
    vname=char(vnames(k));
    fname=char(fnames(k));
    if k==1 % Dealing with the mask (using SST)
        %
        disp(['==========================']);
        disp(['Get the Land Mask tindex = 1']);
        disp(['Get the Land Mask by using SSTK']);
        disp([' '])
        %
        ifile=[char(ecmwf_url),char(catalog(k)),vname,'_',sprintf('%04d',Ymin),sprintf('%02d',Mmin),'.nc'];
        var=oct_extract_ECMWF(ifile,fname,Ymin,Mmin,20,i1min,i1max,i2min,i2max,i3min,i3max,jmin,jmax);
        var(var == min(min(var)))=NaN;
        mask=var-var;
        oct_write_ECMWF_Mask([ECMWF_dir,'land_Y',num2str(Ymin),'M',num2str(Mmin),'.nc'],...
		      'land',lon,lat,mask)
        disp([' '])
        %
    end
    %
    % Loop on the years
    %
    for Y=Ymin:Ymax
        disp(['=========================='])
        disp(['Processing year: ',num2str(Y)])
        disp(['=========================='])
        %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        % Loop on the months
        %
        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
            
            if k<=9   % sst(mask), t2m, sst, u10, v10, q, swh, mwd, pp1d
                
                file =[ECMWF_dir,vname,'_Y',num2str(Y),'M',num2str(M),'.nc'];
                ifile=[char(ecmwf_url),char(catalog(k)),vname,'_',sprintf('%04d',Y),sprintf('%02d',M),'.nc'];
                %
                var=oct_extract_ECMWF(ifile,fname,Y,M,':',i1min,i1max,i2min,i2max,i3min,i3max,jmin,jmax);
                %                
                if (k==1 | k==3)
                 %--> do running daily mean with filter function
                 b=[0.5 1 1 1 0.5]/4.;
                 a=[1 0 0 0 0];
                 var(isnan(var))=9999;
                 var=filtfilt(b,a,var);
                 var(var==9999)=NaN;
                end
                var2=var(5:end-4,:,:);clear var;
                %
                time2=datenum(Y,M,1)-datenum(Yorig,1,1)+0.25*[1:size(var2,1)]-0.25;
                oct_write_ECMWF(file,vname,lon,lat,time2,var2,Yorig);
                clear time;clear var;clear time2;clear var2;
                disp([' ']);
                                        
            else  % str, strd, ssr, tp, ewss, nsss
                
                file=[ECMWF_dir,vname,'_Y',num2str(Y),'M',num2str(M),'.nc'];
                %
                ifile3  = [char(ecmwf_url),char(catalog(k)),vname,'_',sprintf('%04d',Y),sprintf('%02d',M),'_step3.nc'];
                ifile9  = [char(ecmwf_url),char(catalog(k)),vname,'_',sprintf('%04d',Y),sprintf('%02d',M),'_step9.nc'];
                ifile15 = [char(ecmwf_url),char(catalog(k)),vname,'_',sprintf('%04d',Y),sprintf('%02d',M),'_step15.nc'];
                %
                div   = (86400./4.);
                if strcmp(vname,'TP') 
                    div=div/1000.;
                end
                if strcmp(vname,'STR') 
                    div=-1.*div;
                end    
                var3  = oct_extract_ECMWF(ifile3,fname,Y,M,':',i1min,i1max,i2min,i2max,i3min,i3max,jmin,jmax)/div;  %divide by 6h;
                var9  = oct_extract_ECMWF(ifile9,fname,Y,M,':',i1min,i1max,i2min,i2max,i3min,i3max,jmin,jmax)/div;  %divide by 6h;
                var15 = oct_extract_ECMWF(ifile15,fname,Y,M,':',i1min,i1max,i2min,i2max,i3min,i3max,jmin,jmax)/div; %divide by 6h;
                %
                ncid3 = netcdf.open(ifile3, 'NC_NOWRITE');
                time3  = squeeze(netcdf.getVar(ncid3, netcdf.inqVarID(ncid3, 'time')));
                netcdf.close(ncid3);
                %
                var(1:2:(length(time3)-1)*2-1,:,:)=var15(1:end-1,:,:)-var9(1:end-1,:,:);
                var(2:2:(length(time3)-1)*2,:,:)=var9(2:end,:,:)-var3(2:end,:,:);
                %
                var2=var(3:end-4,:,:);
                %
                time2=datenum(Y,M,1)-datenum(Yorig,1,1)+0.25*[1:size(var2,1)]-0.25;
                %
                clear var3; clear var9; clear var15;clear div;         
                %
                oct_write_ECMWF(file,vname,lon,lat,time2,var2,Yorig);
                clear time; clear var;
                disp([' ']);
            end % end if
        end % end loop month
    end % end loop year
end % loop k
%
return