alex-s-gardner / gemb

A 1D column that simulates snow/firn/ice processes and surface-atmosphere mass and energy exchanges

Documentation

GEMB

A 1D column that simulates snow/firn/ice processes and surface-atmosphere mass and energy exchanges.

The Glacier Energy and Mass Balance model (GEMB, the “B” is silent) is a vertical 1-D column model, i.e. no horizontal communication between model nodes, that simulates atmosphere-surface mass and energy exchanges and the internal evolution of snow, firn and ice. The model shares many characteristics with earlier published firn models that also simulate atmosphere-surface exchanges (e.g. Bassford, 2002; Bougamont & Bamber, 2005; Greuell & Konzelmann, 1994). The model is a finite-difference model with tens to hundreds of layers, the thickness of which are managed dynamically. It is forced at its surface with near-surface (2-10 m) estimates of precipitation, air temperature, wind speed, vapor pressure, surface pressure, and downwelling longwave and shortwave radiation fluxes and optional inputs of solar zenith angle, cloud optical thickness and bare ice albedo. At its bottom boundary, the model applies a constant thermal flux. Internally, the model simulates thermal diffusion, shortwave sub-surface penetration, meltwater retention, percolation and refreeze, effective snow grain size, dentricity, and sphericity, and compaction. In this section we detail specific implementation of various processes, and their options, within the model.

Citation Information

If you use the GEMB software, please cite the following:

• Gardner, A. S., Schlegel, N.-J., and Larour, E.: Glacier Energy and Mass Balance (GEMB): a model of firn processes for cryosphere research, Geosci. Model Dev., 16, 2277–2302, https://doi.org/10.5194/gmd-16-2277-2023, 2023.

If you use GEMB model outputs, please use the citation below:

• Schlegel, N.-J., & Gardner, A. (2024). Output from the Glacier Energy and Mass Balance (GEMB v1.0) forced with 3-hourly ERA5 fields and gridded to 10km, Greenland and Antarctica 1979-2023 (1.3-5day) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.10806250.

gemb's Issues

Unreachable section of code in thermo.m

@alex-s-gardner There's a section of code in thermo.m that says:

if false
    [Define constants one way.]
elseif false 
    [Define the constants a different way.] 
elseif false
    [Define the constants yet another way.] 
else 
    [Define the constants.]
end

The section seems to take the form of a switch that would allow the user to select a method of defining constants, but currently none of the false options can ever be reached. Or perhaps this is the result of trying out different methods during the development of the code before settling on a final option?

I'd like to clean this up, so should I:

1. Delete the unused methods of defining the constants, or
2. Rewrite it as a switch and add optional inputs to the thermo function to allow the user to select which method is used. If this, please tell me how I should refer to each option.

dzMin2 inadvertently multipied by zY2 twice

Need to address this issue in the Matlab version of the code.

dzMin2 defines the minimum allowable cell size:
dzMin2(i) = zY2*dzMin2(i-1).. this simply grows the minimum allowable cell size linearly with depth, starting below the equally spaced "surface layer"

dzMax2 defines the maximum allowable cell size:
dzMax2(i)=max(zY2*dzMin2(i-1),dzMin*2.0) which is equivalent to max(dzMin2(i),dzMin*2.0) note change in indexing

This means that below the "surface layer" dzMax2(i) == dzMin2(i) which is not correct and why there is a extra zY2 in line 110

A more appropriate way would have been to define dzMax2 correctly in the first place:
dzMax2(i)=max(zY2*dzMin2(i), dzMin*2.0)
and to drop the extra zY2 on line 110

But this is all in addition to what you noticed: that the same elements of dzMin2 should be deleted on line 89

Originally posted by @alex-s-gardner in #1 (comment)

Can these files be removed from the repo?

There are a few files in the repo that don't appear to serve a purpose, or can't be used due to dependencies. @NJSchlegel Can you offer guidance on which of these files can be removed from the repo?

1. flux_batchRun.PBS
2. modelSettings.asv
3. testGEMBISSM (requires xy2ll, triangle, setflowequation, etc)
4. loadMetData this function wants to either load files like D_1_daily_met_data.mat which don't exist in the repo, or call a metDataCreate.m function, which also doesn't exist in the repo. This one's a bit weird too, because it includes longwaveDn as a subfunction, which seems like it would be a fine standalone function but isn't used by anything else in GEMB as far as I can tell.

Default emissivity disagreement

The Gardner et al., 2023 manuscript says:

The snow–ice is assumed to be a grey body (default emissivity = 0.97) with the incoming longwave radiation absorbed within the top layer and outgoing longwave radiation calculated using the Stefan–Boltzmann constant ( W m−2 K−4).

However, thermo.m currently sets

emissivity = 1.0; 

by default. Does the default emissivity need to be set to 0.97 in thermo.m?

Ambiguous references in turbulentFlux

The turbulentFlux function references Bougamont 2006, but the Gardner et al. 2023 paper references only the following two Bougamont papers:

Bougamont, M. and Bamber, J. L.: A surface mass balance model for the Greenland Ice Sheet, J. Geophys. Res.-Earth Surf., 110, F04018, https://doi.org/10.1029/2005JF000348, 2005. 

Bougamont, M., Bamber, J. L., Ridley, J. K., Gladstone, R. M., Greuell, W., Hanna, E., Payne, A. J., and Rutt, I.: Impact of model physics on estimating the surface mass balance of the Greenland Ice Sheet, Geophys. Res. Lett., 34, L17501, https://doi.org/10.1029/2007GL030700, 2007. 

The function also references Patterson 1998, which I can't find references for in Gardner 2023 or Gardner 2010.

And there's a reference to Wright 1997, but both of the Wright papers referenced in Gardner 2023 are from the 2000's.

There's also a reference to Dingman 2002, which I can't find anywhere.

@alex-s-gardner Can you provide some clarity about which references should be cited for the formulations in turbulentFlux? I'd like to make sure they're cited properly in the function header and documentation.

Unused M0 and M1 calculations in densification.m

There's a section in densification.m where variables M0 and M1 are calculated a few different ways depending on which albedo method is being used. However, immediately after that section, M0 and M1 are calculated again and overwrite the initial versions of those variables, regardless of which albedo method is used.

It seems like this is an error, but I'm not sure the right way to fix it. Should I delete the albedo-dependent calculations that get overwritten, or handle it some other way?

Change input variable names in ISSM version of gridInitialize

To match a recent update to the Matlab version, the input variables to gridInitialize need to be changed as follows:

zTop -> z_top
dzTop -> dz_top
sMax -> z_max
zY -> beta

The reason for the change is to match the variable names presented in Fig 1 of Gardner et al., 2023.

Also, you'll notice that I introduced a new optional output z_center, which is calculated as:

% Optional output: 
if nargout>1
    z_center = -cumsum(dz) + dz/2; 
end

I'm not sure if z_center will be of use to anyone else, but it was helpful for me when creating the gridInitialize documentation to understand what the function does, and perhaps that means it will be useful for others as well.

Check ISSM version for dzMin2 inadvertently multipied by zY2 twice

After fixing this issue in Matlab version, need to check/update ISSM to address the same issue.

#6 (comment)

Grid stretching algorithm in gridInitialize.m

A comment at the bottom of gridInitialize.m reads:

%% ---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------

Is this an issue that needs to be addressed, or can I delete this comment and close it out?

Confused by adThresh

A threshold value adThresh is used by the albedo and densification functions, and it's described here as:

% adThresh
%  Apply below method to all areas with densities below this value,
%  or else apply direct input value, allowing albedo to be altered.
%  Default value is rho water (1023 kg m-3).

Where does this value come from? I'm wondering if it's a mistake, because 1023 is unrealistically high for freshwater.

Potential indexing error in managelayers.

https://github.com/alex-s-gardner/GEMB/blob/7b05ea5bf9dcb03cdb45c8b3fc26a5a54f1e3a39/GEMB/managelayers.m#L105C21-L105C27

@alex-s-gardner I am confused by Line 105 of managelayers and I wonder if it's a mistake. It says:

dzMax2(i)=max(zY2*dzMin2(i-1),dzMin*2.0);

The index referencing dzMin2 is surprising because dzMin2 can have totally different dimensions from the dzMax2 array. That's because a few lines before this loop, a bunch of cells were combined and array elements were deleted with:

% delete combined cells
D = (m <= Delflag+Wtol);
m(D) = []; W(D) = []; dz(D) = []; d(D) = []; T(D) = []; a(D) = [];
re(D) = []; gdn(D) = []; gsp(D) = []; adiff(D)=[]; EI(D)=[]; EW(D)=[];

In Line 105, dzMin2 is the size of the arrays before cells were deleted, whereas dzMax2 is the size of arrays after cells were deleted.

I am wondering if Line 105 is supposed to reference dzMax2 instead of dzMin2? Or if the same elements of dzMin2 should be deleted on line 89 via:

dzMin2(D) = [];

Can you take a look and let me know what you think?

Unused variables in melt.m

@alex-s-gardner There are a handful of variables in melt.m that never get used (e.g., surplusT). Do I have your blessing to delete them?

dzMin2 indexing in managelayers

This issue is to be addressed in the ISSM version of GEMB.

@alex-s-gardner I've made a few edits to managelayers. Most of the edits have been for efficiency and readability, and have not impacted the output of the function. However, the new edit on line 99 does impact the output of managelayers and the GEMB function.

Can you confirm that the change I made on line 99 is correct or if I should revert it?

Originally posted by @chadagreene in #1 (comment)

densification calculation comment/warning

A comment at the top of the densification function reads:

% THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEEMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPENSATED FOR BY TRACES OF SNOW???]

Is this still a concern, or can I remove the comment?

EI and EW overwritten in melt.m

In the melt.m function, it looks like the variables EI and EW are managed by managelayers, and then immediately overwritten. Here's the section of code I'm talking about:

% Manage the layering to match the user defined requirements
[d, T, dz, W, mAdd, dz_add, addE, a, adiff, m, EI, EW, re, gdn, gsp] = ...
    managelayers(T, d, dz, W, a, adiff, m, EI, EW, dzMin, zMax, zMin, re, gdn, gsp, zTop, zY, CI, LF, CtoK);

%% CHECK FOR MASS AND ENERGY CONSERVATION

% Calculate final mass [kg] and energy [J]
sumER = Rsum * (LF + CtoK * CI);
EI    = m .* T * CI;
EW    = W .* (LF + CtoK * CI);

Is that intentional? My assumption would be that either:

1. The managelayers function should handle doing everything that needs to be done to EI and EW, or
2. EI and EW should be calculated fresh.

It's odd that EI and EW are fed into managelayers, then taken as outputs from managelayers, then immediately overwritten, so I'm wondering if there's a mistake.

An analytical solution for gridInitialize

Currently, gridInitialize creates a grid using a while loop that's kind of hard to understand. What are your thoughts on creating the grid analytically, like this:

function dz = gridInitialize_cag(zTop, dzTop, sMax, zY)
% gridInitialize sets up the initial grid spacing and total grid depth.  
% 
%% Syntax 
% 
%  dz = gridInitialize(zTop, dzTop, sMax, zY)
% 
%% Description
% 
% dz = gridInitialize(zTop, dzTop, sMax, zY) creates a 1D grid structure
% dz, given the inputs: 
% 
%  * zTop: 
%  * dzTop: [m] 
%  * sMax: 
%  * zY: 
% 
% The grid structure is set as constant grid length 'dzTop' for the top
% 'zTop' meters of the model grid. Below 'zTop' the grid length increases
% linearly with depth.
% 
%% Example 
% 
%  zTop  = 10; 
%  dzTop = 0.05; 
%  sMax  = 250; 
%  zY    = 1.1;  
%  dz    = gridInitialize(zTop, dzTop, sMax, zY);
%  
%  plot(dz,'o-') 
%  xlabel 'index number'
%  ylabel 'dz' 
%  box off
%  axis tight
% 
%  z_center = -cumsum(dz) + dz/2; 
%  plot(dz,z_center,'o-')
%  xlabel 'dz (meters)' 
%  ylabel 'depth (meters)' 
%  box off
%  axis([0 max(dz) -sMax 0])

%% Error checks: 

% Calculate number of top grid points:
gpTop = zTop/dzTop;

% Check to see if the top grid cell structure length (dzTop) goes evenly 
% into specified top structure depth (zTop)
assert(mod(gpTop,1)==0,['Top grid cell structure length does not go evenly into ' ...
        'specified top structure depth, adjust dzTop or zTop.'])

% Make sure top grid cell structure length (dzTop) is greater than 5 cm
if dzTop < 0.05
    warning('Initial top grid cell length (dzTop) is < 0.05 m.')
end

%% Create grid: 
% Find the number of indices N where dzTop * zY^N = sMax - zTop, using 
% the following formulation:   
% 
% The integral of dzTop * zY^N {wrt N} = dzTop * zY^N / log(zY)
% 
% and dzTop * zY^N / log(zY) = sMax - zTop 
% 
% where N = log((sMax-zTop) * log(zY)/dzTop)/log(zY).

N = floor(log((sMax-zTop) * log(zY)/dzTop)/log(zY)); 

% Define dzB for the bottom of the column: 
dzB = dzTop * zY .^ (1:N)'; 

% Define dzT for the top of the column: 
dzT = ones(gpTop,1) * dzTop;

% Concatenate the top and bottom: 
dz = [dzT;
      dzB]; 

%% ---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------
% See https://github.com/alex-s-gardner/GEMB/issues/8. 

end

A comparison of the function above vs the original shows the two methods agree within numerical noise:

% Inputs: 
zTop  = 10; 
dzTop = 0.05; 
sMax  = 250; 
zY    = 1.1;  

% Calculate grid: 
dz     = gridInitialize(zTop, dzTop, sMax, zY);
dz_cag = gridInitialize_cag(zTop, dzTop, sMax, zY);

figure
subplot(2,1,1)
plot(dz,'o-') 
hold on
plot(dz_cag,'x-')
xlabel 'index number'
ylabel 'dz' 
box off
axis tight
legend('original','cag','location','northwest')

subplot(2,1,2) 
plot(dz_cag-dz)
box off
axis tight
xlabel 'index number'
ylabel 'difference' 

Philosophically, I prefer the analytical solution, and it will be easier to maintain or translate into other languages compared to the while-loop version, but making the change would produce disagreement with the C++ version on the order of 1e-15.

@alex-s-gardner @NJSchlegel What should our approach be in this case and in others like it?

Logical of a sum of a logical in melt.m

@alex-s-gardner I found an error check in melt.m that checks the logical of the sum of a logical. On line 258 it reads:

if sum(W < 0.0-Wtol)

That works, but it's somewhat unconventional, and I wonder if it's an error. I would expect either:

if any(W < 0.0-Wtol)

or

if sum(W) < (0.0-Wtol)

The same check appears again on line 302.

GEMB input variable names and default values

In reference to Issue #19, I'd like to consider renaming the variables that input to the GEMB function and setting their default values as described below. What are your thoughts on these potential changes?

Current Name	Recommended Name	Rec. Default Value	Units	Description
dateN	time		datenum	Time array corresponding to input precipitation, air_temperature, etc.
P0	precipitation		kg m^-2	Precipitation time series.
Ta0	air_temperature		K	Near-surface (2 m) air temperature time series.
V0	wind_speed		m s^-1	Near-surface (10 m) wind speed time series.
dsw0	downward_shortwave_radiation		W m^-3	Downward shortwave radiation flux time series.
dlw0	downward_longwave_radiation		W m^-2	Downward longwave radiation flux time series.
eAir0	vapor_pressure		Pa	Screen-level vapor pressure time series. (Is this a req'd input, or can it be calculated by GEMB?)
pAir0	air_pressure		Pa	Screen-level air pressure time series.
isrestart	initial_state	[delete]	structure	Restart flag can be handled by delcaring initial_state = restart_structure.
S.Adiffini	diffuse_albedo	diffuse_albedo = albedo_sza50 when albedo_method = "gardner_2009", diffuse_albedo = snow_albedo otherwise.	fraction	Surface albedo for diffuse radiation
S.adThresh	albedo_density_threshold	Inf	kg m^-3	Density threshold below which albedo_method is applied.
S.aIce	ice_albedo	0.48	fraction	Ice albedo.
S.aIdx	albedo_method	"gardner_2009"	n/a	Method of calculating albedo and subsurface absorption.
S.Aini	[delete]	[delete]	fraction	Albedo on restart (delete bc snow_albedo value will be overridden by values in restart_structure)
S.aSnow	snow_albedo	0.85	fraction	New snow albedo.
S.aValue	[not necessary]	[delete]	fraction	Recommend setting albedo to snow_albedo value when albedo_method = "constant".
S.C	snow_accumulation_mean	?	kg m^-2 yr^-1	Mean annual snow accumulation.
S.cciceValue	carbon_concentration_ice	0	ppmw	Concentration of light-absorbing carbon for ice.
S.ccsnowValue	carbon_concentration_snow	0	ppmw	Concentration of light-absorbing carbon for snow.
S.cldFrac	cloud_fraction	0.1	fraction	Cloud amount described as n in Greuell & Konzelmann, 1994 https://doi.org/10.1016/0921-8181(94)90010-8
S.cotValue	cloud_optical_thickness	?	unitless?	Cloud optical thickness.
S.denIdx	densification_method	"arthern_2010"	n/a	Method of calculating densification of snow/firn.
S.Dini	[delete]	[delete]	kg m^-3	Density (delete bc ice_density will be overridden by values in restart_structure)
S.dsnowIdx	snow_density_method	"constant"	n/a	Definition of fresh snow density.
S.dswdiffrf	downward_shortwave_radiation_diffuse	0	W m^-2	Downward diffusive shortwave radiation flux
S.dulwrfValue	?	0	W m^-3	Needs a more descriptive variable name. This variable gets renamed deltaULW in thermo.m.
S.Dzini	?	?	m	The code implies Dzini defines dz, but the S2A1D_000001.mat file contains only values of 0.05 for Dzini, which corresponds to z_top.
S.dzMin	dz_min	dz_top/2	m	Minimum allowable grid spacing.
S.dzTop	dz_top	0.05	m	Spacing of the upper portion of the model grid.
S.ECini	evaporation_condensation	0	kg m^-2	Surface evaporation (negative) or condensation (positive).
S.eIdx	emissivity_method	"constant"	n/a	Method of calculating emissivity.
S.Gdnini	grain_dendricity	1	unitless	Grain dendricity
S.Gspini	grain_sphericity	0.5	unitless	Grain sphericity
S.isdeltaLWup	[delete]	false	boolean	Not necessary as an input, because it can be determined by the value of
S.K	constant_K	days/K	7	Timescale constant from Table 2 of Bougamont et al., 2005 https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005JF000348.
S.outputDir	[not necessary]	[delete]	n/a	Will be included in output_filename
S.outputFreq	output_frequency	"monthly"	n/a
S.Reini	grain_radius	2.5	mm	Effective radius of snow grains.
S.runID	output_filename	none	n/a	Name of output file. Must have either .h5 or .mat extension (e.g., "myfile.h5" or "myfile.mat")
S.Sizeini	?	?	?	Needs a more descriptive variable name.
S.spinUp	spinup_cycles	0	n/a	Number of cycles of met data run before output.
S.swIdx	shortwave_penetration	true	boolean	If true, shortwave penetration is calculated based on albedo_method. Otherwise, all shortwave radiation is absorbed by the top grid cell.
S.szaValue	solar_zenith	0	degree	Solar zenith angle.
S.t0dry	albedo_timescale_dry	30	days	Timescale for dry snow albedo decay from Table 2 of Bougamont et al., 2005
S.t0wet	albedo_timescale_wet	15	days	Timescale for wet snow albedo decay from Table 2 of Bougamont et al., 2005
S.tcIdx	conductivity_method	"sturm_1997"	n/a	Method of calculating thermal conductivity.
S.teThresh	emissivity_radius_threshold	10	mm	Grain radius threshold above which emissivity_method is applied.
S.teValue	emissivity	0.97	unitless	Outward longwave radiation thermal emissivity forcing at every element.
S.ThermoDeltaTScaling	thermo_scaling_factor	1/11	unitless	Scaling factor to multiply the thermal diffusion timestep (delta t).
S.Tini	snow_temperature	mean(air_temperature)	K	Snow temperature
S.Tmean	air_temperature_mean	mean(air_temperature)	K	Mean annual surface temperature.
S.Vmean	wind_speed_mean	0	m s^-1	Mean annual wind speed.
S.Wini	water_content	zeros(size(time))	kg m^-2	Water content
S.zMax	z_max	250	m	Maximum thickness of the total column.
S.zMin	z_min	ceil(z_max/2 /10)*10	m	Minimum model depth.
S.zTop	z_top	10	m	Thickness of the upper portion of the model grid, in which grid spacing is constant.
S.zY	beta	1.025	unitless	Grid cell stretching parameter for the lower portion of the model grid, in which grid length increases linearly with depth.
[not a current option]	ice_density	910	kg m^-3	Density of ice.
[not a current option]	add_cells	1000	integer	Currently defined within GEMB as S.addCells: number of addtional vertical levels.
[not a current option]	snow_density	350	kg m^-3	Density of snow used if snow_density_method = "constant".
[not a current option]	print_status	true	boolean	Print the current status of the calculation to the command window.

Restructure the inputs to GEMB

1. The current implementation of GEMB requires the user to manually define >50 constants, variables, and preferences when calling the GEMB function, and I think most of them can be set as defaults that can be overriden according to the user's preference.
2. Some variable names within GEMB are inconsistent between functions and many names are unintuitive and hard to decipher.
3. @alex-s-gardner has been advocating for the time series inputs to the GEMB function to be in table format, and I think that's a good idea because it will force standardization of inputs, eliminate the chance of users inputting variables in the wrong order, and make it easier for code maintainers to follow what's happening inside each function.

To address the points above, I recommend rewriting the input handling of GEMB using arguments, which will allow GEMB to be called in any of the following forms:

Here's simplest-case example:

out = GEMB(input_table)

Above, input_table is a table containing the time series of minimum required variables to run GEMB (air temperature, precipiation, etc).

Overriding the defaults is easy. Here's the same call as above, but overriding the default ice density:

out = GEMB(input_table, ice_density = 917)

Input parsing with arguments makes it easy for the user to manually override multiple defaults:

out = GEMB(input_table, ...
           ice_density = 917, ...
           snow_density = 340, ...
           densification_method = "ligtenberg_2011", ...
           snow_albedo = 0.87)

The ability to restart a run from a previous state could be included like this or something similar:

out = GEMB(input_table, initial_state = restart_structure)

And to save the output, just define an output_filename and the extension will define whether it's a .mat or scientific data format.

GEMB(input_table, output_filename = "myfile.nc") 

Notice that in addition to setting default values, the form I'm suggesting also include renaming variables to be more intuitive and replacing the method-switching index values with descriptive handles like densification_method = "ligtenberg_2011" (which would be defined in the current implementation by S.denIdx = 6).

What do you think about the idea of restructuring the inputs to GEMB as I'm showing above, and propagating those changes into the functions that are called by GEMB?

Change GEMB output to NetCDF

Currently, the GEMB function saves its output as in .mat format. The .mat format is wonderfully easy to work with in Matlab, and there are packages in most other languages that can read .mat, but I think this may be a good time to move away from the proprietary format.

1. Can we change the output of GEMB to NetCDF format?
2. If yes, which variables should be included in the NetCDF? Would we want to include all of the variables listed in READMEout.txt? Including all of the variables within the structure S?

I'm hoping we can reduce the number of variables that are saved to "as much as necessary and nothing more".

Make output filename an optional input to the GEMB function

Currently, the only way to define the GEMB output data filename is by specifying something like:

S.runPfx = 'S2A1D2';

before calling

GEMB(P0, Ta0, V0, dateN, dlw0, dsw0, eAir0, pAir0, S, isrestart)

This approach buries the filename within S, the runPfx handle isn't very intuitive if you're looking for something like S.outputFilename, and it's hard to find or know how to change the filename (it's currently placed 150 lines away from the GEMB call in MASTER_RUN.m).

There's an additional layer of confusion because in the example above, S2A1D2 is not the complete name of the output file, because S2A1D2 is later automagically converted to S2A1D2_000001 inside a function called combineStrucData_GEMB.m, and then that's used by the GEMB function to create S2A1D2_000001.mat.

The GEMB function also overwrites any files that might have the same automagically created name, without warning.

Because GEMB doesn't output anything to the workspace, as a user I get this feeling that I'm giving GEMB a strangely cobbled together collection of inputs, the function does something for five or ten minutes, and then if all goes well, I might end up with a .mat file somewhere on my computer, by a name I didn't explicitly ask for, and no actual data in the workspace.

I propose changing the GEMB call to look something like this:

outputs = GEMB(P0, Ta0, V0, dateN, dlw0, dsw0, eAir0, pAir0, S, isrestart, outputFile="myfile.nc")

The idea is that

1. GEMB should directly provide the outputs to the workspace if the user wants them, which would provide a more direct conceptual link from function input to output, and
2. only save the data to a file if the user explicitly asks for it (and include a guardrail to protect against overwriting).