clc
clear all
close all

L = 5e3;

nx  = 40;
% nzs = [8,10,12,16,20,24,32,40,48,64,80,96,128];
nzs = [32,40,48,64,80,96,128];

choices_flow_law = {
  'Glen'
  };

choices_surf_BC_scheme = {
  'twopoint_onesided'
  'threepoint_onesided'
%   'ghostpoint'
  };

% for nzi = 1: length( nzs)
for nzi = 2: length( nzs)
  
  nz = nzs( nzi);

  Deltax = L / nx;
  Deltay = L / nx;
  zeta = linspace( 0,1,nz);
  
  for fi = 1: length( choices_flow_law)
%   for fi = 1
    for si = 1: length( choices_surf_BC_scheme)
%     for si = 3

      choice_flow_law = choices_flow_law{ fi};
      choice_surf_BC_scheme = choices_surf_BC_scheme{ si};
      run_name = ['ISMIP_HOM_C_5_' choice_flow_law '_' choice_surf_BC_scheme '_nz' num2str( nz)];
      
      choice_base_BC_scheme = choice_surf_BC_scheme;
      
      tstart = tic;
      [grid,u,v] = solve_BPA_ISMIP_HOM_C( L, Deltax, Deltay, zeta, choice_flow_law, choice_surf_BC_scheme, choice_base_BC_scheme, nzs);
      dt = toc( tstart);
      disp(['Finished run "' run_name '" in ' num2str( dt) ' seconds'])
      
      filename = [run_name '.mat'];
      save( filename,'grid','u','v');

    end
  end

end