Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/launch_sasf.sh
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/launch_sasf.sh	(revision 11962)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/launch_sasf.sh	(revision 11996)
@@ -40,5 +40,5 @@
 rsync -avP ${FORC_DIR}/Station_PAPA_50N-145W*.nc ${WORK_DIR}/
 
-for CASE in "ECMWF" "COARE3p6" "NCAR"; do
+for CASE in "ECMWF-noskin" "COARE3p6-noskin" "ECMWF" "COARE3p6" "NCAR"; do
 
     echo ; echo
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg	(revision 11996)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg	(revision 11996)
@@ -0,0 +1,284 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  Configuration namelist : overwrite default values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!             STATION_ASF configuration   
+!!======================================================================
+!!              ***  Domain & Run management namelists  ***           !!
+!!                                                                    !!
+!!   namrun       parameters of the run
+!!   namdom       space and time domain
+!!   namcfg       parameters of the configuration                       (default: user defined GYRE)
+!!   namwad       Wetting and drying                                    (default: OFF)
+!!   namtsd       data: temperature & salinity                          (default: OFF)
+!!   namcrs       coarsened grid (for outputs and/or TOP)               (ln_crs =T)
+!!   namc1d       1D configuration options                              ("key_c1d")
+!!   namc1d_dyndmp 1D newtonian damping applied on currents             ("key_c1d")
+!!   namc1d_uvd   1D data (currents)                                    ("key_c1d")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namusr_def    !   STATION_ASF user defined namelist
+!-----------------------------------------------------------------------
+   rn_dept1    =  1.0   ! depth (m) at which the prescribed SST is taken  == depth of first T point! (IMPORTANT for warm-layer impact, the deeper, the more pronounced!)
+   !                    ! [PAPA station = 1 m]
+/
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   nn_no       =       0   !  Assimilation cycle index
+   cn_exp      =  'STATION_ASF-COARE3p6-noskin'  !  experience name
+   nn_it000    =    1   !  first time step
+   nn_itend    =    26280   !  last  time step (std 5840)
+   nn_date0    =    20160101   !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
+   nn_time0    =       0   !  initial time of day in hhmm
+   nn_leapy    =       0   !  Leap year calendar (1) or not (0)
+   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
+      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
+      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
+      !                          !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart
+      !                          !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart
+      cn_ocerst_in    = ''   !  suffix of ocean restart name (input)
+      cn_ocerst_indir = ''         !  directory from which to read input ocean restarts
+      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output)
+      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts
+   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model
+   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
+   nn_stock    =    26280   ! 1year @ dt=3600 s / frequency of creation of a restart file (modulo referenced to 1)
+   nn_write    =    26280   ! 1year @ dt=3600 s / frequency of write in the output file   (modulo referenced to nn_it000)
+   ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%)
+   ln_cfmeta   = .false.   !  output additional data to netCDF files required for compliance with the CF metadata standard
+   ln_clobber  = .true.    !  clobber (overwrite) an existing file
+   nn_chunksz  =       0   !  chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
+   ln_xios_read = .false.  !  use XIOS to read restart file (only for a single file restart)
+   nn_wxios = 0      !  use XIOS to write restart file 0 - no, 1 - single file output, 2 - multiple file output
+/
+!-----------------------------------------------------------------------
+&namdom        !   time and space domain / STATION_ASF
+!-----------------------------------------------------------------------
+   rn_rdt      =  3600.    !  time step for the dynamics and tracer => frequency of atmospheric forcing
+   ln_meshmask = .false.
+   ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namc1d        !   1D configuration options / STATION_ASF               ("key_c1d" default: PAPA station)
+!-----------------------------------------------------------------------
+   rn_lat1d    =      50.1 !  Column latitude    [PAPA station]
+   rn_lon1d    =    215.1  !  Column longitude   [PAPA station]
+/
+!-----------------------------------------------------------------------
+&namc1d_dyndmp !   U & V newtonian damping                              ("key_c1d" default: OFF)
+!-----------------------------------------------------------------------
+   ln_dyndmp   =  .false.  !  add a damping term (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namc1d_uvd    !   data: U & V currents                                 ("key_c1d" default: OFF)
+!-----------------------------------------------------------------------
+   !                       !  =T read U-V fields for:
+   ln_uvd_init   = .false.       !  ocean initialisation
+   ln_uvd_dyndmp = .false.       !  U-V restoring
+
+   cn_dir      = './'      !  root directory for the U-V data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_ucur     = 'ucurrent'              ,         -1.       ,'u_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Ume'   , ''
+   sn_vcur     = 'vcurrent'              ,         -1.       ,'v_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Vme'   , ''
+/
+
+!!======================================================================
+!!            ***  Surface Boundary Condition namelists  ***          !!
+!!                                                                    !!
+!!                           STATION_ASF case                         !!
+!!                                                                    !!
+!!   namsbc          surface boundary condition manager                 (default: NO selection)
+!!   namsbc_flx      flux               formulation                     (ln_flx     =T)
+!!   namsbc_blk      Bulk formulae formulation                          (ln_blk     =T)
+!!   namsbc_cpl      CouPLed            formulation                     ("key_oasis3" )
+!!   namsbc_sas      Stand-Alone Surface module                         (SAS_SRC  only)
+!!   namsbc_iif      Ice-IF: use observed ice cover                     (nn_ice = 1   )
+!!   namsbc_ssr      sea surface restoring term (for T and/or S)        (ln_ssr     =T)
+!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T)
+!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T)
+!!   namsbc_wave     external fields from wave model                    (ln_wave    =T)
+!!   namberg         iceberg floats                                     (ln_icebergs=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
+!-----------------------------------------------------------------------
+   nn_fsbc     = 1         !!LOLO!  frequency of SBC module call
+   ln_blk      = .true.    !  Bulk formulation                          (T => fill namsbc_blk )
+/
+!-----------------------------------------------------------------------
+&namsbc_blk    !   namsbc_blk  generic Bulk formula / STATION_ASF       (ln_blk =T)
+!-----------------------------------------------------------------------
+   !                    !  bulk algorithm :
+   ln_NCAR      = .false.   ! "NCAR"      algorithm   (Large and Yeager 2008)
+   ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
+   ln_COARE_3p6 = .true.    ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
+      rn_zqt    =  2.5      !  Air temperature & humidity reference height (m) [PAPA station]
+      rn_zu     =  4.       !  Wind vector reference height (m)                [PAPA station]
+      !
+      ! Skin is ONLY available in ECMWF and COARE algorithms:
+      ln_skin_cs = .false.  !  use the cool-skin parameterization => set nn_fsbc=1 and ln_dm2dc=.true.!
+      ln_skin_wl = .false.  !  use the warm-layer        "        => set nn_fsbc=1 and ln_dm2dc=.true.!
+      !
+      ln_humi_sph = .false. !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true.
+      ln_humi_dpt = .false. !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true.
+      ln_humi_rlh = .true.  !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true.
+   !
+   cn_dir = './'  !  root directory for the bulk data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !       weights filename               ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
+   sn_wndi     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'u_air',    .true.     , .false. , 'yearly'  , ''   ,''   , ''
+   sn_wndj     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'v_air',    .true.     , .false. , 'yearly'  , ''   ,''   , ''
+   sn_qsr      = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'rad_sw',   .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_qlw      = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'rad_lw',   .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_tair     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 't_air',    .true.     , .false. , 'yearly'  , ''  , ''   , ''
+   sn_humi     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'rh_air',   .true.     , .false. , 'yearly'  , ''  , ''   , ''
+   sn_prec     = 'Station_PAPA_50N-145W_precip_daily', 24., 'precip',   .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_snow     = 'Station_PAPA_50N-145W_precip_daily', 24., 'snow',     .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_slp      = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'slp',      .true.     , .false. , 'yearly'  , ''  , ''   , ''
+   sn_tdif     = 'Station_PAPA_50N-145W_atm_hourly', 24.  , 'taudif'  , .false.    , .false. , 'yearly'  , ''  , ''   , ''
+
+/
+!-----------------------------------------------------------------------
+&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_sas    !   Stand-Alone Surface boundary condition / STATION_ASF
+!-----------------------------------------------------------------------
+   l_sasread   = .true.    !  =T Read in file ;  =F set all to 0. (see sbcssm)
+      ln_3d_uve   = .false.   !  specify whether we are supplying a 3D u,v and e3 field
+      ln_read_frq = .false.   !  specify whether we must read frq or not
+
+   cn_dir      = './'     !  root directory for the ocean data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_usp  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'ssu'     ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_vsp  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'ssv'     ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_tem  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'sst',        .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_sal  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'sss',        .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_ssh  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'ssh',        .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_e3t  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'e3t_m'   ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_frq  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'frq_m'   ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+/
+!!======================================================================
+!!               ***  Lateral boundary condition  ***                 !!
+!!                                                                    !!
+!!   namlbc        lateral momentum boundary condition                  (default: NO selection)
+!!   namagrif      agrif nested grid   (read by child model only)       ("key_agrif")
+!!   nam_tide      Tidal forcing                                        (default: OFF)
+!!   nambdy        Unstructured open boundaries                         (default: OFF)
+!!   nambdy_dta    Unstructured open boundaries - external data         (see  nambdy)
+!!   nambdy_tide   tidal forcing at open boundaries                     (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+   rn_shlat    =    0.     !  free slip
+/
+!!======================================================================
+!!                ***  Top/Bottom boundary condition  ***             !!
+!!                                                                    !!
+!!   namdrg        top/bottom drag coefficient                          (default: NO selection)
+!!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T)
+!!   namdrg_bot    bottom friction                                      (ln_OFF=F)
+!!   nambbc        bottom temperature boundary condition                (default: OFF)
+!!   nambbl        bottom boundary layer scheme                         (default: OFF)
+!!======================================================================
+!
+!!======================================================================
+!!                        Tracer (T-S) namelists                      !!
+!!                                                                    !!
+!!   nameos        equation of state                                    (default: NO selection)
+!!   namtra_adv    advection scheme                                     (default: NO selection)
+!!   namtra_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namtra_mle    mixed layer eddy param. (Fox-Kemper param.)          (default: OFF)
+!!   namtra_eiv    eddy induced velocity param.                         (default: OFF)
+!!   namtra_dmp    T & S newtonian damping                              (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nameos        !   ocean Equation Of Seawater                           (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_eos80    = .true.          !  = Use EOS80
+/
+!!======================================================================
+!!                      ***  Dynamics namelists  ***                  !!
+!!                                                                    !!
+!!   nam_vvl       vertical coordinate options                          (default: z-star)
+!!   namdyn_adv    formulation of the momentum advection                (default: NO selection)
+!!   namdyn_vor    advection scheme                                     (default: NO selection)
+!!   namdyn_hpg    hydrostatic pressure gradient                        (default: NO selection)
+!!   namdyn_spg    surface pressure gradient                            (default: NO selection)
+!!   namdyn_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namdta_dyn    offline TOP: dynamics read in files                  (OFF_SRC only)
+!!======================================================================
+!
+!!======================================================================
+!!                     vertical physics namelists                     !!
+!!                                                                    !!
+!!    namzdf        vertical physics manager                            (default: NO selection)
+!!    namzdf_ric    richardson number vertical mixing                   (ln_zdfric=T)
+!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
+!!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
+!!    namzdf_osm    OSM vertical diffusion                              (ln_zdfosm=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
+!!======================================================================
+!
+!!======================================================================
+!!                  ***  Diagnostics namelists  ***                   !!
+!!                                                                    !!
+!!   namtrd       dynamics and/or tracer trends                         (default: OFF)
+!!   namptr       Poleward Transport Diagnostics                        (default: OFF)
+!!   namhsb       Heat and salt budgets                                 (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namflo       float parameters                                      (default: OFF)
+!!   nam_diaharm  Harmonic analysis of tidal constituents               (default: OFF)
+!!   nam_diadct   transports through some sections                      (default: OFF)
+!!   nam_diatmb   Top Middle Bottom Output                              (default: OFF)
+!!   nam_dia25h   25h Mean Output                                       (default: OFF)
+!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4")
+!!======================================================================
+!
+!!======================================================================
+!!               ***  Observation & Assimilation  ***                 !!
+!!                                                                    !!
+!!   namobs       observation and model comparison                      (default: OFF)
+!!   nam_asminc   assimilation increments                               ('key_asminc')
+!!======================================================================
+!
+!!======================================================================
+!!                  ***  Miscellaneous namelists  ***                 !!
+!!                                                                    !!
+!!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi")
+!!   namctl            Control prints                                   (default: OFF)
+!!   namsto            Stochastic parametrization of EOS                (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints                                       (default: OFF)
+!-----------------------------------------------------------------------
+ln_ctl = .false.
+/
+!-----------------------------------------------------------------------
+&namsto        ! Stochastic parametrization of EOS                      (default: OFF)
+!-----------------------------------------------------------------------
+/
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg	(revision 11996)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg	(revision 11996)
@@ -0,0 +1,284 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  Configuration namelist : overwrite default values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!             STATION_ASF configuration   
+!!======================================================================
+!!              ***  Domain & Run management namelists  ***           !!
+!!                                                                    !!
+!!   namrun       parameters of the run
+!!   namdom       space and time domain
+!!   namcfg       parameters of the configuration                       (default: user defined GYRE)
+!!   namwad       Wetting and drying                                    (default: OFF)
+!!   namtsd       data: temperature & salinity                          (default: OFF)
+!!   namcrs       coarsened grid (for outputs and/or TOP)               (ln_crs =T)
+!!   namc1d       1D configuration options                              ("key_c1d")
+!!   namc1d_dyndmp 1D newtonian damping applied on currents             ("key_c1d")
+!!   namc1d_uvd   1D data (currents)                                    ("key_c1d")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namusr_def    !   STATION_ASF user defined namelist
+!-----------------------------------------------------------------------
+   rn_dept1    =  1.0   ! depth (m) at which the prescribed SST is taken  == depth of first T point! (IMPORTANT for warm-layer impact, the deeper, the more pronounced!)
+   !                    ! [PAPA station = 1 m]
+/
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   nn_no       =       0   !  Assimilation cycle index
+   cn_exp      =  'STATION_ASF-ECMWF-noskin'  !  experience name
+   nn_it000    =    1   !  first time step
+   nn_itend    =    26280   !  last  time step (std 5840)
+   nn_date0    =    20160101   !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
+   nn_time0    =       0   !  initial time of day in hhmm
+   nn_leapy    =       0   !  Leap year calendar (1) or not (0)
+   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
+      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
+      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
+      !                          !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart
+      !                          !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart
+      cn_ocerst_in    = ''   !  suffix of ocean restart name (input)
+      cn_ocerst_indir = ''         !  directory from which to read input ocean restarts
+      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output)
+      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts
+   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model
+   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
+   nn_stock    =    26280   ! 1year @ dt=3600 s / frequency of creation of a restart file (modulo referenced to 1)
+   nn_write    =    26280   ! 1year @ dt=3600 s / frequency of write in the output file   (modulo referenced to nn_it000)
+   ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%)
+   ln_cfmeta   = .false.   !  output additional data to netCDF files required for compliance with the CF metadata standard
+   ln_clobber  = .true.    !  clobber (overwrite) an existing file
+   nn_chunksz  =       0   !  chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
+   ln_xios_read = .false.  !  use XIOS to read restart file (only for a single file restart)
+   nn_wxios = 0      !  use XIOS to write restart file 0 - no, 1 - single file output, 2 - multiple file output
+/
+!-----------------------------------------------------------------------
+&namdom        !   time and space domain / STATION_ASF
+!-----------------------------------------------------------------------
+   rn_rdt      =  3600.    !  time step for the dynamics and tracer => frequency of atmospheric forcing
+   ln_meshmask = .false.
+   ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namc1d        !   1D configuration options / STATION_ASF               ("key_c1d" default: PAPA station)
+!-----------------------------------------------------------------------
+   rn_lat1d    =      50.1 !  Column latitude    [PAPA station]
+   rn_lon1d    =    215.1  !  Column longitude   [PAPA station]
+/
+!-----------------------------------------------------------------------
+&namc1d_dyndmp !   U & V newtonian damping                              ("key_c1d" default: OFF)
+!-----------------------------------------------------------------------
+   ln_dyndmp   =  .false.  !  add a damping term (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namc1d_uvd    !   data: U & V currents                                 ("key_c1d" default: OFF)
+!-----------------------------------------------------------------------
+   !                       !  =T read U-V fields for:
+   ln_uvd_init   = .false.       !  ocean initialisation
+   ln_uvd_dyndmp = .false.       !  U-V restoring
+
+   cn_dir      = './'      !  root directory for the U-V data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_ucur     = 'ucurrent'              ,         -1.       ,'u_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Ume'   , ''
+   sn_vcur     = 'vcurrent'              ,         -1.       ,'v_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Vme'   , ''
+/
+
+!!======================================================================
+!!            ***  Surface Boundary Condition namelists  ***          !!
+!!                                                                    !!
+!!                           STATION_ASF case                         !!
+!!                                                                    !!
+!!   namsbc          surface boundary condition manager                 (default: NO selection)
+!!   namsbc_flx      flux               formulation                     (ln_flx     =T)
+!!   namsbc_blk      Bulk formulae formulation                          (ln_blk     =T)
+!!   namsbc_cpl      CouPLed            formulation                     ("key_oasis3" )
+!!   namsbc_sas      Stand-Alone Surface module                         (SAS_SRC  only)
+!!   namsbc_iif      Ice-IF: use observed ice cover                     (nn_ice = 1   )
+!!   namsbc_ssr      sea surface restoring term (for T and/or S)        (ln_ssr     =T)
+!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T)
+!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T)
+!!   namsbc_wave     external fields from wave model                    (ln_wave    =T)
+!!   namberg         iceberg floats                                     (ln_icebergs=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
+!-----------------------------------------------------------------------
+   nn_fsbc     = 1         !!LOLO!  frequency of SBC module call
+   ln_blk      = .true.    !  Bulk formulation                          (T => fill namsbc_blk )
+/
+!-----------------------------------------------------------------------
+&namsbc_blk    !   namsbc_blk  generic Bulk formula / STATION_ASF       (ln_blk =T)
+!-----------------------------------------------------------------------
+   !                    !  bulk algorithm :
+   ln_NCAR      = .false.   ! "NCAR"      algorithm   (Large and Yeager 2008)
+   ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
+   ln_COARE_3p6 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   ln_ECMWF     = .true.    ! "ECMWF"     algorithm   (IFS cycle 31)
+      rn_zqt    =  2.5      !  Air temperature & humidity reference height (m) [PAPA station]
+      rn_zu     =  4.       !  Wind vector reference height (m)                [PAPA station]
+      !
+      ! Skin is ONLY available in ECMWF and COARE algorithms:
+      ln_skin_cs = .false.  !  use the cool-skin parameterization => set nn_fsbc=1 and ln_dm2dc=.true.!
+      ln_skin_wl = .false.  !  use the warm-layer        "        => set nn_fsbc=1 and ln_dm2dc=.true.!
+      !
+      ln_humi_sph = .false. !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true.
+      ln_humi_dpt = .false. !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true.
+      ln_humi_rlh = .true.  !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true.
+   !
+   cn_dir = './'  !  root directory for the bulk data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !       weights filename               ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
+   sn_wndi     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'u_air',    .true.     , .false. , 'yearly'  , ''   ,''   , ''
+   sn_wndj     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'v_air',    .true.     , .false. , 'yearly'  , ''   ,''   , ''
+   sn_qsr      = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'rad_sw',   .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_qlw      = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'rad_lw',   .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_tair     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 't_air',    .true.     , .false. , 'yearly'  , ''  , ''   , ''
+   sn_humi     = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'rh_air',   .true.     , .false. , 'yearly'  , ''  , ''   , ''
+   sn_prec     = 'Station_PAPA_50N-145W_precip_daily', 24., 'precip',   .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_snow     = 'Station_PAPA_50N-145W_precip_daily', 24., 'snow',     .false.    , .false. , 'yearly'  , ''  , ''   , ''
+   sn_slp      = 'Station_PAPA_50N-145W_atm_hourly',  1.  , 'slp',      .true.     , .false. , 'yearly'  , ''  , ''   , ''
+   sn_tdif     = 'Station_PAPA_50N-145W_atm_hourly', 24.  , 'taudif'  , .false.    , .false. , 'yearly'  , ''  , ''   , ''
+
+/
+!-----------------------------------------------------------------------
+&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_sas    !   Stand-Alone Surface boundary condition / STATION_ASF
+!-----------------------------------------------------------------------
+   l_sasread   = .true.    !  =T Read in file ;  =F set all to 0. (see sbcssm)
+      ln_3d_uve   = .false.   !  specify whether we are supplying a 3D u,v and e3 field
+      ln_read_frq = .false.   !  specify whether we must read frq or not
+
+   cn_dir      = './'     !  root directory for the ocean data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_usp  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'ssu'     ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_vsp  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'ssv'     ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_tem  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'sst',        .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_sal  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'sss',        .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_ssh  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'ssh',        .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_e3t  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'e3t_m'   ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_frq  = 'Station_PAPA_50N-145W_oce_hourly', 1.         , 'frq_m'   ,   .true.    , .false. , 'yearly'  ,    ''            ,    ''    ,    ''
+/
+!!======================================================================
+!!               ***  Lateral boundary condition  ***                 !!
+!!                                                                    !!
+!!   namlbc        lateral momentum boundary condition                  (default: NO selection)
+!!   namagrif      agrif nested grid   (read by child model only)       ("key_agrif")
+!!   nam_tide      Tidal forcing                                        (default: OFF)
+!!   nambdy        Unstructured open boundaries                         (default: OFF)
+!!   nambdy_dta    Unstructured open boundaries - external data         (see  nambdy)
+!!   nambdy_tide   tidal forcing at open boundaries                     (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+   rn_shlat    =    0.     !  free slip
+/
+!!======================================================================
+!!                ***  Top/Bottom boundary condition  ***             !!
+!!                                                                    !!
+!!   namdrg        top/bottom drag coefficient                          (default: NO selection)
+!!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T)
+!!   namdrg_bot    bottom friction                                      (ln_OFF=F)
+!!   nambbc        bottom temperature boundary condition                (default: OFF)
+!!   nambbl        bottom boundary layer scheme                         (default: OFF)
+!!======================================================================
+!
+!!======================================================================
+!!                        Tracer (T-S) namelists                      !!
+!!                                                                    !!
+!!   nameos        equation of state                                    (default: NO selection)
+!!   namtra_adv    advection scheme                                     (default: NO selection)
+!!   namtra_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namtra_mle    mixed layer eddy param. (Fox-Kemper param.)          (default: OFF)
+!!   namtra_eiv    eddy induced velocity param.                         (default: OFF)
+!!   namtra_dmp    T & S newtonian damping                              (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nameos        !   ocean Equation Of Seawater                           (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_eos80    = .true.          !  = Use EOS80
+/
+!!======================================================================
+!!                      ***  Dynamics namelists  ***                  !!
+!!                                                                    !!
+!!   nam_vvl       vertical coordinate options                          (default: z-star)
+!!   namdyn_adv    formulation of the momentum advection                (default: NO selection)
+!!   namdyn_vor    advection scheme                                     (default: NO selection)
+!!   namdyn_hpg    hydrostatic pressure gradient                        (default: NO selection)
+!!   namdyn_spg    surface pressure gradient                            (default: NO selection)
+!!   namdyn_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namdta_dyn    offline TOP: dynamics read in files                  (OFF_SRC only)
+!!======================================================================
+!
+!!======================================================================
+!!                     vertical physics namelists                     !!
+!!                                                                    !!
+!!    namzdf        vertical physics manager                            (default: NO selection)
+!!    namzdf_ric    richardson number vertical mixing                   (ln_zdfric=T)
+!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
+!!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
+!!    namzdf_osm    OSM vertical diffusion                              (ln_zdfosm=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
+!!======================================================================
+!
+!!======================================================================
+!!                  ***  Diagnostics namelists  ***                   !!
+!!                                                                    !!
+!!   namtrd       dynamics and/or tracer trends                         (default: OFF)
+!!   namptr       Poleward Transport Diagnostics                        (default: OFF)
+!!   namhsb       Heat and salt budgets                                 (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namflo       float parameters                                      (default: OFF)
+!!   nam_diaharm  Harmonic analysis of tidal constituents               (default: OFF)
+!!   nam_diadct   transports through some sections                      (default: OFF)
+!!   nam_diatmb   Top Middle Bottom Output                              (default: OFF)
+!!   nam_dia25h   25h Mean Output                                       (default: OFF)
+!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4")
+!!======================================================================
+!
+!!======================================================================
+!!               ***  Observation & Assimilation  ***                 !!
+!!                                                                    !!
+!!   namobs       observation and model comparison                      (default: OFF)
+!!   nam_asminc   assimilation increments                               ('key_asminc')
+!!======================================================================
+!
+!!======================================================================
+!!                  ***  Miscellaneous namelists  ***                 !!
+!!                                                                    !!
+!!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi")
+!!   namctl            Control prints                                   (default: OFF)
+!!   namsto            Stochastic parametrization of EOS                (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints                                       (default: OFF)
+!-----------------------------------------------------------------------
+ln_ctl = .false.
+/
+!-----------------------------------------------------------------------
+&namsto        ! Stochastic parametrization of EOS                      (default: OFF)
+!-----------------------------------------------------------------------
+/
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/plot_station_asf.py
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/plot_station_asf.py	(revision 11962)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/plot_station_asf.py	(revision 11996)
@@ -5,5 +5,5 @@
 
 import sys
-#from os import path as path
+from os import path as path
 #from string import replace
 import math
@@ -26,5 +26,5 @@
 sys.setdefaultencoding('utf8')
 
-cy1     = '2016' ; # First year 
+cy1     = '2016' ; # First year
 cy2     = '2018' ; # Last year
 
@@ -40,15 +40,26 @@
 
 L_ALGOS = [ 'COARE3p6' , 'ECMWF'   , 'NCAR' ]
+l_xtrns = [ '-noskin'  , '-noskin' ,  ''    ] ; # string to add to algo name (L_ALGOS) to get version without skin params turned on
 l_color = [  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot
 l_width = [     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot
 l_style = [    '-'     ,   '-'     , '--'   ] ; # line-style
 
-L_VNEM  = [ 'qla'  , 'qsb'   , 'qt'   , 'qlw' , 'taum'  , 'dt_skin'  ]
-L_VARO  = [  'Qlat'    ,   'Qsen'    ,  'Qnet'    ,  'Qlw'    , 'Tau'   , 'dT_skin'  ] ; # name of variable on figure
-L_VARL  = [r'$Q_{lat}$',r'$Q_{sens}$',r'$Q_{net}$',r'$Q_{lw}$',r'$|\tau|$',r'$\Delta T_{skin}$' ] ; # name of variable in latex mode
-L_VUNT  = [  r'$W/m^2$'   ,  r'$W/m^2$'    ,  r'$W/m^2$'   ,  r'$W/m^2$'  , r'$N/m^2$' ,   'K'     ]
-L_VMAX  = [    75.     ,    75.      ,   800.     ,    25.    ,   1.2   ,   -0.7     ]
-L_VMIN  = [  -250.     ,  -125.      ,  -400.     ,  -150.    ,   0.    ,    0.7     ]
-L_ANOM  = [  True      ,   True      ,   True     ,   True    ,  True   ,  False     ]
+L_VNEM  = [   'qla'     ,     'qsb'     ,     'qt'     ,   'qlw'     ,  'taum'     ,    'dt_skin'         ]
+L_VARO  = [   'Qlat'    ,    'Qsen'     ,     'Qnet'   ,   'Qlw'     ,  'Tau'      ,    'dT_skin'         ] ; # name of variable on figure
+L_VARL  = [ r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$|\tau|$' , r'$\Delta T_{skin}$' ] ; # name of variable in latex mode
+L_VUNT  = [ r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$N/m^2$'  ,      'K'             ]
+L_VMAX  = [     75.     ,     75.       ,    800.      ,     25.     ,    1.2      ,      -0.7            ]
+L_VMIN  = [   -250.     ,   -125.       ,   -400.      ,   -150.     ,    0.       ,       0.7            ]
+L_ANOM  = [   True      ,    True       ,    True      ,    True     ,   True      ,      False           ]
+
+#L_VNEM  = [   'qlw'      ]
+#L_VARO  = [   'Qlw'      ] ; # name of variable on figure
+#L_VARL  = [ r'$Q_{lw}$'  ] ; # name of variable in latex mode
+#L_VUNT  = [ r'$W/m^2$'   ]
+#L_VMAX  = [     25.      ]
+#L_VMIN  = [   -150.      ]
+#L_ANOM  = [    True      ]
+
+
 
 nb_algos = len(L_ALGOS) ; print(nb_algos)
@@ -56,16 +67,33 @@
 # Getting arguments:
 narg = len(sys.argv)
-if narg != 3:
-    print 'Usage: '+sys.argv[0]+' <DIR_OUT_SASF> <year>'; sys.exit(0)
+if narg != 2:
+    print 'Usage: '+sys.argv[0]+' <DIR_OUT_SASF>'; sys.exit(0)
 cdir_data = sys.argv[1]
-cyyyy     = sys.argv[2]
-
-
-
-# Files to be read:
-cf_in = []
-for ja in range(nb_algos): cf_in.append(cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc')
+
+
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+# Populating and checking existence of files to be read
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+def chck4f(cf):
+    cmesg = 'ERROR: File '+cf+' does not exist !!!'
+    if not path.exists(cf): print cmesg ; sys.exit(0)
+
+###cf_in = nmp.empty((), dtype="S10")
+cf_in = [] ; cf_in_ns = []
+for ja in range(nb_algos):
+    cfi = cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc'
+    chck4f(cfi)
+    cf_in.append(cfi)
+# Same but without skin params:
+for ja in range(nb_algos):
+    cfi = cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+l_xtrns[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc'
+    chck4f(cfi)
+    cf_in_ns.append(cfi)
 print('Files we are goin to use:')
 for ja in range(nb_algos): print(cf_in[ja])
+print('   --- same without cool-skin/warm-layer:')
+for ja in range(nb_algos): print(cf_in_ns[ja])
+#-----------------------------------------------------------------
 
 
@@ -97,68 +125,134 @@
 
 
-
-
-
-for jv in range(nb_var):
-    print('\n *** Treating variable: '+L_VARO[jv]+' !')
-
-    for ja in range(nb_algos):
-        id_in = Dataset(cf_in[ja])
-        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
-        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
-        id_in.close()    
-        #print('    => '+L_ALGOS[ja]+'  => ok!')
-
-
-    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
-
-    ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
-
-    ax1.set_xticks(vtime[::xticks_d])
-    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
-    plt.xticks(rotation='60')
-
-    for ja in range(nb_algos):
-        plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
-
-    ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
-    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
-    
-    ax1.grid(color='k', linestyle='-', linewidth=0.3)
-    plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
-    ax1.annotate(cvar_lnm, xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
-    plt.savefig(L_VARO[jv]+'.'+fig_ext, dpi=int(rDPI), transparent=False)
-    plt.close(jv)
-
-
-
-    if L_ANOM[jv]:
-
-        for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
-
-        # Want a symetric y-range that makes sense for the anomaly we're looking at:
-        rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
-        rmax = max( abs(rmax) , abs(rmin) )
-        romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
-        rmlt = 10.**(int(romagn)) / 2.
-        yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax)
-        #print 'yrng = ', yrng ;  #sys.exit(0)
-        
-        fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
+for ctest in ['skin','noskin']:
+
+    for jv in range(nb_var):
+        print('\n *** Treating variable: '+L_VARO[jv]+' ('+ctest+') !')
+
+        for ja in range(nb_algos):
+            #
+            if ctest == 'skin':   id_in = Dataset(cf_in[ja])
+            if ctest == 'noskin': id_in = Dataset(cf_in_ns[ja])
+            xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
+            if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
+            id_in.close()
+
+        fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
+
         ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
-        
+
         ax1.set_xticks(vtime[::xticks_d])
         ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
         plt.xticks(rotation='60')
-        
+
         for ja in range(nb_algos):
-            plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
-
-        ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
+            plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
+
+        ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
         plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
+
         ax1.grid(color='k', linestyle='-', linewidth=0.3)
         plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
-        ax1.annotate('Anomaly of '+cvar_lnm, xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
-        plt.savefig(L_VARO[jv]+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False)
-        plt.close(10+jv)
-
+        ax1.annotate(cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
+        plt.savefig(L_VARO[jv]+'_'+ctest+'.'+fig_ext, dpi=int(rDPI), transparent=False)
+        plt.close(jv)
+
+
+
+        if L_ANOM[jv]:
+
+            for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
+
+            if nmp.sum(xFa[:,:]) == 0.0:
+                print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
+                print('          ==> skipping anomaly plot...')
+
+            else:
+
+                # Want a symetric y-range that makes sense for the anomaly we're looking at:
+                rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
+                rmax = max( abs(rmax) , abs(rmin) )
+                romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
+                rmlt = 10.**(int(romagn)) / 2.
+                yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax)
+                #print 'yrng = ', yrng ;  #sys.exit(0)
+
+                fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
+                ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
+
+                ax1.set_xticks(vtime[::xticks_d])
+                ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
+                plt.xticks(rotation='60')
+
+                for ja in range(nb_algos):
+                    plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
+
+                ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
+                plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
+                ax1.grid(color='k', linestyle='-', linewidth=0.3)
+                plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
+                ax1.annotate('Anomaly of '+cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
+                plt.savefig(L_VARO[jv]+'_'+ctest+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False)
+                plt.close(10+jv)
+
+
+
+
+# Difference skin vs noskin:
+xFns = nmp.zeros((nbr,nb_algos))
+
+for jv in range(nb_var-1):
+    print('\n *** Treating variable: '+L_VARO[jv]+' ('+ctest+') !')
+
+    for ja in range(nb_algos-1):
+        id_in = Dataset(cf_in[ja])
+        xF[:,ja]   = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
+        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
+        id_in.close()
+        #
+        id_in = Dataset(cf_in_ns[ja])
+        xFns[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
+        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
+        id_in.close()
+
+        xFa[:,ja] = xF[:,ja] - xFns[:,ja]  ; # difference!
+
+
+    # Want a symetric y-range that makes sense for the anomaly we're looking at:
+    rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
+    rmax = max( abs(rmax) , abs(rmin) )
+    romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
+    rmlt = 10.**(int(romagn)) / 2.
+    yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax)
+    print 'yrng = ', yrng ;  #sys.exit(0)
+
+
+
+
+    for ja in range(nb_algos-1):
+
+        calgo = L_ALGOS[ja]
+
+        if nmp.sum(xFa[:,ja]) == 0.0:
+            print('     Well! Seems that for variable '+L_VARO[jv]+', and algo '+calgo+', skin param has no impact')
+            print('          ==> skipping difference plot...')
+
+        else:
+
+            fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
+            ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
+
+            ax1.set_xticks(vtime[::xticks_d])
+            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
+            plt.xticks(rotation='60')
+
+            plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=None, zorder=10+ja)
+
+            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
+            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
+
+            ax1.grid(color='k', linestyle='-', linewidth=0.3)
+            #plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
+            ax1.annotate(cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
+            plt.savefig('diff_skin-noskin_'+L_VARO[jv]+'_'+calgo+'_'+ctest+'.'+fig_ext, dpi=int(rDPI), transparent=False)
+            plt.close(jv)
