qp_driver.py

import sys, glob, os
import argparse
from ipdb import set_trace as mybreak

#class DS_attributes(object):
#  def __init__(self, prefix, output_freq='yearly'):
#    self.prefix = prefix
#    self.output_freq = output_freq
#  return

# --- initialization / default values for config file
run = ''
runname = ''
path_data = ''

# --- path to quickplots
path_quickplots = '../../all_qps/'

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

# --- decide which set of plots to do
do_ocean_plots = True
do_atmosphere_plots = False
do_hamocc_plots = False

# --- grid information
gname = ''
lev   = ''
gname_atm = ''
lev_atm = ''

# --- path to interpolation files
path_grid        = ''
path_grid_atm    = ''
path_ckdtree     = ''
path_ckdtree_atm = 'auto'

# --- grid files and reference data
fpath_tgrid         = '' 
fpath_tgrid_atm     = '' 
fpath_ref_data_oce  = '' 
fpath_ref_data_ham  = ''
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
ave_freq = 12

# --- xarray usage
xr_chunks = None
load_xarray_dset = False
save_data = False
path_nc = './'

# --- information for re-gridding
sec_name_30w   = '30W_300pts'
sec_name_170w   = '170W_300pts'
rgrid_name     = 'global_0.3'
rgrid_name_atm = 'global_1.0_era'

verbose = False
do_write_final_config = False

# --- list containing figures which should not be plotted
red_list = []
# --- list containing figures which should be plotted 
# (if empty all figures will be plotted)
green_list = []

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

For debugging:
%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:
%run qp_driver.py /path/to/config_file.py --no_plots --tave_int=1610-01-01,1620-01-01

Argument list:
--------------
"""

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

# --- necessary arguments
parser.add_argument('fpath_config', metavar='fpath_config', type=str,
                    help='path to quickplot configure file')
# --- optional arguments
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')
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)')
parser.add_argument('--path_quickplots', metavar='path_quickplots', type=str, default='none',
                    help='path where the quickplots website and figures are storred')
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')
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')
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')
iopts = parser.parse_args()

print('--------------------------------------------------------------------------------')
print('Input arguments:')
print('--------------------------------------------------------------------------------')
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:
  matplotlib.use('Agg')
import shutil
import datetime
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import json
#sys.path.append('/home/mpim/m300602/pyicon')
import pyicon as pyic
import pyicon.quickplots as pyicqp
import cartopy
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from qp_cmaps import PyicCmaps

#cm_wbgyr = PyicCmaps().WhiteBlueGreenYellowRed

# ------ r2b4
#exec(open("../../config_qp/conf-icon_08-nib002.py").read())
#exec(open("../../config_qp/conf-mac-icon_08-nib0002.py").read())

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

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

# ------ r2b6
#exec(open("../../config_qp/conf-icon_08-nib0004.py").read())
#exec(open("../../config_qp/conf-icon_08-nib0006.py").read())
#exec(open("../../config_qp/conf-jkr0042.py").read())
#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())

# --- overwrite variables if given by argument parsing
if iopts.path_data!='none':
  path_data = iopts.path_data
if not path_data.endswith('/'):
  path_data += '/'
if iopts.run!='none':
  #path_data = path_data.replace(run, iopts.run)
  run = iopts.run
if iopts.path_quickplots!='none':
  path_quickplots = iopts.path_quickplots
if not iopts.tave_int=='none':
  tave_int = iopts.tave_int.split(',')
  tave_ints = [tave_int]

path_quickplots = os.path.abspath(path_quickplots) + '/'

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}\'
fpath_ref_data_ham  = \'{fpath_ref_data_ham}\'
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}
save_data = {save_data}
path_nc = \'{path_nc}\'

# --- information for re-gridding
sec_name_30w   = \'{sec_name_30w}\'
sec_name_170w   = \'{sec_name_170w}\'
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}
"""

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()

# -------------------------------------------------------------------------------- 
# Settings
# -------------------------------------------------------------------------------- 
projection = 'PlateCarree'
#projection = 'none'

# --- structure of web page
fig_names = []
if do_ocean_plots:
  fig_names += ['sec:Overview']
  fig_names += ['tab_overview']
  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']
  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',]
if do_atmosphere_plots:
  fig_names += ['ts_tas_gmean', 'ts_radtop_gmean']
  fig_names += ['ts_rsdt_gmean', 'ts_rsut_gmean', 'ts_rlut_gmean', 'ts_prec_gmean', 'ts_evap_gmean', 'ts_fwfoce_gmean']
  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']
  fig_names += ['sec:Above surface']
  fig_names += ['atm_geoh_500', 'atm_temp_850']
  fig_names += ['sec:Atmosphere zonal averages']
  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']
  fig_names += ['atm_cloud_cover_zave', 'atm_cloud_water_zave', 'atm_cloud_ice_zave', 'atm_cloud_water_ice_zave', 'atm_psi']
#fig_names += ['sec:TKE and IDEMIX']
#fig_names += ['tke30w', 'iwe30w', 'kv30w']
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']
  fig_names += ['ts_global_surface_nitrate', 'ts_WC_denit', 'ts_sed_denit', 'ts_n2fix', 'ts_global_opal_prod', 'ts_global_caco3_prod']
  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']
  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']
  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']
  fig_names += ['alk_srf_bias_ham', 'alk_1000_bias_ham', 'dic_srf_bias_ham', 'dic_1000_bias_ham']
  fig_names += ['sec:Hamocc sections']
  fig_names += ['po4_30w', 'po4_170w', 'no3_30w', 'no3_170w', 'si_30w', 'si_170w', 'talk_30w', 'talk_170w']
  fig_names += ['sec: Hamocc zonal averages']
  fig_names += ['dic_gzave', 'dic_azave', 'dic_ipzave', 'o2_gzave', 'o2_azave', 'o2_ipzave']
  fig_names += ['sec:Hamocc sediment']
  fig_names += ['sedflux_prcaca', 'sedflux_prorca', 'sedflux_produs', 'sedflux_silpro']


plist = fig_names
fig_names = []
for pitem in plist:
  if not pitem.startswith('sec:') and not pitem in red_list:
    fig_names += [pitem]

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

# --- for debugging
if iopts.debug:
  print('XXXXXXXXXXXXXXXXX Debugging mode! XXXXXXXXXXXXXXX')
  fig_names = []
  #fig_names += ['temp30w', 'salt30w', 'dens30w']
  #fig_names += ['atm_psi']
  fig_names += ['ts_tas_gmean']
  #fig_names += ['sst']
  #fig_names += ['ts_amoc']
  #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',]
  #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',] 
  #fig_names += ['ice_concentration_nh']
  #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']
  #fig_names += ['atm_geoh_500', 'atm_temp_850']
  #fig_names += ['arctic_budgets']
  #fig_names += ['passage_transports', 'tab_passage_transports']
  #fig_names += ['tab_passage_transports']
  #fig_names += ['ts_amoc', 'tab_overview']
  #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']
  #fig_names += ['tab_overview']

fig_names = np.array(fig_names)

# -------------------------------------------------------------------------------- 
# Function to save figures
# -------------------------------------------------------------------------------- 
close_figs = True
def save_fig(title, path_pics, fig_name, FigInf=dict()):
  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)
  if close_figs:
    plt.close('all')
  return
plt.close('all')

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

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

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

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)

# -------------------------------------------------------------------------------- 
# Load all necessary data sets (can take a while)
# -------------------------------------------------------------------------------- 
if do_ocean_plots and not iopts.no_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,
                 fpath_fx     = fpath_fx,
                 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,
                 load_xarray_dset       = load_xarray_dset,
                 xr_chunks              = xr_chunks,
                 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_moc = '%s%s_%s.nc' % (run, oce_moc, tstep)
  print('Dataset %s' % (fname_moc))
  IcD_moc = pyic.IconData(
                 fname        = fname_moc,
                 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_moc.depthc = IcD.depthc
  IcD_moc.depthi = IcD.depthi

  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_ice = '%s%s_%s.nc' % (run, oce_ice, tstep)
  print('Dataset %s' % (fname_ice))
  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,
                 verbose                = verbose,
                )

  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,
                 verbose                = verbose,
                )

if do_atmosphere_plots and not iopts.no_plots:
  fname_atm_2d = '%s%s_%s.nc' % (run, atm_2d, tstep)
  print('Dataset %s' % (fname_atm_2d))
  IcD_atm2d = pyic.IconData(
                 fname        = fname_atm_2d,
                 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   = True,
                 load_rectangular_grid  = True,
                 calc_coeff             = True,
                 verbose                = verbose,
                 time_mode    = 'float2date',
                 model_type   = 'atm',
                )

  fname_atm_3d = '%s%s_%s.nc' % (run, atm_3d, tstep)
  print('Dataset %s' % (fname_atm_3d))
  IcD_atm3d = pyic.IconData(
                 fname        = fname_atm_3d,
                 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  = True,
                 calc_coeff             = False,
                 verbose                = verbose,
                 time_mode    = 'float2date',
                 model_type   = 'atm',
                )
  fpath_ckdtree_atm = IcD_atm3d.rgrid_fpath_dict[rgrid_name_atm]

  fname_atm_mon = '%s%s_%s.nc' % (run, atm_mon, tstep)
  print('Dataset %s' % (fname_atm_mon))
  IcD_atm_mon = pyic.IconData(
                 fname        = fname_atm_mon,
                 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,
                 verbose                = verbose,
                 #time_mode    = 'float2date',
                 model_type   = 'atm',
                )
  
  if do_ocean_plots==False:
    IcD_monthly = IcD_atm_mon

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,
                   fpath_fx     = fpath_fx,
                   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,
                 fpath_fx     = fpath_fx,
                 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

  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

  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,
                )
print('Done reading datasets')


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}')
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}')
print(f'--------------------------------------------------------------------------------')

# -------------------------------------------------------------------------------- 
# timing
# -------------------------------------------------------------------------------- 

for tave_int in tave_ints:
  t1 = tave_int[0].replace(' ', '')
  t2 = tave_int[1].replace(' ', '')

  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
  # -------------------------------------------------------------------------------- 
  # --- path to module
  path_qp_driver = os.path.dirname(pyicqp.__file__)+'/'

  # --- make directory for this current time average
  if runname=='':
    dirname = f'qp-{run}'
  else:
    dirname = f'qp-{runname}-{run}'
  path_qp_sim = f'{path_quickplots}/{dirname}/'
  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
  if not os.path.exists(path_pics):
    os.makedirs(path_pics)
  
  # --- copy css style file
  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)
  
  # --- initialize web page
  if runname=='':
    title_str = run
  else:
    title_str = '%s | %s' % (runname, run)
  qp = pyicqp.QuickPlotWebsite(
    title=title_str, 
    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',
    fpath_html=path_qp_tav+'qp_index.html'
    )
    
  if not iopts.no_plots:
    # -------------------------------------------------------------------------------- 
    # specify time averaging indices
    # -------------------------------------------------------------------------------- 
    mask_int = (IcD_monthly.times>=t1) & (IcD_monthly.times<=t2)
    months = IcD_monthly.times.astype('datetime64[M]').astype(int) % 12 + 1
    it_ave_months = np.where( mask_int )[0]
    # 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]
    if it_ave_mar.size==0:
      it_ave_mar = it_ave_months
    if it_ave_sep.size==0:
      it_ave_sep = it_ave_months
    #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
    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])

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

    # ================================================================================ 
    # start with plotting
    # ================================================================================ 

    # -------------------------------------------------------------------------------- 
    # upper ocean
    # -------------------------------------------------------------------------------- 
    fname = '%s%s_%s.nc' % (run, oce_def, tstep)
    Ddict_global = dict(
      xlim=[-180.,180.], ylim=[-90.,90.],
      rgrid_name=rgrid_name,
      path_ckdtree=path_ckdtree,
      projection=projection,
                )
    
    # ---
    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,
                               IcD=IcD_monthly,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
      save_fig('Mixed layer depth March', path_pics, fig_name, FigInf)
    
    # ---
    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,
                               IcD=IcD_monthly,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
      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',
                               IcD=IcD,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
      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',
                               IcD=IcD,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
      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',
                               IcD=IcD,
                               save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                               **Ddict_global)
      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,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )
      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)
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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 concentration NH', path_pics, fig_name, FigInf)

    # ---
    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)
      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',
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      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',
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      pyic.hplot_base(IcD_monthly, IaV, clim=[0,1], cincr=0.05, cmap='RdYlBu_r',
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      pyic.hplot_base(IcD_monthly, IaV, clim=[0,1], cincr=0.05, cmap='RdYlBu_r',
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[60.,90.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[-90., -50.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      pyic.hplot_base(IcD_monthly, IaV, cmap='RdYlBu_r',
                      clim=[0,1], clevs=np.array([0,1.,5,10,15,20,30,40,50,60,70,80,85,90,95,99,100])/100.,
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[-90., -50.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',
                      )

      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)
      save_fig('Sea ice concentration SH', path_pics, fig_name, FigInf)

    # ---
    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)
      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',
                      projection='PlateCarree', xlim=[-180.,180.], ylim=[-90., -50.],
                      ax=ax, cax=cax,
                      crs_features=False, do_plot_settings=False, do_write_data_range=True,
                      save_data=save_data, fpath_nc=path_nc+fig_name+'.nc',