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
