%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create and fill CROCO clim and bry files with OGCM SODA data.
%
% The on-line reference to SODA is at
% http://iridl.ldeo.columbia.edu./SOURCES/.CARTON-GIESE/.SODA/
% The on-line reference to ECCO is at
% http://ecco.jpl.nasa.gov/cgi-bin/nph-dods/datasets/
%
%  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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%oct_start % to be used in batch mode %
clear all
close all
%%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%
%
% Common parameters
%
crocotools_param

%
%  SODA DODS URL
%
% SODA_2.2.4/ [ C20R-2 1871-2008 / POP2.1 ]
%%url='http://iridl.ldeo.columbia.edu/SOURCES/.CARTON-GIESE/.SODA/.v2p2p4' ;
url='http://apdrc.soest.hawaii.edu:80/dods/public_data/SODA/soda_pop2.2.4' ;
%
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, 'SODA')
  error('Error: OGCM must be "SODA" 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...')
  oct_download_SODA(Ymin,Ymax,Mmin,Mmax,lonmin,lonmax,latmin,latmax,...
                OGCM_dir,OGCM_prefix,url,Yorig)
end
%
%------------------------------------------------------------------------------------
%
% 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(['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);
  oct_nc_add_globatt(ininame,Yorig,Mmin,Dmin,Hmin,Min_min,Smin,OGCM);
  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)
  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(sprintf(Mth_format,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=30; % monthly files (SODA..)
        itolap_a=1; itolap_p=1;
        itolap_tot=itolap_a + itolap_p;
        disp(['Reduced overlap for monthly SODA 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];
        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,M)),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