Skip to content

Commit

Permalink
Merge pull request #14 from tingeman/master
Browse files Browse the repository at this point in the history
Pull tingemans changes to main CryoGrid repository post-ITCH2021
  • Loading branch information
tingeman authored Dec 10, 2021
2 parents 7956f24 + a06d6a5 commit 543ed0c
Show file tree
Hide file tree
Showing 20 changed files with 1,247 additions and 280 deletions.
136 changes: 136 additions & 0 deletions source/IO/FORCING/FORCING_Tair_Tlb_fixed.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
%========================================================================
% CryoGrid FORCING class FORCING_Tair_Tlb_fixed
% simple model forcing for GROUND classes depending only on an upper
% boundary temperature (Tair) forcing, and a fixed lower boundary
% temperature.
% The forcing data must be stored in a Matlab “.mat” file which contains
% a struct FORCING with field “DATA”, which contain the time series of the actual
% forcing data, e.g. FORCING.data.Tair contains the time series of air temperatures.
% Have a look at the existing forcing files in the folder “forcing” and prepare
% new forcing files in the same way.
%
% The mandatory forcing variables are:
% Tair: Air temperature (Tair, in degree Celsius)
% t_span: Timestamp (as Matlab datenum / increment 1 corresponds to one day)
%
% T_lb is specified in the parameter input file, not in the forcing mat file.
%
% T. Ingeman-Nielsen, S. Westermann, December 2021
%========================================================================

classdef FORCING_Tair_Tlb_fixed < matlab.mixin.Copyable

properties
forcing_index
DATA % forcing data time series
TEMP % forcing data interpolated to a timestep
PARA % parameters
STATUS
CONST
end


methods

function forcing = provide_PARA(forcing)
% INITIALIZE_PARA Initializes PARA structure, setting the variables in PARA.

forcing.PARA.filename = []; % filename of Matlab file containing forcing data
forcing.PARA.forcing_path = []; % location (path) of forcing files
forcing.PARA.start_time = []; % start time of the simulations (must be within the range of data in forcing file)
forcing.PARA.end_time = []; % end time of the simulations (must be within the range of data in forcing file)
forcing.PARA.T_lb = []; % Fixed temperature at lower boundary
end

function forcing = provide_CONST(forcing)

end

function forcing = provide_STATVAR(forcing)

end

function forcing = initialize_excel(forcing)

end


function forcing = finalize_init(forcing, tile)
% FINALIZE_INIT Performs all additional property
% initializations and modifications. Checks for some (but not
% all) data validity.

temp = load([forcing.PARA.forcing_path forcing.PARA.filename], 'FORCING');

forcing.DATA.Tair = temp.FORCING.data.Tair;
forcing.DATA.timeForcing = temp.FORCING.data.t_span;

if std(forcing.DATA.timeForcing(2:end,1)-forcing.DATA.timeForcing(1:end-1,1))>=1e-10
error('timestamp of forcing data is not in regular intervals -> check, fix and restart')
else
forcing.STATUS = 1;
end

% handle start time, if specified
if isempty(forcing.PARA.start_time) || ~ischar(forcing.PARA.start_time)
forcing.PARA.start_time = forcing.DATA.timeForcing(1,1);
else
forcing.PARA.start_time = datenum(forcing.PARA.start_time(1,1), forcing.PARA.start_time(2,1), forcing.PARA.start_time(3,1));
end

% handle end time, if specified
if isempty(forcing.PARA.end_time) || isnan(forcing.PARA.end_time(1,1))
forcing.PARA.end_time = floor(forcing.DATA.timeForcing(end,1));
else
forcing.PARA.end_time = datenum(forcing.PARA.end_time(1,1), forcing.PARA.end_time(2,1),forcing.PARA.end_time(3,1));
end

%initialize TEMP
forcing.TEMP.Tair = 0;
forcing.TEMP.T_lb = forcing.PARA.T_lb;

end

function forcing = interpolate_forcing(forcing, tile)
% Interpolate forcing data to timestep tile.t
t = tile.t;
times = forcing.DATA.timeForcing;
posit = floor((t-times(1,1))./(times(2,1)-times(1,1)))+1;

forcing.TEMP.Tair = forcing.lin_interp(t, posit, times, forcing.DATA.Tair);
forcing.TEMP.T_lb = forcing.PARA.T_lb;
forcing.TEMP.t = t;
end

function xls_out = write_excel(forcing)
% XLS_OUT Is a cell array corresponding to the class-specific content of the parameter excel file (refer to function write_controlsheet).
error('This function is not implemented/updated for this specific class')
xls_out = {'FORCING','index',NaN,NaN;'FORCING_seb',1,NaN,NaN;NaN,NaN,NaN,NaN;'filename',NaN,NaN,NaN;'start_time',NaN,NaN,'provide in format dd.mm.yyyy; if left empty, the first timestamp of the forcing data set will be used';'end_time',NaN,NaN,'provide in format dd.mm.yyyy; if left empty, the last timestamp of the forcing data set will be used';'rain_fraction',1,'[-]','rainfall in forcing file multiplied by this number';'snow_fraction',1,'[-]','snowfall in forcing file multiplied by this number';'latitude',NaN,'[degree]','geographical coordinates';'longitude',NaN,'[degree]',NaN;'altitude',NaN,'[m]','a.s.l.';'domain_depth',100,'[m]','should match a GRID point, model domain extends to this depth';'heatFlux_lb',0.0500000000000000,'[W/m2]','geothermal heat flux';'airT_height',2,'[m]','height of air temperature';'FORCING_END',NaN,NaN,NaN};
end

function fig = plot(forcing)
TT = timetable(datetime(forcing.DATA.timeForcing,'ConvertFrom','datenum'), ...
forcing.DATA.Tair, ...
'VariableNames', {'Tair'});
TT.Properties.VariableUnits = {'degC'};
stackedplot(TT);
%hAx=gca;
%hAx.TickLabelFormat='mm-yyyy';
%datetick('x','mm-yyyy');
end

end


methods(Static)
function value = lin_interp(t, posit, times, data)
% t is the current time
% times is the vector of times at which forcing data are available
% data is the vector of forcing data (on parameter)
% posit is an index into the time vector to the largest specified time step before current time
value = data(posit,1) + (data(posit+1,1) - data(posit,1)).*(t-times(posit,1))./(times(2,1)-times(1,1));
end
end


end
3 changes: 3 additions & 0 deletions source/IO/FORCING/FORCING_tair_precip.m
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,9 @@

forcing.PARA.filename = []; % filename of Matlab file containing forcing data
forcing.PARA.forcing_path = []; % location (path) of forcing files
forcing.PARA.latitude = []; %
forcing.PARA.longitude = []; %
forcing.PARA.altitude = []; %
forcing.PARA.start_time = []; % start time of the simulations (must be within the range of data in forcing file)
forcing.PARA.end_time = []; % end time of the simulations (must be within the range of data in forcing file)
forcing.PARA.rain_fraction = []; % rainfall fraction assumed in sumulations (rainfall from the forcing data file is multiplied by this parameter)
Expand Down
115 changes: 115 additions & 0 deletions source/IO/FORCING/FORCING_ubT_lbT_fixed.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
%========================================================================
% CryoGrid FORCING class FORCING_ubT_lbT_fixed
% simple model forcing for GROUND classes depending only on fixed upper
% and lower boundary temperatures
%
% The mandatory userdefined parameters are:
% start_time
% end_time
% time_step
% T_ub Upper boundary temperature.
% T_lb Lower boundary temperature.
%
% T. Ingeman-Nielsen, December 2021
%========================================================================


classdef FORCING_ubT_lbT_fixed < matlab.mixin.Copyable

properties
forcing_index
DATA % forcing data time series
TEMP % forcing data interpolated to a timestep
PARA % parameters
STATUS
CONST
end


methods

function forcing = provide_PARA(forcing)
% INITIALIZE_PARA Initializes PARA structure, setting the variables in PARA.

forcing.PARA.start_time = []; % start time of the simulations (must be within the range of data in forcing file)
forcing.PARA.end_time = []; % end time of the simulations (must be within the range of data in forcing file)
forcing.PARA.time_step = []; % time step of timeseries [days] (e.g. 0.5 days)
forcing.PARA.T_ub = []; %
forcing.PARA.T_lb = []; %
end

function forcing = provide_CONST(forcing)

end

function forcing = provide_STATVAR(forcing)

end

function forcing = initialize_excel(forcing)

end


function forcing = finalize_init(forcing, tile)
% FINALIZE_INIT Performs all additional property
% initializations and modifications. Checks for some (but not
% all) data validity.

forcing.STATUS = 1;

% handle start time
forcing.PARA.start_time = datenum(forcing.PARA.start_time(1,1), forcing.PARA.start_time(2,1), forcing.PARA.start_time(3,1));

% handle end time
forcing.PARA.end_time = datenum(forcing.PARA.end_time(1,1), forcing.PARA.end_time(2,1),forcing.PARA.end_time(3,1));

% generate time sequence
forcing.DATA.timeForcing = forcing.PARA.start_time:forcing.PARA.time_step:forcing.PARA.end_time;
forcing.PARA.time_zero = datenum(datetime(year(datetime(forcing.PARA.start_time,'ConvertFrom','datenum')),1,1,0,0,0));

%initialize TEMP
forcing.TEMP.snowfall=0;
forcing.TEMP.rainfall=0;
forcing.TEMP.Tair=0;
forcing.TEMP.T_ub=0;
forcing.TEMP.T_lb=0;

end


function forcing = interpolate_forcing(forcing, tile)
% Interpolate forcing data to timestep tile.t
t = tile.t;

forcing.TEMP.T_ub = forcing.PARA.T_ub;
forcing.TEMP.T_lb = forcing.PARA.T_lb;
forcing.TEMP.Tair = forcing.TEMP.T_ub;
forcing.TEMP.t = t;
end

function xls_out = write_excel(forcing)
% XLS_OUT Is a cell array corresponding to the class-specific content of the parameter excel file (refer to function write_controlsheet).
error('This function is not implemented/updated for this specific class')
xls_out = {'FORCING','index',NaN,NaN;'FORCING_seb',1,NaN,NaN;NaN,NaN,NaN,NaN;'filename',NaN,NaN,NaN;'start_time',NaN,NaN,'provide in format dd.mm.yyyy; if left empty, the first timestamp of the forcing data set will be used';'end_time',NaN,NaN,'provide in format dd.mm.yyyy; if left empty, the last timestamp of the forcing data set will be used';'rain_fraction',1,'[-]','rainfall in forcing file multiplied by this number';'snow_fraction',1,'[-]','snowfall in forcing file multiplied by this number';'latitude',NaN,'[degree]','geographical coordinates';'longitude',NaN,'[degree]',NaN;'altitude',NaN,'[m]','a.s.l.';'domain_depth',100,'[m]','should match a GRID point, model domain extends to this depth';'heatFlux_lb',0.0500000000000000,'[W/m2]','geothermal heat flux';'airT_height',2,'[m]','height of air temperature';'FORCING_END',NaN,NaN,NaN};
end

function fig = plot(forcing)
error('This function is not implemented/updated for this specific class')
end

end


methods(Static)
function value = lin_interp(t, posit, times, data)
% t is the current time
% times is the vector of times at which forcing data are available
% data is the vector of forcing data (on parameter)
% posit is an index into the time vector to the largest specified time step before current time
value = data(posit,1) + (data(posit+1,1) - data(posit,1)).*(t-times(posit,1))./(times(2,1)-times(1,1));
end
end


end
14 changes: 11 additions & 3 deletions source/IO/OUT/OUT_all_tagged.m
Original file line number Diff line number Diff line change
Expand Up @@ -74,15 +74,15 @@
% forcing: instance of FORCING class
forcing = tile.FORCING;

out.OUTPUT_TIME = forcing.PARA.start_time + out.PARA.output_timestep;
out.OUTPUT_TIME = forcing.PARA.start_time; % + out.PARA.output_timestep;
if isempty(out.PARA.save_interval) || isnan(out.PARA.save_interval)
out.SAVE_TIME = forcing.PARA.end_time;
else
out.SAVE_TIME = min(forcing.PARA.end_time, datenum([out.PARA.save_date num2str(str2num(datestr(forcing.PARA.start_time,'yyyy')) + out.PARA.save_interval)], 'dd.mm.yyyy'));
end

out.TEMP = struct();

out.TEMP.first_step = true;
end

%---------------time integration-------------
Expand All @@ -106,7 +106,15 @@
% It is time to collect output
% Store the current state of the model in the out structure.

disp([datestr(t)])
if ~out.TEMP.first_step
fprintf(repmat('\b',1,19))
else
fprintf('\n')
out.TEMP.first_step = false;
end
fprintf(datestr(t, 'yyyy-mm-dd HH:MM:SS'))
% disp(datestr(t))

out.TIMESTAMP=[out.TIMESTAMP t];

CURRENT =TOP.NEXT;
Expand Down
12 changes: 11 additions & 1 deletion source/IO/PROVIDER/PROVIDER.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
classdef PROVIDER < PROVIDER_EXCEL & PROVIDER_EXCEL3D & PROVIDER_MAT & PROVIDER_YAML
classdef PROVIDER < PROVIDER_EXCEL & PROVIDER_EXCEL3D & PROVIDER_MAT & PROVIDER_YAML & PROVIDER_YAML3D

methods

Expand All @@ -13,6 +13,8 @@
provider = assign_paths_mat(provider, run_name, result_path);
elseif strcmp(init_format, 'YAML')
provider = assign_paths_yaml(provider, run_name, result_path, constant_file);
elseif strcmp(init_format, 'YAML3D')
provider = assign_paths_yaml3d(provider, run_name, result_path, constant_file);
end
end

Expand All @@ -24,6 +26,8 @@
provider = read_const_excel3d(provider);
elseif strcmp(provider.PARA.init_format, 'YAML')
provider = read_const_yaml(provider);
elseif strcmp(provider.PARA.init_format, 'YAML3D')
provider = read_const_yaml3d(provider);
end
end

Expand All @@ -35,6 +39,8 @@
provider = read_parameters_excel3d(provider);
elseif strcmp(provider.PARA.init_format, 'YAML')
provider = read_parameters_yaml(provider);
elseif strcmp(provider.PARA.init_format, 'YAML3D')
provider = read_parameters_yaml3d(provider);
end
end

Expand All @@ -44,6 +50,8 @@
provider = update_parameter_file_excel3d(provider, param_file_number);
elseif strcmp(provider.PARA.init_format, 'YAML')
provider = update_parameter_file_yaml(provider, param_file_number);
elseif strcmp(provider.PARA.init_format, 'YAML3D')
provider = update_parameter_file_yaml3d(provider, param_file_number);
end
end

Expand All @@ -53,6 +61,8 @@
provider = update_run_name_excel3d(provider, worker_number);
elseif strcmp(provider.PARA.init_format, 'YAML')
provider = update_run_name_yaml(provider, worker_number);
elseif strcmp(provider.PARA.init_format, 'YAML3D')
provider = update_run_name_yaml3d(provider, worker_number);
end
end

Expand Down
3 changes: 1 addition & 2 deletions source/IO/PROVIDER/PROVIDER_EXCEL3D.m
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,7 @@
mkdir([provider.PARA.result_path provider.PARA.run_name]);
end
end




end
end
Expand Down
2 changes: 2 additions & 0 deletions source/IO/PROVIDER/PROVIDER_YAML.m
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,8 @@
else
warning(['Unrecognized compound parameter format. Parameter "' fieldname '" in class "' class(new_class) '" not populated.'])
end
elseif ismatrix(contents) && isempty(contents)
new_class.PARA.(fieldname) = [];
elseif isempty(contents)
new_class.PARA.(fieldname) = NaN;
else
Expand Down
Loading

0 comments on commit 543ed0c

Please sign in to comment.