diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3c4c2e6..495409e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
 Unreleased is the current development version, which is currently lying in `main` branch.
 
+## [v0.1.7]
+
+- global_mean can plot also values where observation are missing through the `--addnan` option (#95)
+- Verbosity level reduced in both global mean and performance indices (#95)
+- New colorscale `vlag` for global_mean (#95)
+
 ## [v0.1.6]
 
 - Remove support for python 3.8
diff --git a/docs/sphinx/source/globalmean.rst b/docs/sphinx/source/globalmean.rst
index 4612982..94a4359 100644
--- a/docs/sphinx/source/globalmean.rst
+++ b/docs/sphinx/source/globalmean.rst
@@ -70,6 +70,10 @@ Optional Arguments
 
    Specify the model name, overriding the configuration specified in ``config.yml``.
 
+.. option:: --addnan
+
+   Activate to plot also in the heatmap also fields which does not have a comparison against observations. Default is False.
+
 .. option:: -v LOGLEVEL, --loglevel LOGLEVEL
 
    Define the level of logging. The default is 'warning'.
diff --git a/ecmean/__init__.py b/ecmean/__init__.py
index 50e24ef..a028bfc 100644
--- a/ecmean/__init__.py
+++ b/ecmean/__init__.py
@@ -1,6 +1,6 @@
 """ECmean4 module"""
 
-__version__ = '0.1.6'
+__version__ = '0.1.7'
 
 # functions to be accessible everywhere
 from ecmean.libs.diagnostic import Diagnostic
diff --git a/ecmean/global_mean.py b/ecmean/global_mean.py
index b915ae0..191bdff 100755
--- a/ecmean/global_mean.py
+++ b/ecmean/global_mean.py
@@ -139,7 +139,7 @@ def gm_worker(util, ref, face, diag, varmean, vartrend, varlist):
                             if len(avg) == len(diag.years_joined):
                                 trend[season][region] = np.polyfit(diag.years_joined, avg, 1)[0]
                         if season == 'ALL' and region == 'Global':
-                            print('Average:', var, season, region, result[season][region])
+                            loggy.info('Average: %s %s %s %s', var, season, region, result[season][region])
 
         # nested dictionary, to be redifend as a dict to remove lambdas
         varmean[var] = result
@@ -151,6 +151,7 @@ def global_mean(exp, year1, year2,
                 loglevel='WARNING',
                 numproc=1,
                 interface=None, model=None, ensemble='r1i1p1f1',
+                addnan=False,
                 silent=None, trend=None, line=None,
                 outputdir=None, xdataset=None):
     """The main ECmean4 global mean function
@@ -164,6 +165,7 @@ def global_mean(exp, year1, year2,
     :param interface: interface file to be used, optional (default as specifified in config file)
     :param model: model to be analyzed, optional (default as specifified in config file)
     :param ensemble: ensemble member to be analyzed, optional (default as 'r1i1p1f1')
+    :param add_nan: add to the final plots also fields which cannot be compared against observations
     :param silent: do not print anything to std output, optional
     :param trend: compute yearly trends, optional
     :param line: appends also single line to a table, optional
@@ -207,6 +209,7 @@ def global_mean(exp, year1, year2,
     # get file info
     inifiles = get_inifiles(face, diag)
 
+
     # add missing unit definition
     units_extra_definition()
 
@@ -240,7 +243,7 @@ def global_mean(exp, year1, year2,
     toc = time()
 
     # evaluate tic-toc time  of execution
-    loggy.warning('Analysis done in {:.4f} seconds'.format(toc - tic))
+    loggy.info('Analysis done in {:.4f} seconds'.format(toc - tic))
 
     tic = time()
 
@@ -261,11 +264,15 @@ def global_mean(exp, year1, year2,
         # extract from yaml table for obs mean and standard deviation
         mmm = init_mydict(diag.seasons, diag.regions)
         sss = init_mydict(diag.seasons, diag.regions)
+        # if we have all the obs/std available
         if isinstance(gamma['obs'], dict):
             for season in diag.seasons:
                 for region in diag.regions:
                     mmm[season][region] = gamma['obs'][season][region]['mean']
                     sss[season][region] = gamma['obs'][season][region]['std']
+        # if only global observation is available
+        else:
+            mmm['ALL']['Global'] = gamma['obs']
         obsmean[gamma['longname']] = mmm
         obsstd[gamma['longname']] = sss
 
@@ -324,7 +331,8 @@ def global_mean(exp, year1, year2,
         f'global_mean_{diag.expname}_{diag.modelname}_r1i1p1f1_{diag.year1}_{diag.year2}.pdf'
     loggy.info('Figure file is: %s', mapfile)
 
-    heatmap_comparison_gm(data_table, mean_table, std_table, diag, mapfile)
+    heatmap_comparison_gm(data_table, mean_table, std_table,
+                          diag, mapfile, addnan=diag.addnan)
 
     # Print appending one line to table (for tuning)
     if diag.ftable:
@@ -333,7 +341,8 @@ def global_mean(exp, year1, year2,
 
     toc = time()
     # evaluate tic-toc time of postprocessing
-    loggy.warning(f"Postproc done in {toc - tic:.4f} seconds")
+    loggy.info(f"Postproc done in {toc - tic:.4f} seconds")
+    print('ECmean4 Global Mean succesfully computed!')
 
 
 
@@ -351,8 +360,6 @@ def gm_entry_point():
                 loglevel=args.loglevel,
                 interface=args.interface, config=args.config,
                 model=args.model, ensemble=args.ensemble,
+                addnan=args.addnan,
                 outputdir=args.outputdir)
 
-#if __name__ == "__main__":
-#
-#    gm_entry_point()
diff --git a/ecmean/libs/areas.py b/ecmean/libs/areas.py
index 0a451c3..41fd686 100644
--- a/ecmean/libs/areas.py
+++ b/ecmean/libs/areas.py
@@ -195,7 +195,7 @@ def area_cell(xfield, gridtype=None, formula='triangles'):
     outfield = xr.DataArray(area, dims=area_dims, coords=xfield.coords, name='area')
 
     # check the total area
-    loggy.info('Total Earth Surface: %s Km2',
+    loggy.debug('Total Earth Surface: %s Km2',
                  str(outfield.sum().values / 10**6))
 
     return outfield
diff --git a/ecmean/libs/diagnostic.py b/ecmean/libs/diagnostic.py
index 888ff4b..62c684f 100644
--- a/ecmean/libs/diagnostic.py
+++ b/ecmean/libs/diagnostic.py
@@ -41,6 +41,7 @@ def __init__(self, args):
         self.interface = getattr(args, 'interface', '')
         self.resolution = getattr(args, 'resolution', '')
         self.ensemble = getattr(args, 'ensemble', 'r1i1p1f1')
+        self.addnan = getattr(args, 'addnan', False)
         self.funcname = args.funcname.split(".")[1]
         self.version = __version__
         if self.year1 == self.year2:
diff --git a/ecmean/libs/files.py b/ecmean/libs/files.py
index 0b67dac..c5b9ec7 100644
--- a/ecmean/libs/files.py
+++ b/ecmean/libs/files.py
@@ -85,7 +85,7 @@ def var_is_there(flist, var, face):
         if not varunit:
             varunit = units_avail.get(var_req[0])
             if len(var_req) > 1:
-                loggy.info('%s is a derived var, assuming unit '
+                loggy.debug('%s is a derived var, assuming unit '
                              'as the first of its term', var)
 
     else:
@@ -207,7 +207,7 @@ def _filter_filename_by_year(template, filenames, year):
     if not filternames and len(filenames) > 0:
         loggy.warning('Data for year %s has not been found!', str(year))
 
-    loggy.info('Filtered filenames: %s', filternames)
+    loggy.debug('Filtered filenames: %s', filternames)
     return filternames
 
 
@@ -230,7 +230,7 @@ def _create_filepath(cmorname, face, diag):
         Path(face['model']['basedir']) / \
         Path(face['filetype'][filetype]['dir']) / \
         Path(face['filetype'][filetype]['filename'])
-    loggy.info('Filepath: %s', filepath)
+    loggy.debug('Filepath: %s', filepath)
 
     return filepath
 
@@ -265,7 +265,7 @@ def make_input_filename(cmorname, face, diag):
                 filename1 = filename1 + fname
             filename = filename + filename1
         filename = list(dict.fromkeys(filename))  # Filter unique ones
-        loggy.info("Filenames: %s", filename)
+        loggy.debug("Filenames: %s", filename)
         return filename
 
 
diff --git a/ecmean/libs/parser.py b/ecmean/libs/parser.py
index d3454c3..a5f9bb8 100644
--- a/ecmean/libs/parser.py
+++ b/ecmean/libs/parser.py
@@ -34,6 +34,8 @@ def parse_arguments(args, script):
                         help='variant label (ripf number for cmor)')
     parser.add_argument('-s', '--silent', action='store_true',
                         help='do not print anything to std output')
+    parser.add_argument('--addnan', action='store_true',
+                        help='provide figures also where observations are missing')
     parser.add_argument('-v', '--loglevel', type=str, default='WARNING',
                         help='define the level of logging.')
     parser.add_argument('-o', '--outputdir', type=str,
diff --git a/ecmean/libs/plotting.py b/ecmean/libs/plotting.py
index e9be213..1004d87 100644
--- a/ecmean/libs/plotting.py
+++ b/ecmean/libs/plotting.py
@@ -12,10 +12,11 @@
 import matplotlib.pyplot as plt
 from matplotlib.colors import ListedColormap
 import seaborn as sns
+import numpy as np
 
 
 
-def heatmap_comparison_pi(relative_table, diag, filemap):
+def heatmap_comparison_pi(relative_table, diag, filemap, size_model = 14):
     """Function to produce a heatmap - seaborn based - for Performance Indices
     based on CMIP6 ratio"""
 
@@ -42,7 +43,8 @@ def heatmap_comparison_pi(relative_table, diag, filemap):
     chart = sns.heatmap(myfield, norm=divnorm, cmap=pal,
                         cbar_kws={"ticks": tictoc, 'label': title},
                         ax=axs, annot=True, linewidth=0.5, fmt='.2f',
-                        annot_kws={'fontsize': 14, 'fontweight': 'bold'})
+                        annot_kws={'fontsize': size_model, 'fontweight': 'bold'})
+    
     chart = chart.set_facecolor('whitesmoke')
     axs.set_title(f'{title} {diag.modelname} {diag.expname} {diag.year1} {diag.year2}', fontsize=25)
     axs.vlines(list(range(sss, tot + sss, sss)), ymin=-1, ymax=len(myfield.index), colors='k')
@@ -60,37 +62,54 @@ def heatmap_comparison_pi(relative_table, diag, filemap):
     plt.close()
 
 
-def heatmap_comparison_gm(data_table, mean_table, std_table, diag, filemap):
+def heatmap_comparison_gm(data_table, mean_table, std_table, diag, filemap, addnan = True,
+                          size_model = 14, size_obs = 8):
     """Function to produce a heatmap - seaborn based - for Global Mean
     based on season-averaged standard deviation ratio"""
 
     # define array
     ratio = (data_table - mean_table) / std_table
-    mask = ratio[('ALL', 'Global')].notna()
+    if addnan:
+        mask = data_table[('ALL', 'Global')].notna()
+    else:
+        mask = ratio[('ALL', 'Global')].notna()
     clean = ratio[mask]
 
     # for dimension of plots
     xfig = len(clean.columns)
     yfig = len(clean.index)
-    _, axs = plt.subplots(1, 1, sharey=True, tight_layout=True, figsize=(xfig+5, yfig+2))
+    fig, axs = plt.subplots(1, 1, sharey=True, tight_layout=True, figsize=(xfig+5, yfig+2))
 
     title = 'GLOBAL MEAN'
 
     # set color range and palette
     thr = 10
     tictoc = range(-thr, thr + 1)
-    pal = ListedColormap(sns.color_palette("seismic", n_colors=21))
+    pal = ListedColormap(sns.color_palette("vlag", n_colors=21))
     tot = len(clean.columns)
     sss = len(set([tup[1] for tup in clean.columns]))
 
     chart = sns.heatmap(clean, annot=data_table[mask], vmin=-thr - 0.5, vmax=thr + 0.5, center=0,
-                        annot_kws={'va': 'bottom', 'fontsize': 14},
-                        cbar_kws={'ticks': tictoc,
+                        annot_kws={'va': 'bottom', 'fontsize': size_model},
+                        cbar_kws={'ticks': tictoc, "shrink": .5,
                                   'label': 'Model Bias \n (standard deviation of interannual variability from observations)'},
                         fmt='.2f', cmap=pal)
+    if addnan:
+        empty = np.where(clean.isna(), 0, np.nan)
+        empty = np.where(data_table[mask]==0, np.nan, empty)
+        chart = sns.heatmap(empty, annot=data_table[mask], fmt='.2f',
+                            vmin=-thr - 0.5, vmax=thr + 0.5, center=0,
+                            annot_kws={'va': 'bottom', 'fontsize': size_model, 'color':'dimgrey'}, cbar=False,
+                            cmap=ListedColormap(['none']))
     chart = sns.heatmap(clean, annot=mean_table[mask], vmin=-thr - 0.5, vmax=thr + 0.5, center=0,
-                        annot_kws={'va': 'top', 'fontsize': 8, 'fontstyle': 'italic'},
+                        annot_kws={'va': 'top', 'fontsize': size_obs, 'fontstyle': 'italic'},
                         fmt='.2f', cmap=pal, cbar=False)
+    if addnan:
+        empty = np.where(clean.isna(), 0, np.nan)
+        empty = np.where(mean_table[mask].isna(), np.nan, empty)
+        chart = sns.heatmap(empty, annot=mean_table[mask], vmin=-thr - 0.5, vmax=thr + 0.5, center=0,
+                        annot_kws={'va': 'top', 'fontsize': size_obs, 'fontstyle': 'italic', 'color':'dimgrey'},
+                        fmt='.2f', cmap=ListedColormap(['none']), cbar=False)
 
     chart = chart.set_facecolor('whitesmoke')
     axs.set_title(f'{title} {diag.modelname} {diag.expname} {diag.year1} {diag.year2}', fontsize=25)
diff --git a/ecmean/libs/support.py b/ecmean/libs/support.py
index 57f3cf3..f9608b3 100644
--- a/ecmean/libs/support.py
+++ b/ecmean/libs/support.py
@@ -31,7 +31,7 @@ def __init__(self, component, atmdict, ocedict, areas=True, remap=False, targetg
         """Class for masks, areas and interpolation (xESMF-based)
         for both atmospheric and oceanic component"""
 
-        loggy.info('Running with xesmf version %s', xe.__version__)
+        loggy.debug('Running with xesmf version %s', xe.__version__)
 
         # define the basics
         self.atmareafile = inifiles_priority(atmdict)
diff --git a/ecmean/libs/units.py b/ecmean/libs/units.py
index eacf4c8..f785814 100644
--- a/ecmean/libs/units.py
+++ b/ecmean/libs/units.py
@@ -54,7 +54,7 @@ def __init__(self, var=None, org_units=None, tgt_units=None,
         if face and clim:
             self.init_ecmean(clim, face)
 
-        loggy.info(vars(self))
+        loggy.debug(vars(self))
         # if units are defined and convert called
         if convert and self.org_units and self.tgt_units:
 
@@ -101,12 +101,12 @@ def units_converter(self):
         It will not work if BOTH factor and offset are required"""
 
         units_relation = (self.org_units / self.tgt_units).to_base_units()
-        loggy.info('Original vs Target units relation is %s', units_relation)
+        loggy.debug('Original vs Target units relation is %s', units_relation)
 
         if units_relation.units == units('dimensionless'):
 
             if units_relation.magnitude != 1:
-                loggy.info('Unit conversion required...')
+                loggy.debug('Unit conversion required...')
                 offset_standard = 0 * self.org_units
                 offset = offset_standard.to(self.tgt_units).magnitude
                 if offset == 0:
@@ -119,16 +119,16 @@ def units_converter(self):
                 factor = 1.
 
         elif units_relation.units == units('kg / m^3'):
-            loggy.info("Assuming this as a water flux! Am I correct?")
-            loggy.info("Dividing by water density...")
+            loggy.debug("Assuming this as a water flux! Am I correct?")
+            loggy.debug("Dividing by water density...")
             density_water = units('kg / m^3') * 1000
             offset = 0.
             factor_standard = 1 * self.org_units
             factor = (factor_standard / density_water).to(self.tgt_units).magnitude
 
         elif units_relation.units == units('s'):
-            loggy.info("Assuming this is a cumulated flux...")
-            loggy.info("Dividing by cumulation time (expressed in seconds)...")
+            loggy.debug("Assuming this is a cumulated flux...")
+            loggy.debug("Dividing by cumulation time (expressed in seconds)...")
             if self.cumulation_time:
                 cumtime = units('s') * self.cumulation_time
                 offset = 0.
@@ -145,8 +145,8 @@ def units_converter(self):
         if self.org_direction != self.tgt_direction:
             factor = -1. * factor
 
-        loggy.info('Offset is %s', str(offset))
-        loggy.info('Factor is %s', str(factor))
+        loggy.debug('Offset is %s', str(offset))
+        loggy.debug('Factor is %s', str(factor))
         return offset, factor
 
 
@@ -159,6 +159,7 @@ def units_extra_definition():
     # https://github.com/Unidata/MetPy/issues/2884
     units._on_redefinition = 'ignore'
     units.define('fraction = [] = frac')
+    units.define('Fraction = [] = frac')
     units.define('psu = 1e-3 frac')
     units.define('PSU = 1e-3 frac')
     units.define('million = 1e6 = M')
diff --git a/ecmean/performance_indices.py b/ecmean/performance_indices.py
index abe14e6..e5c6367 100755
--- a/ecmean/performance_indices.py
+++ b/ecmean/performance_indices.py
@@ -99,7 +99,7 @@ def pi_worker(util, piclim, face, diag, field_3d, varstat, varlist):
                 # mean over time and fixing of the units
                 for season in diag.seasons:
 
-                    loggy.info(season)
+                    loggy.debug(season)
 
                     # copy of the full field
                     tmean = outfield.copy(deep=True)
@@ -202,7 +202,7 @@ def pi_worker(util, piclim, face, diag, field_3d, varstat, varlist):
 
                         # diagnostic
                         if region == 'Global':
-                            print('PI for', region, season, var, result[season][region])
+                            logging.info('PI for', region, season, var, result[season][region])
 
         # nested dictionary, to be redifend as a dict to remove lambdas
         varstat[var] = result
@@ -291,7 +291,7 @@ def performance_indices(exp, year1, year2,
     processes = []
 
     toc = time()
-    loggy.warning('Preproc in {:.4f} seconds'.format(toc - tic))
+    loggy.info('Preproc in {:.4f} seconds'.format(toc - tic))
     tic = time()
 
   
@@ -313,7 +313,7 @@ def performance_indices(exp, year1, year2,
 
     toc = time()
     # evaluate tic-toc time  of execution
-    loggy.warning('Done in {:.4f} seconds'.format(toc - tic) + ' with ' + str(diag.numproc) + ' processors')
+    loggy.info('Done in {:.4f} seconds'.format(toc - tic) + ' with ' + str(diag.numproc) + ' processors')
 
     tic = time()
 
@@ -349,7 +349,7 @@ def performance_indices(exp, year1, year2,
 
         # convert output dictionary to pandas dataframe
         data_table = dict_to_dataframe(plotted)
-        loggy.info(data_table)
+        loggy.debug(data_table)
 
 
         # relative pi with re-ordering of rows
@@ -361,7 +361,7 @@ def performance_indices(exp, year1, year2,
         # reordering columns
         lll = [(x, y) for x in diag.seasons for y in diag.regions]
         cmip6_table = cmip6_table[lll]
-        loggy.info(cmip6_table)
+        loggy.debug(cmip6_table)
 
 
         # call the heatmap routine for a plot
@@ -372,7 +372,8 @@ def performance_indices(exp, year1, year2,
 
     toc = time()
     # evaluate tic-toc time of postprocessing
-    loggy.warning('Postproc done in {:.4f} seconds'.format(toc - tic))
+    loggy.info('Postproc done in {:.4f} seconds'.format(toc - tic))
+    print('ECmean4 Performance Indices succesfully computed!')
 
 
 def pi_entry_point():
@@ -391,7 +392,3 @@ def pi_entry_point():
                         model=args.model, ensemble=args.ensemble,
                         outputdir=args.outputdir)
 
-
-#if __name__ == '__main__':
-#
-#    sys.exit(pi_entry_point())
diff --git a/tests/test_global_mean.py b/tests/test_global_mean.py
index 5c21bf8..b1a8f78 100644
--- a/tests/test_global_mean.py
+++ b/tests/test_global_mean.py
@@ -61,7 +61,7 @@ def test_global_mean_amip():
     if os.path.isfile(file1):
         os.remove(file1)
     global_mean(exp='amip', year1=1990, year2=1990, numproc=1, config='tests/config.yml',
-                line=True)
+                line=True, addnan=True)
     
     data1 = load_gm_txt_files(file1)
     data2 = load_gm_txt_files(file2)