-
Notifications
You must be signed in to change notification settings - Fork 11
/
Copy pathCryoGrid3_xice_mpi.m
executable file
·407 lines (337 loc) · 25.2 KB
/
CryoGrid3_xice_mpi.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
% ------------------------------------------------------------------------------
% CryoGRID3
% main script for running the model
%
% Development by: S. Westermann and M. Langer 2015
%
% Development of the parallelized version by: L. Martin and J. Nitzbon 2017-2019
%
% ------------------------------------------------------------------------------
delete( gcp( 'nocreate' ) );
add_modules; %adds required modules
number_of_realizations=3;
if number_of_realizations>1 && isempty( gcp('nocreate') )
parpool(number_of_realizations);
end
spmd
index=labindex; % number identifying the process; change this to e.g. 1 for single realization (non-parallel) run
%---------------define input parameters------------------------------------
% here you provide the ground stratigraphy
% z w/i m o type porosity
% default stratigraphy used in publication:
PARA.soil.layer_properties=[ 0.0 0.80 0.05 0.15 1 0.80 ;...
0.1 0.80 0.15 0.05 1 0.80 ;...
0.4 0.80 0.15 0.05 1 0.55 ;...
3.0 0.50 0.50 0.00 1 0.50 ;...
10.0 0.30 0.70 0.00 1 0.30 ];
% soil stratigraphy
% column 1: start depth of layer (first layer must start with 0) - each layer extends until the beginning of the next layer, the last layer extends until the end of the model domain
% column 2: volumetric water+ice content
% column 3: volumetric mineral content
% column 4: volumetric organic content
% column 5: code for soil type: 1: sand, 2: silt
% column 6: natural porosity - should be the same as 1-mineral-organic if no ground subsidence/thermokarst occurs
%------ model parameters --------------------------------------------------
% parameters related to soil
PARA.soil.albedo=0.2; % albedo snow-free surface
PARA.soil.epsilon=0.97; % emissvity snow-free surface
PARA.soil.z0=1e-3; % roughness length [m] snow-free surface
PARA.soil.rs=50; % surface resistance against evapotransiration [m^-1] snow-free surface
PARA.soil.Qgeo=0.05; % geothermal heat flux [W/m2]
PARA.soil.kh_bedrock=3.0; % thermal conductivity of the mineral soil fraction [W/mK]
% parameters related to hydrology scheme
PARA.soil.fieldCapacity=0.5; % water holding capacity of the soil - this must be adapted to fit the upperlost layers!!
PARA.soil.evaporationDepth=0.10; % depth to which evaporation occurs - place on grid cell boundaries
PARA.soil.rootDepth=0.20; % depth affected by transpiration - place on grid cell boundaries
PARA.soil.ratioET=0.5; % 1: only transpiration; 0: only evaporation, values in between must be made dependent on LAI, etc.
PARA.soil.externalWaterFlux=0.0; % external water flux / drainage in [m/day]
PARA.soil.convectiveDomain=[]; % soil domain where air convection due to buoyancy is possible -> start and end [m] - if empty no convection is possible
PARA.soil.mobileWaterDomain=[0 10.0]; % soil domain where water from excess ice melt is mobile -> start and end [m] - if empty water is not mobile
PARA.soil.relative_maxWater=10.0; % depth at which a water table will form [m] - above excess water is removed, below it pools up
PARA.soil.hydraulic_conductivity = 1e-5;% subsurface saturated hydraulic conductivity assumed for lateral water fluxes [m/s]
PARA.soil.hydraulic_conductivity_subs = 1e-5;
PARA.soil.hydraulic_conductivity_surf = 1e-5;
PARA.soil.infiltration_limit_depth=10.0;% maxiumum depth [m] from the surface to which infiltration occurse
PARA = loadSoilTypes( PARA ); % load the soil types ( silt, sand, water body )
% parameters related to snow
PARA.snow.max_albedo=0.85; % albedo of fresh snow
PARA.snow.min_albedo=0.5; % albedo of old snow
PARA.snow.epsilon=0.99; % surface emissivity snow
PARA.snow.z0=5e-4; % roughness length surface [m]
PARA.snow.rs=0.0; % surface resistance -> should be 0 for snow
PARA.snow.rho_snow=250; % density in [kg/m3]
PARA.snow.tau_1=86400.0; % time constants of snow albedo change (according to ECMWF reanalysis) [sec]
PARA.snow.tau_a=0.008; % [per day]
PARA.snow.tau_f=0.24; % [per day]
PARA.snow.relative_maxSnow=[2.0]; % maximum snow depth that can be reached [m] - excess snow is removed in the model - if empty, no snow threshold
PARA.snow.extinction=25.0; % light extinction coefficient of snow [1/m]
% parameters related to water body on top of soil domain
PARA.water.albedo=0.07; % albedo water (parameterization after Wayne and Burt (1954) in surfaceCondition.m)
PARA.water.epsilon=0.99; % surface emissivity water
PARA.water.rs=0.0; % surface resistance -> should be 0 for water
PARA.water.z0=5e-4; % roughness length surface [m] % value for summer / vegetation
PARA.ice.albedo=0.20; % albedo ice / Lei et al. (2011) shows a range of 0.1 to 0.35
PARA.ice.epsilon=0.98; % surface emissivity snow
PARA.ice.rs=0.0; % surface resistance -> should be 0 for ice
PARA.ice.z0=5e-4; % roughness length surface [m] % value for snow
PARA.ice.extinction=4.5; % [m^-1] light extinction coefficient of ice / Lei et al. (2011) shows a range of 1 to 5 m^-1
% parameters related to lateral sediemnt transport / erosion
% set hillslope diffusivities and critical angle
PARA.soil.weight_diffusion = 1.0; % weighting factor for diffusive sediment transport
PARA.soil.weight_advection = 1.0; % weighting factor for advective sediment transport
PARA.soil.hillslope_diffusivity_land = 3e-10; % in [m^2/sec] approx. 0.01 m^2/yr, reference: [ Kessler et al. 2012, JGR ]
PARA.soil.hillslope_diffusivity_water = 3e-8; % in [m^2/sec] approx 1.0 m^2/yr, reference: [ Kessler et al. 2012, JGR ]
PARA.soil.critical_hillslope_angle = pi/4.0; % critical hillslope angle for which advective fluxes diverge. 45° for unfrozen conditions [ Kessler et al. 2012, JGR ]
% technical parameters
PARA.technical.z=2.0; % height of input air temperature above ground in [m] - assumed constant even when snow depth increases
PARA.technical.SWEperCell=0.005; % SWE per grid cell in [m] - determines size of snow grid cells
PARA.technical.maxSWE=0.4; % in [m] SWE
PARA.technical.arraySizeT=5002; % number of values in the look-up tables for conductivity and capacity
PARA.technical.starttime=datenum(2000,10,1); % starttime of the simulation - if empty start from first value of time series
PARA.technical.endtime=datenum(2001,10,1); % endtime of the simulation - if empty end at last value of time series
PARA.technical.minTimestep=0.1 ./ 3600 ./ 24; % smallest possible time step in [days] - here 0.1 seconds
PARA.technical.maxTimestep=300 ./ 3600 ./ 24; % largest possible time step in [days] - here 300 seconds
PARA.technical.targetDeltaE=1e5; % maximum energy change of a grid cell between time steps in [J/m3] %1e5 corresponds to heating of pure water by 0.025 K
PARA.technical.outputTimestep= 3 ./ 24.0; % output time step in [days] - here three hours
PARA.technical.syncTimeStep = 3 ./ 24.0; % lateral exchange time step in [days] - here three hours
PARA.technical.saveDate='01.01.'; % date of year when output file is written - no effect if "saveInterval" is empty
PARA.technical.saveInterval=[1]; % interval [years] in which output files are written - if empty the entire time series is written - minimum is 1 year
PARA.technical.waterCellSize=0.02; % default size of a newly added water cell when water ponds below water table [m]
% subsurface grid
PARA.technical.subsurfaceGrid = [[0:0.02:2], [2.05:0.05:4.0], [4.1:0.1:10], [10.2:0.2:20], [21:1:30], [35:5:50], [60:10:100], [200:100:1000]]'; % the subsurface K-grid in [m]
% parameters related to the specific location
PARA.location.area=1.0; % area represented by the run [m^2] (here a dummy value of 1.0 is set which is overwritten for individual tiles)
PARA.location.initial_altitude=20.0; % altitude in [m a.s.l.]
% dynamic auxiliary varaibles (stored in the PARA struct for technical reasons)
PARA.location.surface_altitude=PARA.location.initial_altitude; % refers to the surface including water body and snow
PARA.location.altitude = PARA.location.initial_altitude; % refers to the terrain surface, including water but excluding snow; used to generate pressure forcing based on barometric altitude formula, if pressure forcing is not given
PARA.location.soil_altitude = PARA.location.initial_altitude; % refers to the soil surface, excluding water body and snow
PARA.location.infiltration_altitude = nan; % defined at runtime
PARA.location.water_table_altitude = nan; % defined at runtime
PARA.soil.infiltration_limit_altitude=PARA.location.initial_altitude-PARA.soil.infiltration_limit_depth; % absolute altitude to which infiltration occurs, updated when ground subsides
PARA.location.bottomBucketSoilcTIndex = nan; % defined at runtime
% common thresholds
PARA.location.absolute_maxWater_altitude = PARA.location.altitude + PARA.soil.relative_maxWater;
if isempty( PARA.snow.relative_maxSnow )
PARA.location.absolute_maxSnow_altitude = [];
else
PARA.location.absolute_maxSnow_altitude = [ PARA.location.altitude + PARA.snow.relative_maxSnow ];
end
%initial temperature profile -> first column depth [m] -> second column temperature [degree C]
PARA.Tinitial = [ -2 5 ;...
0 0 ;...
2 -2 ;...
5 -7 ;...
10 -9 ;...
25 -9 ;...
100 -8 ;...
1100 10.2 ]; % the geothermal gradient for Qgeo=0.05W/m^2 and K=2.746W/Km is about 18.2 K/km
PARA = loadConstants( PARA ); % load natural constants and thermal properties of soil constituents into the PARA struct
%FORCING data mat-file
PARA.forcing.filename='samoylov_ERA_obs_fitted_1979_2014_spinup_extended2044.mat'; %must be in subfolder "forcing" and follow the conventions for CryoGrid 3 forcing files
PARA.forcing.rain_fraction=1; % scaling factor applied to the entire snowfall forcing data
PARA.forcing.snow_fraction=1; % scaling factor applied to the entire snowfall forcing data
PARA.forcing.snow_scaling=1.0; % scaling factor for incoming snowfall of individual tile, used to emulate lateral snow redistribution
% switches for modules
PARA.modules.infiltration=1; % true if infiltration into unfrozen ground occurs
PARA.modules.xice=1; % true if thaw subsicdence is enabled
PARA.modules.lateral=1; % true if adjacent realizations are run (this does not require actual lateral fluxes)
if PARA.modules.lateral
% switches for lateral processes
PARA.modules.exchange_heat = 1;
PARA.modules.exchange_water = 1;
PARA.modules.exchange_snow = 1;
PARA.modules.exchange_sediment = 1;
%---------overwrites variables for each realization--------------------
% this function must define everything that is realization-specific or dependent of all realizations
PARA = get_parallel_variables( PARA );
end
% ------make output directory (name depends on parameters) ----------------
saveDir = './runs';
run_number = sprintf( [ 'TESTRUN_' datestr( PARA.technical.starttime, 'yyyymm' ) '-' datestr(PARA.technical.endtime, 'yyyymm' ) ] );
mkdir([ saveDir '/' run_number]);
%--------------------------------------------------------------------------
%-----------do not modify from here onwards--------------------------------
%--------------------------------------------------------------------------
[FORCING, success]=load_forcing_from_file(PARA); % load FORCING mat-file
if ~success
warning('A problem with the Forcing occured.');
end
PARA = initializeParameters(PARA, FORCING); %set start time, etc.
%----------------create and initialize the grids --------------------------
GRID=makeGrids(PARA); %create all grids
GRID=createStratigraphy(PARA,GRID); %interpolate input stratigraphy to the soil grid
%----- initializie excess ground ice --------------------------------------
[GRID,PARA] = initializeExcessIce(GRID,PARA);
%----- initializie soil thermal properties --------------------------------
GRID = initializeSoilThermalProperties(GRID, PARA);
%------ initializie snow properties----------------------------------------
GRID = initializeSnow(GRID);
%---- initialize the surface energy balance struct ------------------------
SEB = initializeSEB();
%---- initialize the water body module ------------------------------------
GRID = initializeLAKE(GRID);
%---- initialize temperature profile --------------------------------------
T = inititializeTemperatureProfile_simple(GRID, PARA, FORCING);
%---- modification for infiltration
wc=GRID.soil.cT_water;
GRID.soil.water2pool=0;
%---- modification for lateral erosion
GRID.soil.residualSediment = 0;
GRID.soil.residualOrganic = 0;
GRID.soil.residualMineral = 0;
%---- preallocate temporary arrays for capacity and conductivity-----------
[c_cTgrid, k_cTgrid, k_Kgrid, lwc_cTgrid] = initializeConductivityCapacity(T,wc, GRID, PARA); % this is basically the same as "getThermalProperties" during integration, but without interpolation to K grid
%---- energy and water balance initialization -----------------------------
BALANCE = initializeBALANCE(T, wc, c_cTgrid, lwc_cTgrid, GRID, PARA);
%__________________________________________________________________________
%-------- provide arrays for data storage ---------------------------------
[t, TEMPORARY] = generateTemporary(T, PARA);
OUT = generateOUT();
iSaveSettings( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_settings.mat'] , FORCING, PARA, GRID)
disp('initialization successful');
%% ________________________________________________________________________
% Time Integration Routine I
% I
%_________________________________________________________________________I
%try
while t<PARA.technical.endtime
%------ interpolate forcing data to time t ----------------------------
FORCING = interpolateForcingData(t, FORCING);
%------determine the thermal properties of the model domains ----------
[c_cTgrid, k_cTgrid, k_Kgrid, lwc_cTgrid] = getThermalPropertiesInfiltration(T, wc, c_cTgrid, k_cTgrid, k_Kgrid, lwc_cTgrid, GRID, PARA);
%------- water and energy balance calculations ------------------------
BALANCE = updateBALANCE(T, wc, c_cTgrid, lwc_cTgrid, BALANCE, GRID, PARA);
%------ surface energy balance module ---------------------------------
%set surface conditions (albedo, roughness length, etc.)
[PARA, GRID] = surfaceCondition(GRID, PARA, T);
%calculate the surface energy balance
[SEB, dwc_dt] = surfaceEnergyBalanceInfiltration(T, wc, FORCING, GRID, PARA, SEB);
%------ soil module --------------------------------------------------
%calculate heat conduction
SEB = heatConduction(T, k_Kgrid, GRID, PARA, SEB);
%------ sum up heat fluxes --------------------------------------------
SEB.dE_dt = SEB.dE_dt_cond + SEB.dE_dt_SEB;
%------ determine optimal timestep ------------------------------------
% account for min and max timesteps specified, max. energy change per grid cell and the CFT stability criterion.
% energy change due to advection of heat through water fluxes is still excluded.
% timestep in [days]
timestep = min( [ max( [ min( [ 0.5 * nanmin( GRID.general.K_delta.^2 .* c_cTgrid ./ k_cTgrid ./ (GRID.soil.cT_domain + GRID.snow.cT_domain ) ) ./ (24.*3600), ...
PARA.technical.targetDeltaE .* nanmin( abs(GRID.general.K_delta ./ SEB.dE_dt ) ) ./ (24.*3600), ...
PARA.technical.maxTimestep ] ), ...
PARA.technical.minTimestep ] ), ...
TEMPORARY.outputTime-t ] );
if PARA.modules.lateral
timestep = min( [timestep, TEMPORARY.syncTime-t] );
end
% give a warning when timestep required by CFL criterion is below the minimum timestep specified
if timestep > 0.5 * min( GRID.general.K_delta.^2 .* c_cTgrid ./ k_cTgrid ./ (GRID.soil.cT_domain + GRID.snow.cT_domain) ) / (24.*3600)
warning( 'numerical stability not guaranteed' );
end
%------ update T array ------------------------------------------------
% account for vertical heat fluxes from ground heat flux and heat conduction
T = T + SEB.dE_dt./c_cTgrid./GRID.general.K_delta.*timestep.*24.*3600;
% set grid cells in air to air temperature
T(GRID.air.cT_domain)=FORCING.i.Tair;
%------- water body module --------------------------------------------
T = mixingWaterBody(T, GRID);
%------- snow cover module --------------------------------------------
[T, GRID, PARA, SEB, BALANCE] = CryoGridSnow(T, GRID, FORCING, SEB, PARA, c_cTgrid, timestep, BALANCE);
[GRID, T, BALANCE] = updateGRID_snow(T, GRID, PARA, BALANCE);
%------- infiltration module-------------------------------------------
if PARA.modules.infiltration
[wc, GRID, BALANCE] = CryoGridInfiltration(T, wc, dwc_dt, timestep, GRID, PARA, FORCING, BALANCE);
end
%------- excess ice module --------------------------------------------
if PARA.modules.xice && ~PARA.modules.infiltration
warning( 'energy and water balances are not correct for this combination of modules');
[GRID, PARA] = excessGroundIce(T, GRID, PARA);
wc = wc( end-sum(GRID.soil.cT_domain)+1 : end ); % assure wc has correct length
elseif PARA.modules.xice && PARA.modules.infiltration
[GRID, PARA, wc, meltwaterGroundIce, TEMPORARY] = excessGroundIceInfiltration(T, wc, GRID, PARA, TEMPORARY);
GRID = updateGRID_excessiceInfiltration(meltwaterGroundIce, GRID);
end
%------- update Lstar for next time step ------------------------------
SEB = L_star(FORCING, PARA, SEB);
%------- update auxiliary state variables
PARA.location.altitude = getAltitude( PARA, GRID );
PARA.location.surface_altitude = getSurfaceAltitude( PARA, GRID );
PARA.location.soil_altitude = getSoilAltitude( PARA, GRID );
[PARA.location.infiltration_altitude, PARA.location.bottomBucketSoilcTIndex] = getInfiltrationAltitude( PARA, GRID, T);
[PARA.location.water_table_altitude] = getWaterTableAltitudeFC(T, wc, GRID, PARA);
PARA.soil.infiltration_limit_altitude = PARA.location.soil_altitude - PARA.soil.infiltration_limit_depth;
%------- update threshold variables if no lateral exchange processes occur, otherwise updated at sync time
if ~PARA.modules.lateral
PARA.location.absolute_maxWater_altitude = PARA.location.altitude + PARA.soil.relative_maxWater;
if isempty( PARA.snow.relative_maxSnow )
PARA.location.absolute_maxSnow_altitude = [];
else
PARA.location.absolute_maxSnow_altitude = PARA.location.altitude + PARA.snow.relative_maxSnow;
end
end
%------- lateral exchange module --------------------------------------
% all functions called in this block should go into
% /modules/cryoGridLateral
% calling PARA.ensemble is only allowed here
if PARA.modules.lateral
if t==TEMPORARY.syncTime %communication between workers
fprintf('\n\t\t\tCryoGridLateral: sync - start (Worker %1.0f)\n', labindex);
% update auxiliary variables and common thresholds
labBarrier();
[PARA] = updateAuxiliaryVariablesAndCommonThresholds(T, wc, GRID, PARA);
% HEAT exchange module
if PARA.modules.exchange_heat
[ T, TEMPORARY, BALANCE ] = CryoGridLateralHeat( PARA, GRID, BALANCE, TEMPORARY, T, k_cTgrid, c_cTgrid );
end
% WATER exchange module
if PARA.modules.exchange_water
[wc, GRID, BALANCE] = CryoGridLateralWater( PARA, GRID, BALANCE, T, wc);
end
% SNOW exchange module
if PARA.modules.exchange_snow
PARA = CryoGridLateralSnow( PARA, GRID );
end
% lateral EROSION module // SEDIMENT exchange module
if PARA.modules.exchange_sediment
[wc, GRID, TEMPORARY] = CryoGridLateralErosion(PARA, GRID, wc, T, TEMPORARY);
end
% update auxiliary variables and common thresholds
labBarrier();
[PARA] = updateAuxiliaryVariablesAndCommonThresholds(T, wc, GRID, PARA) ;
% determine next sync time
TEMPORARY.syncTime=round((TEMPORARY.syncTime + PARA.technical.syncTimeStep)./PARA.technical.syncTimeStep).*PARA.technical.syncTimeStep;
fprintf('\t\t\tsync - done\n');
end
end
%------- water balance calculations -----------------------------------
% rainfall
BALANCE.water.dp_rain = BALANCE.water.dp_rain + FORCING.i.rainfall.*timestep; %sum up rainfall in [mm] per output interval
% snowfall
BALANCE.water.dp_snow = BALANCE.water.dp_snow + FORCING.i.snowfall.*timestep; %sum up snowfall in [mm] SWE per output interval
%------- next time step -----------------------------------------------
t=t+timestep;
%---------- sum up + OUTPUT -------------------------------------------
[TEMPORARY, OUT, BALANCE] = sum_up_output_store(t, T, wc, lwc_cTgrid(GRID.soil.cT_domain), timestep, TEMPORARY, BALANCE, PARA, GRID, SEB, OUT, FORCING, saveDir, run_number);
end
% save final state and output at t=endtime
DIAG = diagnose_output_yearly( OUT, PARA, GRID, FORCING );
iSaveDIAG( [ saveDir '/' run_number '/' run_number '_realization' num2str(labindex) '_diagnostics' datestr(t,'yyyy') '.mat' ], DIAG);
iSaveOUT( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_output' datestr(t,'yyyy') '.mat'], OUT)
iSaveState( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_finalState' datestr(t,'yyyy') '.mat'], T, wc, t, SEB, PARA, GRID)
iPlotAltitudes( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_altitudes' datestr(t,'yyyy') '.png'], OUT, PARA );
iPlotTemperature( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_temperature' datestr(t,'yyyy') '.png'], OUT, PARA, GRID);
iPlotWaterContent( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_waterContent' datestr(t,'yyyy') '.png'], OUT, PARA, GRID );
%catch
% fprintf('Catched exception. Saved final state under subscript CRASH.\n');
% save final state and output at crash
% iSaveOUT( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_output' datestr(t,'yyyy') 'CRASH.mat'], OUT)
% iSaveState( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_finalState' datestr(t,'yyyy') 'CRASH.mat'], T, wc, t, SEB, PARA, GRID)
%iPlotAltitudes( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_altitudes' datestr(t,'yyyy') 'CRASH.png'], OUT, PARA );
%iPlotTemperature( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_temperature' datestr(t,'yyyy') 'CRASH.png'], OUT, PARA, GRID);
%iPlotWaterContent( [ saveDir '/' run_number '/' run_number '_realization' num2str(index) '_waterContent' datestr(t,'yyyy') 'CRASH.png'], OUT, PARA, GRID );
%end
end
if number_of_realizations>1
delete(gcp('nocreate'))
end
fprintf('Done.\n');