Index: /NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref	(revision 13373)
+++ /NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref	(revision 13374)
@@ -616,4 +616,5 @@
    !
    ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data)
+      ln_icb_grd           = .false. ! ground icb when icb bottom level hit oce bottom level (need ln_M2016 to be activated)
    !
    cn_dir      = './'      !  root directory for the calving data location
Index: /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90
===================================================================
--- /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90	(revision 13373)
+++ /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90	(revision 13374)
@@ -129,5 +129,5 @@
    LOGICAL , PUBLIC ::   ln_time_average_weight          !: Time average the weight on the ocean    !!gm I don't understand that !
    REAL(wp), PUBLIC ::   rn_speed_limit                  !: CFL speed limit for a berg
-   LOGICAL , PUBLIC ::   ln_M2016                        !: use Nacho's Merino 2016 work
+   LOGICAL , PUBLIC ::   ln_M2016, ln_icb_grd            !: use Nacho's Merino 2016 work
    !
    ! restart
Index: /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90
===================================================================
--- /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90	(revision 13373)
+++ /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90	(revision 13374)
@@ -98,6 +98,6 @@
          zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
          !
-         CALL icb_ground( zxi2, zxi1, zu1,   &
-            &             zyj2, zyj1, zv1, ll_bounced )
+         CALL icb_ground( berg, zxi2, zxi1, zu1,   &
+            &                   zyj2, zyj1, zv1, ll_bounced )
 
          !                                         !**   A2 = A(X2,V2)
@@ -114,6 +114,6 @@
          zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
          !
-         CALL icb_ground( zxi3, zxi1, zu3,   &
-            &                zyj3, zyj1, zv3, ll_bounced )
+         CALL icb_ground( berg, zxi3, zxi1, zu3,   &
+            &                   zyj3, zyj1, zv3, ll_bounced )
 
          !                                         !**   A3 = A(X3,V3)
@@ -130,6 +130,6 @@
          zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
 
-         CALL icb_ground( zxi4, zxi1, zu4,   &
-            &             zyj4, zyj1, zv4, ll_bounced )
+         CALL icb_ground( berg, zxi4, zxi1, zu4,   &
+            &                   zyj4, zyj1, zv4, ll_bounced )
 
          !                                         !**   A4 = A(X4,V4)
@@ -149,6 +149,6 @@
          zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )
 
-         CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
-            &             zyj_n, zyj1, zvvel_n, ll_bounced )
+         CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   &
+            &                   zyj_n, zyj1, zvvel_n, ll_bounced )
 
          pt%uvel = zuvel_n                        !** save in berg structure
@@ -164,6 +164,6 @@
 
 
-   SUBROUTINE icb_ground( pi, pi0, pu,   &
-      &                   pj, pj0, pv, ld_bounced )
+   SUBROUTINE icb_ground( berg, pi, pi0, pu,   &
+      &                         pj, pj0, pv, ld_bounced )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_ground  ***
@@ -174,4 +174,6 @@
       !!                NB two possibilities available one of which is hard-coded here
       !!----------------------------------------------------------------------
+      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
+      !
       REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
       REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
@@ -181,5 +183,9 @@
       INTEGER  ::   ii, ii0
       INTEGER  ::   ij, ij0
+      INTEGER  ::   ikb
       INTEGER  ::   ibounce_method
+      !
+      REAL(wp) :: zD 
+      REAL(wp), DIMENSION(jpk) :: ze3t
       !!----------------------------------------------------------------------
       !
@@ -198,15 +204,21 @@
       !
       ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
-      !IF ( ln_icb_ground ) THEN
-      !   ! interpol needed data
-      !   CALL icb_utl_interp( pxi, pyj, pe3t=ze3t )   ! 3d velocities
-      !  
-      !   !compute bottom level
-      !   CALL icb_utl_getkb( ikb, ze3t, zD )
-      !
-      !   IF(  tmask(ii,ij,ikb)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
-      !ELSE
-      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
-      !END IF
+      IF ( ln_M2016 .AND. ln_icb_grd ) THEN
+         !
+         ! draught (keel depth)
+         zD = rho_berg_1_oce * berg%current_point%thickness
+         !
+         ! interpol needed data
+         CALL icb_utl_interp( pi, pj, pe3t=ze3t )
+         ! 
+         !compute bottom level
+         CALL icb_utl_getkb( ikb, ze3t, zD )
+         !
+         ! berg reach a new t-cell, but an ocean one
+         IF(  tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN
+         !
+      ELSE
+         IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
+      END IF
       !
       ! From here, berg have reach land: treat grounding/bouncing
Index: /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90
===================================================================
--- /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90	(revision 13373)
+++ /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90	(revision 13374)
@@ -415,5 +415,6 @@
          &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   &
          &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      , ln_M2016             ,   &
-         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out
+         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out        ,   &
+         &              ln_icb_grd
       !!----------------------------------------------------------------------
 
@@ -494,4 +495,7 @@
             &                    'bits_erosion_fraction = ', rn_bits_erosion_fraction
 
+         WRITE(numout,*) '   Use icb module modification from Merino et al. (2016) : ln_M2016 = ', ln_M2016
+         WRITE(numout,*) '       ground icebergs if icb bottom lvl hit the oce bottom level : ln_icb_grd = ', ln_icb_grd
+
          WRITE(numout,*) '   Shift of sea-ice concentration in erosion flux modulation ',   &
             &                    '(0<sicn_shift<1)    rn_sicn_shift  = ', rn_sicn_shift
Index: /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90
===================================================================
--- /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90	(revision 13373)
+++ /NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90	(revision 13374)
@@ -164,4 +164,5 @@
       REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask
       REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm
+      REAL(wp), DIMENSION(4,jpk) :: zw1d
       INTEGER                :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner
       INTEGER                :: iiTp, iiTm, ijTp, ijTm
@@ -225,8 +226,22 @@
       ! 3d interpolation
       IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN
-         puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zwU )
-         pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zwV )
+         ! no need to mask as 0 is a valid data for land
+         zw1d(1,:) = zwU(1) ; zw1d(2,:) = zwU(2) ; zw1d(3,:) = zwU(3) ; zw1d(4,:) = zwU(4) ;
+         puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zw1d )
+
+         zw1d(1,:) = zwV(1) ; zw1d(2,:) = zwV(2) ; zw1d(3,:) = zwV(3) ; zw1d(4,:) = zwV(4) ;
+         pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zw1d )
       END IF
-      IF ( PRESENT(ptoce) ) ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zwT * zmskT )
+
+      IF ( PRESENT(ptoce) ) THEN
+         ! for temperature we need to mask the weight properly
+         ! no need of extra halo as it is a T point variable
+         zw1d(1,:) = tmask(iiT  ,ijT  ,:) * zwT(1) * zmskT(1)
+         zw1d(2,:) = tmask(iiT+1,ijT  ,:) * zwT(2) * zmskT(2)
+         zw1d(3,:) = tmask(iiT  ,ijT+1,:) * zwT(3) * zmskT(3)
+         zw1d(4,:) = tmask(iiT+1,ijT+1,:) * zwT(4) * zmskT(4)
+         ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zw1d )
+      END IF
+      !
       IF ( PRESENT(pe3t)  ) pe3t(:)  = e3t_e(iiT,ijT,:)    ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results
       !
@@ -300,6 +315,7 @@
          IF ( ierr > 0 ) THEN
             WRITE(numout,*) 'bottom left corner T point out of bound'
-            WRITE(numout,*) kii, mig( 1 ), mig(jpi)
-            WRITE(numout,*) kij, mjg( 1 ), mjg(jpj)
+            WRITE(numout,*) pi, kii, mig( 1 ), mig(jpi)
+            WRITE(numout,*) pj, kij, mjg( 1 ), mjg(jpj)
+            WRITE(numout,*) pmsk
             CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)')
          END IF
@@ -381,9 +397,10 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp), DIMENSION(4)                   , INTENT(in) ::   pw        ! weight
+      REAL(wp), DIMENSION(4,jpk)               , INTENT(in) ::   pw        ! weight
       INTEGER ,                                  INTENT(in) ::   pii, pij  ! bottom left corner
       REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h
       !
       REAL(wp), DIMENSION(4,jpk) :: zdat ! input data
+      INTEGER :: jk
       !!----------------------------------------------------------------------
       !
@@ -395,5 +412,8 @@
       !
       ! compute interpolated value
-      icb_utl_bilin_3d_h(:) = ( zdat(1,:)*pw(1) + zdat(2,:)*pw(2) + zdat(3,:)*pw(3) + zdat(4,:)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4)) 
+      DO jk=1,jpk
+         icb_utl_bilin_3d_h(jk) =   ( zdat(1,jk)*pw(1,jk) + zdat(2,jk)*pw(2,jk) + zdat(3,jk)*pw(3,jk) + zdat(4,jk)*pw(4,jk) ) &
+            &                     /   MAX(1.e-20, pw(1,jk)+pw(2,jk)+pw(3,jk)+pw(4,jk)) 
+      END DO
       !
    END FUNCTION icb_utl_bilin_3d_h
