qp_driver.py 163 KB
Newer Older
1
import sys, glob, os
2
import argparse
3
from ipdb import set_trace as mybreak
4

5
6
7
8
9
10
#class DS_attributes(object):
#  def __init__(self, prefix, output_freq='yearly'):
#    self.prefix = prefix
#    self.output_freq = output_freq
#  return

11
12
# --- initialization / default values for config file
run = ''
13
runname = ''
14
15
16
path_data = ''

# --- path to quickplots
17
path_quickplots = '../../all_qps/'
18
19

# --- set this to True if the simulation is still running
20
omit_last_file = False
21

22
23
24
25
# --- decide which set of plots to do
do_ocean_plots = True
do_atmosphere_plots = False
do_hamocc_plots = False
26

27
28
29
30
31
# --- grid information
gname = ''
lev   = ''
gname_atm = ''
lev_atm = ''
32

33
34
35
36
# --- path to interpolation files
path_grid        = ''
path_grid_atm    = ''
path_ckdtree     = ''
37
38
path_ckdtree_atm = 'auto'

39
40
41
42
# --- grid files and reference data
fpath_tgrid         = '' 
fpath_tgrid_atm     = '' 
fpath_ref_data_oce  = '' 
43
fpath_ref_data_ham  = ''
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
fpath_ref_data_atm  = '' 
fpath_ref_data_atm = '/mnt/lustre01/work/mh0033/m300602/icon/era/pyicon_prepare_era.nc'
fpath_fx            = '' 

# --- nc file prefixes
oce_def = ''
oce_moc = '_MOC'
oce_mon = '_oceanMonitor'
oce_ice = oce_def
oce_monthly = oce_def

atm_2d      = '' 
atm_3d      = '' 
atm_mon     = '' 

# --- time average information (can be overwritten by qp_driver call)
tave_ints = [['1950-02-01', '1952-01-01']]
# --- decide which data files to take for time series plots
tstep     = '????????????????'
# --- set this to 12 for yearly averages in timeseries plots, set to 0 for no averaging
64
ave_freq = 12
65

66
67
68
# --- xarray usage
xr_chunks = None
load_xarray_dset = False
69
70
save_data = False
path_nc = './'
71

72
73
# --- information for re-gridding
sec_name_30w   = '30W_300pts'
74
sec_name_170w   = '170W_300pts'
75
76
rgrid_name     = 'global_0.3'
rgrid_name_atm = 'global_1.0_era'
77

78
verbose = False
79
do_write_final_config = False
80

81
82
# --- list containing figures which should not be plotted
red_list = []
83
84
# --- list containing figures which should be plotted 
# (if empty all figures will be plotted)
85
86
green_list = []

87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
help_text = """
Driver for pyicon quickplots.

Usage notes:
------------

In ipython / Jupyter:
%run qp_driver.py /path/to/config_file.py

From command line:
python qp_driver.py /path/to/config_file.py

In batch mode (slurm):
python qp_driver.py --batch /path/to/config_file.py

102
For debugging:
103
104
105
%run qp_driver.py /path/to/config_file.py --debug --tave_int=1610-01-01,1620-01-01

Just creating the web page without producing any new figures:
106
107
%run qp_driver.py /path/to/config_file.py --no_plots --tave_int=1610-01-01,1620-01-01

108
109
110
111
112
113
114
Argument list:
--------------
"""

# --- read input arguments
parser = argparse.ArgumentParser(description=help_text, formatter_class=argparse.RawTextHelpFormatter)

115
116
117
118
# --- necessary arguments
parser.add_argument('fpath_config', metavar='fpath_config', type=str,
                    help='path to quickplot configure file')
# --- optional arguments
119
120
121
122
123
124
125
126
parser.add_argument('--batch', type=bool, default=False, 
                    help='suppress displaying figures, required for batch mode (same as --slurm)')
parser.add_argument('--slurm', default=False, 
                    action='store_true', #const=False,
                    help='suppress displaying figures, required for batch mode (same as --batch=True)')
parser.add_argument('--no_plots', default=False, 
                    action='store_true', #const=False,
                    help='do not make any plots')
127
128
129
parser.add_argument('--debug', default=False, 
                    action='store_true', #const=False,
                    help='only limitted number of plots are made (specified in qp_driver debugging section)')
130
parser.add_argument('--path_quickplots', metavar='path_quickplots', type=str, default='none',
131
                    help='path where the quickplots website and figures are storred')
132
133
parser.add_argument('--tave_int', metavar='tave_int', type=str, default='none',
                    help='specify time averaging interval e.g. --tave_int=1600-01-01, 1700-12-31')
134
135
parser.add_argument('--run', metavar='run', type=str, default='none',
                    help='name of simulation; if specified, it is used to alter path_data by replacing "run" from config file')
136
137
parser.add_argument('--path_data', metavar='path_data', type=str, default='none',
                    help='path of simulation data; if specified, it is used to replace path_data')
138
139
140
141
iopts = parser.parse_args()

print('--------------------------------------------------------------------------------')
print('Input arguments:')
142
print('--------------------------------------------------------------------------------')
143
144
145
146
147
148
149
150
151
152
153
154
for var in iopts.__dict__.keys():
  print(var, ' = ', getattr(iopts, var))
print('--------------------------------------------------------------------------------')

fpath_config = iopts.fpath_config

if iopts.batch or iopts.slurm:
  print('using batch mode')

# --- continue loading modules
import matplotlib
if iopts.batch or iopts.slurm:
155
  matplotlib.use('Agg')
156
157
158
159
160
161
import shutil
import datetime
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import json
162
#sys.path.append('/home/mpim/m300602/pyicon')
163
import pyicon as pyic
164
import pyicon.quickplots as pyicqp
165
166
167
import cartopy
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
Nils Brüggemann's avatar
Nils Brüggemann committed
168
from qp_cmaps import PyicCmaps
169

170
#cm_wbgyr = PyicCmaps().WhiteBlueGreenYellowRed
Nils Brüggemann's avatar
Nils Brüggemann committed
171

172
# ------ r2b4
173
174
#exec(open("../../config_qp/conf-icon_08-nib002.py").read())
#exec(open("../../config_qp/conf-mac-icon_08-nib0002.py").read())
175
176

# ------ r2b6
177
#exec(open("../../config_qp/conf-icon_08-nib003.py").read())
178
# does not work so far. diff. grid file
179
#exec(open("../../config_qp/conf-ocean_omip_long_r2b6_19224-YVF.py").read()) 
180
181

# ------ r2b8
182
#exec(open("../../config_qp/conf-ocean_era51h_r2b8_19074-AMK.py").read())
183
#exec(open("../../config_qp/conf-exp.ocean_ncep6h_r2b8_20096-FZA.py").read())
184
185

# ------ r2b6
Nils Brüggemann's avatar
Nils Brüggemann committed
186
#exec(open("../../config_qp/conf-icon_08-nib0004.py").read())
187
#exec(open("../../config_qp/conf-icon_08-nib0006.py").read())
188
#exec(open("../../config_qp/conf-jkr0042.py").read())
189
190
191
192
193
#exec(open("../../config_qp/conf-slo1284.py").read())

if not os.path.isfile(fpath_config):
  raise ValueError("::: Error: Config file %s does not exist! :::" % (fpath_config))
exec(open(fpath_config).read())
Nils Brüggemann's avatar
Nils Brüggemann committed
194

195
# --- overwrite variables if given by argument parsing
196
197
198
199
if iopts.path_data!='none':
  path_data = iopts.path_data
if not path_data.endswith('/'):
  path_data += '/'
200
if iopts.run!='none':
201
  #path_data = path_data.replace(run, iopts.run)
202
  run = iopts.run
203
204
if iopts.path_quickplots!='none':
  path_quickplots = iopts.path_quickplots
205
206
207
208
if not iopts.tave_int=='none':
  tave_int = iopts.tave_int.split(',')
  tave_ints = [tave_int]

209
path_quickplots = os.path.abspath(path_quickplots) + '/'
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

config_file = f"""
run = \'{run}\'
runname = \'{runname}\'
path_data = \'{path_data}\'

# --- path to quickplots
path_quickplots = \'{path_quickplots}\'

# --- set this to True if the simulation is still running
omit_last_file = {omit_last_file}

# --- decide which set of plots to do
do_ocean_plots      = {do_ocean_plots}
do_atmosphere_plots = {do_atmosphere_plots}
do_hamocc_plots     = {do_hamocc_plots}

# --- grid information
gname     = \'{gname}\'
lev       = \'{lev}\'
gname_atm = \'{gname_atm}\'
lev_atm   = \'{lev_atm}\'

# --- path to interpolation files
path_grid        = \'{path_grid}\'
path_grid_atm    = \'{path_grid_atm}\'
path_ckdtree     = \'{path_ckdtree}\'
path_ckdtree_atm = \'{path_ckdtree_atm}\'

# --- grid files and reference data
fpath_tgrid         = \'{fpath_tgrid}\'
fpath_tgrid_atm     = \'{fpath_tgrid_atm}\'
fpath_ref_data_oce  = \'{fpath_ref_data_oce}\'
243
fpath_ref_data_ham  = \'{fpath_ref_data_ham}\'
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
fpath_ref_data_atm  = \'{fpath_ref_data_atm}\'
fpath_fx            = \'{fpath_fx}\'

# --- nc file prefixes
oce_def     = \'{oce_def}\'
oce_moc     = \'{oce_moc}\'
oce_mon     = \'{oce_mon}\'
oce_ice     = \'{oce_ice}\'
oce_monthly = \'{oce_monthly}\'

atm_2d      = \'{atm_2d}\'
atm_3d      = \'{atm_3d}\'
atm_mon     = \'{atm_mon}\'

# --- time average information (can be overwritten by qp_driver call)
tave_ints = {tave_ints}
# --- decide which data files to take for time series plots
tstep     = \'{tstep}\'
# --- set this to 12 for yearly averages in timeseries plots, set to 0 for no averaging
ave_freq = {ave_freq}

# --- xarray usage
xr_chunks = {xr_chunks}
load_xarray_dset = {load_xarray_dset}
268
269
save_data = {save_data}
path_nc = \'{path_nc}\'
270
271
272

# --- information for re-gridding
sec_name_30w   = \'{sec_name_30w}\'
273
sec_name_170w   = \'{sec_name_170w}\'
274
275
276
277
278
279
280
281
282
283
284
285
rgrid_name     = \'{rgrid_name}\'
rgrid_name_atm = \'{rgrid_name_atm}\'

verbose = {verbose}
do_write_final_config = {do_write_final_config}

# --- list containing figures which should not be plotted
red_list = {red_list}
# --- list containing figures which should be plotted 
# (if empty all figures will be plotted)
green_list = {green_list}
"""
286

287
288
289
290
291
292
293
294
295
296
297
298
fpath_config_full = os.path.abspath(fpath_config[:-3]+'_full.py')
print('--------------------------------------------------------------------------------')
print(f'Configurations for pyicon quickplots:')
if do_write_final_config: 
  print(f'Written to:\n{fpath_config_full}')
print('--------------------------------------------------------------------------------')
print(config_file)
print('--------------------------------------------------------------------------------')
if do_write_final_config: 
  f = open(fpath_config_full, 'w')
  f.write(config_file)
  f.close()
299

300
301
302
303
304
305
# -------------------------------------------------------------------------------- 
# Settings
# -------------------------------------------------------------------------------- 
projection = 'PlateCarree'
#projection = 'none'

306
# --- structure of web page
307
fig_names = []
308
if do_ocean_plots:
309
310
  fig_names += ['sec:Overview']
  fig_names += ['tab_overview']
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
  fig_names += ['sec:Upper ocean']
  fig_names += ['ssh', 'ssh_variance', 'sst', 'sss', 'mld_mar', 'mld_sep'] 
  fig_names += ['sec:Ice']
  fig_names += ['ice_concentration_nh', 'ice_thickness_nh', 'snow_thickness_nh',] 
  fig_names += ['ice_concentration_sh', 'ice_thickness_sh', 'snow_thickness_sh',]
  fig_names += ['sec:Sections']
  fig_names += ['temp30w', 'salt30w', 'dens30w']
  fig_names += ['sec:Zonal averages']
  fig_names += ['temp_gzave', 'temp_azave', 'temp_ipzave']
  fig_names += ['salt_gzave', 'salt_azave', 'salt_ipzave']
  fig_names += ['sec:Transports']
  fig_names += ['bstr', 'passage_transports', 'tab_passage_transports']
  fig_names += ['arctic_budgets']
  fig_names += ['amoc', 'pmoc', 'gmoc']
  fig_names += ['ke_100m', 'ke_2000m']
  fig_names += ['heat_flux', 'freshwater_flux']
  fig_names += ['sec:Biases']
  fig_names += ['sst_bias', 'temp_bias_gzave', 'temp_bias_azave', 'temp_bias_ipzave']
  fig_names += ['sss_bias', 'salt_bias_gzave', 'salt_bias_azave', 'salt_bias_ipzave']
  fig_names += ['sec:Time series']
331
  fig_names += ['ts_amoc', 'ts_heat_content', 'ts_ssh', 'ts_sst', 'ts_sss', 'ts_hfl', 'ts_wfl', 'ts_ice_volume_nh', 'ts_ice_volume_sh', 'ts_ice_extent_nh', 'ts_ice_extent_sh',]
332
333
if do_atmosphere_plots:
  fig_names += ['ts_tas_gmean', 'ts_radtop_gmean']
334
  fig_names += ['ts_rsdt_gmean', 'ts_rsut_gmean', 'ts_rlut_gmean', 'ts_prec_gmean', 'ts_evap_gmean', 'ts_fwfoce_gmean']
335
336
337
338
339
340
341
342
343
344
345
346
  fig_names += ['sec:Surface fluxes']
  fig_names += ['atm_zonal_wind_stress', 'atm_meridional_wind_stress']
  fig_names += ['atm_curl_tau', 'atm_wek']
  fig_names += ['sec:Bias surface fluxes']
  fig_names += ['atm_tauu_bias', 'atm_tauv_bias']
  fig_names += ['sec:Atmosphere surface']
  fig_names += ['atm_2m_temp','atm_surface_temp','atm_sea_level_pressure',]
  fig_names += ['atm_column_water_vapour', 'atm_total_precipitation', 'atm_total_cloud_cover', 'atm_p_e', 'atm_10m_wind']
  fig_names += ['sec:Bias atmosphere surface']
  fig_names += ['atm_tas_bias']
  fig_names += ['atm_prw_bias']
  fig_names += ['atm_psl_bias']
347
348
  fig_names += ['sec:Above surface']
  fig_names += ['atm_geoh_500', 'atm_temp_850']
349
  fig_names += ['sec:Atmosphere zonal averages']
350
351
352
353
  fig_names += ['atm_temp_zave', 'atm_temp_zave_bias', 'atm_logv_temp_zave', 'atm_logv_temp_zave_bias']
  fig_names += ['atm_u_zave', 'atm_u_zave_bias', 'atm_logv_u_zave', 'atm_logv_u_zave_bias']
  fig_names += ['atm_v_zave', 'atm_v_zave_bias', 'atm_logv_v_zave', 'atm_logv_v_zave_bias']
  fig_names += ['atm_rel_hum_zave']
354
  fig_names += ['atm_cloud_cover_zave', 'atm_cloud_water_zave', 'atm_cloud_ice_zave', 'atm_cloud_water_ice_zave', 'atm_psi']
355
#fig_names += ['sec:TKE and IDEMIX']
356
#fig_names += ['tke30w', 'iwe30w', 'kv30w']
357
358
359
360
if do_hamocc_plots:
  fig_names += ['sec:Hamocc time series']
  fig_names += ['ts_global_npp', 'ts_global_nppcya', 'ts_global_zoograzing', 'ts_global_netco2flux']
  fig_names += ['ts_global_surface_alk', 'ts_global_surface_dic', 'ts_global_surface_sil', 'ts_global_surface_phos']
361
  fig_names += ['ts_global_surface_nitrate', 'ts_WC_denit', 'ts_sed_denit', 'ts_n2fix', 'ts_global_opal_prod', 'ts_global_caco3_prod']
362
363
364
  fig_names += ['ts_global_OMexp90', 'ts_global_calcexp90', 'ts_global_opalexp90']
  fig_names += ['sec:Hamocc surface maps']
  fig_names += ['srf_phyp', 'srf_zoop', 'srf_cya', 'srf_silicate', 'srf_nitrate', 'srf_phosphate']
365
366
  fig_names += ['srf_iron', 'srf_alk', 'srf_dic', 'srf_pH', 'srf_co2flux'] # srf_hion
  fig_names += ['srf_NP_ratio', 'age_100m', 'age_500m', 'age_2000m']
367
368
369
  fig_names += ['sec:Hamocc bias plots']
  fig_names += ['po4_srf_bias_ham', 'po4_1000_bias_ham', 'no3_srf_bias_ham', 'no3_1000_bias_ham']
  fig_names += ['si_srf_bias_ham', 'si_1000_bias_ham', 'o2_srf_bias_ham', 'o2_1000_bias_ham']
370
  fig_names += ['alk_srf_bias_ham', 'alk_1000_bias_ham', 'dic_srf_bias_ham', 'dic_1000_bias_ham']
371
  fig_names += ['sec:Hamocc sections']
372
373
  fig_names += ['po4_30w', 'po4_170w', 'no3_30w', 'no3_170w', 'si_30w', 'si_170w', 'talk_30w', 'talk_170w']
  fig_names += ['sec: Hamocc zonal averages']
374
  fig_names += ['dic_gzave', 'dic_azave', 'dic_ipzave', 'o2_gzave', 'o2_azave', 'o2_ipzave']
375
376
377
  fig_names += ['sec:Hamocc sediment']
  fig_names += ['sedflux_prcaca', 'sedflux_prorca', 'sedflux_produs', 'sedflux_silpro']

378

379
380
381
plist = fig_names
fig_names = []
for pitem in plist:
382
  if not pitem.startswith('sec:') and not pitem in red_list:
383
384
    fig_names += [pitem]

385
386
387
388
# --- use green list if it hast entries
if not len(green_list)==0:
  fig_names = green_list

389
# --- for debugging
390
if iopts.debug:
391
  print('XXXXXXXXXXXXXXXXX Debugging mode! XXXXXXXXXXXXXXX')
392
  fig_names = []
393
394
  #fig_names += ['temp30w', 'salt30w', 'dens30w']
  #fig_names += ['atm_psi']
395
  fig_names += ['ts_tas_gmean']
396
  #fig_names += ['sst']
397
  #fig_names += ['ts_amoc']
398
  #fig_names += ['ts_amoc', 'ts_ssh', 'ts_sst', 'ts_sss', 'ts_hfl', 'ts_wfl', 'ts_ice_volume_nh', 'ts_ice_volume_sh', 'ts_ice_extent_nh', 'ts_ice_extent_sh',]
399
400
401
402
403
404
405
406
  #fig_names += ['mld_mar', 'mld_sep']
  #fig_names = ['temp_bias_gzave']
  #fig_names = ['sss']
  #fig_names = ['ssh_variance']
  #fig_names += ['amoc']
  #fig_names += ['sst_bias', 'temp_bias_gzave', 'temp_bias_azave', 'temp_bias_ipzave']
  #fig_names += ['sss_bias', 'salt_bias_gzave', 'salt_bias_azave', 'salt_bias_ipzave']
  #fig_names += ['ice_concentration_nh', 'ice_thickness_nh', 'snow_thickness_nh',] 
407
  #fig_names += ['ice_concentration_nh']
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
  #fig_names += ['ice_concentration_sh', 'ice_thickness_sh', 'snow_thickness_sh',]
  #fig_names += ['bstr']
  #fig_names += ['ke_100m', 'ke_2000m']
  #fig_names += ['atm_sea_level_pressure']
  #fig_names += ['atm_2m_temp','atm_sea_level_pressure',]
  #fig_names += ['atm_zonal_wind_stress',]
  #fig_names += ['atm_curl_tau', 'atm_wek']
  #fig_names += ['atm_column_water_vapour', 'atm_total_precipitation', 'atm_total_cloud_cover']
  #fig_names += ['atm_tas_bias']
  #fig_names += ['atm_temp_zave', 'atm_u_zave', 'atm_v_zave', 'atm_rel_hum_zave']
  #fig_names += ['atm_logv_temp_zave']
  #fig_names += ['atm_cloud_cover_zave', 'atm_cloud_water_zave', 'atm_cloud_ice_zave', 'atm_cloud_water_ice_zave']
  #fig_names += ['atm_cloud_water_zave', 'atm_cloud_ice_zave', 'atm_cloud_water_ice_zave']
  #fig_names += ['vort']
  #fig_names += ['np_zonal_wind_stress']
  #fig_names += ['amoc', 'pmoc', 'gmoc']
  #fig_names += ['temp_gzave', 'temp_azave', 'temp_ipzave']
  #fig_names += ['salt30w', 'temp30w']
  #fig_names += ['ts_ice_extent_sh']
  #fig_names += ['atm_temp_zave_bias', 'atm_logv_temp_zave_bias', 'atm_logv_temp_zave', 'atm_temp_zave']
  #fig_names += ['atm_temp_zave', 'atm_temp_zave_bias', 'atm_logv_temp_zave', 'atm_logv_temp_zave_bias']
  #fig_names += ['atm_u_zave', 'atm_u_zave_bias', 'atm_logv_u_zave', 'atm_logv_u_zave_bias']
  #fig_names += ['atm_v_zave', 'atm_v_zave_bias', 'atm_logv_v_zave', 'atm_logv_v_zave_bias']
431
432
433
434
  #fig_names += ['atm_geoh_500', 'atm_temp_850']
  #fig_names += ['arctic_budgets']
  #fig_names += ['passage_transports', 'tab_passage_transports']
  #fig_names += ['tab_passage_transports']
435
  #fig_names += ['ts_amoc', 'tab_overview']
436
437
  #fig_names += ['ts_radtop_gmean', 'tab_overview']
  #fig_names += ['ts_rsdt_gmean', 'ts_rsut_gmean', 'ts_rlut_gmean', 'ts_prec_gmean', 'ts_evap_gmean', 'ts_fwfoce_gmean']
438
  #fig_names += ['tab_overview']
439

Nils Brüggemann's avatar
Nils Brüggemann committed
440
fig_names = np.array(fig_names)
441
442
443
444
445

# -------------------------------------------------------------------------------- 
# Function to save figures
# -------------------------------------------------------------------------------- 
close_figs = True
Nils Brüggemann's avatar
Nils Brüggemann committed
446
def save_fig(title, path_pics, fig_name, FigInf=dict()):
447
448
449
450
451
452
453
  FigInf['name']  = fig_name+'.png'
  FigInf['fpath'] = path_pics+FigInf['name']
  FigInf['title'] = title
  print('Saving {:<40} {:<40}'.format(FigInf['fpath'], FigInf['title']))
  plt.savefig(FigInf['fpath'], dpi=300)
  with open(path_pics+fig_name+'.json', 'w') as fj:
    json.dump(FigInf, fj, sort_keys=True, indent=4)
454
455
  if close_figs:
    plt.close('all')
456
  return
457
plt.close('all')
458

459
460
461
462
463
464
465
466
467
468
469
def save_tab(text, title, path_pics, fig_name, FigInf=dict()):
  FigInf['name']  = fig_name+'.html'
  FigInf['fpath'] = path_pics+FigInf['name']
  FigInf['title'] = title
  print('Saving {:<40} {:<40}'.format(FigInf['fpath'], FigInf['title']))
  with open(FigInf['fpath'], 'w') as f:
    f.write(text)
  with open(path_pics+fig_name+'.json', 'w') as fj:
    json.dump(FigInf, fj, sort_keys=True, indent=4)
  return

Nils Brüggemann's avatar
Nils Brüggemann committed
470
471
472
473
474
475
476
def indfind(array, vals):
  vals = np.array(vals)
  inds = np.zeros(vals.size, dtype=int)
  for nn, val in enumerate(vals):
    inds[nn] = np.argmin(np.abs(array-val))
  return inds

477
478
479
480
481
482
483
def load_era3d_var(fpath, var, isort):
  f = Dataset(fpath, 'r')
  data = f.variables[var][:].mean(axis=0)
  f.close()
  data = data[:,:,isort]
  return data

484
485
486
487
488
489
490
if load_xarray_dset and not xr_chunks is None:
  print('Setting up LocalCluster:')
  from dask.distributed import Client, LocalCluster
  cluster = LocalCluster(interface='ib0')
  client = Client(cluster)
  print(client)

491
492
493
# -------------------------------------------------------------------------------- 
# Load all necessary data sets (can take a while)
# -------------------------------------------------------------------------------- 
494
if do_ocean_plots and not iopts.no_plots:
Nils Brüggemann's avatar
Nils Brüggemann committed
495
  fname = '%s%s_%s.nc' % (run, oce_def, tstep)
496
497
498
499
500
  print('Dataset %s' % (fname))
  IcD = pyic.IconData(
                 fname        = fname,
                 path_data    = path_data,
                 path_grid    = path_grid,
501
                 path_ckdtree = path_ckdtree,
502
                 fpath_fx     = fpath_fx,
503
504
505
                 gname        = gname,
                 lev          = lev,
                 rgrid_name   = rgrid_name,
506
507
508
509
510
511
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = True,
                 load_triangular_grid   = True, # needed for bstr
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
512
513
                 load_xarray_dset       = load_xarray_dset,
                 xr_chunks              = xr_chunks,
514
                 verbose                = verbose,
515

516
517
                )
  fpath_ckdtree = IcD.rgrid_fpath_dict[rgrid_name]
518
  [k100, k500, k800, k1000, k2000, k3000] = indfind(IcD.depthc, [100., 500., 800., 1000., 2000., 3000.])
519
  
Nils Brüggemann's avatar
Nils Brüggemann committed
520
  fname_moc = '%s%s_%s.nc' % (run, oce_moc, tstep)
Nils Brüggemann's avatar
Nils Brüggemann committed
521
  print('Dataset %s' % (fname_moc))
522
523
524
525
  IcD_moc = pyic.IconData(
                 fname        = fname_moc,
                 path_data    = path_data,
                 path_grid    = path_grid,
526
                 path_ckdtree = path_ckdtree,
527
528
529
                 gname        = gname,
                 lev          = lev,
                 rgrid_name   = rgrid_name,
530
531
532
533
534
535
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
536
                 verbose                = verbose,
537
                )
538
539
  IcD_moc.depthc = IcD.depthc
  IcD_moc.depthi = IcD.depthi
540

Nils Brüggemann's avatar
Nils Brüggemann committed
541
542
543
544
545
546
  fname_monthly = '%s%s_%s.nc' % (run, oce_monthly, tstep)
  print('Dataset %s' % (fname_monthly))
  IcD_monthly = pyic.IconData(
                 fname        = fname_monthly,
                 path_data    = path_data,
                 path_grid    = path_grid,
547
                 path_ckdtree = path_ckdtree,
Nils Brüggemann's avatar
Nils Brüggemann committed
548
549
550
                 gname        = gname,
                 lev          = lev,
                 rgrid_name   = rgrid_name,
551
552
553
554
555
556
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
557
                 verbose                = verbose,
Nils Brüggemann's avatar
Nils Brüggemann committed
558
                )
559
  IcD_monthly.wet_c = IcD.wet_c
560

561
  fname_ice = '%s%s_%s.nc' % (run, oce_ice, tstep)
562
  print('Dataset %s' % (fname_ice))
563
564
565
566
567
568
569
570
571
572
573
574
575
  IcD_ice = pyic.IconData(
                 fname        = fname_ice,
                 path_data    = path_data,
                 path_grid    = path_grid,
                 path_ckdtree = path_ckdtree,
                 gname        = gname,
                 lev          = lev,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = False,
                 calc_coeff             = False,
576
                 verbose                = verbose,
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
                )

  fname_mon = '%s%s_%s.nc' % (run, oce_mon, tstep)
  print('Dataset %s' % (fname_mon))
  IcD_mon = pyic.IconData(
                 fname        = fname_mon,
                 path_data    = path_data,
                 path_grid    = path_grid,
                 path_ckdtree = path_ckdtree,
                 gname        = gname,
                 lev          = lev,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = False,
                 calc_coeff             = False,
594
                 verbose                = verbose,
595
596
                )

597
if do_atmosphere_plots and not iopts.no_plots:
598
599
  fname_atm_2d = '%s%s_%s.nc' % (run, atm_2d, tstep)
  print('Dataset %s' % (fname_atm_2d))
Nils Brüggemann's avatar
Nils Brüggemann committed
600
  IcD_atm2d = pyic.IconData(
601
                 fname        = fname_atm_2d,
Nils Brüggemann's avatar
Nils Brüggemann committed
602
603
                 path_data    = path_data,
                 path_grid    = path_grid_atm,
604
                 path_ckdtree = path_ckdtree_atm,
Nils Brüggemann's avatar
Nils Brüggemann committed
605
606
                 gname        = gname_atm,
                 lev          = lev_atm,
607
                 rgrid_name   = rgrid_name_atm,
608
609
610
611
612
613
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = True,
                 load_rectangular_grid  = True,
                 calc_coeff             = True,
614
                 verbose                = verbose,
615
616
                 time_mode    = 'float2date',
                 model_type   = 'atm',
Nils Brüggemann's avatar
Nils Brüggemann committed
617
                )
618

619
620
  fname_atm_3d = '%s%s_%s.nc' % (run, atm_3d, tstep)
  print('Dataset %s' % (fname_atm_3d))
Nils Brüggemann's avatar
Nils Brüggemann committed
621
  IcD_atm3d = pyic.IconData(
622
                 fname        = fname_atm_3d,
Nils Brüggemann's avatar
Nils Brüggemann committed
623
624
                 path_data    = path_data,
                 path_grid    = path_grid_atm,
625
                 path_ckdtree = path_ckdtree_atm,
Nils Brüggemann's avatar
Nils Brüggemann committed
626
627
                 gname        = gname_atm,
                 lev          = lev_atm,
628
                 rgrid_name   = rgrid_name_atm,
629
630
631
632
633
634
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
635
                 verbose                = verbose,
636
637
                 time_mode    = 'float2date',
                 model_type   = 'atm',
Nils Brüggemann's avatar
Nils Brüggemann committed
638
                )
639
  fpath_ckdtree_atm = IcD_atm3d.rgrid_fpath_dict[rgrid_name_atm]
640

641
642
  fname_atm_mon = '%s%s_%s.nc' % (run, atm_mon, tstep)
  print('Dataset %s' % (fname_atm_mon))
643
  IcD_atm_mon = pyic.IconData(
644
                 fname        = fname_atm_mon,
645
646
647
648
649
650
651
652
653
654
655
656
                 path_data    = path_data,
                 path_grid    = path_grid_atm,
                 path_ckdtree = path_ckdtree_atm,
                 gname        = gname_atm,
                 lev          = lev_atm,
                 rgrid_name   = rgrid_name_atm,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = False,
                 calc_coeff             = False,
657
                 verbose                = verbose,
658
659
660
                 #time_mode    = 'float2date',
                 model_type   = 'atm',
                )
Nils Brüggemann's avatar
Nils Brüggemann committed
661
  
662
663
  if do_ocean_plots==False:
    IcD_monthly = IcD_atm_mon
664
665
666
667
668
669
670
671
672
673

if do_hamocc_plots and not iopts.no_plots:
  if not do_ocean_plots:
    fname = '%s%s_%s.nc' % (run, oce_def, tstep)
    print('Dataset %s' % (fname))
    IcD = pyic.IconData(
                   fname        = fname,
                   path_data    = path_data,
                   path_grid    = path_grid,
                   path_ckdtree = path_ckdtree,
674
                   fpath_fx     = fpath_fx,
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
                   gname        = gname,
                   lev          = lev,
                   rgrid_name   = rgrid_name,
                   do_triangulation       = False,
                   omit_last_file         = omit_last_file,
                   load_vertical_grid     = True,
                   load_triangular_grid   = True, # needed for bstr
                   load_rectangular_grid  = True,
                   calc_coeff             = False,
                   verbose                = verbose,
                  )
    fpath_ckdtree = IcD.rgrid_fpath_dict[rgrid_name]
    [k100, k500, k800, k1000, k2000, k3000] = indfind(IcD.depthc, [100., 500., 800., 1000., 2000., 3000.])

    fname_monthly = '%s%s_%s.nc' % (run, oce_monthly, tstep)
    print('Dataset %s' % (fname_monthly))
    IcD_monthly = pyic.IconData(
                   fname        = fname_monthly,
                   path_data    = path_data,
                   path_grid    = path_grid,
                   path_ckdtree = path_ckdtree,
                   gname        = gname,
                   lev          = lev,
                   rgrid_name   = rgrid_name,
                   do_triangulation       = False,
                   omit_last_file         = omit_last_file,
                   load_vertical_grid     = False,
                   load_triangular_grid   = False,
                   load_rectangular_grid  = True,
                   calc_coeff             = False,
                   verbose                = verbose,
                  )
    IcD_monthly.wet_c = IcD.wet_c

  fname_ham_inv = '%s%s_%s.nc' % (run, ham_inv, tstep)
  print('Dataset %s' % (fname_ham_inv))
  IcD_ham_inv = pyic.IconData(
                 fname        = fname_ham_inv,
                 path_data    = path_data,
                 path_grid    = path_grid,
                 path_ckdtree = path_ckdtree,
716
                 fpath_fx     = fpath_fx,
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
                 gname        = gname,
                 lev          = lev,
                 rgrid_name   = rgrid_name,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = True,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
                 verbose                = verbose,
                )
  IcD_ham_inv.wet_c = IcD.wet_c

  fname_ham_2d = '%s%s_%s.nc' % (run, ham_2d, tstep)
  print('Dataset %s' % (fname_ham_2d))
  IcD_ham_2d = pyic.IconData(
                 fname        = fname_ham_2d,
                 path_data    = path_data,
                 path_grid    = path_grid,
                 path_ckdtree = path_ckdtree,
                 gname        = gname,
                 lev          = lev,
                 rgrid_name   = rgrid_name,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
                 verbose                = verbose,
                )
  IcD_ham_2d.wet_c = IcD.wet_c

750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
  fname_ham_sed = '%s%s_%s.nc' % (run, ham_sed, tstep)
  print('Dataset %s' % (fname_ham_sed))
  IcD_ham_sed = pyic.IconData(
                 fname        = fname_ham_sed,
                 path_data    = path_data,
                 path_grid    = path_grid,
                 path_ckdtree = path_ckdtree,
                 gname        = gname,
                 lev          = lev,
                 rgrid_name   = rgrid_name,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = True,
                 calc_coeff             = False,
                 verbose                = verbose,
                )
  IcD_ham_sed.wet_c = IcD.wet_c

770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
  fname_ham_mon = '%s%s_%s.nc' % (run, ham_mon, tstep)
  print('Dataset %s' % (fname_ham_mon))
  IcD_ham_mon = pyic.IconData(
                 fname        = fname_ham_mon,
                 path_data    = path_data,
                 path_grid    = path_grid,
                 path_ckdtree = path_ckdtree,
                 gname        = gname,
                 lev          = lev,
                 do_triangulation       = False,
                 omit_last_file         = omit_last_file,
                 load_vertical_grid     = False,
                 load_triangular_grid   = False,
                 load_rectangular_grid  = False,
                 calc_coeff             = False,
                 verbose                = verbose,
                )
Nils Brüggemann's avatar
Nils Brüggemann committed
787
print('Done reading datasets')
788

789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806

print(f'--------------------------------------------------------------------------------')
if do_ocean_plots and not iopts.no_plots:
  print(f'--- ocean:')
  print(f'fpath_fx          = {IcD.fpath_fx}')
  print(f'fpath_tgrid       = {IcD.fpath_tgrid}')
  print(f'path_data         = {path_data}')
  print(f'gname             = {gname}')
  print(f'lev               = {lev}')
  print(f'rgrid_name        = {rgrid_name}')
if do_atmosphere_plots and not iopts.no_plots:
  print(f'--- atmosphere:')
  #print(f'fpath_fx_atm      = {IcD_atm3d.fpath_fx}')
  print(f'fpath_tgrid_atm   = {IcD_atm3d.fpath_tgrid}')
  print(f'path_data         = {path_data}')
  print(f'gname_atm         = {gname_atm}')
  print(f'lev_atm           = {lev_atm}')
  print(f'rgrid_name_atm    = {rgrid_name_atm}')
807
808
809
810
811
812
813
814
if do_hamocc_plots and not iopts.no_plots:
  print(f'--- hamocc:')
  print(f'fpath_fx          = {IcD.fpath_fx}')
  print(f'fpath_tgrid       = {IcD.fpath_tgrid}')
  print(f'path_data         = {path_data}')
  print(f'gname             = {gname}')
  print(f'lev               = {lev}')
  print(f'rgrid_name        = {rgrid_name}')
815
816
print(f'--------------------------------------------------------------------------------')

817
# -------------------------------------------------------------------------------- 
Nils Brüggemann's avatar
Nils Brüggemann committed
818
# timing
819
# -------------------------------------------------------------------------------- 
820

Nils Brüggemann's avatar
Nils Brüggemann committed
821
for tave_int in tave_ints:
822
823
  t1 = tave_int[0].replace(' ', '')
  t2 = tave_int[1].replace(' ', '')
824

Nils Brüggemann's avatar
Nils Brüggemann committed
825
826
827
828
829
830
831
832
833
834
835
836
  if isinstance(t1,str) and t1=='auto':
    t1 = IcD.times[0]
  else:
    t1 = np.datetime64(t1)
  if isinstance(t2,str) and t2=='auto':
    t2 = IcD.times[-1]
  else:
    t2 = np.datetime64(t2)

  # -------------------------------------------------------------------------------- 
  # making new directories and copying files
  # -------------------------------------------------------------------------------- 
837
838
839
  # --- path to module
  path_qp_driver = os.path.dirname(pyicqp.__file__)+'/'

840
  # --- make directory for this current time average
Nils Brüggemann's avatar
Nils Brüggemann committed
841
842
843
844
  if runname=='':
    dirname = f'qp-{run}'
  else:
    dirname = f'qp-{runname}-{run}'
845
  path_qp_sim = f'{path_quickplots}/{dirname}/'
846
847
848
849
850
851
852
  path_qp_tav = f'{path_quickplots}/{dirname}/tave_{t1}-{t2}/'
  if not os.path.exists(path_qp_tav):
    os.makedirs(path_qp_tav)
  
  # --- make directory for pictures
  rpath_pics = 'pics/'
  path_pics = path_qp_tav+rpath_pics
Nils Brüggemann's avatar
Nils Brüggemann committed
853
854
855
  if not os.path.exists(path_pics):
    os.makedirs(path_pics)
  
856
  # --- copy css style file
857
858
859
860
861
862
863
864
865
  shutil.copyfile(path_qp_driver+'qp_css.css', path_qp_tav+'qp_css.css')
  
  # --- backup qp_driver (this script)
  fname_this_script = __file__.split('/')[-1]
  shutil.copyfile(path_qp_driver+fname_this_script, path_qp_tav+'bcp_'+fname_this_script)

  # --- backup config file
  fname_config = fpath_config.split('/')[-1]
  shutil.copyfile(fpath_config, path_qp_tav+'bcp_'+fname_config)
Nils Brüggemann's avatar
Nils Brüggemann committed
866
  
867
  # --- initialize web page
Nils Brüggemann's avatar
Nils Brüggemann committed
868
  if runname=='':
869
    title_str = run
Nils Brüggemann's avatar
Nils Brüggemann committed
870
  else:
871
    title_str = '%s | %s' % (runname, run)
Nils Brüggemann's avatar
Nils Brüggemann committed
872
  qp = pyicqp.QuickPlotWebsite(
873
    title=title_str, 
Nils Brüggemann's avatar
Nils Brüggemann committed
874
875
876
877
878
    author=os.environ.get('USER'),
    date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S"),
    path_data=path_data,
    info=f'time average from {t1} to {t2}',
    fpath_css='./qp_css.css',
879
    fpath_html=path_qp_tav+'qp_index.html'
Nils Brüggemann's avatar
Nils Brüggemann committed
880
    )
881
882
883
884
885
    
  if not iopts.no_plots:
    # -------------------------------------------------------------------------------- 
    # specify time averaging indices
    # -------------------------------------------------------------------------------- 
886
887
    mask_int = (IcD_monthly.times>=t1) & (IcD_monthly.times<=t2)
    months = IcD_monthly.times.astype('datetime64[M]').astype(int) % 12 + 1
888
    it_ave_months = np.where( mask_int )[0]
889
890
891
892
    # Note: In ICON the date of an average is set to the end of the averaging interval
    #       Thus, the 4th month corresponds to March and the 10th to September
    it_ave_mar = np.where( mask_int & (months==4)  )[0]
    it_ave_sep = np.where( mask_int & (months==10) )[0]
893
894
895
896
    if it_ave_mar.size==0:
      it_ave_mar = it_ave_months
    if it_ave_sep.size==0:
      it_ave_sep = it_ave_months
897
898
    #it_ave_mar = it_ave[2::12] # this only workes if tave_int start with Feb. E.g.: 1610-02-01,1620-01-01
    #it_ave_sep = it_ave[8::12] # this only workes if tave_int start with Feb. E.g.: 1610-02-01,1620-01-01
899
900
901
902
903
    it_ave_years = (IcD.times>=t1) & (IcD.times<=t2)
    print('tpoints for year average from yearly data:  ', IcD.times[it_ave_years])
    print('tpoints for year average from monthly data: ', IcD_monthly.times[it_ave_months])
    print('tpoints for Mar average:                    ', IcD_monthly.times[it_ave_mar])
    print('tpoints for Sep avarage:                    ', IcD_monthly.times[it_ave_sep])
904

905
906
907
    if mask_int.sum()==0:
      raise ValueError(f'::: Error: Cannot find any data in {path_data} for time period from {t1} unitl {t2}! :::')

908
909
910
911
912
913
914
    # ================================================================================ 
    # start with plotting
    # ================================================================================ 

    # -------------------------------------------------------------------------------- 
    # upper ocean
    # -------------------------------------------------------------------------------- 
Nils Brüggemann's avatar
Nils Brüggemann committed
915
    fname = '%s%s_%s.nc' % (run, oce_def, tstep)
916
917
918
919
920
921
    Ddict_global = dict(
      xlim=[-180.,180.], ylim=[-90.,90.],
      rgrid_name=rgrid_name,
      path_ckdtree=path_ckdtree,
      projection=projection,
                )
Nils Brüggemann's avatar
Nils Brüggemann committed
922
    
923
924
925
926
927
928
929
930
931
    # ---
    fig_name = 'mld_mar'
    if fig_name in fig_names:
      FigInf = pyicqp.qp_hplot(fpath=path_data+fname, var='mld', depth=0,
                               it_ave=it_ave_mar,
                               title='mixed layer depth March [m]',
                               #clim=[0,5000.], cincr=250., cmap=PyicCmaps().WhiteBlueGreenYellowRed,
                               clim=[0,5000.], cincr=250., cmap='RdYlBu_r',
                               do_mask=True,
932
933
934
                               IcD=IcD_monthly,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
935
      save_fig('Mixed layer depth March', path_pics, fig_name, FigInf)
Nils Brüggemann's avatar
Nils Brüggemann committed
936
    
937
938
939
940
941
942
943
944
    # ---
    fig_name = 'mld_sep'
    if fig_name in fig_names:
      FigInf = pyicqp.qp_hplot(fpath=path_data+fname, var='mld', depth=0,
                               it_ave=it_ave_sep,
                               title='mixed layer depth September [m]',
                               clim=[0,5000.], cincr=250., cmap='RdYlBu_r',
                               do_mask=True,
945
946
947
                               IcD=IcD_monthly,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
948
949
950
951
952
953
954
955
      save_fig('Mixed layer depth September', path_pics, fig_name, FigInf)

    # ---
    fig_name = 'sst'
    if fig_name in fig_names:
      FigInf = pyicqp.qp_hplot(fpath=path_data+fname, var='to', depth=0, it=0,
                               t1=t1, t2=t2,
                               clim=[-2.,30.], cincr=2.0, cmap='cmo.thermal',
956
957
958
                               IcD=IcD,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
959
960
961
962
963
964
965
966
967
      save_fig('SST', path_pics, fig_name, FigInf)
    
    # ---
    fig_name = 'sss'
    if fig_name in fig_names:
      FigInf = pyicqp.qp_hplot(fpath=path_data+fname, var='so', depth=0, it=0,
                               t1=t1, t2=t2,
                               #clim=[32.,37], cincr=0.25, cmap='cmo.haline',
                               clim=[25.,40.], clevs=[25, 28, 30, 32, 32.5, 33, 33.5, 34, 34.5, 35, 35.5, 36, 37, 38, 40], cmap='cmo.haline',
968
969
970
                               IcD=IcD,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
971
972
973
974
975
976
977
978
      save_fig('SSS', path_pics, fig_name, FigInf)
    
    # ---
    fig_name = 'ssh'
    if fig_name in fig_names:
      FigInf = pyicqp.qp_hplot(fpath=path_data+fname, var='zos', depth=0, it=0,
                               t1=t1, t2=t2,
                               clim=2, cincr=0.2, cmap='RdBu_r',
979
980
981
                               IcD=IcD,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
      save_fig('SSH', path_pics, fig_name, FigInf)

    # ---
    fig_name = 'ssh_variance'
    if fig_name in fig_names:
      zos, it_ave          = pyic.time_average(IcD, 'zos',        t1=t1, t2=t2)
      zos_square, it_ave   = pyic.time_average(IcD, 'zos_square', t1=t1, t2=t2)
      zos_var = np.sqrt(zos_square-zos**2)
      IaV = pyic.IconVariable('zos_var', 'm', 'ssh variance')
      IaV.data = zos_var
      IaV.interp_to_rectgrid(fpath_ckdtree)
      pyic.hplot_base(IcD, IaV, clim=[-2,0], cincr=0.2, cmap='RdYlBu_r',
                      projection=projection, xlim=[-180.,180.], ylim=[-90.,90.],
                      logplot=True,
                      title='log$_{10}$(ssh variance) [m]', do_write_data_range=True,
997
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
                      )
      save_fig('SSH variance', path_pics, fig_name)
    
    # ---
    fig_name = 'ice_concentration_nh'
    if fig_name in fig_names:

      conc_mar, it_ave = pyic.time_average(IcD_monthly, 'conc', it_ave=it_ave_mar, iz='all')
      hi_mar, it_ave   = pyic.time_average(IcD_monthly, 'hi', it_ave=it_ave_mar, iz='all')
      hs_mar, it_ave   = pyic.time_average(IcD_monthly, 'hs', it_ave=it_ave_mar, iz='all')
      hiconc_mar = (hi_mar*conc_mar)[0,:]
      hsconc_mar = (hs_mar*conc_mar)[0,:]

      conc_sep, it_ave = pyic.time_average(IcD_monthly, 'conc', it_ave=it_ave_sep, iz='all')
      hi_sep, it_ave   = pyic.time_average(IcD_monthly, 'hi', it_ave=it_ave_sep, iz='all')
      hs_sep, it_ave   = pyic.time_average(IcD_monthly, 'hs', it_ave=it_ave_sep, iz='all')
      hiconc_sep = (hi_sep*conc_sep)[0,:]
      hsconc_sep = (hs_sep*conc_sep)[0,:]

    # ---
    fig_name = 'ice_concentration_nh'
    if fig_name in fig_names:
      hca, hcb = pyic.arrange_axes(2,1, plot_cb=True, asp=1., fig_size_fac=2.,
                   sharex=True, sharey=True, xlabel="", ylabel="",
                   projection=ccrs.NorthPolarStereo(),
                                  )
      ii=-1

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('conc_mar', '', 'sea ice concentration March')
      IaV.data = conc_mar[0,:]
      IaV.interp_to_rectgrid(fpath_ckdtree)
1030
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
1031
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
1032
1033
1034
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1035
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1036
1037
1038
1039
1040
1041
                      )

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('conc_sep', 'm', 'sea ice concentration September')
      IaV.data = conc_sep[0,:]
      IaV.interp_to_rectgrid(fpath_ckdtree)
1042
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
1043
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
1044
1045
1046
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1047
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1048
1049
1050
1051
1052
1053
1054
1055
1056
                      )

      for ax in hca:
        ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())
        ax.gridlines()
        ax.add_feature(cartopy.feature.LAND)
        ax.coastlines()

      FigInf = dict(long_name=IaV.long_name)
1057
      save_fig('Sea ice concentration NH', path_pics, fig_name, FigInf)
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071

    # ---
    fig_name = 'ice_thickness_nh'
    if fig_name in fig_names:
      hca, hcb = pyic.arrange_axes(2,1, plot_cb=True, asp=1., fig_size_fac=2.,
                   sharex=True, sharey=True, xlabel="", ylabel="",
                   projection=ccrs.NorthPolarStereo(),
                                  )
      ii=-1

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('hiconc_mar', 'm', 'ice equiv. thickness March')
      IaV.data = hiconc_mar
      IaV.interp_to_rectgrid(fpath_ckdtree)
1072
1073
      pyic.hplot_base(IcD_monthly, IaV, 
                      clim=[0,6], clevs=[0, 0.01, 0.1, 0.2, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6], cmap='RdYlBu_r',
1074
1075
1076
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1077
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1078
1079
1080
1081
1082
1083
                      )

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('hiconc_sep', 'm', 'ice equiv. thickness September')
      IaV.data = hiconc_sep
      IaV.interp_to_rectgrid(fpath_ckdtree)
1084
1085
      pyic.hplot_base(IcD_monthly, IaV, 
                      clim=[0,6], clevs=[0, 0.01, 0.1, 0.2, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6], cmap='RdYlBu_r',
1086
1087
1088
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1089
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
                      )

      for ax in hca:
        ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())
        ax.gridlines()
        ax.add_feature(cartopy.feature.LAND)
        ax.coastlines()

      FigInf = dict(long_name=IaV.long_name)
      save_fig('Sea ice equiv. thickness NH', path_pics, fig_name, FigInf)

    # ---
    fig_name = 'snow_thickness_nh'
    if fig_name in fig_names:
      hca, hcb = pyic.arrange_axes(2,1, plot_cb=True, asp=1., fig_size_fac=2.,
                   sharex=True, sharey=True, xlabel="", ylabel="",
                   projection=ccrs.NorthPolarStereo(),
                                  )
      ii=-1

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('hsconc_mar', 'm', 'snow equiv. thickness March')
      IaV.data = hsconc_mar
      IaV.interp_to_rectgrid(fpath_ckdtree)
1114
      pyic.hplot_base(IcD_monthly, IaV, clim=[0,1], cincr=0.05, cmap='RdYlBu_r',
1115
1116
1117
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1118
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1119
1120
1121
1122
1123
1124
                      )

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('hsconc_sep', 'm', 'snow equiv. thickness September')
      IaV.data = hsconc_sep
      IaV.interp_to_rectgrid(fpath_ckdtree)
1125
      pyic.hplot_base(IcD_monthly, IaV, clim=[0,1], cincr=0.05, cmap='RdYlBu_r',
1126
1127
1128
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1129
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
                      )

      for ax in hca:
        ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())
        ax.gridlines()
        ax.add_feature(cartopy.feature.LAND)
        ax.coastlines()

      FigInf = dict(long_name=IaV.long_name)
      save_fig('Snow equiv. thickness NH', path_pics, fig_name, FigInf)

    # ---
    fig_name = 'ice_concentration_sh'
    if fig_name in fig_names:
      hca, hcb = pyic.arrange_axes(2,1, plot_cb=True, asp=1., fig_size_fac=2.,
                   sharex=True, sharey=True, xlabel="", ylabel="",
                   projection=ccrs.SouthPolarStereo(),
                                  )
      ii=-1

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('conc_mar', '', 'sea ice concentration March')
      IaV.data = conc_mar[0,:]
      IaV.interp_to_rectgrid(fpath_ckdtree)
1154
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
1155
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
1156
1157
1158
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[-90., -50.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1159
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1160
1161
1162
1163
1164
1165
                      )

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('conc_sep', 'm', 'sea ice concentration September')
      IaV.data = conc_sep[0,:]
      IaV.interp_to_rectgrid(fpath_ckdtree)
1166
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
1167
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
1168
1169
1170
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[-90., -50.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1171
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
1172
1173
1174
1175
1176
1177
1178
1179
1180
                      )

      for ax in hca:
        ax.set_extent([-180, 180, -90., -50.], ccrs.PlateCarree())
        ax.gridlines()
        ax.add_feature(cartopy.feature.LAND)
        ax.coastlines()

      FigInf = dict(long_name=IaV.long_name)
1181
      save_fig('Sea ice concentration SH', path_pics, fig_name, FigInf)
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195

    # ---
    fig_name = 'ice_thickness_sh'
    if fig_name in fig_names:
      hca, hcb = pyic.arrange_axes(2,1, plot_cb=True, asp=1., fig_size_fac=2.,
                   sharex=True, sharey=True, xlabel="", ylabel="",
                   projection=ccrs.SouthPolarStereo(),
                                  )
      ii=-1

      ii+=1; ax=hca[ii]; cax=hcb[ii]
      IaV = pyic.IconVariable('hiconc_mar', 'm', 'ice equiv. thickness March')
      IaV.data = hiconc_mar
      IaV.interp_to_rectgrid(fpath_ckdtree)
1196
      pyic.hplot_base(IcD_monthly, IaV, clim=[0,6], clevs=[0.1, 0.2, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6], cmap='RdYlBu_r',
1197
1198
1199
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[-90., -50.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
1200
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',