Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SHARED/namelist_ref	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SHARED/namelist_ref	(revision 14031)
@@ -629,5 +629,8 @@
    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)
+      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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icb_oce.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icb_oce.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icb_oce.F90	(revision 14031)
@@ -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, sss    ! properties of iceberg environment 
+      REAL(wp) ::   ssu, ssv, ui, vi, ua, va, ssh_x, ssh_y, sst, sss, cn, hi  ! properties of iceberg environment 
       REAL(wp) ::   mass_of_bits, heat_density
+      INTEGER  ::   kb                                                   ! icb bottom level
    END TYPE point
 
@@ -85,9 +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 ::   ff_e, tt_e, fr_e, ss_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssu_e, ssv_e
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   sst_e, sss_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
+   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
@@ -117,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
@@ -124,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, ln_icb_grd            !: use Nacho's Merino 2016 work
    !
    ! restart
@@ -135,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
@@ -171,6 +178,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) ,                            &
@@ -178,7 +185,7 @@
          &      hi_e(0:jpi+1,0:jpj+1) ,                            &
 #endif
-         &      ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1)  ,   &
-         &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   &
-         &      ss_e(0:jpi+1,0:jpj+1) ,                            & 
+         &      fr_e(0:jpi+1,0:jpj+1) ,                            &
+         &      sst_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,  &
+         &      sss_e(0:jpi+1,0:jpj+1) ,                           & 
          &      first_width(nclasses) , first_length(nclasses) ,   &
          &      src_calving (jpi,jpj) ,                            &
@@ -186,6 +193,12 @@
       icb_alloc = icb_alloc + ill
 
+      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), &
-         &      STAT=ill)
+         &      rlon_e(0:jpi+1,0:jpj+1) , rlat_e(0:jpi+1,0:jpj+1) , ff_e(0:jpi+1,0:jpj+1)   , STAT=ill)
       icb_alloc = icb_alloc + ill
 
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbclv.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbclv.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbclv.F90	(revision 14031)
@@ -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)
@@ -172,9 +172,4 @@
       END DO
       !
-      DO jn = 1, nclasses
-         CALL lbc_lnk( 'icbclv', berg_grid%stored_ice(:,:,jn), 'T', 1._wp )
-      END DO
-      CALL lbc_lnk( 'icbclv', berg_grid%stored_heat, 'T', 1._wp )
-      !
       IF( nn_verbose_level > 0 .AND. icntmax > 1 )   WRITE(numicb,*) 'icb_clv: icnt=', icnt,' on', narea
       !
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbdyn.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbdyn.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbdyn.F90	(revision 14031)
@@ -14,4 +14,5 @@
    USE dom_oce        ! NEMO ocean domain
    USE phycst         ! NEMO physical constants
+   USE in_out_manager                      ! IO parameters
    !
    USE icb_oce        ! define iceberg arrays
@@ -97,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)
@@ -113,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)
@@ -129,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)
@@ -148,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
@@ -156,8 +157,4 @@
          pt%yj   = zyj_n
 
-         ! update actual position
-         pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
-         pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
-
          berg => berg%next                         ! switch to the next berg
          !
@@ -167,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  ***
@@ -177,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
@@ -184,5 +183,9 @@
       INTEGER  ::   ii, ii0
       INTEGER  ::   ij, ij0
+      INTEGER  ::   ikb
       INTEGER  ::   ibounce_method
+      !
+      REAL(wp) :: zD 
+      REAL(wp), DIMENSION(jpk) :: ze3t
       !!----------------------------------------------------------------------
       !
@@ -200,5 +203,23 @@
       ij  = mj1( ij  )
       !
-      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
+      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
+      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
+         ! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0)
+         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
@@ -257,23 +278,27 @@
       REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
       !
-      INTEGER  ::   itloop
-      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss
-      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, pe1, zuo, zui, zua, zssh_x,                     &
-         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss )
+      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
@@ -282,7 +307,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.
@@ -309,4 +334,22 @@
       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 
+         CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb)
+         CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb)
+      ELSE
+         zuo = zssu
+         zvo = zssv
+      END IF
+
       zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
       !
@@ -321,7 +364,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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbini.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbini.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbini.F90	(revision 14031)
@@ -73,13 +73,21 @@
       !
       IF( .NOT. ln_icebergs )   RETURN
-
+      !
       !                          ! allocate gridded fields
       IF( icb_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'icb_alloc : unable to allocate arrays' )
       !
       !                          ! 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   ;   ss_e(:,:) = 0._wp   ;
+      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   ;   sss_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
       hi_e(:,:) = 0._wp   ;
@@ -100,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
@@ -240,7 +254,17 @@
       vmask_e(:,:) = 0._wp   ;   vmask_e(1:jpi,1:jpj) = vmask(:,:,1)
       CALL lbc_lnk_icb( 'icbini', tmask_e, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbini', umask_e, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbini', vmask_e, 'T', +1._wp, 1, 1 )
-      !
+      CALL lbc_lnk_icb( 'icbini', umask_e, 'U', +1._wp, 1, 1 )
+      CALL lbc_lnk_icb( 'icbini', vmask_e, 'V', +1._wp, 1, 1 )
+
+      ! definition of extended lat/lon array needed by icb_bilin_h
+      rlon_e(:,:) = 0._wp     ;  rlon_e(1:jpi,1:jpj) = glamt(:,:) 
+      rlat_e(:,:) = 0._wp     ;  rlat_e(1:jpi,1:jpj) = gphit(:,:)
+      CALL lbc_lnk_icb( 'icbini', rlon_e, 'T', +1._wp, 1, 1 )
+      CALL lbc_lnk_icb( 'icbini', rlat_e, 'T', +1._wp, 1, 1 )
+      !
+      ! definnitionn of extennded ff_f array needed by icb_utl_interp
+      ff_e(:,:) = 0._wp       ;  ff_e(1:jpi,1:jpj) = ff_f(:,:)
+      CALL lbc_lnk_icb( 'icbini', ff_e, 'F', +1._wp, 1, 1 )
+
       ! assign each new iceberg with a unique number constructed from the processor number
       ! and incremented by the total number of processors
@@ -338,6 +362,5 @@
                localpt%xi = REAL( mig(ji), wp )
                localpt%yj = REAL( mjg(jj), wp )
-               localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T' )
-               localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T' )
+               CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon )   
                localpt%mass      = rn_initial_mass     (iberg)
                localpt%thickness = rn_initial_thickness(iberg)
@@ -350,4 +373,5 @@
                localpt%uvel = 0._wp
                localpt%vvel = 0._wp
+               localpt%kb   = 1
                CALL icb_utl_incr()
                localberg%number(:) = num_bergs(:)
@@ -383,6 +407,7 @@
          &              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      ,                          &
-         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out
+         &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      , ln_M2016             ,   &
+         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out        ,   &
+         &              ln_icb_grd
       !!----------------------------------------------------------------------
 
@@ -463,4 +488,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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbstp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbstp.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbstp.F90	(revision 14031)
@@ -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
@@ -70,4 +71,7 @@
       !
       nktberg = kt
+      !
+      !CALL test_icb_utl_getkb
+      !CALL ctl_stop('end test icb')
       !
       IF( nn_test_icebergs < 0 .OR. ln_use_calving ) THEN !* read calving data
@@ -92,5 +96,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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbthm.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbthm.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbthm.F90	(revision 14031)
@@ -49,11 +49,12 @@
       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
+      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) ::   zSSS, zfzpt
-      REAL(wp) ::   zMv, zMe, zMb, zmelt, zdvo, 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, zzMv
       TYPE(iceberg), POINTER ::   this, next
       TYPE(point)  , POINTER ::   pt
@@ -85,7 +86,18 @@
          pt => this%current_point
          nknberg = this%number(1)
-         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   &
-            &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   &
-            &                 pt%sst, pt%cn, pt%hi, zff, pt%sss )
+
+         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, &
+             &                 psss=pt%sss             )
+
+         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,                 &
+               &                 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
          !
          zSST = pt%sst
@@ -95,6 +107,5 @@
          zM   = pt%mass
          zT   = pt%thickness                               ! total thickness
-       ! D   = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth)
-       ! F   = zT - D ! freeboard
+         zD   = rho_berg_1_oce * zT                        ! draught (keel depth)
          zW   = pt%width
          zL   = pt%length
@@ -108,16 +119,48 @@
 
          ! 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))                             ! limit pt%kb by mbkt 
+                                                                     ! => bottom temperature used to fill ztoce(mbkt:jpk)
+            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)
+         ! Buoyant convection at sides (eqn M.A10)
+         IF ( ln_M2016 ) THEN
+            ! averaging along all the iceberg draft
+            zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
+            CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb )
+         ELSE
+            zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
+         END IF
+         !
+         ! Basal turbulent melting     (eqn M.A7 )
          IF ( zSST > zfzpt ) THEN                                                              ! Calculate basal melting only if SST above freezing point  
-            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 )
+            zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
          ELSE
             zMb = 0._wp                                                                        ! No basal melting if SST below freezing point     
          ENDIF
-         zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3)))     , 0._wp ) * z1_rday      ! Wave erosion                (eqn M.A8 )
+         !
+         ! 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
@@ -207,5 +250,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
@@ -224,4 +267,5 @@
 
 !!gm  add a test to avoid over melting ?
+!!pm  I agree, over melting could break conservation (more melt than calving)
 
          IF( zMnew <= 0._wp ) THEN       ! Delete the berg if completely melted
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbtrj.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbtrj.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbtrj.F90	(revision 14031)
@@ -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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbutl.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbutl.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbutl.F90	(revision 14031)
@@ -8,12 +8,16 @@
    !!            -    !                            Removal of mapping from another grid
    !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!           4.2   !  2020-07  (P. Mathiot)     simplification of interpolation routine
+   !!                 !                            and add Nacho Merino work
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   icb_utl_interp   :
-   !!   icb_utl_bilin    :
-   !!   icb_utl_bilin_e  :
+   !!   icb_utl_pos      : compute bottom left corner indice, weight and mask
+   !!   icb_utl_bilin_h  : interpolation field to icb position
+   !!   icb_utl_bilin_e  : interpolation of scale factor to icb position
    !!----------------------------------------------------------------------
    USE par_oce                             ! ocean parameters
+   USE oce,    ONLY: ts, uu, vv
    USE dom_oce                             ! ocean domain
    USE in_out_manager                      ! IO parameters
@@ -31,8 +35,14 @@
    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   test_icb_utl_getkb    ! routine called in icbdyn and icbthm modules
+   PUBLIC   icb_utl_zavg          ! routine called in icbdyn and icbthm modules
    PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
-   PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules
-   PUBLIC   icb_utl_bilin_x       ! routine called in icbdyn module
+   PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module
    PUBLIC   icb_utl_add           ! routine called in icbini.F90, icbclv, icblbc and icbrst modules
    PUBLIC   icb_utl_delete        ! routine called in icblbc, icbthm modules
@@ -54,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE icb_utl_copy()
+   SUBROUTINE icb_utl_copy( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_copy  ***
@@ -62,29 +72,31 @@
       !! ** 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)
+      !
+      ! 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(:,:)
+      sss_e(1:jpi,1:jpj) = sss_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
       ff_e(1:jpi,1:jpj) = ff_f (:,:) 
-      tt_e(1:jpi,1:jpj) = sst_m(:,:)
-      ss_e(1:jpi,1:jpj) = sss_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 )
-      CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +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 )
-      CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', ss_e, 'T', +1._wp, 1, 1 )
+      !
+      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 )
 #if defined key_si3
       hi_e(1:jpi, 1:jpj) = hm_i (:,:)  
@@ -96,18 +108,35 @@
       ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1)
       !
-      CALL lbc_lnk_icb( 'icbutl', hi_e , 'T', +1._wp, 1, 1 )
       CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 )
       CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 )
 #else
-      ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)
+      ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)         
 #endif
-      CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 )
+      !
+      ! (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, pe1, puo, pui, pua, pssh_i,   &
-      &                       pj, pe2, pvo, pvi, pva, pssh_j,   &
-      &                       psst, pcn, phi, pff, psss        )
+   SUBROUTINE icb_utl_interp( pi, pj, pe1 , pssu, pui, pua, pssh_i,         &
+      &                               pe2 , pssv, pvi, pva, pssh_j,         &
+      &                               psst, psss, pcn, phi, pff   ,         &
+      &                               plon, plat, ptoce, puoce, pvoce, pe3t )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_interp  ***
@@ -127,49 +156,100 @@
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential
-      REAL(wp), INTENT(  out) ::   pe1, pe2                       ! i- and j scale factors
-      REAL(wp), INTENT(  out) ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds
-      REAL(wp), INTENT(  out) ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
-      REAL(wp), INTENT(  out) ::   psst, pcn, phi, pff, psss      ! SST, ice concentration, ice thickness, Coriolis, SSS
-      !
+      REAL(wp), INTENT(  out), OPTIONAL ::   pe1, pe2                       ! i- and j scale factors
+      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, psss, pcn, phi, pff      ! SST, SSS, 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
+      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
       REAL(wp) ::   zcd, zmod       ! local scalars
       !!----------------------------------------------------------------------
-
-      pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors
-      pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
-      !
-      puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U', .false.  )    ! ocean velocities
-      pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false.  )
-      psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true.   )    ! SST
-      psss = icb_utl_bilin_h( ss_e, pi, pj, 'T', .true.   )    ! SSS
-      pcn  = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true.   )    ! ice concentration
-      pff  = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false.  )    ! Coriolis parameter
-      !
-      pua  = icb_utl_bilin_h( ua_e, pi, pj, 'U', .true.   )    ! 10m wind
-      pva  = icb_utl_bilin_h( va_e, pi, pj, 'V', .true.   )    ! here (ua,va) are stress => rough conversion from stress to speed
-      zcd  = 1.22_wp * 1.5e-3_wp                               ! air density * drag coefficient 
-      zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  )
-      pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0
-      pva  = pva * zmod
-
+      !
+      ! get position, weight and mask 
+      CALL icb_utl_pos( pi, pj, 'T', iiT, ijT, zwT, zmskT )
+      CALL icb_utl_pos( pi, pj, 'U', iiU, ijU, zwU, zmskU )
+      CALL icb_utl_pos( pi, pj, 'V', iiV, ijV, zwV, zmskV )
+      CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF )
+      !
+      ! metrics and coordinates
+      IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors
+      IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
+      IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true.  )
+      IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. )
+      !
+      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(psss) ) psss = icb_utl_bilin_h( sss_e, iiT, ijT, zwT * zmskT, .false. ) ! sss
+      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
+         pua  = icb_utl_bilin_h( ua_e, iiU, ijU, zwU * zmskU, .false. ) ! 10m wind
+         pva  = icb_utl_bilin_h( va_e, iiV, ijV, zwV * zmskV, .false. ) ! here (ua,va) are stress => rough conversion from stress to speed
+         zcd  = 1.22_wp * 1.5e-3_wp                               ! air density * drag coefficient 
+         zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  )
+         pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0
+         pva  = pva * zmod
+      END IF
+      !
 #if defined key_si3
-      pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. )    ! sea-ice velocities
-      pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. )
-      phi = icb_utl_bilin_h( hi_e , pi, pj, 'T', .true.  )    ! ice thickness
+      IF ( PRESENT(pui) ) pui = icb_utl_bilin_h( ui_e , iiU, ijU, zwU        , .false. ) ! sea-ice velocities
+      IF ( PRESENT(pvi) ) pvi = icb_utl_bilin_h( vi_e , iiV, ijV, zwV        , .false. )
+      IF ( PRESENT(phi) ) phi = icb_utl_bilin_h( hi_e , iiT, ijT, zwT * zmskT, .false. ) ! ice thickness
 #else
-      pui = 0._wp
-      pvi = 0._wp
-      phi = 0._wp
+      IF ( PRESENT(pui) ) pui = 0._wp
+      IF ( PRESENT(pvi) ) pvi = 0._wp
+      IF ( PRESENT(phi) ) phi = 0._wp
 #endif
-
+      !
       ! Estimate SSH gradient in i- and j-direction (centred evaluation)
-      pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T', .true. ) -   &
-         &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T', .true. )  ) / ( 0.2_wp * pe1 )
-      pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T', .true. ) -   &
-         &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T', .true. )  ) / ( 0.2_wp * pe2 )
+      IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN
+         CALL icb_utl_pos( pi+0.1, pj    , 'T', iiTp, ijTp, zwTp, zmskTp )
+         CALL icb_utl_pos( pi-0.1, pj    , 'T', iiTm, ijTm, zwTm, zmskTm )
+         !
+         IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )
+         pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
+            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe1 )
+         !
+         CALL icb_utl_pos( pi    , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp )
+         CALL icb_utl_pos( pi    , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm )
+         !
+         IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
+         pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
+            &       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
+         ! 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) ) 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
       !
    END SUBROUTINE icb_utl_interp
 
-
-   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask )
+   SUBROUTINE icb_utl_pos( pi, pj, cd_type, kii, kij, pw, pmsk )
       !!----------------------------------------------------------------------
       !!                  ***  FUNCTION icb_utl_bilin  ***
@@ -182,13 +262,12 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
-      CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
-      LOGICAL                             , INTENT(in) ::   plmask    ! special treatment of mask point
-      !
-      INTEGER  ::   ii, ij   ! local integer
-      REAL(wp) ::   zi, zj   ! local real
-      REAL(wp) :: zw1, zw2, zw3, zw4
-      REAL(wp), DIMENSION(4) :: zmask
+      REAL(wp)              , INTENT(IN)  ::   pi, pj    ! targeted coordinates in (i,j) referential
+      CHARACTER(len=1)      , INTENT(IN)  ::   cd_type   ! point type
+      REAL(wp), DIMENSION(4), INTENT(OUT) ::   pw, pmsk  ! weight and mask
+      INTEGER ,               INTENT(OUT) ::   kii, kij  ! bottom left corner position in local domain
+      !
+      REAL(wp) :: zwi, zwj ! distance to bottom left corner
+      INTEGER  :: ierr 
+      !
       !!----------------------------------------------------------------------
       !
@@ -198,73 +277,85 @@
          ! since we're looking for four T points containing quadrant we're in of 
          ! current T cell
-         ii = MAX(0, INT( pi     ))
-         ij = MAX(0, INT( pj     ))    ! T-point
-         zi = pi - REAL(ii,wp)
-         zj = pj - REAL(ij,wp)
+         kii = MAX(0, INT( pi        ))
+         kij = MAX(0, INT( pj        ))    ! T-point
+         zwi = pi - REAL(kii,wp)
+         zwj = pj - REAL(kij,wp)
       CASE ( 'U' )
-         ii = MAX(0, INT( pi-0.5_wp ))
-         ij = MAX(0, INT( pj     ))    ! U-point
-         zi = pi - 0.5_wp - REAL(ii,wp)
-         zj = pj - REAL(ij,wp)
+         kii = MAX(0, INT( pi-0.5_wp ))
+         kij = MAX(0, INT( pj        ))    ! U-point
+         zwi = pi - 0.5_wp - REAL(kii,wp)
+         zwj = pj - REAL(kij,wp)
       CASE ( 'V' )
-         ii = MAX(0, INT( pi     ))
-         ij = MAX(0, INT( pj-0.5_wp ))    ! V-point
-         zi = pi - REAL(ii,wp)
-         zj = pj - 0.5_wp - REAL(ij,wp)
+         kii = MAX(0, INT( pi        ))
+         kij = MAX(0, INT( pj-0.5_wp ))    ! V-point
+         zwi = pi - REAL(kii,wp)
+         zwj = pj - 0.5_wp - REAL(kij,wp)
       CASE ( 'F' )
-         ii = MAX(0, INT( pi-0.5_wp ))
-         ij = MAX(0, INT( pj-0.5_wp ))    ! F-point
-         zi = pi - 0.5_wp - REAL(ii,wp)
-         zj = pj - 0.5_wp - REAL(ij,wp)
+         kii = MAX(0, INT( pi-0.5_wp ))
+         kij = MAX(0, INT( pj-0.5_wp ))    ! F-point
+         zwi = pi - 0.5_wp - REAL(kii,wp)
+         zwj = pj - 0.5_wp - REAL(kij,wp)
       END SELECT
+      !
+      ! compute weight
+      pw(1) = (1._wp-zwi) * (1._wp-zwj)
+      pw(2) =        zwi  * (1._wp-zwj)
+      pw(3) = (1._wp-zwi) *        zwj
+      pw(4) =        zwi  *        zwj
+      !
+      ! find position in this processor. Prevent near edge problems (see #1389)
+      !
+      IF (TRIM(cd_type) == 'T' ) THEN
+         ierr = 0
+         IF    ( kii <  mig( 1 ) ) THEN   ;  ierr = ierr + 1
+         ELSEIF( kii >= mig(jpi) ) THEN   ;  ierr = ierr + 1
+         ENDIF
+         !
+         IF    ( kij <  mjg( 1 ) ) THEN   ;   ierr = ierr + 1
+         ELSEIF( kij >= mjg(jpj) ) THEN   ;   ierr = ierr + 1
+         ENDIF
+         !
+         IF ( ierr > 0 ) THEN
+            WRITE(numout,*) 'bottom left corner T point out of bound'
+            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
+      END IF
       !
       ! find position in this processor. Prevent near edge problems (see #1389)
       ! (PM) will be useless if extra halo is used in NEMO
       !
-      IF    ( ii <= mig(1)-1 ) THEN   ;   ii = 0
-      ELSEIF( ii  > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                            ;   ii = mi1(ii)
+      IF    ( kii <= mig(1)-1 ) THEN   ;   kii = 0
+      ELSEIF( kii  > mig(jpi) ) THEN   ;   kii = jpi
+      ELSE                             ;   kii = mi1(kii)
       ENDIF
-      IF    ( ij <= mjg(1)-1 ) THEN   ;   ij = 0
-      ELSEIF( ij  > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                            ;   ij = mj1(ij)
+      IF    ( kij <= mjg(1)-1 ) THEN   ;   kij = 0
+      ELSEIF( kij  > mjg(jpj) ) THEN   ;   kij = jpj
+      ELSE                             ;   kij = mj1(kij)
       ENDIF
       !
       ! define mask array 
-      IF (plmask) THEN
-         ! land value is not used in the interpolation
-         SELECT CASE ( cd_type )
-         CASE ( 'T' )
-            zmask = (/tmask_e(ii,ij), tmask_e(ii+1,ij), tmask_e(ii,ij+1), tmask_e(ii+1,ij+1)/)
-         CASE ( 'U' )
-            zmask = (/umask_e(ii,ij), umask_e(ii+1,ij), umask_e(ii,ij+1), umask_e(ii+1,ij+1)/)
-         CASE ( 'V' )
-            zmask = (/vmask_e(ii,ij), vmask_e(ii+1,ij), vmask_e(ii,ij+1), vmask_e(ii+1,ij+1)/)
-         CASE ( 'F' )
-            ! F case only used for coriolis, ff_f is not mask so zmask = 1
-            zmask = 1.
-         END SELECT
-      ELSE
-         ! land value is used during interpolation
-         zmask = 1.
-      END iF
-      !
-      ! compute weight
-      zw1 = zmask(1) * (1._wp-zi) * (1._wp-zj)
-      zw2 = zmask(2) *        zi  * (1._wp-zj)
-      zw3 = zmask(3) * (1._wp-zi) *        zj
-      zw4 = zmask(4) *        zi  *        zj
-      !
-      ! compute interpolated value
-      icb_utl_bilin_h = ( pfld(ii,ij)*zw1 + pfld(ii+1,ij)*zw2 + pfld(ii,ij+1)*zw3 + pfld(ii+1,ij+1)*zw4 ) / MAX(1.e-20, zw1+zw2+zw3+zw4) 
-      !
-   END FUNCTION icb_utl_bilin_h
-
-
-   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
+      ! land value is not used in the interpolation
+      SELECT CASE ( cd_type )
+      CASE ( 'T' )
+         pmsk = (/tmask_e(kii,kij), tmask_e(kii+1,kij), tmask_e(kii,kij+1), tmask_e(kii+1,kij+1)/)
+      CASE ( 'U' )
+         pmsk = (/umask_e(kii,kij), umask_e(kii+1,kij), umask_e(kii,kij+1), umask_e(kii+1,kij+1)/)
+      CASE ( 'V' )
+         pmsk = (/vmask_e(kii,kij), vmask_e(kii+1,kij), vmask_e(kii,kij+1), vmask_e(kii+1,kij+1)/)
+      CASE ( 'F' )
+         ! F case only used for coriolis, ff_f is not mask so zmask = 1
+         pmsk = 1.
+      END SELECT
+   END SUBROUTINE icb_utl_pos
+
+   REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon )
       !!----------------------------------------------------------------------
       !!                  ***  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
@@ -272,63 +363,35 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
-      CHARACTER(len=1)            , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
-      !
-      INTEGER  ::   ii, ij   ! local integer
-      REAL(wp) ::   zi, zj   ! local real
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( cd_type )
-         CASE ( 'T' )
-            ! note that here there is no +0.5 added
-            ! since we're looking for four T points containing quadrant we're in of 
-            ! current T cell
-            ii = MAX(1, INT( pi     ))
-            ij = MAX(1, INT( pj     ))    ! T-point
-            zi = pi - REAL(ii,wp)
-            zj = pj - REAL(ij,wp)
-         CASE ( 'U' )
-            ii = MAX(1, INT( pi-0.5 ))
-            ij = MAX(1, INT( pj     ))    ! U-point
-            zi = pi - 0.5 - REAL(ii,wp)
-            zj = pj - REAL(ij,wp)
-         CASE ( 'V' )
-            ii = MAX(1, INT( pi     ))
-            ij = MAX(1, INT( pj-0.5 ))    ! V-point
-            zi = pi - REAL(ii,wp)
-            zj = pj - 0.5 - REAL(ij,wp)
-         CASE ( 'F' )
-            ii = MAX(1, INT( pi-0.5 ))
-            ij = MAX(1, INT( pj-0.5 ))    ! F-point
-            zi = pi - 0.5 - REAL(ii,wp)
-            zj = pj - 0.5 - REAL(ij,wp)
-      END SELECT
-      !
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                           ;   ii = mi1(ii)
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp), DIMENSION(4)              , INTENT(in) ::   pw        ! weight
+      LOGICAL                             , INTENT(in) ::   pllon     ! input data is a longitude
+      INTEGER ,                             INTENT(in) ::   pii, pij  ! bottom left corner
+      !
+      REAL(wp), DIMENSION(4) :: 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)
+      !
+      IF( pllon .AND. MAXVAL(zdat) - MINVAL(zdat) > 90._wp ) THEN
+         WHERE( zdat < 0._wp ) zdat = zdat + 360._wp
       ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi )   ii = ii-1      
-      IF( ij == jpj )   ij = ij-1
-      !
-      icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
-         &          + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
-      !
-   END FUNCTION icb_utl_bilin
-
-
-   REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj )
-      !!----------------------------------------------------------------------
-      !!                  ***  FUNCTION icb_utl_bilin_x  ***
+      !
+      ! compute interpolated value
+      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
-      !!                Special case for interpolating longitude
+      !!                this version deals with extra halo points
       !!
       !!       !!gm  CAUTION an optional argument should be added to handle
@@ -336,48 +399,26 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
-      !
-      INTEGER                                  ::   ii, ij   ! local integer
-      REAL(wp)                                 ::   zi, zj   ! local real
-      REAL(wp)                                 ::   zret     ! local real
-      REAL(wp), DIMENSION(4)                   ::   z4
-      !!----------------------------------------------------------------------
-      !
-      ! note that here there is no +0.5 added
-      ! since we're looking for four T points containing quadrant we're in of 
-      ! current T cell
-      ii = MAX(1, INT( pi     ))
-      ij = MAX(1, INT( pj     ))    ! T-point
-      zi = pi - REAL(ii,wp)
-      zj = pj - REAL(ij,wp)
-      !
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi )   ii = ii-1      
-      IF( ij == jpj )   ij = ij-1
-      !
-      z4(1) = pfld(ii  ,ij  )
-      z4(2) = pfld(ii+1,ij  )
-      z4(3) = pfld(ii  ,ij+1)
-      z4(4) = pfld(ii+1,ij+1)
-      IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN
-         WHERE( z4 < 0._wp ) z4 = z4 + 360._wp
-      ENDIF
-      !
-      zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
-      IF( zret > 180._wp ) zret = zret - 360._wp
-      icb_utl_bilin_x = zret
-      !
-   END FUNCTION icb_utl_bilin_x
-
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) ::   pfld      ! field to be interpolated
+      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
+      !!----------------------------------------------------------------------
+      !
+      ! 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
+      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
 
    REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
@@ -390,10 +431,9 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pet, peu, pev, pef   ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
-      REAL(wp)                , INTENT(in) ::   pi, pj               ! targeted coordinates in (i,j) referential
-      !
-      INTEGER  ::   ii, ij, icase, ierr   ! local integer
+      REAL(wp)                , INTENT(IN) ::   pi , pj              ! iceberg position
       !
       ! weights corresponding to corner points of a T cell quadrant
       REAL(wp) ::   zi, zj          ! local real
+      INTEGER  ::   ii, ij          ! bottom left corner coordinate in local domain
       !
       ! values at corner points of a T cell quadrant
@@ -402,6 +442,7 @@
       !!----------------------------------------------------------------------
       !
+      ! cannot used iiT because need ii/ij reltaive to global indices not local one
       ii = MAX(1, INT( pi ))   ;   ij = MAX(1, INT( pj ))            ! left bottom T-point (i,j) indices
-
+      ! 
       ! fractional box spacing
       ! 0   <= zi < 0.5, 0   <= zj < 0.5   -->  NW quadrant of current T cell
@@ -413,20 +454,7 @@
       zj = pj - REAL(ij,wp)
 
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      !
-      ierr = 0
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1       ; ierr = ierr + 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi     ; ierr = ierr + 1
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1       ; ierr = ierr + 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj     ; ierr = ierr + 1
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi ) THEN ; ii = ii-1 ; ierr = ierr + 1 ; END IF     
-      IF( ij == jpj ) THEN ; ij = ij-1 ; ierr = ierr + 1 ; END IF
-      !
-      IF ( ierr > 0 ) CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)')
+      ! conversion to local domain (no need to do a sanity check already done in icbpos)
+      ii = mi1(ii)
+      ij = mj1(ij)
       !
       IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
@@ -465,4 +493,51 @@
    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 = pe3(1) ; kb = 2
+      DO WHILE ( zdepw <  pD)
+         zdepw = zdepw + pe3(kb)
+         kb = kb + 1
+      END DO
+      kb = MIN(kb - 1,jpk)
+   END SUBROUTINE
+
+   SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_getkb         ***
+      !!
+      !! ** Purpose :   compute the vertical average of ocean properties affected by icb
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                INTENT(in ) :: kb        ! deepest level affected by icb
+      REAL(wp), DIMENSION(:), INTENT(in ) :: pe3, pdat ! vertical profile
+      REAL(wp),               INTENT(in ) :: pD        ! draft
+      REAL(wp),               INTENT(out) :: pzavg     ! z average
+      !!----------------------------------------------------------------------
+      INTEGER  :: jk
+      REAL(wp) :: zdep
+      !!----------------------------------------------------------------------
+      pzavg = 0.0 ; zdep = 0.0
+      DO jk = 1,kb-1
+         pzavg = pzavg + pe3(jk)*pdat(jk)
+         zdep  = zdep  + pe3(jk)
+      END DO
+      ! if kb is limited by mbkt  => bottom value is used between bathy and icb tail
+      ! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v)
+      pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD
+   END SUBROUTINE
 
    SUBROUTINE icb_utl_add( bergvals, ptvals )
@@ -653,5 +728,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))
@@ -679,5 +754,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) )
@@ -823,4 +898,75 @@
    END FUNCTION icb_utl_heat
 
+   SUBROUTINE test_icb_utl_getkb
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION test_icb_utl_getkb  ***
+      !!
+      !! ** Purpose : Test routine icb_utl_getkb, icb_utl_zavg
+      !! ** Methode : Call each subroutine with specific input data
+      !!              What should be the output is easy to determined and check 
+      !!              if NEMO return the correct answer.
+      !! ** Comments : not called, if needed a CALL test_icb_utl_getkb need to be added in icb_step
+      !!----------------------------------------------------------------------
+      INTEGER :: ikb
+      REAL(wp) :: zD, zout
+      REAL(wp), DIMENSION(jpk) :: ze3, zin
+      WRITE(numout,*) 'Test icb_utl_getkb : '
+      zD = 0.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 8000000.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 80.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 85.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 75.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      WRITE(numout,*) '=================================='
+      WRITE(numout,*) 'Test icb_utl_zavg'
+      zD = 0.0 ; ze3= 20.0 ; zin=1.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      zD = 50.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+      CALL FLUSH(numout)
+
+      zD = 80.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      zD = 80 ; ze3= 20.0 ; zin=1.0 ; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      ikb = 2
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      CALL FLUSH(numout)
+
+   END SUBROUTINE test_icb_utl_getkb
+
    !!======================================================================
 END MODULE icbutl
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcmod.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcmod.F90	(revision 14031)
@@ -475,5 +475,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. 
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/nemogcm.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/nemogcm.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/nemogcm.F90	(revision 14031)
@@ -34,4 +34,5 @@
    USE diu_layers     ! diurnal bulk SST and coolskin
    USE step_diu       ! diurnal bulk SST timestepping (called from here if run offline)
+   USE icb_oce        ! icebergs
    !
    USE prtctl         ! Print control
@@ -393,4 +394,10 @@
       ! ==>
                            CALL icb_init( rn_Dt, nit000)   ! initialise icebergs instance
+
+      ! compatibility check
+      IF( ln_icebergs .AND. ln_M2016 ) THEN
+         IF( lwp ) WRITE(numout,*) '         ==>>>   ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)'
+         CALL ctl_stop('ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)')
+      END IF
       !
       IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/sbcssm.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/sbcssm.F90	(revision 14030)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/sbcssm.F90	(revision 14031)
@@ -21,4 +21,5 @@
    USE zpshde         ! z-coord. with partial steps: horizontal derivatives
    USE closea         ! for ln_closea
+   USE icb_oce        ! for icebergs
    !
    USE in_out_manager ! I/O manager
@@ -226,5 +227,8 @@
          ln_closea = .false.
       ENDIF
-
+      IF( ln_icebergs .AND. ln_M2016 ) THEN
+         IF( lwp ) WRITE(numout,*) '         ==>>>   ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)'
+         CALL ctl_stop('ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)')
+      END IF
       !
       IF( l_sasread ) THEN                       ! store namelist information in an array
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/axis_def_nemo.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/axis_def_nemo.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/axis_def_nemo.xml	(revision 14031)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/axis_def_nemo.xml
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/context_nemo.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/context_nemo.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/context_nemo.xml	(revision 14031)
@@ -0,0 +1,37 @@
+<!--
+ ============================================================================================== 
+    NEMO context
+============================================================================================== 
+-->
+<context id="nemo">
+    <!-- $id$ -->
+    <variable_definition>
+       <!-- Year/Month/Day of time origin for NetCDF files; defaults to 1800-01-01 -->
+       <variable id="ref_year"  type="int"> 1900 </variable>
+       <variable id="ref_month" type="int"> 01 </variable>
+       <variable id="ref_day"   type="int"> 01 </variable>
+       <variable id="rau0"      type="float" > 1026.0 </variable>
+       <variable id="cpocean"   type="float" > 3991.86795711963 </variable>
+       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable>
+       <variable id="rhoic"     type="float" > 917.0 </variable>
+       <variable id="rhosn"     type="float" > 330.0 </variable>
+       <variable id="missval"  type="float" > 1.e20 </variable>	       
+    </variable_definition>
+
+<!-- Fields definition -->
+    <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     -->
+
+<!-- Files definition -->
+    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     -->
+
+<!-- Axis definition -->
+    <axis_definition src="./axis_def_nemo.xml"/>
+ 
+<!-- Domain definition -->
+    <domain_definition src="./domain_def_nemo.xml"/>
+
+<!-- Grids definition -->
+    <grid_definition   src="./grid_def_nemo.xml"/>
+  
+
+</context>
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/domain_def_nemo.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/domain_def_nemo.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/domain_def_nemo.xml	(revision 14031)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/domain_def_nemo.xml
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/field_def_nemo-oce.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/field_def_nemo-oce.xml	(revision 14031)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-oce.xml
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/file_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/file_def_nemo-oce.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/file_def_nemo-oce.xml	(revision 14031)
@@ -0,0 +1,58 @@
+<?xml version="1.0"?>
+    <!-- 
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own files                                         =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+    
+    <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1d" min_digits="4">
+    
+      <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."> <!-- 1 time step files -->
+	<file id="file11" output_freq="1ts" name_suffix="_grid_T" description="ocean T grid variables" >
+	  <field field_ref="empmr"        name="soemp"     />
+	  <field field_ref="berg_buoy_melt"  name="sobergbuoy"  />
+	  <field field_ref="berg_eros_melt"  name="sobergeros"  />
+	  <field field_ref="berg_conv_melt"  name="sobergconv"  />
+        </file>
+      </file_group>
+      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."> <!-- 1h files -->
+	<file id="file1" output_freq="1h" name_suffix="_grid_T" description="ocean T grid variables" >
+	  <field field_ref="toce"         name="votemper"  />
+	  <field field_ref="soce"         name="vosaline"  />
+	  <field field_ref="ssh"          name="sossheig"  />
+	  <field field_ref="empmr"        name="soemp"     />
+	  <field field_ref="berg_buoy_melt"  name="sobergbuoy"  />
+	  <field field_ref="berg_eros_melt"  name="sobergeros"  />
+	  <field field_ref="berg_conv_melt"  name="sobergconv"  />
+        </file>
+	<file id="file2" output_freq="1h" name_suffix="_grid_U" description="ocean U grid variables" >
+          <field field_ref="uoce"         name="vozocrtx"  />
+        </file>
+	<file id="file3" output_freq="1h" name_suffix="_grid_V" description="ocean V grid variables" >
+          <field field_ref="voce"         name="vomecrty"  /> 
+        </file>
+      </file_group>
+
+      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->
+     
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
+      <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->  
+
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
+      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m" output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m" output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+      <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."/> <!-- real yearly files -->
+      <file_group id="2y"  output_freq="2y" output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/grid_def_nemo.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/grid_def_nemo.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/grid_def_nemo.xml	(revision 14031)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/grid_def_nemo.xml
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/iodef.xml
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/iodef.xml	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/iodef.xml	(revision 14031)
@@ -0,0 +1,26 @@
+<?xml version="1.0"?>
+<simulation>
+
+<!-- ============================================================================================ -->
+<!-- XIOS context                                                                                 -->
+<!-- ============================================================================================ -->
+
+  <context id="xios" >
+
+      <variable_definition>
+
+          <variable id="info_level"                type="int">10</variable>
+          <variable id="using_server"              type="bool">false</variable>
+          <variable id="using_oasis"               type="bool">false</variable>
+          <variable id="oasis_codes_id"            type="string" >oceanx</variable>
+
+      </variable_definition>
+  </context>
+
+<!-- ============================================================================================ -->
+<!-- NEMO  CONTEXT add and suppress the components you need                                       -->
+<!-- ============================================================================================ -->
+
+  <context id="nemo" src="./context_nemo.xml"/>       <!--  NEMO       -->
+
+</simulation>
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/namelist_cfg	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/namelist_cfg	(revision 14031)
@@ -0,0 +1,484 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  Configuration namelist : overwrite default values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!                         ISOMIP configuration                       !!
+!!======================================================================
+!!              ***  Domain & Run management namelists  ***           !!
+!!                                                                    !!
+!!   namrun       parameters of the run
+!!   namdom       space and time domain
+!!   namcfg       parameters of the configuration                       (default: user defined GYRE)
+!!   namwad       Wetting and drying                                    (default: OFF)
+!!   namtsd       data: temperature & salinity                          (default: OFF)
+!!   namcrs       coarsened grid (for outputs and/or TOP)               (ln_crs =T)
+!!   namc1d       1D configuration options                              ("key_c1d")
+!!   namc1d_dyndmp 1D newtonian damping applied on currents             ("key_c1d")
+!!   namc1d_uvd   1D data (currents)                                    ("key_c1d")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   cn_exp      =  "ICB" !  experience name
+   nn_it000    =       1   !  first time step
+   nn_itend    =     480   !  last  time step
+   nn_leapy    =       0   !  Leap year calendar (1) or not (0)
+   ln_clobber  = .true.    !  clobber (overwrite) an existing file
+   nn_istate   =       0   !  output the initial state (1) or not (0)
+   nn_stock    =  99999999 !  frequency of creation of a restart file (modulo referenced to 1)
+   nn_write    =      48   !  frequency of write in the output file   (modulo referenced to nn_it000)
+   nn_istate   =       0   !  output the initial state (1) or not (0)
+/
+!-----------------------------------------------------------------------
+&namusr_def    !   ISOMIP user defined namelist  
+!-----------------------------------------------------------------------
+   ln_zps      = .false.   ! z-partial-step coordinate
+   ln_zco      = .true.    ! z-full-step coordinate
+   ln_sco      = .false.   ! s-coordinate
+   rn_dx       =   1000    !  zonal      grid-spacing (degrees)
+   rn_dy       =   1000    !  meridional grid-spacing (degrees)
+   rn_e3       =   10.     !  vertical resolution
+/
+!-----------------------------------------------------------------------
+&namdom        !   time and space domain
+!-----------------------------------------------------------------------
+   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time
+   rn_Dt       =  180.     !  time step for the dynamics (and tracer if nn_acc=0)
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtsd        !    Temperature & Salinity Data  (init/dmp)             (default: OFF)
+!-----------------------------------------------------------------------
+   
+/
+!-----------------------------------------------------------------------
+&namwad        !   Wetting and Drying (WaD)                             (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namcrs        !   coarsened grid (for outputs and/or TOP)              (ln_crs =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namc1d        !   1D configuration options                             ("key_c1d" default: PAPA station)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namc1d_dyndmp !   U & V newtonian damping                              ("key_c1d" default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namc1d_uvd    !   data: U & V currents                                 ("key_c1d" default: OFF)
+!-----------------------------------------------------------------------
+
+/
+
+!!======================================================================
+!!            ***  Surface Boundary Condition namelists  ***          !!
+!!                                                                    !!
+!!   namsbc          surface boundary condition manager                 (default: NO selection)
+!!   namsbc_flx      flux               formulation                     (ln_flx     =T)
+!!   namsbc_blk      Bulk formulae formulation                          (ln_blk     =T)
+!!   namsbc_cpl      CouPLed            formulation                     ("key_oasis3" )
+!!   namsbc_sas      Stand-Alone Surface module                         (SAS_SRC  only)
+!!   namsbc_iif      Ice-IF: use observed ice cover                     (nn_ice = 1   )
+!!   namtra_qsr      penetrative solar radiation                        (ln_traqsr  =T)
+!!   namsbc_ssr      sea surface restoring term (for T and/or S)        (ln_ssr     =T)
+!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T)
+!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T)
+!!   namsbc_wave     external fields from wave model                    (ln_wave    =T)
+!!   namberg         iceberg floats                                     (ln_icebergs=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
+!-----------------------------------------------------------------------
+   nn_fsbc     = 1         !  frequency of SBC module call
+      !                    !  (control sea-ice & iceberg model call)
+   ln_usr      = .true.    !  user defined formulation                  (T => check usrdef_sbc)
+/
+!-----------------------------------------------------------------------
+&namsbc_flx    !   surface boundary condition : flux formulation        (ln_flx =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_blk    !   namsbc_blk  generic Bulk formula                     (ln_blk =T)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&namsbc_sas    !   Stand-Alone Surface module: ocean data               (SAS_SRC  only)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&namsbc_iif    !   Ice-IF : use observed ice cover                      (nn_ice = 1)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
+!-----------------------------------------------------------------------
+   
+/
+!-----------------------------------------------------------------------
+&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&namsbc_rnf    !   runoffs                                              (ln_rnf =T)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&namsbc_apr    !   Atmospheric pressure used as ocean forcing           (ln_apr_dyn =T)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&namisf       !  Top boundary layer (ISF)                               (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model                        (ln_wave=T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namberg       !   iceberg parameters                                   (default: OFF)
+!-----------------------------------------------------------------------
+   ln_icebergs = .true.       ! activate iceberg floats (force =F with "key_agrif")
+   !                          ! diagnostics:
+   ln_bergdia        = .true.        ! Calculate budgets
+   nn_verbose_level  = 0             ! Turn on more verbose output if level > 0
+   nn_verbose_write  = 15            ! Timesteps between verbose messages
+   nn_sample_rate    = 1             ! Timesteps between sampling for trajectory storage
+   !
+   !                          ! iceberg setting:
+   !                                 ! Initial mass required for an iceberg of each class
+   rn_initial_mass   = 8.8e7, 4.1e8, 3.3e9, 1.8e10, 3.8e10, 7.5e10, 1.2e11, 2.2e11, 3.9e11, 7.4e11
+   !                                 ! Proportion of calving mass to apportion to each class
+   rn_distribution   = 0.24, 0.12, 0.15, 0.18, 0.12, 0.07, 0.03, 0.03, 0.03, 0.02
+   !                                 ! Ratio between effective and real iceberg mass (non-dim)
+   !                                 ! i.e. number of icebergs represented at a point
+   rn_mass_scaling   = 2000., 200., 50., 20., 10., 5., 2., 1., 1., 1.
+                                     ! thickness of newly calved bergs (m)
+   rn_initial_thickness     = 40., 67., 133., 175., 250., 250., 250., 250., 250., 100.
+   !
+   rn_rho_bergs            = 1025.    ! Density of icebergs
+   rn_LoW_ratio            = 1.5     ! Initial ratio L/W for newly calved icebergs
+   ln_operator_splitting   = .true.  ! Use first order operator splitting for thermodynamics
+   rn_bits_erosion_fraction = 0.     ! Fraction of erosion melt flux to divert to bergy bits
+   rn_sicn_shift           = 0.      ! Shift of sea-ice concn in erosion flux (0<sicn_shift<1)
+   ln_passive_mode         = .false. ! iceberg - ocean decoupling
+   nn_test_icebergs        =  10     ! Create test icebergs of this class (-1 = no)
+   !                                 ! Put a test iceberg at each gridpoint in box (lon1,lon2,lat1,lat2)
+   rn_test_box             = 1500.0,  39500.0, 9500.0, 10500.0
+   !
+   ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data)
+      ln_icb_grd           = .false.
+/
+
+!!======================================================================
+!!               ***  Lateral boundary condition  ***                 !!
+!!                                                                    !!
+!!   namlbc        lateral momentum boundary condition                  (default: NO selection)
+!!   namagrif      agrif nested grid   (read by child model only)       ("key_agrif")
+!!   nam_tide      Tidal forcing                                        (default: OFF)
+!!   nambdy        Unstructured open boundaries                         (default: OFF)
+!!   nambdy_dta    Unstructured open boundaries - external data         (see  nambdy)
+!!   nambdy_tide   tidal forcing at open boundaries                     (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+   rn_shlat    =    0.     !  free slip
+/
+!-----------------------------------------------------------------------
+&namagrif      !  AGRIF zoom                                            ("key_agrif")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_tide      !   tide parameters                                      (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nambdy        !  unstructured open boundaries                          (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nambdy_dta    !  open boundaries - external data                       (see nam_bdy)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&nambdy_tide   !  tidal forcing at open boundaries                      (default: OFF)
+!-----------------------------------------------------------------------
+/
+
+!!======================================================================
+!!                ***  Top/Bottom boundary condition  ***             !!
+!!                                                                    !!
+!!   namdrg        top/bottom drag coefficient                          (default: NO selection)
+!!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T)
+!!   namdrg_bot    bottom friction                                      (ln_OFF=F)
+!!   nambbc        bottom temperature boundary condition                (default: OFF)
+!!   nambbl        bottom boundary layer scheme                         (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namdrg        !   top/bottom drag coefficient                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_drg_OFF     = .true.    !  free-slip       : Cd = 0                  (F => fill namdrg_bot
+/
+!-----------------------------------------------------------------------
+&namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nambbc        !   bottom temperature boundary condition                (default: OFF)
+!-----------------------------------------------------------------------
+
+/
+!-----------------------------------------------------------------------
+&nambbl        !   bottom boundary layer scheme                         (default: OFF)
+!-----------------------------------------------------------------------
+/
+
+!!======================================================================
+!!                        Tracer (T-S) namelists                      !!
+!!                                                                    !!
+!!   nameos        equation of state                                    (default: NO selection)
+!!   namtra_adv    advection scheme                                     (default: NO selection)
+!!   namtra_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namtra_mle    mixed layer eddy param. (Fox-Kemper param.)          (default: OFF)
+!!   namtra_eiv    eddy induced velocity param.                         (default: OFF)
+!!   namtra_dmp    T & S newtonian damping                              (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nameos        !   ocean Equation Of Seawater                           (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_eos80    = .true.          !  = Use EOS80 equation of state
+/
+!-----------------------------------------------------------------------
+&namtra_adv    !   advection scheme for tracer                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_traadv_OFF = .true. !  No tracer advection
+/
+!-----------------------------------------------------------------------
+&namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_traldf_OFF   = .true.    !  No explicit diffusion
+/
+!-----------------------------------------------------------------------
+&namtra_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtra_eiv    !   eddy induced velocity param.                         (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtra_dmp    !   tracer: T & S newtonian damping                      (default: OFF)
+!-----------------------------------------------------------------------
+/
+
+!!======================================================================
+!!                      ***  Dynamics namelists  ***                  !!
+!!                                                                    !!
+!!   nam_vvl       vertical coordinate options                          (default: z-star)
+!!   namdyn_adv    formulation of the momentum advection                (default: NO selection)
+!!   namdyn_vor    advection scheme                                     (default: NO selection)
+!!   namdyn_hpg    hydrostatic pressure gradient                        (default: NO selection)
+!!   namdyn_spg    surface pressure gradient                            (default: NO selection)
+!!   namdyn_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namdta_dyn    offline TOP: dynamics read in files                  (OFF_SRC only)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nam_vvl       !   vertical coordinate options                          (default: z-star)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdyn_adv    !   formulation of the momentum advection                (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynadv_OFF = .true.  !  linear dynamics (no momentum advection)
+/
+!-----------------------------------------------------------------------
+&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynvor_ene = .true.  !  energy conserving scheme
+/
+!-----------------------------------------------------------------------
+&namdyn_hpg    !   Hydrostatic pressure gradient option                 (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_hpg_sco  = .true.    !  s-coordinate adapted for isf (standard jacobian formulation)
+/
+!-----------------------------------------------------------------------
+&namdyn_spg    !   surface pressure gradient                            (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynspg_ts  = .true.  !  split-explicit free surface
+/
+!-----------------------------------------------------------------------
+&namdyn_ldf    !   lateral diffusion on momentum                        (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynldf_OFF = .true.     !  No operator (i.e. no explicit diffusion)
+/
+!-----------------------------------------------------------------------
+&namdta_dyn    !   offline ocean input files                            (OFF_SRC only)
+!-----------------------------------------------------------------------
+
+/
+
+!!======================================================================
+!!                     vertical physics namelists                     !!
+!!                                                                    !!
+!!    namzdf        vertical physics manager                            (default: NO selection)
+!!    namzdf_ric    richardson number vertical mixing                   (ln_zdfric=T)
+!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
+!!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
+!!    namzdf_osm    OSM vertical diffusion                              (ln_zdfosm=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namzdf        !   vertical physics manager                             (default: NO selection)
+!-----------------------------------------------------------------------
+   !                       ! type of vertical closure (required)
+   ln_zdfcst   = .true.       !  constant mixing
+   !
+   !                       ! convection
+   ln_zdfevd   = .true.       !  enhanced vertical diffusion
+   rn_evd      =  0.1         ! mixing coefficient [m2/s]
+   !
+   !                       ! coefficients
+   rn_avm0     =   1.e-3      !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst)
+   rn_avt0     =   5.e-5      !  vertical eddy diffusivity [m2/s]       (background Kz if ln_zdfcst)
+/
+!-----------------------------------------------------------------------
+&namzdf_ric    !   richardson number dependent vertical diffusion       (ln_zdfric =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_gls    !   GLS vertical diffusion                               (ln_zdfgls =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
+!-----------------------------------------------------------------------
+/
+
+!!======================================================================
+!!                  ***  Diagnostics namelists  ***                   !!
+!!                                                                    !!
+!!   namtrd       dynamics and/or tracer trends                         (default: OFF)
+!!   namptr       Poleward Transport Diagnostics                        (default: OFF)
+!!   namhsb       Heat and salt budgets                                 (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namflo       float parameters                                      (default: OFF)
+!!   nam_diadct   transports through some sections                      (default: OFF)
+!!   nam_dia25h   25h Mean Output                                       (default: OFF)
+!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namtrd        !   trend diagnostics                                    (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namptr        !   Poleward Transport Diagnostic                        (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namhsb        !  Heat and salt budgets                                 (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdiu        !   Cool skin and warm layer models                      (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namflo        !   float parameters                                     ("key_float")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_diaharm   !   Harmonic analysis of tidal constituents              ("key_diaharm")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdct        ! transports through some sections                       ("key_diadct")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_diatmb    !  Top Middle Bottom Output                              (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_dia25h    !  25h Mean Output                                       (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4")
+!-----------------------------------------------------------------------
+/
+
+!!======================================================================
+!!               ***  Observation & Assimilation  ***                 !!
+!!                                                                    !!
+!!   namobs       observation and model comparison                      (default: OFF)
+!!   nam_asminc   assimilation increments                               ('key_asminc')
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namobs        !  observation usage switch                              (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_asminc    !   assimilation increments                              ('key_asminc')
+!-----------------------------------------------------------------------
+/
+
+!!======================================================================
+!!                  ***  Miscellaneous namelists  ***                 !!
+!!                                                                    !!
+!!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi")
+!!   namctl            Control prints                                   (default: OFF)
+!!   namsto            Stochastic parametrization of EOS                (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints                                       (default: OFF)
+!-----------------------------------------------------------------------
+   sn_cfctl%l_oceout  = .TRUE.
+   ln_timing   = .true.           !  timing by routine write out in timing.output file
+/
+!-----------------------------------------------------------------------
+&namsto        ! Stochastic parametrization of EOS                      (default: OFF)
+!-----------------------------------------------------------------------
+/
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/namelist_ref
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/namelist_ref	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/namelist_ref	(revision 14031)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ref
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/run_job.sh
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/run_job.sh	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/EXPREF/run_job.sh	(revision 14031)
@@ -0,0 +1,59 @@
+#!/bin/bash 
+#OAR -n nemo_sette
+#OAR -l core=2,walltime=00:05:00
+#OAR -O nemo_sette.o%jobid%
+#OAR -E nemo_sette.e%jobid%
+#OAR --project tipaccs
+
+. ~/bin/load_intelmodule.bash
+
+  export PBS_O_WORKDIR=$(readlink -f $PBS_O_WORKDIR)
+  export OMP_NUM_THREADS=1
+  cd $PBS_O_WORKDIR
+  export XIO_HOME=$xios_path
+#
+  echo " ";
+  export OMP_NUM_THREADS=1
+  O_PER_NODE=2
+  X_PER_NODE=2
+  OCORES=2
+  XCORES=0
+  if [ $OCORES -le 32 ] ; then O_PER_NODE=$OCORES; fi
+  export SETTE_DIR=/bettik/mathiotp/NEMO/NEMO_dev/branches/2020/tickets_icb_1900/sette
+
+  export EXE_DIR=/bettik/mathiotp/NEMO/NEMO_dev/branches/2020/tickets_icb_1900/tests/ICB/EXP00/
+  ulimit -c unlimited
+  ulimit -s unlimited
+#
+# end of set up
+###############################################################
+#
+# change to the working directory 
+#
+  cd $EXE_DIR
+  echo Directory is `pwd`
+  
+  if [ $XCORES -gt 0 ]; then
+#
+#  Run MPMD case
+#
+     #XIOS will run on a separate node so will run in parallel queue
+     if [ ! -f ./xios_server.exe ] && [ -f ${XIO_HOME}/bin/xios_server.exe ]; then
+        cp ${XIO_HOME}/bin/xios_server.exe .
+     fi
+     if [ ! -f ./xios_server.exe ]; then
+        echo "./xios_server.exe not found"
+        echo "run aborted"
+        exit
+     fi
+       echo time aprun -b -n $XCORES -N $X_PER_NODE ./xios_server.exe : -n $OCORES -N $O_PER_NODE ./nemo
+            time aprun -b -n $XCORES -N $X_PER_NODE ./xios_server.exe : -n $OCORES -N $O_PER_NODE ./nemo
+#
+  else
+#
+# Run SPMD case
+#
+    echo time mpirun -n $OCORES ./nemo
+         time mpirun -n $OCORES ./nemo
+  fi
+#
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/step.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/step.F90	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/step.F90	(revision 14031)
@@ -0,0 +1,175 @@
+MODULE step
+   !!======================================================================
+   !!                       ***  MODULE step  ***
+   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
+   !!======================================================================
+   !! History :  OPA  !  1991-03  (G. Madec)  Original code
+   !!             -   !  1991-11  (G. Madec)
+   !!             -   !  1992-06  (M. Imbard)  add a first output record
+   !!             -   !  1996-04  (G. Madec)  introduction of dynspg
+   !!             -   !  1996-04  (M.A. Foujols)  introduction of passive tracer
+   !!            8.0  !  1997-06  (G. Madec)  new architecture of call
+   !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
+   !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit
+   !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
+   !!             -   !  2004-08  (C. Talandier) New trends organization
+   !!             -   !  2005-01  (C. Ethe) Add the KPP closure scheme
+   !!             -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !!             -   !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!             -   !  2006-07  (S. Masson)  restart using iom
+   !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
+   !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl
+   !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
+   !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
+   !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
+   !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs
+   !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state
+   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
+   !!            3.7  !  2014-10  (G. Madec)  LDF simplication 
+   !!             -   !  2014-12  (G. Madec) remove KPP scheme
+   !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface)
+   !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy)
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
+   !!----------------------------------------------------------------------
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+#else
+   !!----------------------------------------------------------------------
+   !!   stp             : OPA system time-stepping
+   !!----------------------------------------------------------------------
+   USE step_oce         ! time stepping definition modules
+   !
+   USE iom              ! xIOs server
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp   ! called by nemogcm.F90
+
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 13237 2020-07-03 09:12:53Z smasson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_agrif
+   RECURSIVE SUBROUTINE stp( )
+      INTEGER             ::   kstp   ! ocean time-step index
+#else
+   SUBROUTINE stp( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+#endif
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp  ***
+      !!
+      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
+      !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.)
+      !!              - Time stepping of TRC  (passive tracer eqs.)
+      !!
+      !! ** Method  : -1- Update forcings and data
+      !!              -2- Update 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, hdiv,w)
+      !!              -8- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+!!gm kcall can be removed, I guess
+      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
+      !! ---------------------------------------------------------------------
+#if defined key_agrif
+      IF( nstop > 0 ) RETURN   ! avoid to go further if an error was detected during previous time step (child grid)
+      kstp = nit000 + Agrif_Nb_Step()
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
+# if defined key_iomput
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
+# endif
+#endif
+      !
+      IF( ln_timing )   CALL timing_start('stp')
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! model timestep
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF( l_1st_euler ) THEN  
+         ! start or restart with Euler 1st time-step
+         rDt =  rn_Dt   
+         r1_Dt = 1._wp / rDt
+      ENDIF
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom)
+         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
+                             CALL iom_init_closedef
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice)
+
+                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
+
+!!gm : This does not only concern the dynamics ==>>> add a new title
+!!gm2: why ouput restart before AGRIF update?
+!!
+!!jc: That would be better, but see comment above
+!!
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! File manipulation at the end of the first time step
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<                         
+      IF( kstp == nit000 ) THEN                          ! 1st time step only
+                                        CALL iom_close( numror )   ! close input  ocean restart file
+         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
+         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
+      ENDIF
+      !
+#if defined key_iomput
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Finalize contextes if end of simulation or error detected
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<                         
+      IF( kstp == nitend .OR. nstop > 0 ) THEN 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+         IF( lrxios ) CALL iom_context_finalize(      crxios_context         )
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
+      ENDIF
+#endif
+      !
+      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
+         rDt = 2._wp * rn_Dt   
+         r1_Dt = 1._wp / rDt
+         l_1st_euler = .FALSE.      
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('stp')
+      !
+   END SUBROUTINE stp
+   !
+#endif
+   !!======================================================================
+END MODULE step
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_hgr.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_hgr.F90	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_hgr.F90	(revision 14031)
@@ -0,0 +1,122 @@
+MODULE usrdef_hgr
+   !!======================================================================
+   !!                       ***  MODULE usrdef_hgr   ***
+   !!
+   !!                  ===  LOCK_EXCHANGE configuration  ===
+   !!
+   !! User defined :   mesh and Coriolis parameter of a user configuration
+   !!======================================================================
+   !! History :  NEMO  ! 2016-08  (S. Flavoni, G. Madec)   Original code
+   !!                  ! 2017-02  (P. Mathiot, S. Flavoni) Adapt code to ISOMIP case
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_hgr    : initialize the horizontal mesh for ISOMIP configuration
+   !!----------------------------------------------------------------------
+   USE dom_oce  ,  ONLY: nimpp, njmpp       ! ocean space and time domain
+   USE par_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE usrdef_nam, ONLY: rn_dx, rn_dy
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_hgr   ! called by domhgr.F90
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_hgr.F90 12740 2020-04-12 09:03:06Z smasson $ 
+   !! 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
+      !!                ISOMIP configuration
+      !!
+      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 
+      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
+      !!              - define u- & v-surfaces (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) ::   zfact, zlam0, zphi0, zti, zui, zvi, zfi, ztj, zuj, zvj, zfj      ! local scalars
+      !!-------------------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'usr_def_hgr : ISOMIP configuration'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*)
+         WRITE(numout,*) '   ===>>  geographical mesh on the sphere with regular grid-spacing'
+         WRITE(numout,*) '          given by rn_e1deg and rn_e2deg'
+      ENDIF
+      !
+      zlam0 = 0.0
+      zphi0 = 0.0
+      DO_2D (1,1,1,1)
+         zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 )
+         zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp
+         
+         plamt(ji,jj) = zlam0 + rn_dx * zti
+         plamu(ji,jj) = zlam0 + rn_dx * zui
+         plamv(ji,jj) = plamt(ji,jj) 
+         plamf(ji,jj) = plamu(ji,jj) 
+         
+         pphit(ji,jj) = zphi0 + rn_dy * ztj
+         pphiv(ji,jj) = zphi0 + rn_dy * zvj
+         pphiu(ji,jj) = pphit(ji,jj) 
+         pphif(ji,jj) = pphiv(ji,jj) 
+      END_2D
+      !
+      !                       !==  Horizontal scale factors  ==!   (in meters)
+      DO_2D (1,1,1,1)
+         !                       ! e1   (zonal)
+         pe1t(ji,jj) = rn_dx
+         pe1u(ji,jj) = rn_dx
+         pe1v(ji,jj) = rn_dx
+         pe1f(ji,jj) = rn_dx
+         !                       ! e2   (meridional)
+         pe2t(ji,jj) = rn_dy
+         pe2u(ji,jj) = rn_dy
+         pe2v(ji,jj) = rn_dy
+         pe2f(ji,jj) = rn_dy
+      END_2D
+      !                             ! 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 = 1                       ! Coriolis parameter set to 0
+      pff_f(:,:) = 0._wp            ! CAUTION: set to zero to avoid error with some compilers that
+      pff_t(:,:) = 0._wp            !             require an initialization of INTENT(out) arguments
+      !
+   END SUBROUTINE usr_def_hgr
+
+   !!======================================================================
+END MODULE usrdef_hgr
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_istate.F90	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_istate.F90	(revision 14031)
@@ -0,0 +1,73 @@
+MODULE usrdef_istate
+   !!======================================================================
+   !!                     ***  MODULE usrdef_istate   ***
+   !!
+   !!                  ===  ISOMIP configuration  ===
+   !!
+   !! User defined : set the initial state of a user configuration
+   !!======================================================================
+   !! History :  NEMO ! 2016-11 (S. Flavoni)             Original code
+   !!                 ! 2017-02 (P. Mathiot, S. Flavoni) Adapt code to ISOMIP case
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  usr_def_istate : initial state in Temperature and salinity
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean space and time domain
+   USE dom_oce
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_istate   ! called by istate.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_istate.F90 10074 2018-08-28 16:15:49Z nicolasmartin $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+  
+   SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv, pssh )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE usr_def_istate  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers
+      !!                Here ISOMIP configuration 
+      !!
+      !! ** Method  : - set temperature field
+      !!              - set salinity    field
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pdept   ! depth of t-point               [m]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   ptmask  ! t-point ocean mask             [m]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts     ! T & S fields      [Celsius ; g/kg]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pu      ! i-component of the velocity  [m/s] 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pv      ! j-component of the velocity  [m/s] 
+      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out) ::   pssh    ! sea-surface height
+      !
+      INTEGER  ::   jk     ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_istate : ISOMIP configuration, analytical definition of initial state'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   Ocean at rest, with a constant salinity and temperature.      '
+      pu  (:,:,:) = 1._wp * umask(:,:,:)       ! ocean at rest
+      pv  (:,:,:) = 0._wp * vmask(:,:,:) 
+      pu  (:,:,3:jpkm1) = 0._wp * umask(:,:,3:jpkm1)        ! ocean at rest
+      pv  (:,:,3:jpkm1) = 1._wp * vmask(:,:,3:jpkm1)
+      pssh(:,:)   = 0._wp
+      !
+      !                          ! T & S profiles
+      pts(:,:,:,jp_tem) = 2  * ptmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields
+      pts(:,:,:,jp_sal) = 34 * ptmask(:,:,:)
+
+      pts(:,:,3:jpkm1,jp_tem) = 1  * ptmask(:,:,3:jpkm1)          ! ISOMIP configuration : start from constant T+S fields
+      !   
+   END SUBROUTINE usr_def_istate
+
+   !!======================================================================
+END MODULE usrdef_istate
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_nam.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_nam.F90	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_nam.F90	(revision 14031)
@@ -0,0 +1,108 @@
+MODULE usrdef_nam
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_nam  ***
+   !!
+   !!                     ===  ISOMIP configuration  ===
+   !!
+   !! User defined : set the domain characteristics of a user configuration
+   !!======================================================================
+   !! History :  NEMO ! 2016-03  (S. Flavoni, G. Madec)   Original code
+   !!                 ! 2017-02  (P. Mathiot, S. Flavoni) Adapt code to ISOMIP case
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   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            ! i- & j-indices of the local 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
+   USE timing         ! Timing
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_nam   ! called by nemogcm.F90
+
+   !                                         !!* namusr_def namelist *!!
+   REAL(wp), PUBLIC ::   rn_dx, rn_dy         !: horizontal resolution   [degrees]
+   REAL(wp), PUBLIC ::   rn_e3                !: vertical   resolution         [m]
+   
+   REAL(wp), PARAMETER, PUBLIC ::   rbathy = 100._wp   !: depth of the seafloor   [m]
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_nam.F90 12377 2020-02-12 14:39:06Z acc $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_nam( 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 ISOMIP configuration
+      !!
+      !! ** input   : - namusr_def namelist found in namelist_cfg
+      !!----------------------------------------------------------------------
+      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   ! Local integer
+      !!
+      NAMELIST/namusr_def/ ln_zco, ln_zps, ln_sco, rn_dx, rn_dy, rn_e3
+      !!----------------------------------------------------------------------
+      !
+      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' )
+      !
+      IF(lwm)   WRITE( numond, namusr_def )
+      !
+      cd_cfg = 'ISOMIP'           ! name & resolution (not used)
+      kk_cfg = INT( rn_e3 )
+      !
+      ! Global Domain size:  ICB domain is  15° x 10° x 900 m
+      kpi = 40 + 2     ! add 2 for t-point in the east  & west  coasts
+      kpj = 40 + 2     !     -        -           north & south   -
+      kpk = INT( rbathy / rn_e3    ) + 1     ! add 1 for t-point in the seafloor
+      !
+      !                             ! Set the lateral boundary condition of the global domain
+      kperio = 0                    ! ICB configuration : box
+      !
+      !                             ! control print
+      IF(lwp) THEN
+         WRITE(numout,*) '   '
+         WRITE(numout,*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '   Namelist namusr_def : ISOMIP test case'
+         WRITE(numout,*) '      type of vertical coordinate : '
+         WRITE(numout,*) '         z-coordinate flag                     ln_zco   = ', ln_zco
+         WRITE(numout,*) '         z-partial-step coordinate flag        ln_zps   = ', ln_zps
+         WRITE(numout,*) '         s-coordinate flag                     ln_sco   = ', ln_sco
+         WRITE(numout,*) '      resolution'
+         WRITE(numout,*) '         zonal      resolution                 rn_dx    = ', rn_dx   , ' meters '
+         WRITE(numout,*) '         meridional resolution                 rn_dy    = ', rn_dy   , ' meters '
+         WRITE(numout,*) '         vertical   resolution                 rn_e3    = ', rn_e3   , ' meters'
+         WRITE(numout,*) '      ICB domain = 15° x 10° x 900 m'
+         WRITE(numout,*) '         resulting global domain size :        jpiglo   = ', kpi
+         WRITE(numout,*) '                                               jpjglo   = ', kpj
+         WRITE(numout,*) '                                               jpkglo   = ', kpk
+         WRITE(numout,*) '   '
+         WRITE(numout,*) '   Lateral boundary condition of the global domain'
+         WRITE(numout,*) '      ICB : closed basin                    jperio   = ', kperio
+      ENDIF
+      !
+   END SUBROUTINE usr_def_nam
+
+   !!======================================================================
+END MODULE usrdef_nam
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_sbc.F90	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_sbc.F90	(revision 14031)
@@ -0,0 +1,91 @@
+MODULE usrdef_sbc
+   !!======================================================================
+   !!                     ***  MODULE  usrdef_sbc  ***
+   !! 
+   !!                  ===  ISOMIP configuration  ===
+   !!
+   !! User defined :   surface forcing of a user configuration
+   !!======================================================================
+   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)   user defined interface
+   !!                  ! 2017-02  (P. Mathiot, S. Flavoni) adapt code to ISOMIP case
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_sbc    : user defined surface bounday conditions in ISOMIP case
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbc_ice         ! Surface boundary condition: ice fields
+   USE phycst          ! physical constants
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usrdef_sbc_oce      ! routine called in sbcmod module
+   PUBLIC   usrdef_sbc_ice_tau  ! routine called by icestp.F90 for ice dynamics
+   PUBLIC   usrdef_sbc_ice_flx  ! routine called by icestp.F90 for ice thermo
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_sbc.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usrdef_sbc_oce( kt, Kbb )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE usr_def_sbc  ***
+      !!              
+      !! ** Purpose :   provide at each time-step the surface boundary
+      !!              condition, i.e. the momentum, heat and freshwater fluxes.
+      !!
+      !! ** Method  :   all 0 fields, for ISOMIP case
+      !!                CAUTION : never mask the surface stress field !
+      !!
+      !! ** Action  : - set to ZERO all the ocean surface boundary condition, i.e.   
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb  ! ocean time index
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)' usr_sbc : ISOMIP case: NO surface forcing'
+         IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~   utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0'
+         !
+      ENDIF
+      !
+      utau(:,:) = 0._wp
+      vtau(:,:) = 0._wp
+      taum(:,:) = 0._wp
+      wndm(:,:) = 0._wp
+      !
+      emp (:,:) = 0._wp
+      sfx (:,:) = 0._wp
+      qns (:,:) = 0._wp
+      qsr (:,:) = 0._wp
+      !         
+      !
+   END SUBROUTINE usrdef_sbc_oce
+
+   SUBROUTINE usrdef_sbc_ice_tau( kt )
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+   END SUBROUTINE usrdef_sbc_ice_tau
+
+
+   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
+   END SUBROUTINE usrdef_sbc_ice_flx
+
+   !!======================================================================
+END MODULE usrdef_sbc
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_zgr.F90
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_zgr.F90	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/MY_SRC/usrdef_zgr.F90	(revision 14031)
@@ -0,0 +1,116 @@
+MODULE usrdef_zgr
+   !!======================================================================
+   !!                   ***  MODULE  usrdef_zgr  ***
+   !!
+   !!                   ===      ICE_AGRIF case     ===
+   !!
+   !! Ocean domain : user defined vertical coordinate system 
+   !!======================================================================
+   !! History :  4.0  ! 2016-08  (G. Madec, S. Flavoni)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_zgr   : user defined vertical coordinate system (required)
+   !!---------------------------------------------------------------------
+   USE oce            ! ocean variables
+   USE dom_oce        ! ocean domain
+   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
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_zgr   ! called by domzgr.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_zgr.F90 12597 2020-03-25 08:57:21Z smasson $
+   !! 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
+      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
+      !
+      INTEGER  ::   jk, k_dz  ! dummy indices
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_zgr : ICE_AGRIF configuration '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (slab ocean - advection of an ice patch in a biperiodic square box domain)'
+      !
+      !
+      ! type of vertical coordinate  ==>>>   here ICE_AGRIF : slab ocean always
+      ! ---------------------------
+      ld_zco    = .TRUE.       ! z-full-step coordinate
+      ld_zps    = .FALSE.      ! z-partial-step coordinate
+      ld_sco    = .FALSE.      ! s-coordinate
+      ld_isfcav = .FALSE.      ! ISF Ice Shelves Flag
+      !
+      !
+      ! Build the vertical coordinate system
+      ! ------------------------------------
+      !
+      !                       !==  UNmasked meter bathymetry  ==!
+      !
+      ! 
+      k_dz = 10
+      pdepw_1d(1) = 0.0
+      pdept_1d(1) = k_dz / 2.0
+      pe3w_1d = k_dz
+      pe3t_1d = k_dz
+      !
+      DO jk = 2, jpk
+         pdepw_1d(jk) = pdepw_1d(jk-1) +   k_dz
+         pdept_1d(jk) = pdept_1d(jk-1) +   k_dz
+      END DO
+      !                       !==  top masked level bathymetry  ==!  (all coordinates)
+      !
+      ! no ocean cavities : top ocean level is ONE, except over land
+      k_top(:,:) = 1
+      !
+      !                       !==  z-coordinate  ==!   (step-like topography)
+      !                                !* bottom ocean compute from the depth of grid-points
+      jpkm1 = jpk-1
+      k_bot(:,:) = 10         ! here use k_top as a land mask
+      WHERE (( glamt < 10000 .OR. glamt > 30000 ))
+         k_bot(:,:) = 4
+      END WHERE
+      !                                !* horizontally uniform coordinate (reference z-co everywhere)
+      DO jk = 1, jpk
+         pdept(:,:,jk) = pdept_1d(jk)
+         pdepw(:,:,jk) = pdepw_1d(jk)
+         pe3t (:,:,jk) = pe3t_1d (jk)
+         pe3u (:,:,jk) = pe3t_1d (jk)
+         pe3v (:,:,jk) = pe3t_1d (jk)
+         pe3f (:,:,jk) = pe3t_1d (jk)
+         pe3w (:,:,jk) = pe3w_1d (jk)
+         pe3uw(:,:,jk) = pe3w_1d (jk)
+         pe3vw(:,:,jk) = pe3w_1d (jk)
+      END DO
+      !
+   END SUBROUTINE usr_def_zgr
+
+   !!======================================================================
+END MODULE usrdef_zgr
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/README
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/README	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/README	(revision 14031)
@@ -0,0 +1,38 @@
+ICB 
+---
+
+ICB is a very idealized configuration used to test and debug the icb module.
+The configuration is box with a shallow shelf (40m) on the east and west part of the domain with a deep central trough (> 100m).
+ICB are generating using the test capability of the icb model along a E-W line (this can easily be tuned).
+
+In this configuration, ONLY the ICB module is activated. It means that:
+- SI3 is OFF
+- NO atmospheric fluxes. See usr routine in MY_SRC.
+- NO time variation of ocean properties. T, S, U, V, ssh (...) are intialized and never changed in time. This is achieved by simplifying the step.F90 routine in MY_SRC.
+
+Resolution and Initial condition are:
+- Homogeneous S (34)
+- Spatially constant T     set to 2     over the first 2 levels and 1     below
+- Spatially constant (U,V) set to (1,0) over the first 2 levels and (0,1) below
+- e3 is set to 10m
+- iceberg thicness is set to 100 in the class tested (10 level affected).
+- iceberg density is set to 1026 as the ocean reference density => icb tikness = icb draft (easier to evaluate the output)
+
+Results expected are:
+- ln_M2016 = F => icb moving along the W direction only and pilling up on the W side of the domain.
+- ln_M2016 = T => icb on shallow bank move at 45 degrees (NW) and 76 degrees (NNNNW) in the deeper through. 
+- ln_M2016 = T and ln_icb_grd = T => as before but now icb cannot move onto the shallow bank on the W and migrate northward along the shelf break. Some get stuck at specific point because of the grounding law (they always get back to there intial position at every time step).
+
+This set up allow us to easily test various option of the ICB code as:
+- ln_M2016   (use properties in depth to compute dyn and thd)
+- ln_icb_grd (use icb bottom level of the icb (tail) to ground the icb or not)
+
+It can also help us to:
+- reproduce bugs (interpolation issue for example)
+- develop and test new functionalities (new grounding scheme for exemple)
+- evaluate ICB performance alone
+
+The side effect of the implementation is that:
+- the run.stat (if produced) is meaningless as NEMO tracer and dynamics are not computed
+
+
Index: /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/cpp_ICB.fcm
===================================================================
--- /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/cpp_ICB.fcm	(revision 14031)
+++ /NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/tests/ICB/cpp_ICB.fcm	(revision 14031)
@@ -0,0 +1,1 @@
+ bld::tool::fppkeys   key_iomput key_mpp_mpi 
