Skip to content
Snippets Groups Projects
Commit e85e58fe authored by mey's avatar mey
Browse files

added template-file for eros8

parent 1b8181f5
No related branches found
No related tags found
No related merge requests found
function eros_continue(projectpath)
% This function enables the continuation of a eros/floodos run from the
% time of the last output. The inflow timeseries will also be started from where it was stopped.
% The function allows for changing boundary conditions and/or parameters in the continued run.
% USAGE: 1. navigate into the output folder of the run to be continued
% before calling eros_continue
% 2. projectpath is the path 1 step upwards of the folder where
% the grids are stored, i.e. the folder in which the topo
% folder resides
%
%
% After you have started a floodos/eros run ALWAYS save the LEM
% object as LEM.mat to the output folder. If you haven't then you
% need to construct it again using your run template.
load LEM.mat % load previous model setup
T = dir('*.txt');
[~,index] = sortrows({T.datenum}.');
T = T(index);
T = readtable(T(end).name);
ix = find(~strcmp(T{:,end},'-')); % find the index of the last output
T(ix(end)+1:end,:)=[];
tablename = dir('*.txt');
[~,index] = sortrows({tablename.datenum}.');
tablename = tablename(index);
tablename=tablename(end).name;
writetable(T,tablename); % prune the table to the last output
H = dir('*.alt');
W = dir('*.water');
S = dir('*.sed');
[~,index] = sortrows({H.datenum}.');
H = H(index);
W = W(index);
LEM.dem = grd2GRIDobj(H(end).name,LEM.dem);
LEM.water = grd2GRIDobj(W(end).name,LEM.dem);
LEM.dem.name = 'continue';
if ~isempty(S)
S = S(index);
sediment = grd2GRIDobj(S(end).name,dem);
LEM.sed = sediment;
end
cd(projectpath) %<<< path to your project folder
climate = importdata(LEM.climate,'\t');
climate.data = climate.data(:,~isnan(climate.data(1,:))); % in case NaNs appear during import
interupt_state = (length(ix)-1)*LEM.draw;
unpassed = climate.data(:,1)>=interupt_state;
climate.data = climate.data(unpassed,:); % start timeseries from the last output time
fileID = fopen([projectpath,'\Topo\continue.climate'],'w');
fprintf(fileID, ['time flow:relative\n']);
% fprintf(fileID, [climate.colheaders{1}, ' ', climate.colheaders{2}, ' ', climate.colheaders{3},'\n']);
for i = 1:length(climate.data)
% fprintf(fileID, [num2str(climate.data(i,1)),' ',num2str(climate.data(i,2)),' ' num2str(climate.data(i,3)),'\n']);
fprintf(fileID, [num2str(climate.data(i,1)),' ',num2str(climate.data(i,2)),'\n']);
end
fclose(fileID);
LEM.end = LEM.end-interupt_state; % end time reduces by the already computed time
GRIDobj2grd(LEM.dem,['./Topo/',dem.name,'.alt']);
GRIDobj2grd(LEM.rain,['./Topo/',dem.name,'.rain']);
try
GRIDobj2grd(LEM.sed,['./Topo/',dem.name,'.sed']);
catch nosediment
end
GRIDobj2grd(LEM.water,['./Topo/',dem.name,'.water']);
try
GRIDobj2grd(LEM.cs,['./Topo/',dem.name,'.cs']);
catch nocsgrid
end
try
GRIDobj2grd(LEM.uplift,['./Topo/',dem.name,'.uplift']);
catch nouplift
end
writeErosInputs(LEM);
system([LEM.experiment,'.bat'])
\ No newline at end of file
%--------------------------------------------------------------------------
% PREPARE GRIDS
%--------------------------------------------------------------------------
%
% addpath('.\mfiles')
% ALT (elevation model)
clear
dem=GRIDobj('./Topo/cop30DEM_utm32n_hochrhein_carved_filled.tif');
dem = resample(dem,60);dem = crop(dem);
dem.name = 'hochrhein_60m';
% RAIN (sources (>0) and sinks (-1))
rain = GRIDobj('.\Topo\map_wc2.tif'); % MAP after WorldClim2 (mm/yr/m^2)
rain = resample(rain,dem);
rain = rain/1000; % convert to m/yr/m^2
rain = rain/3600/24/365.25; % convert to m/s/m^2
upstream = GRIDobj('.\Topo\rain_boundaries_ms-1.tif');
upstream = resample(upstream,dem);
upstream.Z(isnan(upstream.Z))=0;
inflow_factor = 2.45;
rain = (rain+upstream)*inflow_factor;
% INFLOWS
inflow_rhine = 346*inflow_factor; % (m^3/s)
% inflow_rhine_y = [273:279]; % y-location (column) of inlet
% inflow_rhine_x = ones(7,1)'*3266; % x-location (row) of inlet
inflow_rhine_y = [137:140]; % y-location (column) of inlet
inflow_rhine_x = ones(4,1)'*1632; % x-location (row) of inlet
inflow_rhine = inflow_rhine/length(inflow_rhine_x)/dem.cellsize.^2; % divide by number of inflow cells and by cellsize^2
rain.Z(inflow_rhine_y,inflow_rhine_x(1)) = ones(length(inflow_rhine_x),1)*inflow_rhine;
inflow_aare = 349*inflow_factor; % (m^3/s)
% inflow_aare_x = [1461:1462]; % y-location (column) of inlet
% inflow_aare_y = ones(2,1)'*998; % x-location (row) of inlet
inflow_aare_x = [731]; % y-location (column) of inlet
inflow_aare_y = ones(length(inflow_aare_x),1)'*498; % x-location (row) of inlet
inflow_aare = inflow_aare/length(inflow_aare_x)/dem.cellsize.^2; % divide by number of inflow cells and by cellsize^2
rain.Z(inflow_aare_y,inflow_aare_x(1)) = ones(length(inflow_aare_x),1)*inflow_aare;
inflow_reuss = 140*inflow_factor; % (m^3/s)
% inflow_reuss_x = [1680:1682]; % y-location (column) of inlet
% inflow_reuss_y = ones(3,1)'*998; % x-location (row) of inlet
inflow_reuss_x = [842]; % y-location (column) of inlet
inflow_reuss_y = ones(length(inflow_reuss_x),1)'*498; % x-location (row) of inlet
inflow_reuss = inflow_reuss/length(inflow_reuss_x)/dem.cellsize.^2; % divide by number of inflow cells and by cellsize^2
rain.Z(inflow_reuss_y,inflow_reuss_x(1)) = ones(length(inflow_reuss_x),1)*inflow_reuss;
inflow_limmat = 114*inflow_factor; % (m^3/s)
% inflow_limmat_x = [1851:1852]; % y-location (column) of inlet
% inflow_limmat_y = ones(2,1)'*998; % x-location (row) of inlet
inflow_limmat_x = [926]; % y-location (column) of inlet
inflow_limmat_y = ones(length(inflow_limmat_x),1)'*498; % x-location (row) of inlet
inflow_limmat = inflow_limmat/length(inflow_limmat_x)/dem.cellsize.^2; % divide by number of inflow cells and by cellsize^2
rain.Z(inflow_limmat_y,inflow_limmat_x(1)) = ones(length(inflow_limmat_x),1)*inflow_limmat;
% rain = resample(rain,dem60);
rain.Z(1:end,1)=-1;
% rain.Z(1,1:end)=-1;
% rain.Z(end,1:end)=-1;
% rivers = GRIDobj('.\Topo\rivers2raster.tif');
% INITIAL SEDIMENT CONCENTRATION
cs = rain;
cs.Z=cs.Z*0;
cs.Z(1,:) = 1;
cs.Z(:,end) = 1;
cs.Z(end,:) = 1;
climate = 'Topo\\boundary_1Ma_GT_clear.climate';
% WATER
% water = GRIDobj('.\Topo\HochRhein_WATER+LAKE_1000m.tif');
% UPLIFT
% uplift = GRIDobj('.\Topo\uplift_m_per_s_baselevel_Basel.tif');
uplift = GRIDobj('.\Topo\uplift_elicited_updated.tif')-5.63e-5; % [m/yr] minus the uplift at the outlet to make it the baselevel
uplift = resample(uplift,dem);
uplift.Z(1,:) = uplift.Z(2,:);
uplift.Z(:,1) = uplift.Z(:,2);
% SED (sediment thickness in meters)
sed = GRIDobj('.\Topo\mqu_140715g_utm32n_aug.tif');
sed = resample(sed,dem);
sed.Z(isnan(sed.Z))=0;
LEM.dem = dem;
LEM.rain = rain;
LEM.sed = sed;
LEM.uplift = uplift;
% LEM.water = water;
LEM.cs = cs;
LEM.climate = climate;
GRIDobj2grd(dem,['./Topo/',dem.name,'.alt']);
GRIDobj2grd(rain,['./Topo/',dem.name,'.rain']);
GRIDobj2grd(sed,['./Topo/',dem.name,'.sed']);
GRIDobj2grd(uplift,['./Topo/',dem.name,'.uplift']);
% GRIDobj2grd(water,['./Topo/',dem.name,'.water']);
GRIDobj2grd(cs,['./Topo/',dem.name,'.cs']);
%--------------------------------------------------------------------------
%% DEFINE INPUT PARAMETERS
%--------------------------------------------------------------------------
LEM.experiment = 'hochrhein'; % Project name
LEM.ErosPath = 'D:\\USER\\mey'; % Path to .exe
LEM.outfolder = 'hochrhein\\60m\\eros8\\clear\\'; % folder to store results in
LEM.eros_version = 'eros8.0.131M';
% LEM.inflow = 1060; % [m3s-1]water inflow at source cells
LEM.rainfall = 1; % Sets the precipitation rate per unit surface when multiplied by the rainfall map
LEM.initial_sediment_stock = '0.01'; % % The total "stock" of sediment at the precipiton landing is: input_sediment_concentration*cs_map[i]*Precipiton_volume
LEM.inertia = 0; % refers to inertia term in shallow water equation
LEM.time_unit = 'year';
LEM.begin = 0; LEM.begin_option = 'time'; % start time
LEM.end = 1e+5; LEM.end_option = 'time'; % length of model run
LEM.draw = 1e+3; LEM.draw_option = 'time'; % output interval
LEM.step = '2.5e2:dir'; LEM.step_option = 'volume:adapt';
LEM.stepmin = 0.1e1;
LEM.stepmax = 1e4;
LEM.initphase ='time:init:begin=10:end=1:step=1:log:op=*:tu=0.0001';
% LEM.TU_coefficient = 0.001; % sets the proportion of rain pixels that make up 1 TU
LEM.TU = 'TU:surface=rain:coefficient=0.001';
LEM.flow_model = 'stationary:pow';
LEM.erosion_multiplier = '1000'; % enhance the erosion by this factor for each precipiton
LEM.erosion_multiplier_increment = '3:log10';
LEM.erosion_multiplier_op = '*';
LEM.erosion_multiplier_min = 1000; % smallest erosion_multiplier during calculation
LEM.erosion_multiplier_max = 1000'; % largest erosion_multiplier during calculation;
LEM.erosion_multiplier_default_min = '10000'; % minimum number of defaults that triggers a time step increase
LEM.erosion_multiplier_default_max = 20000'; % maximum number of defaults that triggers a time step decrease
LEM.time_extension = '2e+3';
LEM.uplift_multiplier = 1;
LEM.limiter_dh = '0.1:dir'; % define how excessive must be the topography variations and whether this limits topography or not
LEM.limiter_mode = 'unclip:dir';
%--------------------------------------------------------------------------
% EROSION/DEPOSITION
%--------------------------------------------------------------------------
LEM.erosion_model = 'MPM'; % (stream_power, shear_stress, shear_mpm)
LEM.deposition_model = 'constant'; % need to know whether there are other options!
LEM.stress_model = 'rgqs';
% ALLUVIAL
LEM.fluvial_stress_exponent = 1.5; % exponent in sediment flux eq. (MPM): qs = E(tau-tau_c)^a
LEM.fluvial_sediment_erodability = '2.6e-8'; % [kg-1.5 m-3.5 s-2] E in MPM equation
LEM.fluvial_sediment_threshold = 0.05; % [Pa] critical shear stress (tau_c) in MPM equation
LEM.deposition_length = '60'; % [m] xi in vertical erosion term: edot = qs/xi
% LATERAL EROSION/DEPOSITION
LEM.fluvial_lateral_erosion_coefficient = '1e-2'; % dimensionless coefficient (Eq. 17 in Davy, Croissant, Lague (2017))
LEM.fluvial_lateral_deposition_coefficient = 1e-2;
LEM.lateral_erosion_model = 1;
LEM.lateral_deposition_model = 'constant';
% BEDROCK
LEM.fluvial_basement_erodability = 0.01;
LEM.fluvial_basement_threshold = 0.5;
LEM.outbend_erosion_coefficient = 1.000000;
LEM.inbend_erosion_coefficient = 1.00000;
LEM.poisson_coefficient = 5;
LEM.diffusion_coefficient = 4;
LEM.sediment_grain = 0.0025;
LEM.basement_grain = 0.025;
% LEM.survey_points = '171130';
%--------------------------------------------------------------------------
% FLOW MODEL
%--------------------------------------------------------------------------
LEM.friction_model = 'manning';
LEM.friction_coefficient = 0.025; %
LEM.flow_boundary = 'free';
%--------------------------------------------------------------------------
% OUTPUTS TO WRITE
%--------------------------------------------------------------------------
% options:
% stress, waters, discharge, downward, slope, qs, capacity, sediment, precipiton_flux, stock
LEM.str_write = 'stress:water:discharge:slope:qs:capacity:sediment:precipiton_flux:stock';
writeErosInputs8(LEM);
%% run model
system([LEM.experiment,'.bat'])
\ No newline at end of file
function writeErosInputs8(LEM)
% Writes the neccessary files for running Eros.exe
fileID = fopen([LEM.experiment,'.arg'],'w');
fprintf(fileID, ['dat=dat\\inputs.dat\n']);
fclose(fileID);
% write input file
fileID = fopen('./dat/inputs.dat','w');
fprintf(fileID, ['erosion_model=',LEM.erosion_model,'\n']);
fprintf(fileID, ['deposition_model=',LEM.deposition_model,'\n']);
% ALLUVIAL
fprintf(fileID, ['deposition_length_fluvial=',num2str(LEM.deposition_length),'\n']);
fprintf(fileID, ['sediment_grain=',num2str(LEM.sediment_grain),'\n']);
% Lateral erosion/deposition
fprintf(fileID, ['lateral_erosion_model=',num2str(LEM.lateral_erosion_model),'\n']);
fprintf(fileID, ['lateral_deposition_model=',num2str(LEM.lateral_deposition_model),'\n']);
fprintf(fileID, ['lateral_erosion_coefficient_fluvial=',num2str(LEM.fluvial_lateral_erosion_coefficient),'\n']);
fprintf(fileID, ['lateral_deposition_coefficient_fluvial=',num2str(LEM.fluvial_lateral_deposition_coefficient),'\n']);
try
fprintf(fileID, ['lateral_erosion_outbend=',num2str(LEM.outbend_erosion_coefficient),'\n']);
fprintf(fileID, ['lateral_erosion_inbend=',num2str(LEM.inbend_erosion_coefficient),'\n']);
catch
end
% BEDROCK
fprintf(fileID, ['poisson_coefficient=',num2str(LEM.poisson_coefficient),'\n']);
fprintf(fileID, ['diffusion_coefficient=',num2str(LEM.diffusion_coefficient),'\n']);
fprintf(fileID, ['basement_grain=',num2str(LEM.basement_grain),'\n']);
fprintf(fileID, ['basement_erodibility=',num2str(LEM.fluvial_basement_erodability),'\n']);
try
fprintf(fileID, ['sediment_erodability=',num2str(LEM.fluvial_sediment_erodability),'\n']);
catch
end
% FLOW MODEL
fprintf(fileID, ['flow_model=',LEM.flow_model,'\n']);
fprintf(fileID, ['friction_coefficient=',num2str(LEM.friction_coefficient),'\n']);
fprintf(fileID, ['flow_boundary=',num2str(LEM.flow_boundary),'\n']);
fprintf(fileID, ['stress_model=',num2str(LEM.stress_model),'\n']);
% Boundary conditions
% Topo
fprintf(fileID, ['topo=Topo\\',LEM.dem.name,'.alt\n']);
if isfield(LEM,'rain')
fprintf(fileID, ['rain=Topo\\',LEM.dem.name,'.rain\n']);
end
if isfield(LEM,'water')
fprintf(fileID, ['water=Topo\\',LEM.dem.name,'.water\n']);
end
if isfield(LEM,'sed')
fprintf(fileID, ['sed=Topo\\',LEM.dem.name,'.sed\n']);
end
if isfield(LEM,'uplift')
fprintf(fileID, ['uplift=Topo\\',LEM.dem.name,'.uplift\n']);
end
if isfield(LEM,'cs')
fprintf(fileID, ['cs=Topo\\',LEM.dem.name,'.cs\n']);
end
if isfield(LEM,'climate')
fprintf(fileID, ['climate=',LEM.climate,'\n']);
end
% INFLOW/RAINFALL CONDITIONS
if isfield(LEM,'inflow')
fprintf(fileID, ['inflow=',num2str(LEM.inflow),':dir\n']);
end
if isfield(LEM,'rainfall')
fprintf(fileID, ['rainfall=',num2str(LEM.rainfall),'\n']);
end
try
fprintf(fileID, ['input_sediment_concentration=',num2str(LEM.initial_sediment_stock),'\n']);
catch
end
try
fprintf(fileID, ['model=',LEM.model,'\n']);
catch
end
% TIME
try % some parameters that are available from eros 7.5.92 onwards
fprintf(fileID, ['time:unit=',LEM.time_unit,'\n']);
fprintf(fileID, ['time_extension=',num2str(LEM.time_extension),'\n']);
catch % in case we use an older eros version
end
fprintf(fileID, ['time:end=',num2str(LEM.end),':',LEM.end_option,'\n']);
fprintf(fileID, ['time:draw=',num2str(LEM.draw),':',LEM.draw_option,'\n']);
fprintf(fileID, ['time:step=',num2str(LEM.step),':',LEM.step_option,'\n']);
fprintf(fileID, ['time:step:min=',num2str(LEM.stepmin),':max=',num2str(LEM.stepmax),'\n']);
fprintf(fileID, [LEM.initphase,'\n']);
fprintf(fileID, ['erosion_multiplier=',num2str(LEM.erosion_multiplier),'\n']);
fprintf(fileID, ['erosion_multiplier_op=',num2str(LEM.erosion_multiplier_op),'\n']);
fprintf(fileID, ['erosion_multiplier_min=',num2str(LEM.erosion_multiplier_min),'\n']);
fprintf(fileID, ['erosion_multiplier_max=',num2str(LEM.erosion_multiplier_max),'\n']);
fprintf(fileID, ['erosion_multiplier_increment=',num2str(LEM.erosion_multiplier_increment),'\n']);
fprintf(fileID, ['erosion_multiplier_default_min=',num2str(LEM.erosion_multiplier_default_min),'\n']);
fprintf(fileID, ['erosion_multiplier_default_max=',num2str(LEM.erosion_multiplier_default_max),'\n']);
fprintf(fileID, ['uplift_rate=',num2str(LEM.uplift_multiplier),'\n']);
try
fprintf(fileID, ['Point=',LEM.survey_points,'\n']);
catch
end
% Default management
fprintf(fileID, ['limiter_dh=',num2str(LEM.limiter_dh),'\n']);
fprintf(fileID, ['limiter_mode=',num2str(LEM.limiter_mode),'\n']);
% fprintf(fileID, 'default:model=all:min=20:max=10000:step=4:op=*:log10\n');
% Save parameters
fprintf(fileID, ['write=',LEM.str_write,'\n']);
fprintf(fileID, [LEM.TU,'\n']); % unknown parameter
fprintf(fileID, ['flow_inertia_coefficient=',num2str(LEM.inertia),'\n']); % inertia in shallow water equation
fprintf(fileID, ['friction_model=',LEM.friction_model,'\n']); % floodos mode
fclose(fileID);
% write .bat file
fileID = fopen([LEM.experiment,'.bat'],'w');
fprintf(fileID, '@rem Run Eros program with following arguments\n');
fprintf(fileID, '@rem\n');
fprintf(fileID, '@echo off\n');
fprintf(fileID, ['@set EROS_PROG=..\\bin\\',LEM.eros_version,'.exe\n']);
fprintf(fileID, '@set COMMAND=%%EROS_PROG%% %%*\n');
fprintf(fileID, '@echo on\n');
fprintf(fileID, '@rem\n\n');
fprintf(fileID, 'goto:todo\n\n');
fprintf(fileID, ':not_todo\n\n');
fprintf(fileID, ':todo\n\n\n');
fprintf(fileID, ['start /LOW %%COMMAND%% -dir=',LEM.outfolder,'\\ ',LEM.experiment,'.arg']);
fclose(fileID);
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment