Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11636)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11637)
@@ -153,7 +153,7 @@
       !!             ***  ROUTINE sbc_blk_cswl_alloc ***
       !!-------------------------------------------------------------------
-      !PRINT *, '*** LB: allocating tsk!'
+      !WRITE(numout,*) '*** LB: allocating tsk!'
       ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc )
-      !PRINT *, '*** LB: done!'
+      !WRITE(numout,*) '*** LB: done!'
       CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc )
       IF( sbc_blk_cswl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' )
@@ -510,8 +510,8 @@
          zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
       CASE( np_humi_dpt )
-         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !'
+         IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !'
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       CASE( np_humi_rlh )
-         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !'
+         IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !'
          zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       END SELECT
@@ -529,5 +529,5 @@
 
          CASE( np_COARE_3p0 )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day=', nsec_day
+            IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
             CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
@@ -536,5 +536,5 @@
 
          CASE( np_COARE_3p6 )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
+            IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
             CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
@@ -543,5 +543,5 @@
 
          CASE( np_ECMWF     )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day=', nsec_day
+            IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1)
             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
@@ -576,15 +576,15 @@
 
          CASE( np_COARE_3p0 )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!'
+            IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!'
             CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_COARE_3p6 )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!'
+            IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!'
             CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_ECMWF     )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
+            IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
@@ -854,8 +854,8 @@
          zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
       CASE( np_humi_dpt )
-         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !'
+         IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !'
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       CASE( np_humi_rlh )
-         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !'
+         IF (lwp) WRITE(numout,*) ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !'
          zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       END SELECT
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11636)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11637)
@@ -151,6 +151,4 @@
       !!---------------------------------------------------------------------
 
-      IF (lwp) WRITE(numout,*) 'LOLO: WL_ECMWF@sbcblk_skin_ecmwf.F90: depth for SST bulk takent at', gdept_1d(1), ' m'
-
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -164,6 +162,4 @@
             flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$
             dT_wl = pdT(ji,jj) / ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
-            !PRINT *, 'LOLO/mod_wl_ecmwf.f90: dT_wl2=', dT_wl
-            !PRINT *, ''
 
             zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/README.md
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/README.md	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/README.md	(revision 11637)
@@ -0,0 +1,3 @@
+## Welcome to ```STATION_ASF```, a NEMO setup which mimics an off-shore station dedicated to the measure of Air-Sea Fluxes (ASF).
+
+
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/diawri.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/diawri.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/diawri.F90	(revision 11637)
@@ -0,0 +1,383 @@
+MODULE diawri
+   !!======================================================================
+   !!                     ***  MODULE  diawri  ***
+   !! Ocean diagnostics :  write ocean output files
+   !!=====================================================================
+   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
+   !!            4.0  ! 1991-11  (G. Madec)
+   !!                 ! 1992-06  (M. Imbard)  correction restart file
+   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
+   !!                 ! 1993-03  (M. Imbard)  suppress writibm
+   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
+   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
+   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
+   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
+   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
+   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
+   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_wri       : create the standart output files
+   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE dianam         ! build name of file (routine)
+   !LB:USE diahth         ! thermocline diagnostics
+   USE zdf_oce        ! ocean vertical physics
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   !LB:USE sbc_ice        ! Surface boundary condition: ice fields
+   !LB:USE sbcssr         ! restoring term toward SST/SSS climatology
+   !LB:USE zdfmxl         ! mixed layer
+   !LB:USE zdfddm         ! vertical  physics: double diffusion
+   !
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager ! I/O manager
+   USE iom            ! 
+   USE ioipsl         ! 
+#if defined key_si3
+   USE ice 
+   USE icewri 
+#endif
+   USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_wri                 ! routines called by step.F90
+   PUBLIC   dia_wri_state
+   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
+
+   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
+   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
+   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
+   INTEGER ::   ndex(1)                              ! ???
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/SAS 4.0 , NEMO Consortium (2018)
+   !! $Id: diawri.F90 10425 2018-12-19 21:54:16Z smasson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_iomput
+   !!----------------------------------------------------------------------
+   !!   'key_iomput'                                        use IOM library
+   !!----------------------------------------------------------------------
+   INTEGER FUNCTION dia_wri_alloc()
+      !
+      dia_wri_alloc = 0
+      !
+   END FUNCTION dia_wri_alloc
+
+   
+   SUBROUTINE dia_wri( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_wri  ***
+      !!                   
+      !! ** Purpose :   Standard output of opa: dynamics and tracer fields 
+      !!      NETCDF format is used by default 
+      !!      Standalone surface scheme 
+      !!
+      !! ** Method  :  use iom_put
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!----------------------------------------------------------------------
+      ! 
+      ! Output the initial state and forcings
+      IF( ninist == 1 ) THEN                       
+         CALL dia_wri_state( 'output.init' )
+         ninist = 0
+      ENDIF
+      !
+   END SUBROUTINE dia_wri
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                  use IOIPSL  library
+   !!----------------------------------------------------------------------
+
+   INTEGER FUNCTION dia_wri_alloc()
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( ndex_hT(jpi*jpj), ndex_hU(jpi*jpj), ndex_hV(jpi*jpj), STAT=dia_wri_alloc )
+      CALL mpp_sum( 'diawri', dia_wri_alloc )
+      !
+   END FUNCTION dia_wri_alloc
+
+   
+   SUBROUTINE dia_wri( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_wri  ***
+      !!                   
+      !! ** Purpose :   Standard output of opa: dynamics and tracer fields 
+      !!      NETCDF format is used by default 
+      !!
+      !! ** Method  :   At the beginning of the first time step (nit000), 
+      !!      define all the NETCDF files and fields
+      !!      At each time step call histdef to compute the mean if ncessary
+      !!      Each nwrite time step, output the instantaneous or mean fields
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !
+      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
+      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
+      INTEGER  ::   inum = 11                                ! temporary logical unit
+      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
+      INTEGER  ::   ierr                                     ! error code return from allocation
+      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
+      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
+      !!----------------------------------------------------------------------
+      ! 
+      IF( ln_timing )   CALL timing_start('dia_wri')
+      !
+      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
+         CALL dia_wri_state( 'output.init' )
+         ninist = 0
+      ENDIF
+      !
+      ! 0. Initialisation
+      ! -----------------
+
+      ll_print = .FALSE.                  ! local variable for debugging
+      ll_print = ll_print .AND. lwp
+
+      ! Define frequency of output and means
+      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
+      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
+      ENDIF
+#if defined key_diainstant
+      zsto = nwrite * rdt
+      clop = "inst("//TRIM(clop)//")"
+#else
+      zsto=rdt
+      clop = "ave("//TRIM(clop)//")"
+#endif
+      zout = nwrite * rdt
+      zmax = ( nitend - nit000 + 1 ) * rdt
+
+      ! Define indices of the horizontal output zoom and vertical limit storage
+      iimi = 1      ;      iima = jpi
+      ijmi = 1      ;      ijma = jpj
+      ipk = jpk
+
+      ! define time axis
+      it = kt
+      itmod = kt - nit000 + 1
+
+
+      ! 1. Define NETCDF files and fields at beginning of first time step
+      ! -----------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+
+         ! Define the NETCDF files (one per grid)
+
+         ! Compute julian date from starting date of the run
+         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
+         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
+         IF(lwp)WRITE(numout,*)
+         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
+            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
+         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
+                                 ' limit storage in depth = ', ipk
+
+         ! WRITE root name in date.file for use by postpro
+         IF(lwp) THEN
+            CALL dia_nam( clhstnam, nwrite,' ' )
+            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            WRITE(inum,*) clhstnam
+            CLOSE(inum)
+         ENDIF
+
+         ! Define the T grid FILE ( nid_T )
+
+         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
+         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
+            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
+         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
+            &           "m", ipk, gdept_1d, nz_T, "down" )
+         !                                                            ! Index of ocean points
+         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
+
+         ! Define the U grid FILE ( nid_U )
+
+         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
+         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
+            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
+         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
+            &           "m", ipk, gdept_1d, nz_U, "down" )
+         !                                                            ! Index of ocean points
+         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
+
+         ! Define the V grid FILE ( nid_V )
+
+         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
+         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
+         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
+            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
+         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
+            &          "m", ipk, gdept_1d, nz_V, "down" )
+         !                                                            ! Index of ocean points
+         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
+
+         ! No W grid FILE
+
+         ! Declare all the output fields as NETCDF variables
+
+         !                                                                                      !!! nid_T : 3D
+         CALL histdef( nid_T, "sst_m", "Sea Surface temperature"            , "C"      ,   &  ! sst
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sss_m", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! (sfx)
+             &         jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
+            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+
+         CALL histend( nid_T, snc4chunks=snc4set )
+
+         !                                                                                      !!! nid_U : 3D
+         CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s"   ,         &  ! ssu
+            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
+            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+
+         CALL histend( nid_U, snc4chunks=snc4set )
+
+         !                                                                                      !!! nid_V : 3D
+         CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s",            &  ! ssv_m
+            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
+            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+
+         CALL histend( nid_V, snc4chunks=snc4set )
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
+         IF(ll_print) CALL FLUSH(numout )
+
+      ENDIF
+
+      ! 2. Start writing data
+      ! ---------------------
+
+      ! ndex(1) est utilise ssi l'avant dernier argument est different de 
+      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
+      ! donne le nombre d'elements, et ndex la liste des indices a sortir
+
+      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN 
+         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
+         WRITE(numout,*) '~~~~~~ '
+      ENDIF
+
+      ! Write fields on T grid
+      CALL histwrite( nid_T, "sst_m", it, sst_m, ndim_hT, ndex_hT )   ! sea surface temperature
+      CALL histwrite( nid_T, "sss_m", it, sss_m, ndim_hT, ndex_hT )   ! sea surface salinity
+      CALL histwrite( nid_T, "sowaflup", it, (emp - rnf )  , ndim_hT, ndex_hT )   ! upward water flux
+      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux 
+                                                                                  ! (includes virtual salt flux beneath ice 
+                                                                                  ! in linear free surface case)
+
+      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
+      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
+      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
+      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
+
+         ! Write fields on U grid
+      CALL histwrite( nid_U, "ssu_m"   , it, ssu_m         , ndim_hU, ndex_hU )   ! i-current speed
+      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
+
+         ! Write fields on V grid
+      CALL histwrite( nid_V, "ssv_m"   , it, ssv_m         , ndim_hV, ndex_hV )   ! j-current speed
+      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
+
+      ! 3. Close all files
+      ! ---------------------------------------
+      IF( kt == nitend ) THEN
+         CALL histclo( nid_T )
+         CALL histclo( nid_U )
+         CALL histclo( nid_V )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('dia_wri')
+      !
+   END SUBROUTINE dia_wri
+#endif
+
+   SUBROUTINE dia_wri_state( cdfile_name )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_wri_state  ***
+      !!        
+      !! ** Purpose :   create a NetCDF file named cdfile_name which contains 
+      !!      the instantaneous ocean state and forcing fields.
+      !!        Used to find errors in the initial state or save the last
+      !!      ocean state in case of abnormal end of a simulation
+      !!
+      !! ** Method  :   NetCDF files using ioipsl
+      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
+      !!      File 'output.abort.nc' is created in case of abnormal job end
+      !!----------------------------------------------------------------------
+      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
+      !!
+      INTEGER :: inum
+      !!----------------------------------------------------------------------
+      ! 
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
+      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
+
+#if defined key_si3
+     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
+#else
+     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
+#endif
+
+      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
+      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
+      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
+      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
+      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
+      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
+      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
+      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
+      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
+      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
+      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
+      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
+ 
+#if defined key_si3
+      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
+         CALL ice_wri_state( inum )
+      ENDIF
+#endif
+      !
+      CALL iom_close( inum )
+      ! 
+   END SUBROUTINE dia_wri_state
+
+   !!======================================================================
+END MODULE diawri
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/nemogcm.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/nemogcm.F90	(revision 11637)
@@ -0,0 +1,514 @@
+MODULE nemogcm
+   !!======================================================================
+   !!                       ***  MODULE nemogcm   ***
+   !! StandAlone Surface module : surface fluxes
+   !!======================================================================
+   !! History :  3.6  ! 2011-11  (S. Alderson, G. Madec) original code
+   !!             -   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) nemo_northcomms: setup avoiding MPI communication 
+   !!             -   ! 2014-12  (G. Madec) remove KPP scheme and cross-land advection (cla)
+   !!            4.0  ! 2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   nemo_gcm      : solve ocean dynamics, tracer, biogeochemistry
+   !!   nemo_init     : initialization of the NEMO system
+   !!   nemo_ctl      : initialisation of the contol print
+   !!   nemo_closefile: close remaining open files
+   !!   nemo_alloc    : dynamical allocation
+   !!----------------------------------------------------------------------
+   USE step_oce       ! module used in the ocean time stepping module
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constant                  (par_cst routine)
+   USE domain         ! domain initialization   (dom_init & dom_cfg routines)
+   USE closea         ! treatment of closed seas (for ln_closea)
+   USE usrdef_nam     ! user defined configuration
+   USE daymod         ! calendar
+   USE restart        ! open  restart file
+   !LB:USE step           ! NEMO time-stepping                 (stp     routine)
+   !LB:USE cpl_oasis3     !
+   USE c1d            ! 1D configuration
+   USE step_c1d       ! Time stepping loop for the 1D configuration
+   USE sbcssm         !
+   !LB:USE icbini         ! handle bergs, initialisation
+   !LB:USE icbstp         ! handle bergs, calving, themodynamics and transport
+   !LB:USE bdyini         ! open boundary cond. setting       (bdy_init routine). mandatory for sea-ice
+   !LB:USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine). mandatory for sea-ice
+   !
+   USE lib_mpp        ! distributed memory computing
+   USE mppini         ! shared/distributed memory setting (mpp_init routine)
+   !LB:USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges 
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+#if defined key_iomput
+   USE xios           ! xIOserver
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   nemo_gcm    ! called by model.F90
+   PUBLIC   nemo_init   ! needed by AGRIF
+
+   CHARACTER(lc) ::   cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
+
+   !!----------------------------------------------------------------------
+   !! NEMO/SAS 4.0 , NEMO Consortium (2018)
+   !! $Id: nemogcm.F90 10601 2019-01-28 20:53:01Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE nemo_gcm
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_gcm  ***
+      !!
+      !! ** Purpose :   NEMO solves the primitive equations on an orthogonal
+      !!              curvilinear mesh on the sphere.
+      !!
+      !! ** Method  : - model general initialization
+      !!              - launch the time-stepping (stp routine)
+      !!              - finalize the run by closing files and communications
+      !!
+      !! References : Madec, Delecluse, Imbard, and Levy, 1997:  internal report, IPSL.
+      !!              Madec, 2008, internal report, IPSL.
+      !!----------------------------------------------------------------------
+      INTEGER ::   istp   ! time step index
+      !!----------------------------------------------------------------------
+      !
+      !                            !-----------------------!
+      CALL nemo_init               !==  Initialisations  ==!
+      !                            !-----------------------!
+      ! check that all process are still there... If some process have an error,
+      ! they will never enter in step and other processes will wait until the end of the cpu time!
+      CALL mpp_max( 'nemogcm', nstop )
+
+      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
+
+      !                            !-----------------------!
+      !                            !==   time stepping   ==!
+      !                            !-----------------------!
+      istp = nit000
+      !
+#if defined key_c1d
+      DO WHILE ( istp <= nitend .AND. nstop == 0 )    !==  C1D time-stepping  ==!
+         CALL stp_c1d( istp )
+         istp = istp + 1
+      END DO
+#else
+      !
+      IF( .NOT.ln_diurnal_only ) THEN                 !==  Standard time-stepping  ==!
+         !
+         DO WHILE( istp <= nitend .AND. nstop == 0 )
+#if defined key_mpp_mpi
+            ncom_stp = istp
+            IF ( istp == ( nit000 + 1 ) ) elapsed_time = MPI_Wtime()
+            IF ( istp ==         nitend ) elapsed_time = MPI_Wtime() - elapsed_time
+#endif
+            CALL stp        ( istp ) 
+            istp = istp + 1
+         END DO
+         !
+      ELSE                                            !==  diurnal SST time-steeping only  ==!
+         !
+         DO WHILE( istp <= nitend .AND. nstop == 0 )
+            CALL stp_diurnal( istp )   ! time step only the diurnal SST 
+            istp = istp + 1
+         END DO
+         !
+      ENDIF
+      !
+#endif
+      !
+      !                            !------------------------!
+      !                            !==  finalize the run  ==!
+      !                            !------------------------!
+      IF(lwp) WRITE(numout,cform_aaa)        ! Flag AAAAAAA
+      !
+      IF( nstop /= 0 .AND. lwp ) THEN        ! error print
+         WRITE(numout,cform_err)
+         WRITE(numout,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found'
+         WRITE(numout,*)
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_finalize
+      !
+      CALL nemo_closefile
+      !
+#if defined key_iomput
+                                    CALL xios_finalize  ! end mpp communications with xios
+#else
+      IF( lk_mpp   ) THEN   ;   CALL mppstop( ldfinal = .TRUE. )   ! end mpp communications
+      ENDIF
+#endif
+      !
+      IF(lwm) THEN
+         IF( nstop == 0 ) THEN   ;   STOP 0
+         ELSE                    ;   STOP 999
+         ENDIF
+      ENDIF
+      !
+   END SUBROUTINE nemo_gcm
+
+
+   SUBROUTINE nemo_init
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_init  ***
+      !!
+      !! ** Purpose :   initialization of the NEMO GCM
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji                 ! dummy loop indices
+      INTEGER  ::   ios, ilocal_comm   ! local integers
+      CHARACTER(len=120), DIMENSION(60) ::   cltxt, cltxt2, clnam
+      !!
+      NAMELIST/namctl/ ln_ctl   , sn_cfctl, nn_print, nn_ictls, nn_ictle,   &
+         &             nn_isplt , nn_jsplt, nn_jctls, nn_jctle,             &
+         &             ln_timing, ln_diacfl
+      NAMELIST/namcfg/ ln_read_cfg, cn_domcfg, ln_closea, ln_write_cfg, cn_domcfg_out, ln_use_jattr
+      !!----------------------------------------------------------------------
+      !
+      cltxt  = ''
+      cltxt2 = ''
+      clnam  = ''  
+      cxios_context = 'nemo'
+      !
+      !                             ! Open reference namelist and configuration namelist files
+      CALL ctl_opn( numnam_ref, 'namelist_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+      CALL ctl_opn( numnam_cfg, 'namelist_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+      !
+      REWIND( numnam_ref )              ! Namelist namctl in reference namelist
+      READ  ( numnam_ref, namctl, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namctl in reference namelist', .TRUE. )
+      REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist
+      READ  ( numnam_cfg, namctl, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namctl in configuration namelist', .TRUE. )
+      !
+      REWIND( numnam_ref )              ! Namelist namcfg in reference namelist
+      READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 )
+903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. )
+      REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist
+      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 )
+904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )   
+
+      !                             !--------------------------!
+      !                             !  Set global domain size  !   (control print return in cltxt2)
+      !                             !--------------------------!
+      IF( ln_read_cfg ) THEN              ! Read sizes in domain configuration file
+         CALL domain_cfg ( cltxt2,        cn_cfg, nn_cfg, jpiglo, jpjglo, jpkglo, jperio )
+         !
+      ELSE                                ! user-defined namelist
+         CALL usr_def_nam( cltxt2, clnam, cn_cfg, nn_cfg, jpiglo, jpjglo, jpkglo, jperio )
+      ENDIF
+      !
+      !
+      !                             !--------------------------------------------!
+      !                             !  set communicator & select the local node  !
+      !                             !  NB: mynode also opens output.namelist.dyn !
+      !                             !      on unit number numond on first proc   !
+      !                             !--------------------------------------------!
+#if defined key_iomput
+      IF( Agrif_Root() ) THEN
+            CALL xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )    ! nemo local communicator given by xios
+      ENDIF
+      ! Nodes selection (control print return in cltxt)
+      narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
+#else
+         ilocal_comm = 0                                    ! Nodes selection (control print return in cltxt)
+         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop )
+#endif
+
+      narea = narea + 1                                     ! mynode return the rank of proc (0 --> jpnij -1 )
+
+      IF( sn_cfctl%l_config ) THEN
+         ! Activate finer control of report outputs
+         ! optionally switch off output from selected areas (note this only
+         ! applies to output which does not involve global communications)
+         IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax  ) .OR. &
+           & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) )    &
+           &   CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. )
+      ELSE
+         ! Use ln_ctl to turn on or off all options.
+         CALL nemo_set_cfctl( sn_cfctl, ln_ctl, .TRUE. )
+      ENDIF
+
+      lwm = (narea == 1)                                    ! control of output namelists
+      lwp = (narea == 1) .OR. ln_ctl                        ! control of all listing output print
+
+      IF(lwm) THEN               ! write merged namelists from earlier to output namelist 
+         !                       ! now that the file has been opened in call to mynode. 
+         !                       ! NB: nammpp has already been written in mynode (if lk_mpp_mpi)
+         WRITE( numond, namctl )
+         WRITE( numond, namcfg )
+         IF( .NOT.ln_read_cfg ) THEN
+            DO ji = 1, SIZE(clnam)
+               IF( TRIM(clnam(ji)) /= '' )   WRITE(numond, * ) clnam(ji)     ! namusr_def print
+            END DO
+         ENDIF
+      ENDIF
+
+      IF(lwp) THEN                            ! open listing units
+         !
+         CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
+         !
+         WRITE(numout,*)
+         WRITE(numout,*) '   CNRS - NERC - Met OFFICE - MERCATOR-ocean - INGV - CMCC'
+         WRITE(numout,*) '                       NEMO team'
+         WRITE(numout,*) '            Ocean General Circulation Model'
+         WRITE(numout,*) '                NEMO version 4.0  (2019) '
+         WRITE(numout,*) '                        SASC1D '
+         WRITE(numout,*) "           ._      ._      ._      ._      ._    "
+         WRITE(numout,*) "       _.-._)`\_.-._)`\_.-._)`\_.-._)`\_.-._)`\_ "
+         WRITE(numout,*)
+         WRITE(numout,*) "           o         _,           _,             "
+         WRITE(numout,*) "            o      .' (        .-' /             "
+         WRITE(numout,*) "           o     _/..._'.    .'   /              "
+         WRITE(numout,*) "      (    o .-'`      ` '-./  _.'               "
+         WRITE(numout,*) "       )    ( o)           ;= <_         (       "
+         WRITE(numout,*) "      (      '-.,\\__ __.-;`\   '.        )      "
+         WRITE(numout,*) "       )  )       \) |`\ \)  '.   \      (   (   "
+         WRITE(numout,*) "      (  (           \_/       '-._\      )   )  "
+         WRITE(numout,*) "       )  ) jgs                     `    (   (   "
+         WRITE(numout,*) "     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ "
+         WRITE(numout,*)
+         
+         DO ji = 1, SIZE(cltxt)
+            IF( TRIM(cltxt (ji)) /= '' )   WRITE(numout,*) TRIM(cltxt(ji))    ! control print of mynode
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*)
+         DO ji = 1, SIZE(cltxt2)
+            IF( TRIM(cltxt2(ji)) /= '' )   WRITE(numout,*) TRIM(cltxt2(ji))   ! control print of domain size
+         END DO
+         !
+         WRITE(numout,cform_aaa)                                        ! Flag AAAAAAA
+         !
+      ENDIF
+      ! open /dev/null file to be able to supress output write easily
+      CALL ctl_opn( numnul, '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+      !
+      !                                      ! Domain decomposition
+      CALL mpp_init                          ! MPP
+
+      ! Now we know the dimensions of the grid and numout has been set: we can allocate arrays
+      CALL nemo_alloc()
+
+      !                             !-------------------------------!
+      !                             !  NEMO general initialization  !
+      !                             !-------------------------------!
+
+      CALL nemo_ctl                          ! Control prints
+      !
+      !                                      ! General initialization
+      IF( ln_timing    )   CALL timing_init     ! timing
+      IF( ln_timing    )   CALL timing_start( 'nemo_init')
+      !
+                           CALL     phy_cst         ! Physical constants
+                           CALL     eos_init        ! Equation of state
+!LB:                           CALL dom_init('SAS') ! Domain
+      IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration !LB
+                           CALL     dom_init("OPA") ! Domain
+      IF( ln_ctl       )   CALL prt_ctl_init        ! Print control
+      !PRINT *, 'LOLO/nemogcm.F90: after prt_ctl_init!';  STOP
+      
+                           CALL day_init        ! model calendar (using both namelist and restart infos)
+      IF( ln_rstart )      CALL rst_read_open
+
+      !                                      ! external forcing 
+                           CALL sbc_init        ! Forcings : surface module 
+
+      ! ==> clem: open boundaries init. is mandatory for sea-ice because ice BDY is not decoupled from  
+      !           the environment of ocean BDY. Therefore bdy is called in both OPA and SAS modules. 
+      !           This is not clean and should be changed in the future. 
+                           !CALL bdy_init
+      !
+      IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA
+      !
+      IF( ln_timing    )   CALL timing_stop( 'nemo_init')
+      !
+   END SUBROUTINE nemo_init
+
+
+   SUBROUTINE nemo_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_ctl  ***
+      !!
+      !! ** Purpose :   control print setting
+      !!
+      !! ** Method  : - print namctl information and check some consistencies
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'nemo_ctl: Control prints'
+         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) '   Namelist namctl'
+         WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl
+         WRITE(numout,*) '       finer control over o/p sn_cfctl%l_config  = ', sn_cfctl%l_config
+         WRITE(numout,*) '                              sn_cfctl%l_runstat = ', sn_cfctl%l_runstat
+         WRITE(numout,*) '                              sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat
+         WRITE(numout,*) '                              sn_cfctl%l_oceout  = ', sn_cfctl%l_oceout
+         WRITE(numout,*) '                              sn_cfctl%l_layout  = ', sn_cfctl%l_layout
+         WRITE(numout,*) '                              sn_cfctl%l_mppout  = ', sn_cfctl%l_mppout
+         WRITE(numout,*) '                              sn_cfctl%l_mpptop  = ', sn_cfctl%l_mpptop
+         WRITE(numout,*) '                              sn_cfctl%procmin   = ', sn_cfctl%procmin  
+         WRITE(numout,*) '                              sn_cfctl%procmax   = ', sn_cfctl%procmax  
+         WRITE(numout,*) '                              sn_cfctl%procincr  = ', sn_cfctl%procincr 
+         WRITE(numout,*) '                              sn_cfctl%ptimincr  = ', sn_cfctl%ptimincr 
+         WRITE(numout,*) '      level of print                  nn_print   = ', nn_print
+         WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls
+         WRITE(numout,*) '      End i indice for SUM control    nn_ictle   = ', nn_ictle
+         WRITE(numout,*) '      Start j indice for SUM control  nn_jctls   = ', nn_jctls
+         WRITE(numout,*) '      End j indice for SUM control    nn_jctle   = ', nn_jctle
+         WRITE(numout,*) '      number of proc. following i     nn_isplt   = ', nn_isplt
+         WRITE(numout,*) '      number of proc. following j     nn_jsplt   = ', nn_jsplt
+         WRITE(numout,*) '      timing by routine               ln_timing  = ', ln_timing
+         WRITE(numout,*) '      CFL diagnostics                 ln_diacfl  = ', ln_diacfl
+      ENDIF
+      !
+      nprint    = nn_print          ! convert DOCTOR namelist names into OLD names
+      nictls    = nn_ictls
+      nictle    = nn_ictle
+      njctls    = nn_jctls
+      njctle    = nn_jctle
+      isplt     = nn_isplt
+      jsplt     = nn_jsplt
+
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namcfg'
+         WRITE(numout,*) '      read domain configuration file                ln_read_cfg      = ', ln_read_cfg
+         WRITE(numout,*) '         filename to be read                           cn_domcfg     = ', TRIM(cn_domcfg)
+         WRITE(numout,*) '         keep closed seas in the domain (if exist)     ln_closea     = ', ln_closea
+         WRITE(numout,*) '      create a configuration definition file        ln_write_cfg     = ', ln_write_cfg
+         WRITE(numout,*) '         filename to be written                        cn_domcfg_out = ', TRIM(cn_domcfg_out)
+         WRITE(numout,*) '      use file attribute if exists as i/p j-start   ln_use_jattr     = ', ln_use_jattr
+      ENDIF
+      IF( .NOT.ln_read_cfg )   ln_closea = .false.   ! dealing possible only with a domcfg file
+      !
+      !                             ! Parameter control
+      !
+      IF( ln_ctl ) THEN                 ! sub-domain area indices for the control prints
+         IF( lk_mpp .AND. jpnij > 1 ) THEN
+            isplt = jpni   ;   jsplt = jpnj   ;   ijsplt = jpni*jpnj   ! the domain is forced to the real split domain
+         ELSE
+            IF( isplt == 1 .AND. jsplt == 1  ) THEN
+               CALL ctl_warn( ' - isplt & jsplt are equal to 1',   &
+                  &           ' - the print control will be done over the whole domain' )
+            ENDIF
+            ijsplt = isplt * jsplt            ! total number of processors ijsplt
+         ENDIF
+         IF(lwp) WRITE(numout,*)'          - The total number of processors over which the'
+         IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt
+         !
+         !                              ! indices used for the SUM control
+         IF( nictls+nictle+njctls+njctle == 0 )   THEN    ! print control done over the default area
+            lsp_area = .FALSE.
+         ELSE                                             ! print control done over a specific  area
+            lsp_area = .TRUE.
+            IF( nictls < 1 .OR. nictls > jpiglo )   THEN
+               CALL ctl_warn( '          - nictls must be 1<=nictls>=jpiglo, it is forced to 1' )
+               nictls = 1
+            ENDIF
+            IF( nictle < 1 .OR. nictle > jpiglo )   THEN
+               CALL ctl_warn( '          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo' )
+               nictle = jpiglo
+            ENDIF
+            IF( njctls < 1 .OR. njctls > jpjglo )   THEN
+               CALL ctl_warn( '          - njctls must be 1<=njctls>=jpjglo, it is forced to 1' )
+               njctls = 1
+            ENDIF
+            IF( njctle < 1 .OR. njctle > jpjglo )   THEN
+               CALL ctl_warn( '          - njctle must be 1<=njctle>=jpjglo, it is forced to jpjglo' )
+               njctle = jpjglo
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      IF( 1._wp /= SIGN(1._wp,-0._wp)  )   CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows f2003 standard.',  &
+         &                                                'Compile with key_nosignedzero enabled:',   &
+         &                                                '--> add -Dkey_nosignedzero to the definition of %CPP in your arch file' )
+      !
+   END SUBROUTINE nemo_ctl
+
+
+   SUBROUTINE nemo_closefile
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_closefile  ***
+      !!
+      !! ** Purpose :   Close the files
+      !!----------------------------------------------------------------------
+      !
+      IF( lk_mpp )   CALL mppsync
+      !
+      CALL iom_close                                 ! close all input/output files managed by iom_*
+      !
+      IF( numstp          /= -1 )   CLOSE( numstp          )   ! time-step file
+      IF( numrun          /= -1 )   CLOSE( numrun          )   ! run statistics file
+      IF( numnam_ref      /= -1 )   CLOSE( numnam_ref      )   ! oce reference namelist
+      IF( numnam_cfg      /= -1 )   CLOSE( numnam_cfg      )   ! oce configuration namelist
+      IF( lwm.AND.numond  /= -1 )   CLOSE( numond          )   ! oce output namelist
+      IF( numnam_ice_ref  /= -1 )   CLOSE( numnam_ice_ref  )   ! ice reference namelist
+      IF( numnam_ice_cfg  /= -1 )   CLOSE( numnam_ice_cfg  )   ! ice configuration namelist
+      IF( lwm.AND.numoni  /= -1 )   CLOSE( numoni          )   ! ice output namelist
+      IF( numevo_ice      /= -1 )   CLOSE( numevo_ice      )   ! ice variables (temp. evolution)
+      IF( numout          /=  6 )   CLOSE( numout          )   ! standard model output file
+      !
+      numout = 6                                     ! redefine numout in case it is used after this point...
+      !
+   END SUBROUTINE nemo_closefile
+
+
+   SUBROUTINE nemo_alloc
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_alloc  ***
+      !!
+      !! ** Purpose :   Allocate all the dynamic arrays of the OPA modules
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      USE diawri    , ONLY : dia_wri_alloc
+      USE dom_oce   , ONLY : dom_oce_alloc
+      !USE bdy_oce   , ONLY : ln_bdy, bdy_oce_alloc
+      USE oce       ! mandatory for sea-ice because needed for bdy arrays
+      !
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ierr =        dia_wri_alloc()
+      ierr = ierr + dom_oce_alloc()          ! ocean domain
+      ierr = ierr + oce_alloc    ()          ! (tsn...) needed for agrif and/or SI3 and bdy
+      !ierr = ierr + bdy_oce_alloc()          ! bdy masks (incl. initialization)
+      !
+      CALL mpp_sum( 'nemogcm', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc: unable to allocate standard ocean arrays' )
+      !
+   END SUBROUTINE nemo_alloc
+
+   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_set_cfctl  ***
+      !!
+      !! ** Purpose :   Set elements of the output control structure to setto.
+      !!                for_all should be .false. unless all areas are to be
+      !!                treated identically.
+      !!
+      !! ** Method  :   Note this routine can be used to switch on/off some
+      !!                types of output for selected areas but any output types
+      !!                that involve global communications (e.g. mpp_max, glob_sum)
+      !!                should be protected from selective switching by the
+      !!                for_all argument
+      !!----------------------------------------------------------------------
+      LOGICAL :: setto, for_all
+      TYPE(sn_ctl) :: sn_cfctl
+      !!----------------------------------------------------------------------
+      IF( for_all ) THEN
+         sn_cfctl%l_runstat = setto
+         sn_cfctl%l_trcstat = setto
+      ENDIF
+      sn_cfctl%l_oceout  = setto
+      sn_cfctl%l_layout  = setto
+      sn_cfctl%l_mppout  = setto
+      sn_cfctl%l_mpptop  = setto
+   END SUBROUTINE nemo_set_cfctl
+
+   !!======================================================================
+END MODULE nemogcm
+
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/sbcssm.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/sbcssm.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/sbcssm.F90	(revision 11637)
@@ -0,0 +1,315 @@
+MODULE sbcssm
+   !!======================================================================
+   !!                       ***  MODULE  sbcssm  ***
+   !! Off-line : interpolation of the physical fields
+   !!======================================================================
+   !! History :  3.4  ! 2012-03 (S. Alderson)  original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_ssm_init  : initialization, namelist read, and SAVEs control
+   !!   sbc_ssm       : Interpolation of the fields
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE c1d            ! 1D configuration: lk_c1d
+   USE dom_oce        ! ocean domain: variables
+   !LB:USE zdf_oce        ! ocean vertical physics: variables
+   USE sbc_oce        ! surface module: variables
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state - Brunt Vaisala frequency
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   !LB:USE zpshde         ! z-coord. with partial steps: horizontal derivatives
+   !LB:USE closea         ! for ln_closea
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE lib_mpp        ! distributed memory computing library
+   USE prtctl         ! print control
+   USE fldread        ! read input fields 
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_ssm        ! routine called by step.F90
+   PUBLIC   sbc_ssm_init   ! routine called by sbcmod.F90
+
+   CHARACTER(len=100) ::   cn_dir        ! Root directory for location of ssm files
+   LOGICAL            ::   ln_3d_uve     ! specify whether input velocity data is 3D
+   LOGICAL            ::   ln_read_frq   ! specify whether we must read frq or not
+   
+   LOGICAL            ::   l_sasread     ! Ice intilisation: =T read a file ; =F anaytical initilaistion
+   LOGICAL            ::   l_initdone = .false.
+   INTEGER     ::   nfld_3d
+   INTEGER     ::   nfld_2d
+
+   INTEGER     ::   jf_tem         ! index of temperature
+   INTEGER     ::   jf_sal         ! index of salinity
+   INTEGER     ::   jf_usp         ! index of u velocity component
+   INTEGER     ::   jf_vsp         ! index of v velocity component
+   INTEGER     ::   jf_ssh         ! index of sea surface height
+   INTEGER     ::   jf_e3t         ! index of first T level thickness
+   INTEGER     ::   jf_frq         ! index of fraction of qsr absorbed in the 1st T level
+
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_ssm_3d  ! structure of input fields (file information, fields read)
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_ssm_2d  ! structure of input fields (file information, fields read)
+
+   !!----------------------------------------------------------------------
+   !! NEMO/SAS 4.0 , NEMO Consortium (2018)
+   !! $Id: sbcssm.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_ssm( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_ssm  ***
+      !!
+      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
+      !!               for an off-line simulation using surface processes only
+      !!
+      !! ** Method : calculates the position of data 
+      !!             - interpolates data if needed
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
+      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'sbc_ssm')
+     
+      IF ( l_sasread ) THEN
+         IF( nfld_3d > 0 ) CALL fld_read( kt, 1, sf_ssm_3d )      !==   read data at kt time step   ==!
+         IF( nfld_2d > 0 ) CALL fld_read( kt, 1, sf_ssm_2d )      !==   read data at kt time step   ==!
+         ! 
+         IF( ln_3d_uve ) THEN
+            IF( .NOT. ln_linssh ) THEN
+               e3t_m(:,:) = sf_ssm_3d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor 
+            ELSE
+               e3t_m(:,:) = e3t_0(:,:,1)                                 ! vertical scale factor
+            ENDIF
+            ssu_m(:,:) = sf_ssm_3d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
+            ssv_m(:,:) = sf_ssm_3d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
+         ELSE
+            IF( .NOT. ln_linssh ) THEN
+               e3t_m(:,:) = sf_ssm_2d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor 
+            ELSE
+               e3t_m(:,:) = e3t_0(:,:,1)                                 ! vertical scale factor
+            ENDIF
+            ssu_m(:,:) = sf_ssm_2d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
+            ssv_m(:,:) = sf_ssm_2d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
+         ENDIF
+         !
+         sst_m(:,:) = sf_ssm_2d(jf_tem)%fnow(:,:,1) * tmask(:,:,1)    ! temperature
+         sss_m(:,:) = sf_ssm_2d(jf_sal)%fnow(:,:,1) * tmask(:,:,1)    ! salinity
+         ssh_m(:,:) = sf_ssm_2d(jf_ssh)%fnow(:,:,1) * tmask(:,:,1)    ! sea surface height
+         IF( ln_read_frq ) THEN
+            frq_m(:,:) = sf_ssm_2d(jf_frq)%fnow(:,:,1) * tmask(:,:,1) ! solar penetration
+         ELSE
+            frq_m(:,:) = 1._wp
+         ENDIF
+      ELSE
+         sss_m(:,:) = 35._wp                             ! =35. to obtain a physical value for the freezing point
+         CALL eos_fzp( sss_m(:,:), sst_m(:,:) )          ! sst_m is set at the freezing point
+         ssu_m(:,:) = 0._wp
+         ssv_m(:,:) = 0._wp
+         ssh_m(:,:) = 0._wp
+         IF( .NOT. ln_linssh ) e3t_m(:,:) = e3t_0(:,:,1) !clem: necessary at least for sas2D
+         frq_m(:,:) = 1._wp                              !              - -
+         sshn (:,:) = 0._wp                              !              - -
+      ENDIF
+      
+      IF ( nn_ice == 1 ) THEN
+         tsn(:,:,1,jp_tem) = sst_m(:,:)
+         tsn(:,:,1,jp_sal) = sss_m(:,:)
+         tsb(:,:,1,jp_tem) = sst_m(:,:)
+         tsb(:,:,1,jp_sal) = sss_m(:,:)
+      ENDIF
+      ub (:,:,1) = ssu_m(:,:)
+      vb (:,:,1) = ssv_m(:,:)
+ 
+      IF(ln_ctl) THEN                  ! print control
+         CALL prt_ctl(tab2d_1=sst_m, clinfo1=' sst_m   - : ', mask1=tmask   )
+         CALL prt_ctl(tab2d_1=sss_m, clinfo1=' sss_m   - : ', mask1=tmask   )
+         CALL prt_ctl(tab2d_1=ssu_m, clinfo1=' ssu_m   - : ', mask1=umask   )
+         CALL prt_ctl(tab2d_1=ssv_m, clinfo1=' ssv_m   - : ', mask1=vmask   )
+         CALL prt_ctl(tab2d_1=ssh_m, clinfo1=' ssh_m   - : ', mask1=tmask   )
+         IF( .NOT.ln_linssh )   CALL prt_ctl(tab2d_1=ssh_m, clinfo1=' e3t_m   - : ', mask1=tmask   )
+         IF( ln_read_frq    )   CALL prt_ctl(tab2d_1=frq_m, clinfo1=' frq_m   - : ', mask1=tmask   )
+      ENDIF
+      !
+      IF( l_initdone ) THEN          !   Mean value at each nn_fsbc time-step   !
+         CALL iom_put( 'ssu_m', ssu_m )
+         CALL iom_put( 'ssv_m', ssv_m )
+         CALL iom_put( 'sst_m', sst_m )
+         CALL iom_put( 'sss_m', sss_m )
+         CALL iom_put( 'ssh_m', ssh_m )
+         IF( .NOT.ln_linssh )   CALL iom_put( 'e3t_m', e3t_m )
+         IF( ln_read_frq    )   CALL iom_put( 'frq_m', frq_m )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop( 'sbc_ssm')
+      !
+   END SUBROUTINE sbc_ssm
+
+
+   SUBROUTINE sbc_ssm_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_ssm_init  ***
+      !!
+      !! ** Purpose :   Initialisation of sea surface mean data     
+      !!----------------------------------------------------------------------
+      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
+      INTEGER  :: ifpr                               ! dummy loop indice
+      INTEGER  :: inum, idv, idimv, jpm              ! local integer
+      INTEGER  ::   ios                              ! Local integer output status for namelist read
+      !!
+      CHARACTER(len=100)                     ::  cn_dir       ! Root directory for location of core files
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::  slf_3d       ! array of namelist information on the fields to read
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::  slf_2d       ! array of namelist information on the fields to read
+      TYPE(FLD_N) ::   sn_tem, sn_sal                     ! information about the fields to be read
+      TYPE(FLD_N) ::   sn_usp, sn_vsp
+      TYPE(FLD_N) ::   sn_ssh, sn_e3t, sn_frq
+      !!
+      NAMELIST/namsbc_sas/ l_sasread, cn_dir, ln_3d_uve, ln_read_frq,   &
+         &                 sn_tem, sn_sal, sn_usp, sn_vsp, sn_ssh, sn_e3t, sn_frq
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_rstart .AND. nn_components == jp_iam_sas )   RETURN
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_ssm_init : sea surface mean data initialisation '
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      !
+      REWIND( numnam_ref )              ! Namelist namsbc_sas in reference namelist : Input fields
+      READ  ( numnam_ref, namsbc_sas, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_sas in reference namelist', lwp )
+      REWIND( numnam_cfg )              ! Namelist namsbc_sas in configuration namelist : Input fields
+      READ  ( numnam_cfg, namsbc_sas, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_sas in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_sas )
+      !           
+      IF(lwp) THEN                              ! Control print
+         WRITE(numout,*) '   Namelist namsbc_sas'
+         WRITE(numout,*) '      Initialisation using an input file                                 l_sasread   = ', l_sasread 
+         WRITE(numout,*) '      Are we supplying a 3D u,v and e3 field                             ln_3d_uve   = ', ln_3d_uve
+         WRITE(numout,*) '      Are we reading frq (fraction of qsr absorbed in the 1st T level)   ln_read_frq = ', ln_read_frq
+      ENDIF
+      !
+      !! switch off stuff that isn't sensible with a standalone module
+      !! note that we need sbc_ssm called first in sbc
+      !
+      IF( ln_apr_dyn ) THEN
+         IF( lwp ) WRITE(numout,*) '         ==>>>   No atmospheric gradient needed with StandAlone Surface scheme'
+         ln_apr_dyn = .FALSE.
+      ENDIF
+      IF( ln_rnf ) THEN
+         IF( lwp ) WRITE(numout,*) '         ==>>>   No runoff needed with StandAlone Surface scheme'
+         ln_rnf = .FALSE.
+      ENDIF
+      IF( ln_ssr ) THEN
+         IF( lwp ) WRITE(numout,*) '         ==>>>   No surface relaxation needed with StandAlone Surface scheme'
+         ln_ssr = .FALSE.
+      ENDIF
+      IF( nn_fwb > 0 ) THEN
+         IF( lwp ) WRITE(numout,*) '         ==>>>   No freshwater budget adjustment needed with StandAlone Surface scheme'
+         nn_fwb = 0
+      ENDIF
+      
+      !                  
+      IF( l_sasread ) THEN                       ! store namelist information in an array
+         ! 
+         !! following code is a bit messy, but distinguishes between when u,v are 3d arrays and
+         !! when we have other 3d arrays that we need to read in
+         !! so if a new field is added i.e. jf_new, just give it the next integer in sequence
+         !! for the corresponding dimension (currently if ln_3d_uve is true, 4 for 2d and 3 for 3d,
+         !! alternatively if ln_3d_uve is false, 6 for 2d and 1 for 3d), reset nfld_3d, nfld_2d,
+         !! and the rest of the logic should still work
+         !
+         jf_tem = 1   ;   jf_ssh = 3   ! default 2D fields index
+         jf_sal = 2   ;   jf_frq = 4   !
+         !
+         IF( ln_3d_uve ) THEN
+            jf_usp = 1   ;   jf_vsp = 2   ;   jf_e3t = 3     ! define 3D fields index
+            nfld_3d  = 2 + COUNT( (/.NOT.ln_linssh/) )       ! number of 3D fields to read
+            nfld_2d  = 3 + COUNT( (/ln_read_frq/) )          ! number of 2D fields to read
+         ELSE
+            jf_usp = 4   ;   jf_e3t = 6                      ! update 2D fields index
+            jf_vsp = 5   ;   jf_frq = 6 + COUNT( (/.NOT.ln_linssh/) )
+            !
+            nfld_3d  = 0                                     ! no 3D fields to read
+            nfld_2d  = 5 + COUNT( (/.NOT.ln_linssh/) ) + COUNT( (/ln_read_frq/) )    ! number of 2D fields to read
+         ENDIF
+         !
+         IF( nfld_3d > 0 ) THEN
+            ALLOCATE( slf_3d(nfld_3d), STAT=ierr )         ! set slf structure
+            IF( ierr > 0 ) THEN
+               CALL ctl_stop( 'sbc_ssm_init: unable to allocate slf 3d structure' )   ;   RETURN
+            ENDIF
+            slf_3d(jf_usp) = sn_usp
+            slf_3d(jf_vsp) = sn_vsp
+            IF( .NOT.ln_linssh )   slf_3d(jf_e3t) = sn_e3t
+         ENDIF
+         !
+         IF( nfld_2d > 0 ) THEN
+            ALLOCATE( slf_2d(nfld_2d), STAT=ierr )         ! set slf structure
+            IF( ierr > 0 ) THEN
+               CALL ctl_stop( 'sbc_ssm_init: unable to allocate slf 2d structure' )   ;   RETURN
+            ENDIF
+            slf_2d(jf_tem) = sn_tem   ;   slf_2d(jf_sal) = sn_sal   ;   slf_2d(jf_ssh) = sn_ssh
+            IF( ln_read_frq )   slf_2d(jf_frq) = sn_frq
+            IF( .NOT. ln_3d_uve ) THEN
+               slf_2d(jf_usp) = sn_usp ; slf_2d(jf_vsp) = sn_vsp
+               IF( .NOT.ln_linssh )   slf_2d(jf_e3t) = sn_e3t
+            ENDIF
+         ENDIF
+         !
+         ierr1 = 0    ! default definition if slf_?d(ifpr)%ln_tint = .false. 
+         IF( nfld_3d > 0 ) THEN
+            ALLOCATE( sf_ssm_3d(nfld_3d), STAT=ierr )         ! set sf structure
+            IF( ierr > 0 ) THEN
+               CALL ctl_stop( 'sbc_ssm_init: unable to allocate sf structure' )   ;   RETURN
+            ENDIF
+            DO ifpr = 1, nfld_3d
+                                            ALLOCATE( sf_ssm_3d(ifpr)%fnow(jpi,jpj,jpk)    , STAT=ierr0 )
+               IF( slf_3d(ifpr)%ln_tint )   ALLOCATE( sf_ssm_3d(ifpr)%fdta(jpi,jpj,jpk,2)  , STAT=ierr1 )
+               IF( ierr0 + ierr1 > 0 ) THEN
+                  CALL ctl_stop( 'sbc_ssm_init : unable to allocate sf_ssm_3d array structure' )   ;   RETURN
+               ENDIF
+            END DO
+            !                                         ! fill sf with slf_i and control print
+            CALL fld_fill( sf_ssm_3d, slf_3d, cn_dir, 'sbc_ssm_init', '3D Data in file', 'namsbc_ssm' )
+         ENDIF
+         !
+         IF( nfld_2d > 0 ) THEN
+            ALLOCATE( sf_ssm_2d(nfld_2d), STAT=ierr )         ! set sf structure
+            IF( ierr > 0 ) THEN
+               CALL ctl_stop( 'sbc_ssm_init: unable to allocate sf 2d structure' )   ;   RETURN
+            ENDIF
+            DO ifpr = 1, nfld_2d
+                                            ALLOCATE( sf_ssm_2d(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
+               IF( slf_2d(ifpr)%ln_tint )   ALLOCATE( sf_ssm_2d(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
+               IF( ierr0 + ierr1 > 0 ) THEN
+                  CALL ctl_stop( 'sbc_ssm_init : unable to allocate sf_ssm_2d array structure' )   ;   RETURN
+               ENDIF
+            END DO
+            !
+            CALL fld_fill( sf_ssm_2d, slf_2d, cn_dir, 'sbc_ssm_init', '2D Data in file', 'namsbc_ssm' )
+         ENDIF
+         !
+         IF( nfld_3d > 0 )   DEALLOCATE( slf_3d, STAT=ierr )
+         IF( nfld_2d > 0 )   DEALLOCATE( slf_2d, STAT=ierr )
+         !
+      ENDIF
+      !
+      CALL sbc_ssm( nit000 )   ! need to define ss?_m arrays used in iceistate
+      l_initdone = .TRUE.
+      !
+   END SUBROUTINE sbc_ssm_init
+
+   !!======================================================================
+END MODULE sbcssm
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/step_c1d.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/step_c1d.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/step_c1d.F90	(revision 11637)
@@ -0,0 +1,119 @@
+MODULE step_c1d
+   !!======================================================================
+   !!                       ***  MODULE step_c1d  ***
+   !! Time-stepping    : manager of the ocean, tracer and ice time stepping - c1d case
+   !!======================================================================
+   !! History :   2.0  !  2004-04  (C. Ethe)  adapted from step.F90 for C1D
+   !!             3.0  !  2008-04  (G. Madec)  redo the adaptation to include SBC
+   !!----------------------------------------------------------------------
+#if defined key_c1d
+   !!----------------------------------------------------------------------
+   !!   'key_c1d'                                       1D Configuration
+   !!----------------------------------------------------------------------  
+   !!   stp_c1d        : NEMO system time-stepping in c1d case
+   !!----------------------------------------------------------------------
+   USE step_oce        ! time stepping definition modules 
+   !LB:USE dyncor_c1d      ! Coriolis term (c1d case)         (dyn_cor_1d     )
+   !LB:USE dynnxt          ! time-stepping                    (dyn_nxt routine)
+   !LB:USE dyndmp          ! U & V momentum damping           (dyn_dmp routine)
+   USE restart         ! restart 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_c1d   ! called by nemogcm.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step_c1d.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_c1d( kstp )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_c1d  ***
+      !!                      
+      !! ** Purpose :  - Time stepping of SBC including sea ice (dynamic and thermodynamic eqs.)
+      !!               - Time stepping of OPA (momentum and active tracer eqs.)
+      !!               - Time stepping of TOP (passive tracer eqs.)
+      !! 
+      !! ** Method  : -1- Update forcings and data  
+      !!              -2- Update vertical ocean physics 
+      !!              -3- Compute the t and s trends 
+      !!              -4- Update t and s 
+      !!              -5- Compute the momentum trends
+      !!              -6- Update the horizontal velocity
+      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
+      !!              -8- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+      !
+      INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   indic    ! error indicator if < 0
+      !! ---------------------------------------------------------------------
+
+                             indic = 0                ! reset to no error condition
+      IF( kstp == nit000 )   CALL iom_init( "nemo")   ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1, "nemo" )   ! say to iom that we are at time step kstp
+
+                             CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+
+                             CALL dia_wri( kstp )         ! ocean model: outputs
+
+                             
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !LOLO:                       
+                            ! LOLO: should write a stp_ctl, who checks on flux realism (current SAS version checks ice fields...)
+                            ! CALL stp_ctl( kstp, indic )
+      !IF( indic < 0  )  THEN
+      !                       CALL ctl_stop( 'step: indic < 0' )
+      !                       CALL dia_wri_state( 'output.abort' )
+      !ENDIF
+      !LOLO.
+                             
+      IF( kstp == nit000   ) CALL iom_close( numror )     ! close input  ocean restart file
+
+      !LOLO: from C1D:
+      IF( lrst_oce       )   CALL rst_write( kstp )        ! write output ocean restart file
+      !
+#if defined key_iomput
+      IF( kstp == nitend .OR. indic < 0 )   CALL xios_context_finalize()   ! needed for XIOS
+#endif
+
+      !LOLO: from SAS"
+!#if defined key_iomput
+!      IF( kstp == nitrst ) THEN
+!         IF(.NOT.lwxios) THEN
+!            CALL iom_close( numrow )
+!         ELSE
+!            CALL iom_context_finalize( cwxios_context )
+!         ENDIF
+!         lrst_oce = .FALSE.
+!      ENDIF
+!      IF( kstp == nitend .OR. indic < 0 ) THEN
+!                             CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
+!      ENDIF
+!#endif
+!      !
+!      IF( ln_timing .AND.  kstp == nit000  )   CALL timing_reset
+      !
+   END SUBROUTINE stp_c1d
+
+#else
+   !!----------------------------------------------------------------------                                                                     
+   !!   Default key                                            NO 1D Config                                                                     
+   !!----------------------------------------------------------------------                                                                     
+CONTAINS
+   SUBROUTINE stp_c1d ( kt )      ! dummy routine                                                                                               
+      IMPLICIT NONE
+      INTEGER, INTENT( in ) :: kt
+      WRITE(*,*) 'stp_c1d: You should not have seen this print! error?', kt
+   END SUBROUTINE stp_c1d
+#endif
+   
+   !!======================================================================
+END MODULE step_c1d
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_hgr.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_hgr.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_hgr.F90	(revision 11637)
@@ -0,0 +1,110 @@
+MODULE usrdef_hgr
+   !!======================================================================
+   !!                     ***  MODULE usrdef_hgr   ***
+   !!
+   !!                     ===  STATION_ASF configuration  ===
+   !!
+   !! User defined :   mesh and Coriolis parameter of a user configuration
+   !!======================================================================
+   !! History :  4.0 ! 2016-03  (S. Flavoni)
+   !! History :  4.0 ! 2019-10  (L. Brodeau)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_hgr   : initialize the horizontal mesh 
+   !!----------------------------------------------------------------------
+   USE dom_oce  , ONLY: nimpp, njmpp        ! ocean space and time domain
+   USE c1d      ,  ONLY: rn_lon1d, rn_lat1d ! ocean lon/lat define by namelist
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE usrdef_nam     !
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_hgr   ! called in domhgr.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_hgr.F90 10072 2018-08-28 15:21:50Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
+      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
+      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain not on the sphere)
+      &                    pe1t  , pe1u  , pe1v   , pe1f   ,   &   ! scale factors       (required)
+      &                    pe2t  , pe2u  , pe2v   , pe2f   ,   &   !
+      &                    ke1e2u_v      , pe1e2u , pe1e2v     )   ! u- & v-surfaces (if gridsize reduction is used in strait(s))
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE usr_def_hgr  ***
+      !!
+      !! ** Purpose :   user defined mesh and Coriolis parameter
+      !!
+      !! ** Method  :   set all intent(out) argument to a proper value
+      !!
+      !!                Here STATION_ASF configuration :
+      !!          Rectangular 3x3 domain 
+      !!          - Located at 150E-50N
+      !!          - a constant horizontal resolution  
+      !!
+      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 
+      !!              - define coriolis parameter at f-point if the domain in not on the sphere
+      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
+      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
+      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors                             [m]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors                             [m]
+      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces computed here, =0 otherwise 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
+      !
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0   ! local scalars
+      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme, zfact !   -      -
+      !!-------------------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_hgr : C1D configuration'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   uniform grid spacing WITH Coriolis force'
+      IF(lwp) WRITE(numout,*) ' LOLO => rn_lon1d, rn_lat1d =', rn_lon1d, rn_lat1d
+      !
+      !                       !==  grid point position  ==!   (Special case of 1D config: all lon/lat are the same)
+      !                       ! longitude
+      plamt(:,:) = rn_lon1d
+      plamu(:,:) = rn_lon1d 
+      plamv(:,:) = rn_lon1d
+      plamf(:,:) = rn_lon1d
+      !                       ! latitude
+      pphit(:,:) = rn_lat1d
+      pphiu(:,:) = rn_lat1d
+      pphiv(:,:) = rn_lat1d
+      pphif(:,:) = rn_lat1d
+      !                       !== Horizontal scale factors ==! (in meters)
+      !                     
+      !                                         ! constant grid spacing
+      pe1t(:,:) = 100.  ;   pe2t(:,:) = 100.
+      pe1u(:,:) = 100.  ;   pe2u(:,:) = 100.
+      pe1v(:,:) = 100.  ;   pe2v(:,:) = 100.
+      pe1f(:,:) = 100.  ;   pe2f(:,:) = 100.
+      !
+      !                                         ! NO reduction of grid size in some straits 
+      ke1e2u_v = 0                              !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
+      pe1e2u(:,:) = 0._wp                       !    CAUTION: set to zero to avoid error with some compilers that
+      pe1e2v(:,:) = 0._wp                       !             require an initialization of INTENT(out) arguments
+      !
+      !
+      !                       !==  Coriolis parameter  ==!
+      kff = 0                 !  indicate to compute Coriolis parameter afterward
+      !
+   END SUBROUTINE usr_def_hgr
+
+   !!======================================================================
+END MODULE usrdef_hgr
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_nam.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_nam.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_nam.F90	(revision 11637)
@@ -0,0 +1,105 @@
+MODULE usrdef_nam
+   !!======================================================================
+   !!                     ***  MODULE usrdef_nam   ***
+   !!
+   !!                     ===  STATION_ASF configuration  ===
+   !!
+   !! User defined : set the domain characteristics of a user configuration
+   !!======================================================================
+   !! History :  4.0  ! 2016-03  (S. Flavoni, G. Madec)  Original code
+   !! History :  trunk  ! 2019-10 (L. Brodeau) for STATION_ASF (C1D meets SAS)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_nam   : read user defined namelist and set global domain size
+   !!   usr_def_hgr   : initialize the horizontal mesh 
+   !!----------------------------------------------------------------------
+   USE dom_oce  , ONLY: nimpp, njmpp             ! ocean space and time domain
+   USE dom_oce  , ONLY: ln_zco, ln_zps, ln_sco   ! flag of type of coordinate
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_nam   ! called in nemogcm.F90 module
+
+   !                              !!* namusr_def namelist *!!
+   LOGICAL, PUBLIC ::   ln_bench   ! =T benchmark test with gyre: the gridsize is constant (no need to adjust timestep or viscosity)
+   INTEGER, PUBLIC ::   nn_GYRE    ! 1/nn_GYRE = the resolution chosen in degrees and thus defining the horizontal domain size
+   REAL(wp), PUBLIC::   rn_dept1   ! Depth (m) at which the SST is taken/measured == depth of first T point!
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_nam.F90 10072 2018-08-28 15:21:50Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_nam( ldtxt, ldnam, cd_cfg, kk_cfg, kpi, kpj, kpk, kperio )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read user defined namelist and define the domain size
+      !!
+      !! ** Method  :   read in namusr_def containing all the user specific namelist parameter
+      !!
+      !!                Here C1D configuration
+      !!
+      !! ** input   : - namusr_def namelist found in namelist_cfg
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt, ldnam    ! stored print information
+      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name
+      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution
+      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
+      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
+      !
+      INTEGER ::   ios, ii   ! Local integer
+      !!
+      NAMELIST/namusr_def/ rn_dept1
+      !!----------------------------------------------------------------------
+      !
+      ii = 1
+      !
+      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only)
+      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', .TRUE. )
+      !
+      WRITE( ldnam(:), namusr_def )
+      !
+      cd_cfg = 'STATION_ASF'               ! name & resolution (not used)
+      kk_cfg = 0
+
+      ! Global Domain size: STATION_ASF domain is 3 x 3 grid-points x 75 or vertical levels
+      kpi = 3
+      kpj = 3
+      kpk = 1
+      !
+      !                             ! control print
+      WRITE(ldtxt(ii),*) '   '                                                                            ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'     ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '~~~~~~~~~~~ '                                                                   ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '   Namelist namusr_def : C1 case'                                               ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      type of vertical coordinate : '                                           ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '         z-coordinate flag                     ln_zco = ', ln_zco               ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '         z-partial-step coordinate flag        ln_zps = ', ln_zps               ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '         s-coordinate flag                     ln_sco = ', ln_sco               ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      C1D domain = 3 x 3 x 1 grid-points                '                       ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '         resulting global domain size :        jpiglo = ', kpi                  ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '                                               jpjglo = ', kpj                  ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '                                               jpkglo = ', kpk                  ;   ii = ii + 1
+
+      !
+      !                             ! Set the lateral boundary condition of the global domain
+      kperio =  7                   ! C1D configuration : 3x3 basin with cyclic Est-West and Norht-South condition
+      !
+      WRITE(ldtxt(ii),*) '   Lateral boundary condition of the global domain'                           ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      C1D : closed basin                 jperio = ', kperio                   ;   ii = ii + 1
+      !
+   END SUBROUTINE usr_def_nam
+
+   !!======================================================================
+END MODULE usrdef_nam
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_zgr.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_zgr.F90	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/usrdef_zgr.F90	(revision 11637)
@@ -0,0 +1,90 @@
+MODULE usrdef_zgr
+   !!======================================================================
+   !!                     ***  MODULE usrdef_zgr  ***
+   !!
+   !!                       ===  STATION_ASF case  ===
+   !!
+   !! user defined :  vertical coordinate system of a user configuration
+   !!======================================================================
+   !! History :  4.0  ! 2019-10  (L. Brodeau)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_zgr   : user defined vertical coordinate system (required)
+   !!---------------------------------------------------------------------
+   USE oce            ! ocean variables
+   !USE dom_oce        ! ocean domain
+   !USE depth_e3       ! depth <=> e3
+   USE usrdef_nam     ! User defined : namelist variables
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! distributed memory computing library
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_zgr   ! called by domzgr.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_zgr.F90 10072 2018-08-28 15:21:50Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
+      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
+      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
+      &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors
+      &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      -
+      &                    k_top  , k_bot    )                             ! top & bottom ocean level
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE usr_def_zgr  ***
+      !!
+      !! ** Purpose :   User defined the vertical coordinates
+      !!
+      !!----------------------------------------------------------------------
+      LOGICAL                   , INTENT(  out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags ( read in namusr_def )
+      LOGICAL                   , INTENT(  out) ::   ld_isfcav                   ! under iceshelf cavity flag
+      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth        [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors
+      INTEGER , DIMENSION(:,:)  , INTENT(  out) ::   k_top, k_bot                ! first & last ocean level
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_zgr : LOLO STATION_ASF configuration, setting first level properties.'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !
+
+      ld_zco    = .TRUE.         ! z-coordinate without ocean cavities
+      ld_zps    = .FALSE.
+      ld_sco    = .FALSE.
+      ld_isfcav = .FALSE.
+      
+      pdept_1d(1) = rn_dept1 ! depth (m) at which the SST is taken/measured == depth of first T point!
+      pdepw_1d(1) = 0._wp
+      pe3t_1d(1)  = 2._wp*rn_dept1
+      pe3w_1d(1)  = rn_dept1 ! LB???
+
+      pdept(:,:,:) = rn_dept1
+      pdepw(:,:,:) = 0._wp
+      pe3t(:,:,:) = 2._wp*rn_dept1
+      pe3u(:,:,:) = 2._wp*rn_dept1
+      pe3v(:,:,:) = 2._wp*rn_dept1
+      pe3f(:,:,:) = 2._wp*rn_dept1
+      pe3w(:,:,:)  = rn_dept1  ! LB???
+      pe3uw(:,:,:) = rn_dept1  ! LB???
+      pe3vw(:,:,:) = rn_dept1  ! LB???
+      k_top = 1
+      k_bot = 1
+      !
+   END SUBROUTINE usr_def_zgr
+   !!======================================================================
+END MODULE usrdef_zgr
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/cpp_STATION_ASF.fcm
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/cpp_STATION_ASF.fcm	(revision 11637)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/cpp_STATION_ASF.fcm	(revision 11637)
@@ -0,0 +1,1 @@
+ bld::tool::fppkeys   key_c1d key_mpp_mpi key_iomput
