Index: NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref	(revision 13357)
@@ -614,5 +614,7 @@
    ln_use_calving          = .false. ! Use calving data even when nn_test_icebergs > 0
    rn_speed_limit          = 0.      ! CFL speed limit for a berg
-
+   !
+   ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data)
+   !
    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 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90	(revision 13357)
@@ -57,10 +57,11 @@
    TYPE, PUBLIC ::   point              !: properties of an individual iceberg (position, mass, size, etc...)
       INTEGER  ::   year
-      REAL(wp) ::   xi , yj                                              ! iceberg coordinates in the (i,j) referential (global)
+      REAL(wp) ::   xi , yj , zk                                         ! iceberg coordinates in the (i,j) referential (global) and deepest level affected
       REAL(wp) ::   e1 , e2                                              ! horizontal scale factors at the iceberg position
       REAL(wp) ::   lon, lat, day                                        ! geographic position
       REAL(wp) ::   mass, thickness, width, length, uvel, vvel           ! iceberg physical properties
-      REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi    ! properties of iceberg environment 
+      REAL(wp) ::   ssu, ssv, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi  ! properties of iceberg environment 
       REAL(wp) ::   mass_of_bits, heat_density
+      INTEGER  ::   kb                                                   ! icb bottom level
    END TYPE point
 
@@ -85,15 +86,12 @@
    ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position
    ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tt_e, fr_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssu_e, ssv_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   sst_e, fr_e
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tmask_e, umask_e, vmask_e
    REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   rlon_e, rlat_e, ff_e
-      ! prepration to Nacho work
-      ! tt3d_e
-      ! uo3d_e
-      ! vo3d_e
-      !
+   REAl(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   uoce_e, voce_e, toce_e, e3t_e
+   !
 #if defined key_si3 || defined key_cice
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   hi_e, ui_e, vi_e
@@ -123,4 +121,5 @@
    INTEGER , PUBLIC ::   nn_verbose_write                !: timesteps between verbose messages
    REAL(wp), PUBLIC ::   rn_rho_bergs                    !: Density of icebergs
+   REAL(wp), PUBLIC ::   rho_berg_1_oce                  !: convertion factor (thickness to draft) (rn_rho_bergs/pp_rho_seawater)
    REAL(wp), PUBLIC ::   rn_LoW_ratio                    !: Initial ratio L/W for newly calved icebergs
    REAL(wp), PUBLIC ::   rn_bits_erosion_fraction        !: Fraction of erosion melt flux to divert to bergy bits
@@ -130,4 +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
    !
    ! restart
@@ -141,4 +141,5 @@
    REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_initial_thickness !  Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   src_calving, src_calving_hflx    !: accumulate input ice
+   INTEGER , PUBLIC             , SAVE                     ::   micbkb                           !: deepest level affected by icebergs
    INTEGER , PUBLIC             , SAVE                     ::   numicb                           !: iceberg IO
    INTEGER , PUBLIC             , SAVE, DIMENSION(nkounts) ::   num_bergs                        !: iceberg counter
@@ -180,6 +181,6 @@
       !
       ! expanded arrays for bilinear interpolation
-      ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,    &
-         &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,    &
+      ALLOCATE( ssu_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   &
+         &      ssv_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   &
 #if defined key_si3 || defined key_cice
          &      ui_e(0:jpi+1,0:jpj+1) ,                            &
@@ -188,5 +189,5 @@
 #endif
          &      fr_e(0:jpi+1,0:jpj+1) ,                            &
-         &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   &
+         &      sst_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,  &
          &      first_width(nclasses) , first_length(nclasses) ,   &
          &      src_calving (jpi,jpj) ,                            &
@@ -194,8 +195,9 @@
       icb_alloc = icb_alloc + ill
 
-      ! prepration to Nacho work
-      ! tt3d_e
-      ! uo3d_e
-      ! vo3d_e
+      IF ( ln_M2016 ) THEN
+         ALLOCATE( uoce_e(0:jpi+1,0:jpj+1,jpk), voce_e(0:jpi+1,0:jpj+1,jpk), &
+            &      toce_e(0:jpi+1,0:jpj+1,jpk), e3t_e(0:jpi+1,0:jpj+1,jpk) , STAT=ill )
+         icb_alloc = icb_alloc + ill
+      END IF
       !
       ALLOCATE( tmask_e(0:jpi+1,0:jpj+1), umask_e(0:jpi+1,0:jpj+1), vmask_e(0:jpi+1,0:jpj+1), &
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbclv.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbclv.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbclv.F90	(revision 13357)
@@ -21,5 +21,4 @@
    USE lbclnk         ! NEMO boundary exchanges for gridded data
 
-   USE icb_oce        ! iceberg variables
    USE icbdia         ! iceberg diagnostics
    USE icbutl         ! iceberg utility routines
@@ -142,4 +141,5 @@
                   newpt%mass           = rn_initial_mass     (jn)
                   newpt%thickness      = rn_initial_thickness(jn)
+                  newpt%kb             = 1         ! compute correctly in icbthm if needed        
                   newpt%width          = first_width         (jn)
                   newpt%length         = first_length        (jn)
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90	(revision 13357)
@@ -197,5 +197,16 @@
       ij  = mj1( ij  )
       !
+      ! 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
       !
       ! From here, berg have reach land: treat grounding/bouncing
@@ -254,35 +265,27 @@
       REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
       !
-      INTEGER  ::   itloop
-      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi
-      REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y
+      INTEGER  ::   itloop, ikb, jk
+      REAL(wp) ::   zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi
+      REAL(wp) ::   zvo, zssv, zvi, zva, zvwave, zssh_y
       REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
       REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
-      REAL(wp) ::   z_ocn, z_atm, z_ice
+      REAL(wp) ::   z_ocn, z_atm, z_ice, zdep
       REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
       REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
       REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
+      REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw
       !!----------------------------------------------------------------------
 
       ! Interpolate gridded fields to berg
       nknberg = berg%number(1)
-      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,   &   ! scale factor
-         &                 puo=zuo, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
-         &                 pvo=zvo, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities
-         &                 pssh_i=zssh_x, pssh_j=zssh_y, &   ! ssh gradient
-         &                 phi=zhi, pff=zff )               ! ice thickness and coriolis
-
-      IF (lwp) THEN
-         WRITE(numout,*) 'icbdyn ', pxi, pyj 
-         WRITE(numout,*) pe1, zuo, zui, zua, zssh_x
-         WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y
-         WRITE(numout,*) zhi, zff 
-         WRITE(numout,*) ''
-      END IF
-
+      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,     &   ! scale factor
+         &                 pssu=zssu, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
+         &                 pssv=zssv, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities
+         &                 pssh_i=zssh_x, pssh_j=zssh_y,   &   ! ssh gradient
+         &                 phi=zhi, pff=zff)                   ! ice thickness and coriolis
 
       zM = berg%current_point%mass
       zT = berg%current_point%thickness               ! total thickness
-      zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth)
+      zD = rho_berg_1_oce * zT                        ! draught (keel depth)
       zF = zT - zD                                    ! freeboard
       zW = berg%current_point%width
@@ -291,7 +294,7 @@
       zhi   = MIN( zhi   , zD    )
       zD_hi = MAX( 0._wp, zD-zhi )
-
-      ! Wave radiation
-      zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model
+ 
+     ! Wave radiation
+      zuwave = zua - zssu   ;   zvwave = zva - zssv   ! Use wind speed rel. to ocean for wave model
       zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
       !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
@@ -318,4 +321,28 @@
       IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
 
+      ! lateral velocities
+      ! default ssu and ssv
+      ! ln_M2016: mean velocity along the profile
+      IF ( ln_M2016 ) THEN
+         ! interpol needed data
+         CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )   ! 3d velocities
+        
+         !compute bottom level
+         CALL icb_utl_getkb( ikb, ze3t, zD )
+         
+         ! compute mean velocity 
+         zuo = 0.0 ; zvo = 0.0; zdep = 0.0
+         DO jk = 1, ikb-1
+            zuo = zuo + zuoce(jk) * ze3t(jk)
+            zvo = zvo + zvoce(jk) * ze3t(jk)
+            zdep = zdep + ze3t(jk)
+         END DO
+         zuo = (zuo + zuoce(ikb) * (zD - zdep)) / zD
+         zvo = (zvo + zvoce(ikb) * (zD - zdep)) / zD
+      ELSE
+         zuo = zssu
+         zvo = zssv
+      END IF
+
       zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
       !
@@ -330,7 +357,7 @@
          ! Explicit accelerations
          !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
-         !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
+         !    -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
          !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
-         !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
+         !    -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
          zaxe = -grav * zssh_x + zwave_rad * zuwave
          zaye = -grav * zssh_y + zwave_rad * zvwave
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90	(revision 13357)
@@ -78,13 +78,15 @@
       !
       !                          ! initialised variable with extra haloes to zero
-      uo_e(:,:) = 0._wp   ;   vo_e(:,:) = 0._wp   ;
-      ua_e(:,:) = 0._wp   ;   va_e(:,:) = 0._wp   ;
-      ff_e(:,:) = 0._wp   ;   tt_e(:,:) = 0._wp   ;
-      fr_e(:,:) = 0._wp   ;
-      ! prepration to Nacho work
-      ! e3t_e
-      ! tt3d_e
-      ! uo3d_e
-      ! vo3d_e
+      ssu_e(:,:) = 0._wp   ;   ssv_e(:,:) = 0._wp   ;
+      ua_e(:,:)  = 0._wp   ;   va_e(:,:)  = 0._wp   ;
+      ff_e(:,:)  = 0._wp   ;   sst_e(:,:)  = 0._wp   ;
+      fr_e(:,:)  = 0._wp   ;
+      !
+      IF ( ln_M2016 ) THEN
+         toce_e(:,:,:) = 0._wp
+         uoce_e(:,:,:) = 0._wp
+         voce_e(:,:,:) = 0._wp
+         e3t_e(:,:,:)  = 0._wp
+      END IF
       !
 #if defined key_si3
@@ -106,4 +108,10 @@
       first_width (:) = SQRT(  rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) )  )
       first_length(:) = rn_LoW_ratio * first_width(:)
+      rho_berg_1_oce  = rn_rho_bergs / pp_rho_seawater  ! scale factor used for convertion thickness to draft
+      !
+      ! deepest level affected by icebergs
+      ! can be tuned but the safest is this 
+      ! (with z* and z~ the depth of each level change overtime, so the more robust micbkb is jpk)
+      micbkb = jpk
 
       berg_grid%calving      (:,:)   = 0._wp
@@ -365,4 +373,5 @@
                localpt%uvel = 0._wp
                localpt%vvel = 0._wp
+               localpt%kb   = 1
                CALL icb_utl_incr()
                localberg%number(:) = num_bergs(:)
@@ -398,5 +407,5 @@
          &              rn_bits_erosion_fraction        , rn_sicn_shift       , ln_passive_mode      ,   &
          &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   &
-         &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      ,                          &
+         &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      , ln_M2016             ,   &
          &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out
       !!----------------------------------------------------------------------
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbstp.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbstp.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbstp.F90	(revision 13357)
@@ -52,5 +52,5 @@
 CONTAINS
 
-   SUBROUTINE icb_stp( kt )
+   SUBROUTINE icb_stp( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_stp  ***
@@ -61,4 +61,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       LOGICAL ::   ll_sample_traj, ll_budget, ll_verbose   ! local logical
@@ -92,5 +93,5 @@
       !
       !                                   !* copy nemo forcing arrays into iceberg versions with extra halo
-      CALL icb_utl_copy()                 ! only necessary for variables not on T points
+      CALL icb_utl_copy( Kmm )                 ! only necessary for variables not on T points
       !
       !
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90	(revision 13357)
@@ -48,10 +48,11 @@
       INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg
       !
-      INTEGER  ::   ii, ij
-      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn
-      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv
+      INTEGER  ::   ii, ij, jk, ikb
+      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb
+      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv
       REAL(wp) ::   zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
       REAL(wp) ::   zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
-      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2
+      REAL(wp) ::   zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw
+      REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, ztmp
       TYPE(iceberg), POINTER ::   this, next
       TYPE(point)  , POINTER ::   pt
@@ -83,16 +84,16 @@
          pt => this%current_point
          nknberg = this%number(1)
+
+         CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position
+             &                 pssu=pt%ssu, pua=pt%ua, &   ! oce/atm velocities
+             &                 pssv=pt%ssv, pva=pt%va, &   ! oce/atm velocities
+             &                 psst=pt%sst, pcn=pt%cn )
+
          IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN
             CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2,                 &
-               &                 puo=pt%uo, pui=pt%ui, pua=pt%ua, pssh_i=pt%ssh_x,   &
-               &                 pvo=pt%vo, pvi=pt%vi, pva=pt%va, pssh_j=pt%ssh_y,   &
-               &                 psst=pt%sst, pcn=pt%cn, phi=pt%hi,                  &
-               &                 plat=pt%lat, plon=pt%lon                            )
-         ELSE
-            CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position
-               &                 puo=pt%uo, pua=pt%ua,   &   ! oce/atm velocities
-               &                 pvo=pt%vo, pva=pt%va,   &   ! oce/atm velocities
-               &                 psst=pt%sst, pcn=pt%cn )    ! sst and ice concentration
-               ! preparation Nacho add t3d and uo, vo, to basal
+               &                 pui=pt%ui, pssh_i=pt%ssh_x, &
+               &                 pvi=pt%vi, pssh_j=pt%ssh_y, &
+               &                 phi=pt%hi,                  &
+               &                 plat=pt%lat, plon=pt%lon )
          END IF
          !
@@ -101,5 +102,5 @@
          zM   = pt%mass
          zT   = pt%thickness                               ! total thickness
-       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
+         zD   = rho_berg_1_oce * zT          ! draught (keel depth)
        ! F   = zT - D ! freeboard
          zW   = pt%width
@@ -114,12 +115,47 @@
 
          ! Environment
-         zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 )
-         zdva = SQRT( (pt%ua  -pt%uo)**2 + (pt%va  -pt%vo)**2 )
-         zSs  = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva                ! Sea state      (eqn M.A9)
-
+         ! default sst, ssu and ssv
+         ! ln_M2016: use temp, u and v profile
+         IF ( ln_M2016 ) THEN
+
+            ! load t, u, v and e3 profile at icb position
+            CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )
+            
+            !compute bottom level
+            CALL icb_utl_getkb( pt%kb, ze3t, zD )
+
+            ikb = MIN(pt%kb,mbkt(ii,ij))
+            ztb = ztoce(ikb)                                                ! basal temperature
+            zub = zuoce(ikb)
+            zvb = zvoce(ikb)
+         ELSE
+            ztb = pt%sst
+            zub = pt%ssu
+            zvb = pt%ssv
+         END IF
+
+         zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 )        ! relative basal velocity
+         zdva  = SQRT( (pt%ua  -pt%ssu)**2 + (pt%va  -pt%ssv)**2 )  ! relative wind
+         zSs   = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva              ! Sea state      (eqn M.A9)
+         !
          ! Melt rates in m/s (i.e. division by rday)
-         zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2)                    , 0._wp ) * z1_rday   ! Buoyant convection at sides (eqn M.A10)
-         zMb = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 )
-         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
+         IF ( ln_M2016 ) THEN
+            ! Buoyant convection at sides (eqn M.A10) but averaging along all the iceberg draft
+            ztmp(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
+            zMv = 0.0; zdepw = 0.0
+            DO jk = 1,ikb-1
+               zMv = zMv + ze3t(jk)*ztmp(jk)
+               zdepw = zdepw + ze3t(jk)
+            END DO
+            zMv = (zMv + (zD - zdepw)*ztmp(ikb)) / zD
+         ELSE
+            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
+         END IF
+
+         ! Basal turbulent melting     (eqn M.A7 ) but using the basal temperature
+         zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
+
+         ! Wave erosion                (eqn M.A8 )
+         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday
 
          IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
@@ -209,5 +245,5 @@
 
          ! Rolling
-         zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn       ! draught (keel depth)
+         zDn = rho_berg_1_oce * zTn       ! draught (keel depth)
          IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
             zT  = zTn
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbtrj.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbtrj.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbtrj.F90	(revision 13357)
@@ -40,5 +40,5 @@
    INTEGER ::   numberid, nstepid, nscaling_id
    INTEGER ::   nlonid, nlatid, nxid, nyid, nuvelid, nvvelid, nmassid
-   INTEGER ::   nuoid, nvoid, nuaid, nvaid, nuiid, nviid
+   INTEGER ::   nssuid, nssvid, nuaid, nvaid, nuiid, nviid
    INTEGER ::   nsshxid, nsshyid, nsstid, ncntid, nthkid
    INTEGER ::   nthicknessid, nwidthid, nlengthid
@@ -111,6 +111,6 @@
       iret = NF90_DEF_VAR( ntrajid, 'uvel'          , NF90_DOUBLE, n_dim          , nuvelid          )
       iret = NF90_DEF_VAR( ntrajid, 'vvel'          , NF90_DOUBLE, n_dim          , nvvelid          )
-      iret = NF90_DEF_VAR( ntrajid, 'uto'           , NF90_DOUBLE, n_dim          , nuoid            )
-      iret = NF90_DEF_VAR( ntrajid, 'vto'           , NF90_DOUBLE, n_dim          , nvoid            )
+      iret = NF90_DEF_VAR( ntrajid, 'ssu'           , NF90_DOUBLE, n_dim          , nssuid           )
+      iret = NF90_DEF_VAR( ntrajid, 'ssv'           , NF90_DOUBLE, n_dim          , nssvid           )
       iret = NF90_DEF_VAR( ntrajid, 'uta'           , NF90_DOUBLE, n_dim          , nuaid            )
       iret = NF90_DEF_VAR( ntrajid, 'vta'           , NF90_DOUBLE, n_dim          , nvaid            )
@@ -148,8 +148,8 @@
       iret = NF90_PUT_ATT( ntrajid, nvvelid         , 'long_name', 'meridional velocity' )
       iret = NF90_PUT_ATT( ntrajid, nvvelid         , 'units'    , 'm/s' )
-      iret = NF90_PUT_ATT( ntrajid, nuoid           , 'long_name', 'ocean u component' )
-      iret = NF90_PUT_ATT( ntrajid, nuoid           , 'units'    , 'm/s' )
-      iret = NF90_PUT_ATT( ntrajid, nvoid           , 'long_name', 'ocean v component' )
-      iret = NF90_PUT_ATT( ntrajid, nvoid           , 'units'    , 'm/s' )
+      iret = NF90_PUT_ATT( ntrajid, nssuid          , 'long_name', 'ocean u component' )
+      iret = NF90_PUT_ATT( ntrajid, nssuid          , 'units'    , 'm/s' )
+      iret = NF90_PUT_ATT( ntrajid, nssvid          , 'long_name', 'ocean v component' )
+      iret = NF90_PUT_ATT( ntrajid, nssvid          , 'units'    , 'm/s' )
       iret = NF90_PUT_ATT( ntrajid, nuaid           , 'long_name', 'atmosphere u component' )
       iret = NF90_PUT_ATT( ntrajid, nuaid           , 'units'    , 'm/s' )
@@ -231,6 +231,6 @@
          iret = NF90_PUT_VAR( ntrajid, nuvelid         , pt%uvel          , (/ jn /) )
          iret = NF90_PUT_VAR( ntrajid, nvvelid         , pt%vvel          , (/ jn /) )
-         iret = NF90_PUT_VAR( ntrajid, nuoid           , pt%uo            , (/ jn /) )
-         iret = NF90_PUT_VAR( ntrajid, nvoid           , pt%vo            , (/ jn /) )
+         iret = NF90_PUT_VAR( ntrajid, nssuid          , pt%ssu           , (/ jn /) )
+         iret = NF90_PUT_VAR( ntrajid, nssvid          , pt%ssv           , (/ jn /) )
          iret = NF90_PUT_VAR( ntrajid, nuaid           , pt%ua            , (/ jn /) )
          iret = NF90_PUT_VAR( ntrajid, nvaid           , pt%va            , (/ jn /) )
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90	(revision 13357)
@@ -19,4 +19,5 @@
    !!----------------------------------------------------------------------
    USE par_oce                             ! ocean parameters
+   USE oce,    ONLY: ts, uu, vv
    USE dom_oce                             ! ocean domain
    USE in_out_manager                      ! IO parameters
@@ -34,5 +35,10 @@
    PRIVATE
 
+   INTERFACE icb_utl_bilin_h
+      MODULE PROCEDURE icb_utl_bilin_2d_h, icb_utl_bilin_3d_h
+   END INTERFACE
+
    PUBLIC   icb_utl_copy          ! routine called in icbstp module
+   PUBLIC   icb_utl_getkb         ! routine called in icbdyn and icbthm modules
    PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
    PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module
@@ -56,5 +62,5 @@
 CONTAINS
 
-   SUBROUTINE icb_utl_copy()
+   SUBROUTINE icb_utl_copy( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_copy  ***
@@ -64,21 +70,27 @@
       !! ** Method  : - blah blah
       !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1) :: ztmp
 #if defined key_si3
       REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m    !    ocean surface (ssh_m) if ice is not embedded
       !                                              !    ocean surface in leads if ice is embedded   
 #endif
+      INTEGER :: jk   ! vertical loop index
+      INTEGER :: Kmm  ! ocean time levelindex
+      !
       ! copy nemo forcing arrays into iceberg versions with extra halo
       ! only necessary for variables not on T points
       ! and ssh which is used to calculate gradients
-
-      uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
-      vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
-      tt_e(1:jpi,1:jpj) = sst_m(:,:)
-      fr_e(1:jpi,1:jpj) = fr_i (:,:)
-      ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-      va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-      !
-      CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 )
+      !
+      ! surface forcing
+      !
+      ssu_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
+      ssv_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
+      sst_e(1:jpi,1:jpj) = sst_m(:,:)
+      fr_e (1:jpi,1:jpj) = fr_i (:,:)
+      ua_e (1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
+      va_e (1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
+      !
+      CALL lbc_lnk_icb( 'icbutl', ssu_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( 'icbutl', ssv_e, 'V', -1._wp, 1, 1 )
       CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 )
       CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 )
@@ -98,10 +110,29 @@
 #endif
       !
+      ! (PM) could be improve with a 3d lbclnk gathering both variables
+      ! should be done once extra haloe generalised
+      IF ( ln_M2016 ) THEN
+         DO jk = 1,jpk
+            ! uoce
+            ztmp(1:jpi,1:jpj) = uu(:,:,jk,Kmm)
+            CALL lbc_lnk_icb( 'icbutl', ztmp, 'U', -1._wp, 1, 1 )
+            uoce_e(:,:,jk) = ztmp(:,:)
+            !
+            ! voce
+            ztmp(1:jpi,1:jpj) = vv(:,:,jk,Kmm)
+            CALL lbc_lnk_icb( 'icbutl', ztmp, 'V', -1._wp, 1, 1 )
+            voce_e(:,:,jk) = ztmp(:,:)
+         END DO
+         toce_e(1:jpi,1:jpj,1:jpk) = ts(:,:,:,1,Kmm)
+         e3t_e (1:jpi,1:jpj,1:jpk) = e3t(:,:,:,Kmm)
+      END IF
+      !
    END SUBROUTINE icb_utl_copy
 
 
-   SUBROUTINE icb_utl_interp( pi, pj, pe1, puo, pui, pua, pssh_i,   &
-      &                               pe2, pvo, pvi, pva, pssh_j,   &
-      &                               psst, pcn, phi, pff, plon, plat)
+   SUBROUTINE icb_utl_interp( pi, pj, pe1, pssu, pui, pua, pssh_i,     &
+      &                               pe2, pssv, pvi, pva, pssh_j,     &
+      &                               psst, pcn, phi, pff, plon, plat, &
+      &                               ptoce, puoce, pvoce, pe3t       )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_interp  ***
@@ -122,8 +153,9 @@
       REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential
       REAL(wp), INTENT(  out), OPTIONAL ::   pe1, pe2                       ! i- and j scale factors
-      REAL(wp), INTENT(  out), OPTIONAL ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds
+      REAL(wp), INTENT(  out), OPTIONAL ::   pssu, pssv, pui, pvi, pua, pva   ! ocean, ice and wind speeds
       REAL(wp), INTENT(  out), OPTIONAL ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
       REAL(wp), INTENT(  out), OPTIONAL ::   psst, pcn, phi, pff            ! SST, ice concentration, ice thickness, Coriolis
       REAL(wp), INTENT(  out), OPTIONAL ::   plat, plon                     ! position
+      REAL(wp), DIMENSION(jpk), INTENT(  out), OPTIONAL ::   ptoce, puoce, pvoce, pe3t   ! 3D variables
       !
       REAL(wp), DIMENSION(4) :: zwT  , zwU  , zwV  , zwF   ! interpolation weight
@@ -147,9 +179,9 @@
       IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. )
       !
-      IF ( PRESENT(puo ) ) puo  = icb_utl_bilin_h( uo_e, iiU, ijU, zwU        , .false. )         ! ocean velocities
-      IF ( PRESENT(pvo ) ) pvo  = icb_utl_bilin_h( vo_e, iiV, ijV, zwV        , .false. )         !
-      IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( tt_e, iiT, ijT, zwT * zmskT, .false. ) ! sst
-      IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e, iiT, ijT, zwT * zmskT, .false. ) ! ice concentration
-      IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e, iiF, ijF, zwF        , .false. )         ! Coriolis parameter
+      IF ( PRESENT(pssu) ) pssu = icb_utl_bilin_h( ssu_e, iiU, ijU, zwU        , .false. )         ! ocean velocities
+      IF ( PRESENT(pssv) ) pssv = icb_utl_bilin_h( ssv_e, iiV, ijV, zwV        , .false. )         !
+      IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( sst_e, iiT, ijT, zwT * zmskT, .false. ) ! sst
+      IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration
+      IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e , iiF, ijF, zwF        , .false. )         ! Coriolis parameter
       !
       IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN
@@ -188,4 +220,12 @@
             &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe2 )
       END IF
+      !
+      ! 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 )
+      END IF
+      IF ( PRESENT(ptoce) ) ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zwT * zmskT )
+      IF ( PRESENT(pe3t)  ) pe3t(:)  = e3t_e(iiT,ijT,:)    ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results
       !
    END SUBROUTINE icb_utl_interp
@@ -291,5 +331,5 @@
    END SUBROUTINE icb_utl_pos
 
-   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pii, pij, pw, pllon )
+   REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon )
       !!----------------------------------------------------------------------
       !!                  ***  FUNCTION icb_utl_bilin  ***
@@ -321,9 +361,39 @@
       !
       ! compute interpolated value
-      icb_utl_bilin_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)) 
-      !
-      IF( pllon .AND. icb_utl_bilin_h > 180._wp ) icb_utl_bilin_h = icb_utl_bilin_h - 360._wp
-      !
-   END FUNCTION icb_utl_bilin_h
+      icb_utl_bilin_2d_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)) 
+      !
+      IF( pllon .AND. icb_utl_bilin_2d_h > 180._wp ) icb_utl_bilin_2d_h = icb_utl_bilin_2d_h - 360._wp
+      !
+   END FUNCTION icb_utl_bilin_2d_h
+
+   FUNCTION icb_utl_bilin_3d_h( pfld, pii, pij, pw )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION icb_utl_bilin  ***
+      !!
+      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
+      !!                this version deals with extra halo points
+      !!
+      !!       !!gm  CAUTION an optional argument should be added to handle
+      !!             the slip/no-slip conditions  ==>>> to be done later
+      !!
+      !!----------------------------------------------------------------------
+      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
+      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
+      !!----------------------------------------------------------------------
+      !
+      ! data
+      zdat(1,:) = pfld(pii  ,pij  ,:)
+      zdat(2,:) = pfld(pii+1,pij  ,:)
+      zdat(3,:) = pfld(pii  ,pij+1,:)
+      zdat(4,:) = pfld(pii+1,pij+1,:)
+      !
+      ! 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)) 
+      !
+   END FUNCTION icb_utl_bilin_3d_h
 
    REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
@@ -405,4 +475,26 @@
    END FUNCTION icb_utl_bilin_e
 
+   SUBROUTINE icb_utl_getkb( kb, pe3, pD )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_getkb         ***
+      !!
+      !! ** Purpose :   compute the latest level affected by icb
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(out)               :: kb
+      REAL(wp), DIMENSION(:), INTENT(in) :: pe3
+      REAL(wp),               INTENT(in) :: pD
+      !!
+      INTEGER  :: jk
+      REAL(wp) :: zdepw
+      !!
+      zdepw = 0.0
+      kb = 1
+      DO WHILE ( zdepw <  pD)
+         zdepw = zdepw + pe3(kb)
+         kb = kb + 1
+      END DO
+      kb = kb - 1
+   END SUBROUTINE
 
    SUBROUTINE icb_utl_add( bergvals, ptvals )
@@ -593,5 +685,5 @@
       WRITE(numicb, 9200) kt, berg%number(1), &
                    pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
-                   pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
+                   pt%ssu, pt%ssv, pt%ua, pt%va, pt%ui, pt%vi
       CALL flush( numicb )
  9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
@@ -619,5 +711,5 @@
          WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
          WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
-            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
+            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','ssu,ssv','ua,va','ui,vi'
       ENDIF
       DO WHILE( ASSOCIATED(this) )
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcmod.F90	(revision 13356)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcmod.F90	(revision 13357)
@@ -456,5 +456,5 @@
 
       IF( ln_icebergs    )   THEN
-                                     CALL icb_stp( kt )           ! compute icebergs
+                                     CALL icb_stp( kt, Kmm )           ! compute icebergs
          ! Icebergs do not melt over the haloes. 
          ! So emp values over the haloes are no more consistent with the inner domain values. 
