Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca025_like
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca025_like	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca025_like	(revision 10179)
@@ -53,4 +53,5 @@
 !-----------------------------------------------------------------------
   ln_usr      = .true.    !  user defined formulation                  (T => check usrdef_sbc)
+  ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr)
 /
 
@@ -83,4 +84,11 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
+!-----------------------------------------------------------------------
+   !                       !  type of penetration                        (default: NO selection)
+   ln_qsr_rgb  = .true.      !  RGB light penetration (Red-Green-Blue)
+   nn_chldta   =      0      !  RGB : Chl data (=1) or cst value (=0)
+/
 !-----------------------------------------------------------------------
 &namdrg        !   top/bottom drag coefficient                          (default: NO selection)
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca12_like
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca12_like	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca12_like	(revision 10179)
@@ -53,4 +53,5 @@
 !-----------------------------------------------------------------------
   ln_usr      = .true.    !  user defined formulation                  (T => check usrdef_sbc)
+  ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr)
 /
 
@@ -83,4 +84,11 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
+!-----------------------------------------------------------------------
+   !                       !  type of penetration                        (default: NO selection)
+   ln_qsr_rgb  = .true.      !  RGB light penetration (Red-Green-Blue)
+   nn_chldta   =      0      !  RGB : Chl data (=1) or cst value (=0)
+/
 !-----------------------------------------------------------------------
 &namdrg        !   top/bottom drag coefficient                          (default: NO selection)
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca1_like
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca1_like	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_cfg_orca1_like	(revision 10179)
@@ -53,4 +53,5 @@
 !-----------------------------------------------------------------------
   ln_usr      = .true.    !  user defined formulation                  (T => check usrdef_sbc)
+  ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr)
 /
 
@@ -83,4 +84,11 @@
 !!======================================================================
 !
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
+!-----------------------------------------------------------------------
+   !                       !  type of penetration                        (default: NO selection)
+   ln_qsr_rgb  = .true.      !  RGB light penetration (Red-Green-Blue)
+   nn_chldta   =      0      !  RGB : Chl data (=1) or cst value (=0)
+/
 !-----------------------------------------------------------------------
 &namdrg        !   top/bottom drag coefficient                          (default: NO selection)
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_ice_cfg	(revision 10179)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_ice_cfg	(revision 10179)
@@ -0,0 +1,84 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! SI3 configuration namelist: Overwrites SHARED/namelist_ice_ref
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface forcing                (namforcing)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (namthd_pnd)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!
+!------------------------------------------------------------------------------
+&nampar         !   Generic parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namitd         !   Ice discretization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn         !   Ice dynamics
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn_rdgrft  !   Ice ridging/rafting
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn_rhg     !   Ice rheology
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn_adv     !   Ice advection
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namforcing     !   Ice surface forcing
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd         !   Ice thermodynamics
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_zdf     !   Ice heat diffusion
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_da      !   Ice lateral melting
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_do      !   Ice growth in open water
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_sal     !   Ice salinity
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
+!------------------------------------------------------------------------------
+/
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_ice_ref
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_ice_ref	(revision 10179)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/EXPREF/namelist_ice_ref	(revision 10179)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ice_ref
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/diawri.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/diawri.F90	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/diawri.F90	(revision 10179)
@@ -408,4 +408,10 @@
    SUBROUTINE dia_wri( kt )
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
+      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
+         CALL dia_wri_state( 'output.init', kt )
+         ninist = 0
+      ENDIF
+
    END SUBROUTINE dia_wri
 
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/usrdef_istate.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/usrdef_istate.F90	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/usrdef_istate.F90	(revision 10179)
@@ -53,5 +53,6 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
-      INTEGER :: ji, jj, jk
+      REAL(wp) ::   zfact
+      INTEGER  ::   ji, jj, jk
       !!----------------------------------------------------------------------
       !
@@ -60,22 +61,24 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   '
       !
-      ! define unique value on each point
+      ! define unique value on each point. z2d ranging from 0.05 to -0.05
       DO jj = 1, jpj
          DO ji = 1, jpi
-            z2d(ji,jj) = ( 0.5_wp - REAL( nimpp + ji - 1 + ( njmpp + jj - 2 ) * jpiglo) / REAL ( jpiglo * jpjglo ) ) * 0.1_wp
+            z2d(ji,jj) = 0.1 * ( 0.5 - REAL( nimpp + ji - 1 + ( njmpp + jj - 2 ) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) )
          ENDDO
       ENDDO
       !
       ! sea level:
-      pssh(:,:) = z2d(:,:)
+      pssh(:,:) = z2d(:,:)                                                ! +/- 0.05 m
       !
       DO jk = 1, jpk
-         ! temperature:
-         pts(:,:,jk,jp_tem) = 10._wp + z2d(:,:) * 2._wp
+         zfact = REAL(jk-1,wp) / REAL(jpk-1,wp)   ! 0 to 1 to add a basic stratification
+         ! temperature choosen to lead to 20% ice
+         pts(:,:,jk,jp_tem) = 2._wp - 0.1_wp * zfact + z2d(:,:) * 100._wp ! 2 to 1.9 +/- 5 degG
+         WHERE ( pts(:,:,jk,jp_tem) < -1.5_wp ) pts(:,:,jk,jp_tem) = -1.5_wp + z2d(:,:) * 0.2_wp  
          ! salinity:  
-         pts(:,:,jk,jp_sal) = 35._wp + z2d(:,:)
+         pts(:,:,jk,jp_sal) = 30._wp + 1._wp * zfact + z2d(:,:)           ! 30 to 31 +/- 0.05 psu
          ! velocities:
-         pu(:,:,jk) = z2d(:,:) * 0.1_wp
-         pv(:,:,jk) = z2d(:,:) * 0.01_wp
+         pu(:,:,jk) = z2d(:,:) * 0.1_wp                                   ! +/- 0.005  m/s
+         pv(:,:,jk) = z2d(:,:) * 0.01_wp                                  ! +/- 0.0005 m/s
       ENDDO
       !
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/usrdef_sbc.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/usrdef_sbc.F90	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/MY_SRC/usrdef_sbc.F90	(revision 10179)
@@ -13,6 +13,17 @@
    !!   usr_def_sbc    : user defined surface bounday conditions in BENCH case
    !!----------------------------------------------------------------------
+   USE par_oce        ! ocean space and time domain
+   USE dom_oce        
+   USE oce             ! ocean dynamics and tracers
    USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbc_ice         ! Surface boundary condition: ocean fields
    USE in_out_manager  ! I/O manager
+   USE phycst          ! physical constants
+   USE lib_mpp         ! MPP library
+   USE lbclnk          ! lateral boundary conditions - mpp exchanges
+
+#if defined key_si3
+   USE ice, ONLY       : at_i_b, a_i_b
+#endif
 
    IMPLICIT NONE
@@ -61,5 +72,11 @@
          qns (:,:) = 0._wp
          qsr (:,:) = 0._wp
-         !         
+         !
+         utau_b(:,:) = 0._wp 
+         vtau_b(:,:) = 0._wp
+         emp_b (:,:) = 0._wp
+         sfx_b (:,:) = 0._wp
+         qns_b (:,:) = 0._wp
+         !
       ENDIF
       
@@ -67,10 +84,94 @@
    END SUBROUTINE usrdef_sbc_oce
 
+   
    SUBROUTINE usrdef_sbc_ice_tau( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE usrdef_sbc_ice_tau  ***
+      !!
+      !! ** Purpose :   provide the surface boundary (momentum) condition over
+      !sea-ice
+      !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
+      INTEGER  ::   ji, jj
+      !!---------------------------------------------------------------------
+#if defined key_si3
+      IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : BENCH case: constant stress forcing'
+      !
+      ! define unique value on each point. z2d ranging from 0.05 to -0.05
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            z2d(ji,jj) = 0.1 * ( 0.5 - REAL( nimpp + ji - 1 + ( njmpp + jj - 2 ) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) )
+         ENDDO
+      ENDDO
+      utau_ice(:,:) = 0.1_wp +  z2d(:,:)
+      vtau_ice(:,:) = 0.1_wp +  z2d(:,:)
+
+      CALL lbc_lnk_multi( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+#endif
+      !
    END SUBROUTINE usrdef_sbc_ice_tau
 
-   SUBROUTINE usrdef_sbc_ice_flx( kt )
+   
+   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE usrdef_sbc_ice_flx  ***
+      !!
+      !! ** Purpose :   provide the surface boundary (flux) condition over
+      !sea-ice
+      !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
+      !!
+      REAL(wp) ::   zfr1, zfr2                 ! local variables
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsnw   ! snw distribution after wind blowing
+      !!---------------------------------------------------------------------
+      !
+#if defined key_si3
+      IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : BENCH case: NO flux forcing'
+      !
+      ! ocean variables (renaming)
+      emp_oce (:,:)   = 0._wp   ! uniform value for freshwater budget (E-P)
+      qsr_oce (:,:)   = 0._wp   ! uniform value for     solar radiation
+      qns_oce (:,:)   = 0._wp   ! uniform value for non-solar radiation
+
+      ! ice variables
+      alb_ice (:,:,:) = 0.7_wp  ! useless
+      qsr_ice (:,:,:) = 0._wp   ! uniform value for     solar radiation
+      qns_ice (:,:,:) = 0._wp   ! uniform value for non-solar radiation
+      sprecip (:,:)   = 0._wp   ! uniform value for snow precip
+      evap_ice(:,:,:) = 0._wp   ! uniform value for sublimation
+
+      ! ice fields deduced from above
+      zsnw(:,:) = 1._wp
+      !!CALL lim_thd_snwblow( at_i_b, zsnw )  ! snow distribution over ice after
+      !wind blowing 
+      emp_ice  (:,:)   = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw(:,:)
+      emp_oce  (:,:)   = emp_oce(:,:) - sprecip(:,:) * (1._wp - zsnw(:,:) )
+      qevap_ice(:,:,:) =   0._wp
+      qprec_ice(:,:)   =   rhos * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) !  in J/m3
+      qemp_oce (:,:)   = - emp_oce(:,:) * sst_m(:,:) * rcp
+      qemp_ice (:,:)   =   sprecip(:,:) * zsnw * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! solid precip (only)
+
+      ! total fluxes
+      emp_tot (:,:) = emp_ice  + emp_oce
+      qns_tot (:,:) = at_i_b(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
+      qsr_tot (:,:) = at_i_b(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
+
+      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
+      !
+      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
+      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
+         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
+      ELSEWHERE                                                         ! zero when hs>0
+         qtr_ice_top(:,:,:) = 0._wp 
+      END WHERE
+#endif
+
    END SUBROUTINE usrdef_sbc_ice_flx
 
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/cpp_BENCH.fcm
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/cpp_BENCH.fcm	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/BENCH/cpp_BENCH.fcm	(revision 10179)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys  key_mpp_mpi  
+ bld::tool::fppkeys  key_mpp_mpi  key_si3
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/demo_cfgs.txt
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/demo_cfgs.txt	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/demo_cfgs.txt	(revision 10179)
@@ -6,3 +6,3 @@
 VORTEX OCE NST
 WAD OCE
-BENCH OCE
+BENCH OCE ICE
