Index: branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 8933)
@@ -45,4 +45,8 @@
 !-----------------------------------------------------------------------
 &namtra_qsr    !   penetrative solar radiation
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_cpl
 !-----------------------------------------------------------------------
 /
Index: branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8933)
@@ -187,4 +187,5 @@
    ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 )
    ln_mixcpl   = .false.   !  forced-coupled mixed formulation          ( requires key_oasis3 )
+   ln_meto_cpl = .false.   !  Met-Office coupling specifics
    nn_components = 0       !  configuration of the opa-sas OASIS coupling
                            !  =0 no opa-sas OASIS coupling: default single executable config.
@@ -272,32 +273,41 @@
    sn_snd_crt    =   'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'
    sn_snd_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
-   sn_snd_crtw   =       'none'                 ,    'no'    ,     ''      ,         ''           , 'U,V'
-   sn_snd_ifrac  =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
-   sn_snd_wlev   =       'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_crtw   =   'none'                 ,    'no'    ,     ''      ,         ''           , 'U,V'
+   sn_snd_ifrac  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_wlev   =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_cond   =   'weighted ice'         ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_thick1 =   'ice and snow'         ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_mpnd   =   'weighted ice'         ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_sstfrz =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_ttilyr =   'weighted ice'         ,    'no'    ,     ''      ,         ''           ,   ''
 ! receive
-   sn_rcv_w10m   =   'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_taumod =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_tau    =   'oce only'             ,    'no'    , 'cartesian' , 'eastward-northward',  'U,V'
-   sn_rcv_dqnsdt =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_qsr    =   'oce and ice'          ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_qns    =   'oce and ice'          ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_emp    =   'conservative'         ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_rnf    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_cal    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_hsig   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_iceflx =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_sdrfx  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_wstrf  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_w10m   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_taumod =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_tau    =   'oce only'             ,    'no'    , 'cartesian' , 'eastward-northward' ,  'U,V'
+   sn_rcv_dqnsdt =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_qsr    =   'oce and ice'          ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_qns    =   'oce and ice'          ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_emp    =   'conservative'         ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_rnf    =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_cal    =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_hsig   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_iceflx =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_mslp   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_phioc  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_sdrfx  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_sdrfy  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wper   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wnum   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wstrf  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wdrag  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_ts_ice =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_isf    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_icb    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
 !
    nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentialy sending/receiving data
    ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models
    !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8933)
@@ -313,4 +313,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d     !: maximum ice concentration 2d array
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice        !: transmitted solar radiation under ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice         !: temperature of the first layer (Jules coupling) [K]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice        !: effective conductivity at the top of ice/snow (Jules coupling) [W.m-2.K-1]
 
    !!----------------------------------------------------------------------
@@ -325,5 +327,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_su      !: Sea-Ice Surface Temperature (K)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   s_i       !: Sea-Ice Bulk salinity (ppt)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sv_i     !: Sea-Ice Bulk salinity times volume per area (ppt.m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sv_i      !: Sea-Ice Bulk salinity times volume per area (ppt.m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (s)
@@ -364,6 +366,6 @@
    !! * Old values of global variables
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b  !: snow and ice volumes/thickness
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b        !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b    !: snow and ice volumes/thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b         !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                         !: snow heat content
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                         !: ice temperatures
@@ -449,15 +451,15 @@
       ! * Ice global state variables
       ii = ii + 1
-      ALLOCATE( ftr_ice(jpi,jpj,jpl) ,                                                 &
-         &      h_i   (jpi,jpj,jpl) , a_i   (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,     &
-         &      v_s    (jpi,jpj,jpl) , h_s  (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,     &
-         &      s_i   (jpi,jpj,jpl) , sv_i (jpi,jpj,jpl) , o_i   (jpi,jpj,jpl) ,     &
-         &      oa_i   (jpi,jpj,jpl) , bv_i  (jpi,jpj,jpl) ,  STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                       &
-         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     &
-         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     &
-         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) ,     &
-         &      om_i (jpi,jpj) , tau_icebfr(jpi,jpj)                              , STAT=ierr(ii) )
+      ALLOCATE( ftr_ice(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , &
+         &      h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) ,     &
+         &      v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , t_su(jpi,jpj,jpl) ,     &
+         &      s_i (jpi,jpj,jpl) , sv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) ,     &
+         &      oa_i(jpi,jpj,jpl) , bv_i(jpi,jpj,jpl) ,  STAT=ierr(ii) )
+      ii = ii + 1
+      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                      &
+         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,     &
+         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , bvm_i(jpi,jpj) ,     &
+         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s (jpi,jpj) ,     &
+         &      om_i (jpi,jpj) , tau_icebfr(jpi,jpj)                             , STAT=ierr(ii) )
       ii = ii + 1
       ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
@@ -466,6 +468,5 @@
 
       ii = ii + 1
-      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , &
-         &      h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
+      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
       ii = ii + 1
       ALLOCATE( at_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
@@ -473,7 +474,7 @@
       ! * Old values of global variables
       ii = ii + 1
-      ALLOCATE( v_s_b  (jpi,jpj,jpl) , v_i_b  (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   &
-         &      a_i_b  (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b (jpi,jpj,nlay_i,jpl) , e_s_b (jpi,jpj,nlay_s,jpl) ,   &
-         &      oa_i_b (jpi,jpj,jpl)                                                     , STAT=ierr(ii) )
+      ALLOCATE( v_s_b  (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   &
+         &      a_i_b  (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   &
+         &      oa_i_b (jpi,jpj,jpl)                                                   , STAT=ierr(ii) )
       ii = ii + 1
       ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
@@ -491,7 +492,5 @@
       ! * SIMIP diagnostics
       ii = ii + 1
-      ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
-                diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
-                STAT = ierr(ii) )
+      ALLOCATE( t_si (jpi,jpj,jpl) , tm_si(jpi,jpj) , diag_fc_bo(jpi,jpj) , diag_fc_su(jpi,jpj) , STAT = ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90	(revision 8933)
@@ -38,6 +38,8 @@
    
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qml_ice_1d     !: heat available for snow / ice surface melting [W/m2] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qcn_ice_1d     !: heat available for snow / ice surface sublimation¬ [W/m2] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_tr_1d  !: solar flux transmitted below the ice surface¬ [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qcn_ice_1d     !: heat available for snow / ice surface sublimation [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_tr_1d  !: solar flux transmitted below the ice surface [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t1_ice_1d      !: temperature of the 1st layer (Jules coupling) [K]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   cnd_ice_1d     !: conductivity at the top of ice/snow (Jules coupling) [W/K/m2]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d
@@ -176,11 +178,11 @@
 
       ii = 1
-      ALLOCATE( nptidx   (jpij) ,   &
-         &      qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   &
-         &      qns_ice_1d(jpij)  ,   &
-         &      qml_ice_1d(jpij), qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , &
-         &      t_bo_1d   (jpij) ,                                         &
-         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    & 
-         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,hfx_dyn_1d(jpij) ,    &
+      ALLOCATE( nptidx    (jpij) ,   &
+         &      qlead_1d  (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   &
+         &      qns_ice_1d(jpij) ,   &
+         &      qml_ice_1d(jpij) , qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , &
+         &      cnd_ice_1d(jpij) , t1_ice_1d(jpij)  , t_bo_1d   (jpij) ,      &
+         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) , hfx_bog_1d(jpij) ,    & 
+         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) , hfx_dyn_1d(jpij) ,    &
          &      rn_amax_1d(jpij) ,                                         &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8933)
@@ -354,5 +354,7 @@
          !
          tn_ice (:,:,:) = t_su (:,:,:)
-         !
+         t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu
+         cnd_ice(:,:,:) = 0._wp           ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
+
          ! Melt pond volume and fraction
          IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
@@ -388,5 +390,9 @@
             END DO
          END DO
-         !
+
+         tn_ice (:,:,:) = t_i (:,:,1,:)
+         t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu
+         cnd_ice(:,:,:) = 0._wp           ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
+         
          a_ip(:,:,:)      = 0._wp
          v_ip(:,:,:)      = 0._wp
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8933)
@@ -19,5 +19,5 @@
    USE ice            ! sea-ice variables
    USE dom_oce        ! ocean domain
-   USE sbc_oce , ONLY : nn_fsbc
+   USE sbc_oce , ONLY : nn_fsbc, ln_meto_cpl
    USE icectl
    !
@@ -129,14 +129,11 @@
       CALL iom_rstput( iter, nitrst, numriw, 'a_i' , a_i  )
       CALL iom_rstput( iter, nitrst, numriw, 't_su', t_su )
-      !
       ! Melt ponds
       CALL iom_rstput( iter, nitrst, numriw, 'a_ip', a_ip )
       CALL iom_rstput( iter, nitrst, numriw, 'v_ip', v_ip )
-      !
 !!gm dangerous !!!!!  ===>>>> better reading writing all snow layers !
       ! Snow enthalpy (1st snow layer only)
       z3d = e_s(:,:,1,:)
       CALL iom_rstput( iter, nitrst, numriw, 'tempt_sl1' , z3d )
-      !
       ! Ice enthalpy (all ice layers)
       DO jk = 1, nlay_i 
@@ -146,9 +143,14 @@
          CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
       END DO
-      !
       ! ice velocity
       CALL iom_rstput( iter, nitrst, numriw, 'u_ice', u_ice ) ! u_ice
       CALL iom_rstput( iter, nitrst, numriw, 'v_ice', v_ice ) ! v_ice
-      
+      ! fields needed for Jules coupling
+      IF (ln_meto_cpl) THEN
+         CALL iom_rstput( iter, nitrst, numriw, 'cnd_ice', cnd_ice )
+         CALL iom_rstput( iter, nitrst, numriw, 't1_ice' , t1_ice  )
+      ENDIF
+      !
+
       ! close restart file
       ! ------------------
@@ -207,5 +209,4 @@
       CALL iom_get( numrir, jpdom_autoglo, 'a_i' , a_i  )
       CALL iom_get( numrir, jpdom_autoglo, 't_su', t_su )
-      !
       ! Melt ponds
       id1 = iom_varid( numrir, 'a_ip_htc01' , ldstop = .FALSE. )
@@ -218,10 +219,8 @@
          v_ip(:,:,:) = 0._wp
       ENDIF
-      !
 !!gm dangerous !!!!!  ===>>>> better reading writing all snow layers !
       ! Snow enthalpy (1st snow layer only)
       CALL iom_get( numrir, jpdom_autoglo, 'tempt_sl1' , z3d )
       e_s(:,:,1,:) = z3d
-      !
       ! Ice enthalpy (all ice layers)
       DO jk = 1, nlay_i 
@@ -231,8 +230,12 @@
          e_i(:,:,jk,:) = z3d(:,:,:)
       END DO
-      !
       ! ice velocity
       CALL iom_get( numrir, jpdom_autoglo, 'u_ice', u_ice )
       CALL iom_get( numrir, jpdom_autoglo, 'v_ice', v_ice )
+      ! fields needed for Jules coupling
+      IF (ln_meto_cpl) THEN
+         CALL iom_get( numrir, jpdom_autoglo, 'cnd_ice', cnd_ice )
+         CALL iom_get( numrir, jpdom_autoglo, 't1_ice' , t1_ice  )
+      ENDIF
 
    END SUBROUTINE ice_rst_read
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8933)
@@ -319,4 +319,8 @@
          IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
       ENDIF
+      !
+	  IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
+	     CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
+	  ENDIF
 !     IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
 !        CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
@@ -410,7 +414,8 @@
 
       ! SIMIP diagnostics
-      diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp
-
-      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
+      diag_fc_bo(:,:) = 0._wp ; diag_fc_su(:,:) = 0._wp
+
+      tau_icebfr(:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
+      cnd_ice   (:,:,:) = 0._wp   ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
       
    END SUBROUTINE diag_set0
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8933)
@@ -527,5 +527,7 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d(1:npti), qns_ice(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), ftr_ice_1d(1:npti), ftr_ice(:,:,kl) )
-         !
+         ! effective conductivity and 1st layer temperature (for Jules coupling)
+         CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
          ! SIMIP diagnostics         
          CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d      (1:npti), t_si(:,:,kl) )
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8933)
@@ -840,20 +840,20 @@
 
       IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF
-      
+         
          CALL ice_thd_enmelt       
-
+         
          !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
          DO ji = 1, npti
             zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
                &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
-               
+            
             IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN
                
-                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-                     zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
-                  ELSE                          ! case T_su = 0degC
-                     zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
-                  ENDIF
-            
+               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
+                  zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
+               ELSE                          ! case T_su = 0degC
+                  zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
+               ENDIF
+               
             ELSE ! RCV CASE
             
@@ -892,11 +892,24 @@
       ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
       !---------------------------------------------------------------------------------------
+      ! effective conductivity and 1st layer temperature (Jules coupling)
+      DO ji = 1, npti
+         IF (h_s_1d(ji) > 0.1 ) THEN 
+             cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
+         ELSE
+             IF (h_i_1d(ji) > 0.1 ) THEN
+                 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
+             ELSE
+                 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1
+             ENDIF
+         ENDIF
+         t1_ice_1d (ji) = ( isnow(ji) * t_s_1d  (ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d  (ji,1) )
+      END DO
       !
       IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode
          !
          ! Restore temperatures to their initial values
-         t_s_1d(1:npti,:)        = ztsold (1:npti,:)
-         t_i_1d(1:npti,:)        = ztiold (1:npti,:)
-         qcn_ice_1d(1:npti)      = fc_su(1:npti)
+         t_s_1d    (1:npti,:) = ztsold(1:npti,:)
+         t_i_1d    (1:npti,:) = ztiold(1:npti,:)
+         qcn_ice_1d(1:npti)   = fc_su (1:npti)
          !  
       ENDIF
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 8933)
@@ -70,4 +70,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: wind speed module at T-point                 [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice       !: sea ice surface skin temperature (on categories)
 #endif
 
@@ -123,13 +125,13 @@
 
 #if defined key_lim3
-      ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
-         &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
-         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl)  , alb_ice (jpi,jpj,jpl) ,   &
-         &      qml_ice(jpi,jpj,jpl)  , qcn_ice(jpi,jpj,jpl)   , qsr_ice_tr(jpi,jpj,jpl),  &
-         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   &
-         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
-         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
-         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce (jpi,jpj)  ,   &
-         &      emp_ice(jpi,jpj)      ,  STAT= ierr(2) )
+      ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice  (jpi,jpj,jpl) ,     &
+         &      qla_ice (jpi,jpj,jpl) , dqla_ice (jpi,jpj,jpl) ,     &
+         &      dqns_ice(jpi,jpj,jpl) , tn_ice   (jpi,jpj,jpl) , alb_ice   (jpi,jpj,jpl) ,   &
+         &      qml_ice (jpi,jpj,jpl) , qcn_ice  (jpi,jpj,jpl) , qsr_ice_tr(jpi,jpj,jpl) ,   &
+         &      utau_ice(jpi,jpj)     , vtau_ice (jpi,jpj)     , wndm_ice  (jpi,jpj)     ,   &
+         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice (jpi,jpj)     ,   &
+         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce  (jpi,jpj)     ,   &
+         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce   (jpi,jpj)     ,   &
+         &      emp_ice (jpi,jpj)     , tsfc_ice (jpi,jpj,jpl) , sstfrz    (jpi,jpj)     , STAT= ierr(2) )
 #endif
 
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8933)
@@ -41,4 +41,5 @@
    LOGICAL , PUBLIC ::   ln_cpl         !: ocean-atmosphere coupled formulation
    LOGICAL , PUBLIC ::   ln_mixcpl      !: ocean-atmosphere forced-coupled mixed formulation
+   LOGICAL , PUBLIC ::   ln_meto_cpl    !: Met Office coupling formulation, with surface exchange carried out in atmosphere
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8933)
@@ -114,6 +114,7 @@
    INTEGER, PARAMETER ::   jpr_isf    = 52
    INTEGER, PARAMETER ::   jpr_icb    = 53
-
-   INTEGER, PARAMETER ::   jprcv      = 53   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_ts_ice = 54   ! Sea ice surface temp
+
+   INTEGER, PARAMETER ::   jprcv      = 55   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -149,5 +150,14 @@
    INTEGER, PARAMETER ::   jps_ocyw   = 31   ! currents on grid 2
    INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level 
-   INTEGER, PARAMETER ::   jpsnd      = 32   ! total number of fields sent 
+   INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes)
+   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area
+   INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness
+   INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity
+   INTEGER, PARAMETER ::   jps_sstfrz = 37   ! sea surface freezing temperature
+   INTEGER, PARAMETER ::   jps_ttilyr = 38   ! sea ice top layer temp
+
+   INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent 
+
+   INTEGER :: nn_cats_cpl   ! number of sea ice categories over which the coupling is carried out
 
    !                                  !!** namelist namsbc_cpl **
@@ -160,7 +170,9 @@
    END TYPE FLD_C
    !                                   ! Send to the atmosphere  
-   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                        
+   TYPE(FLD_C) ::   sn_snd_temp  , sn_snd_alb , sn_snd_thick, sn_snd_crt   , sn_snd_co2,  &
+      &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr         
    !                                   ! Received from the atmosphere
-   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
+   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr,  &
+      &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf                              
    ! Send to waves 
@@ -183,5 +195,5 @@
    REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0) 
 
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) ::   nrcvinfo           ! OASIS info argument
 
    !! Substitution
@@ -236,11 +248,11 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      & 
-         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
-         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   & 
-         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   &
-         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &
-         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp  ,   &
-         &                  sn_rcv_icb , sn_rcv_isf
+      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   & 
+         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   & 
+         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   & 
+         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_wstrf,   &
+         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   &
+         &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   &
+         &                  sn_rcv_icb   , sn_rcv_isf   , nn_cats_cpl  
 
       !!---------------------------------------------------------------------
@@ -249,4 +261,13 @@
       !
       CALL wrk_alloc( jpi,jpj,   zacs, zaos )
+
+      IF( ln_meto_cpl ) THEN
+         tsfc_ice(:,:,:) = 0.0
+         a_ip    (:,:,:) = 0.0
+         v_ip    (:,:,:) = 0.0
+         t1_ice  (:,:,:) = rt0
+         cnd_ice (:,:,:) = 0.0
+         sstfrz  (:,:)   = 0.0
+      ENDIF
 
       ! ================================ !
@@ -294,6 +315,8 @@
          WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
+         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
+         WRITE(numout,*)'      top ice layer temperature       = ', TRIM(sn_snd_ttilyr%cldes), ' (', TRIM(sn_snd_ttilyr%clcat), ')'
          WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
          WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
@@ -304,4 +327,7 @@
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
+         WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes  ), ' (', TRIM(sn_snd_cond%clcat  ), ')'
+         WRITE(numout,*)'      meltponds fraction and depth    = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat  ), ')'
+         WRITE(numout,*)'      sea surface freezing temp       = ', TRIM(sn_snd_sstfrz%cldes), ' (', TRIM(sn_snd_sstfrz%clcat), ')'
          WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')' 
          WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')' 
@@ -312,4 +338,5 @@
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
       ENDIF
 
@@ -501,6 +528,9 @@
       !
       ! non solar sensitivity mandatory for LIM ice model
-      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 3 .AND. nn_components /= jp_iam_sas ) &
-         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
+      IF (.NOT. ln_meto_cpl) THEN
+         IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 3 .AND. nn_components /= jp_iam_sas )  &
+            &   CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
+      ENDIF
+
       ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
       IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
@@ -547,4 +577,13 @@
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
+      !                                                      ! ----------------------------- !
+
+      !!!!! To get NEMO4-LIM working at Met Office
+      srcv(jpr_ts_ice)%clname = 'OTsfIce'
+      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' )   srcv(jpr_ts_ice)%laction = .TRUE.
+      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' )   srcv(jpr_ts_ice)%nct     = nn_cats_cpl
+      IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl
+      !!!!!
+
       !                                                      ! ------------------------- !
       !                                                      !      Wave breaking        !    
@@ -707,11 +746,12 @@
       !                                                      !    Surface temperature    !
       !                                                      ! ------------------------- !
-      ssnd(jps_toce)%clname = 'O_SSTSST'
-      ssnd(jps_tice)%clname = 'O_TepIce'
-      ssnd(jps_tmix)%clname = 'O_TepMix'
+      ssnd(jps_toce)%clname   = 'O_SSTSST'
+      ssnd(jps_tice)%clname   = 'O_TepIce'
+      ssnd(jps_ttilyr)%clname = 'O_TtiLyr'
+      ssnd(jps_tmix)%clname   = 'O_TepMix'
       SELECT CASE( TRIM( sn_snd_temp%cldes ) )
       CASE( 'none'                                 )       ! nothing to do
       CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
-      CASE( 'oce and ice' , 'weighted oce and ice' )
+      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice' )
          ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
          IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
@@ -744,12 +784,18 @@
       !                                                      !  Ice fraction & Thickness ! 
       !                                                      ! ------------------------- !
-      ssnd(jps_fice)%clname = 'OIceFrc'
+      ssnd(jps_fice)%clname  = 'OIceFrc'
       ssnd(jps_ficet)%clname = 'OIceFrcT' 
-      ssnd(jps_hice)%clname = 'OIceTck'
-      ssnd(jps_hsnw)%clname = 'OSnwTck'
+      ssnd(jps_hice)%clname  = 'OIceTck'
+      ssnd(jps_a_p)%clname   = 'OPndFrc'
+      ssnd(jps_ht_p)%clname  = 'OPndTck'
+      ssnd(jps_hsnw)%clname  = 'OSnwTck'
+      ssnd(jps_fice1)%clname = 'OIceFrd'
       IF( k_ice /= 0 ) THEN
          ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
+         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used
+                                                          ! producing atmos-to-ice fluxes
 ! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
-         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
+         IF ( TRIM( sn_snd_thick%clcat  ) == 'yes' ) ssnd(jps_fice)%nct = jpl
+         IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = jpl
       ENDIF
       
@@ -769,4 +815,38 @@
       END SELECT
 
+      !                                                      ! ------------------------- ! 
+      !                                                      ! Ice Meltponds             ! 
+      !                                                      ! ------------------------- ! 
+
+
+      !!!!! Getting NEMO4-LIM to work at Met Office
+      ssnd(jps_a_p)%clname = 'OPndFrc'    
+      ssnd(jps_ht_p)%clname = 'OPndTck'    
+      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) ) 
+      CASE ( 'none' ) 
+         ssnd(jps_a_p)%laction = .FALSE. 
+         ssnd(jps_ht_p)%laction = .FALSE. 
+      CASE ( 'ice only' )  
+         ssnd(jps_a_p)%laction = .TRUE. 
+         ssnd(jps_ht_p)%laction = .TRUE. 
+         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
+            ssnd(jps_a_p)%nct = nn_cats_cpl 
+            ssnd(jps_ht_p)%nct = nn_cats_cpl 
+         ELSE 
+            IF ( jpl > 1 ) THEN 
+               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' ) 
+            ENDIF 
+         ENDIF 
+      CASE ( 'weighted ice' )  
+         ssnd(jps_a_p)%laction = .TRUE. 
+         ssnd(jps_ht_p)%laction = .TRUE. 
+         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
+            ssnd(jps_a_p)%nct = nn_cats_cpl  
+            ssnd(jps_ht_p)%nct = nn_cats_cpl  
+         ENDIF 
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes; '//sn_snd_mpnd%cldes ) 
+      END SELECT 
+      !!!!!
+ 
       !                                                      ! ------------------------- !
       !                                                      !      Surface current      !
@@ -818,5 +898,59 @@
       !                                                      ! ------------------------- !
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
-
+      ! 
+      
+      !!!!! Getting NEMO4-LIM to work at the Met Office
+      !                                                      ! ------------------------- ! 
+      !                                                      ! Sea surface freezing temp ! 
+      !                                                      ! ------------------------- ! 
+      ssnd(jps_sstfrz)%clname = 'O_SSTFrz' ; IF( TRIM(sn_snd_sstfrz%cldes) == 'coupled' )  ssnd(jps_sstfrz)%laction = .TRUE. 
+      !!!!!
+
+      ! 
+      !                                                      ! ------------------------- ! 
+      !                                                      !    Ice conductivity       ! 
+      !                                                      ! ------------------------- ! 
+      ! Note that ultimately we will move to passing an ocean effective conductivity as well so there 
+      ! will be some changes to the parts of the code which currently relate only to ice conductivity 
+
+      ssnd(jps_ttilyr )%clname = 'O_TtiLyr' 
+      SELECT CASE ( TRIM( sn_snd_ttilyr%cldes ) ) 
+      CASE ( 'none' ) 
+         ssnd(jps_ttilyr)%laction = .FALSE. 
+      CASE ( 'ice only' ) 
+         ssnd(jps_ttilyr)%laction = .TRUE. 
+         IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN 
+            ssnd(jps_ttilyr)%nct = nn_cats_cpl 
+         ELSE 
+            IF ( nn_cats_cpl > 1 ) THEN 
+               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_ttilyr%cldes if not exchanging category fields' ) 
+            ENDIF 
+         ENDIF 
+      CASE ( 'weighted ice' ) 
+         ssnd(jps_ttilyr)%laction = .TRUE. 
+         IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl 
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_ttilyr%cldes;'//sn_snd_ttilyr%cldes ) 
+      END SELECT 
+
+      ssnd(jps_kice )%clname = 'OIceKn' 
+      SELECT CASE ( TRIM( sn_snd_cond%cldes ) ) 
+      CASE ( 'none' ) 
+         ssnd(jps_kice)%laction = .FALSE. 
+      CASE ( 'ice only' ) 
+         ssnd(jps_kice)%laction = .TRUE. 
+         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN 
+            ssnd(jps_kice)%nct = nn_cats_cpl 
+         ELSE 
+            IF ( nn_cats_cpl > 1 ) THEN 
+               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' ) 
+            ENDIF 
+         ENDIF 
+      CASE ( 'weighted ice' ) 
+         ssnd(jps_kice)%laction = .TRUE. 
+         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl 
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes;'//sn_snd_cond%cldes ) 
+      END SELECT 
+      ! 
+	    
       !                                                      ! ------------------------- ! 
       !                                                      !     Sea surface height    ! 
@@ -1126,4 +1260,25 @@
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
       ! 
+
+      !!!!! Getting NEMO4-LIM to work at the Met Office
+      !  ! Sea ice surface skin temp: 
+      IF( srcv(jpr_ts_ice)%laction ) THEN 
+        DO jn = 1, jpl 
+          DO jj = 1, jpj 
+            DO ji = 1, jpi 
+              IF (frcv(jpr_ts_ice)%z3(ji,jj,jn) > 0.0) THEN 
+                tsfc_ice(ji,jj,jn) = 0.0 
+              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jn) < -60.0) THEN 
+                tsfc_ice(ji,jj,jn) = -60.0 
+              ELSE 
+                tsfc_ice(ji,jj,jn) = frcv(jpr_ts_ice)%z3(ji,jj,jn) 
+              ENDIF 
+            END DO 
+          END DO 
+        END DO 
+      ENDIF 
+      !!!!!
+
+
       !                                                      ! ========================= ! 
       !                                                      ! Mean Sea Level Pressure   !   (taum) 
@@ -1584,7 +1739,7 @@
       REAL(wp) ::   ztri         ! local scalar
       REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zevap_ice, zdevap_ice
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice    !!gm , zfrqsr_tr_i
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice    !!gm , zfrqsr_tr_i
       !!----------------------------------------------------------------------
       !
@@ -1628,5 +1783,9 @@
 
       ! --- evaporation over ice (kg/m2/s) --- !
-      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
+      DO jl=1,jpl
+         IF (sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl)
+         ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF
+      ENDDO
+
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
       ! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
@@ -1656,15 +1815,15 @@
          tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
          DO jl = 1, jpl
-            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
-            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
+            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:,jl) * zmsk(:,:)
+            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:)    * zmsk(:,:)
          END DO
       ELSE
-         emp_tot(:,:) = zemp_tot(:,:)
-         emp_ice(:,:) = zemp_ice(:,:)
-         emp_oce(:,:) = zemp_oce(:,:)     
-         sprecip(:,:) = zsprecip(:,:)
-         tprecip(:,:) = ztprecip(:,:)
+         emp_tot (:,:)   = zemp_tot (:,:)
+         emp_ice (:,:)   = zemp_ice (:,:)
+         emp_oce (:,:)   = zemp_oce (:,:)     
+         sprecip (:,:)   = zsprecip (:,:)
+         tprecip (:,:)   = ztprecip (:,:)
+         evap_ice(:,:,:) = zevap_ice(:,:,:)
          DO jl = 1, jpl
-            evap_ice (:,:,jl) = zevap_ice (:,:)
             devap_ice(:,:,jl) = zdevap_ice(:,:)
          END DO
@@ -1934,12 +2093,15 @@
          dqns_ice(:,:,:) = zdqns_ice(:,:,:)
       ENDIF
-
-      !                                                      ! ========================= !
-      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
-      !                                                      ! ========================= !
-      CASE ('coupled')
-         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
-         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
-      END SELECT
+      
+      IF( ln_meto_cpl ) THEN
+         !                                                      ! ========================= !
+         SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
+         !                                                      ! ========================= !
+         CASE ('coupled')
+            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i(:,:,:)
+            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i(:,:,:)
+         END SELECT
+      ENDIF
+
 
 #if defined key_lim3
@@ -1950,32 +2112,4 @@
       CASE( np_jules_OFF    )       !==  No Jules coupler  ==!
          !
-!!gm         ! former coding was    
-!!gm         ! Coupled case: since cloud cover is not received from atmosphere 
-!!gm         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-!!gm         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-!!gm         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-!!gm         
-!!gm         ! to retrieve that coding, we needed to access h_i & h_s from here
-!!gm         ! we could even retrieve cloud fraction from the coupler
-!!gm         !
-!!gm         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter
-!!gm         !
-!!gm         DO jl = 1, jpl
-!!gm            DO jj = 1 , jpj
-!!gm               DO ji = 1, jpi
-!!gm                  !              !--- surface transmission parameter (Grenfell Maykut 77) --- !
-!!gm                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice 
-!!gm                  !
-!!gm                  !              ! --- influence of snow and thin ice --- !
-!!gm                  IF ( phs(ji,jj,jl) >= 0.0_wp )   zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque
-!!gm                  IF ( phi(ji,jj,jl) <= 0.1_wp )   zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation
-!!gm               END DO
-!!gm            END DO
-!!gm         END DO
-!!gm         !
-!!gm         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation 
-!!gm         !
-!!gm better coding of the above calculation:
-         !
          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
          ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
@@ -1984,5 +2118,4 @@
          WHERE( phs(:,:,:) >= 0.0_wp )   qsr_ice_tr(:,:,:) = 0._wp            ! snow fully opaque
          WHERE( phi(:,:,:) <= 0.1_wp )   qsr_ice_tr(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
-!!gm end
          !     
       CASE( np_jules_ACTIVE )       !==  Jules coupler is active  ==!
@@ -2064,4 +2197,15 @@
                CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
                END SELECT
+            CASE( 'oce and weighted ice')    ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0  
+               SELECT CASE( sn_snd_temp%clcat ) 
+               CASE( 'yes' )    
+                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
+               CASE( 'no' ) 
+                  ztmp3(:,:,:) = 0.0 
+                  DO jl=1,jpl 
+                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl) 
+                  ENDDO 
+               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 
+               END SELECT 
             CASE( 'mixed oce-ice'        )   
                ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
@@ -2076,4 +2220,18 @@
          IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
       ENDIF
+
+      !!!!! Getting NEMO4-LIM working at Met Office
+      ! Top layer ice temperature
+      IF( ssnd(jps_ttilyr)%laction) THEN
+         SELECT CASE( sn_snd_ttilyr%cldes)
+         CASE ('weighted ice')
+            ztmp3(:,:,1:jpl) = t1_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
+         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ttilyr%cldes' )
+         END SELECT
+         IF( ssnd(jps_ttilyr)%laction )   CALL cpl_snd( jps_ttilyr, isec, ztmp3, info )
+      ENDIF
+      !!!!!
+
+
       !                                                      ! ------------------------- !
       !                                                      !           Albedo          !
@@ -2134,4 +2292,13 @@
          END SELECT
          IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
+      ENDIF
+
+      IF( ssnd(jps_fice1)%laction ) THEN
+         SELECT CASE( sn_snd_thick1%clcat )
+         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
+         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
+         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' )
+         END SELECT
+         CALL cpl_snd( jps_fice1, isec, ztmp3, info )
       ENDIF
       
@@ -2179,4 +2346,50 @@
          IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
       ENDIF
+
+      ! NEMO4 - Jules coupling - Met Office
+      ! Send meltpond fields  
+      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
+         SELECT CASE( sn_snd_mpnd%cldes)  
+         CASE( 'ice only' )  
+            SELECT CASE( sn_snd_mpnd%clcat )  
+            CASE( 'yes' )  
+               ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl)
+               ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)  
+            CASE( 'no' )  
+               ztmp3(:,:,:) = 0.0  
+               ztmp4(:,:,:) = 0.0  
+               DO jl=1,jpl  
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip(:,:,jpl)  
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl) 
+               ENDDO  
+            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' )  
+            END SELECT  
+         CASE( 'default' )    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%cldes' )     
+         END SELECT  
+         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )     
+         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )     
+         ! 
+         ! Send ice effective conductivity 
+         SELECT CASE( sn_snd_cond%cldes) 
+         CASE( 'weighted ice' )    
+            SELECT CASE( sn_snd_cond%clcat ) 
+            CASE( 'yes' )    
+	               ztmp3(:,:,1:jpl) =  cnd_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
+            CASE( 'no' ) 
+               ztmp3(:,:,:) = 0.0 
+               DO jl=1,jpl 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + cnd_ice(:,:,jl) * a_i(:,:,jl) 
+               ENDDO 
+            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' ) 
+            END SELECT 
+         CASE( 'ice only' )    
+           ztmp3(:,:,1:jpl) = cnd_ice(:,:,1:jpl) 
+         END SELECT 
+         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info ) 
+      ENDIF 
+      !   
+      !!!!!
+
+
       !                                                      ! ------------------------- !
       !                                                      !  CO2 flux from PISCES     ! 
@@ -2500,4 +2713,9 @@
       IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
 
+      ! NEMO4 - Jules coupling - Met Office
+      CALL eos_fzp(tsn(:,:,1,jp_sal), sstfrz)
+      ztmp1(:,:) = sstfrz(:,:) + rt0
+      IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info)
+
       CALL wrk_dealloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
       CALL wrk_dealloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8933)
@@ -92,5 +92,5 @@
       NAMELIST/namsbc/ nn_fsbc  ,                                                    &
          &             ln_usr   , ln_flx   , ln_blk       ,                          &
-         &             ln_cpl   , ln_mixcpl, nn_components,                          &
+         &             ln_cpl   , ln_mixcpl, ln_meto_cpl  , nn_components,           &
          &             nn_ice   , ln_ice_embd,                                       &
          &             ln_traqsr, ln_dm2dc ,                                         &
@@ -138,4 +138,5 @@
          WRITE(numout,*) '         mixed forced-coupled     formulation       ln_mixcpl     = ', ln_mixcpl
 !!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist 
+         WRITE(numout,*) '         Met Office coupling specifics              ln_meto_cpl   = ', ln_meto_cpl
          WRITE(numout,*) '         OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis
          WRITE(numout,*) '         components of your executable              nn_components = ', nn_components
