Index: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8878)
+++ branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8879)
@@ -29,12 +29,8 @@
    USE ice            ! ice variables
 #endif
-#if defined key_lim2
-   USE par_ice_2      ! ice parameters
-   USE ice_2          ! ice variables
-#endif
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      ! 
    USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
-   USE albedo         ! 
+   USE ocealb         ! 
    USE eosbn2         ! 
    USE sbcrnf, ONLY : l_rnfcpl
@@ -44,5 +40,5 @@
 #endif
 #if defined key_lim3
-   USE limthd_dh      ! for CALL lim_thd_snwblow
+   USE icethd_dh      ! for CALL ice_thd_snwblow
 #endif
    !
@@ -58,8 +54,8 @@
 
    PUBLIC   sbc_cpl_init      ! routine called by sbcmod.F90
-   PUBLIC   sbc_cpl_rcv       ! routine called by sbc_ice_lim(_2).F90
+   PUBLIC   sbc_cpl_rcv       ! routine called by icestp.F90
    PUBLIC   sbc_cpl_snd       ! routine called by step.F90
-   PUBLIC   sbc_cpl_ice_tau   ! routine called by sbc_ice_lim(_2).F90
-   PUBLIC   sbc_cpl_ice_flx   ! routine called by sbc_ice_lim(_2).F90
+   PUBLIC   sbc_cpl_ice_tau   ! routine called by icestp.F90
+   PUBLIC   sbc_cpl_ice_flx   ! routine called by icestp.F90
    PUBLIC   sbc_cpl_alloc     ! routine called in sbcice_cice.F90
 
@@ -117,6 +113,8 @@
    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 ::   jpr_rcv    = 55   
+
+   INTEGER, PARAMETER ::   jprcv      = 55   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -152,5 +150,13 @@
    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 **
@@ -163,7 +169,7 @@
    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 
@@ -181,5 +187,5 @@
    TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
@@ -205,7 +211,7 @@
       ierr(:) = 0
       !
-      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
+      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
       
-#if ! defined key_lim3 && ! defined key_lim2 && ! defined key_cice
+#if ! defined key_lim3 && ! defined key_cice
       ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
 #endif
@@ -242,11 +248,13 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
+      LOGICAL ::   ln_iceshelf_init_atmos
       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_thick1,  sn_snd_cond, sn_snd_mpnd, sn_snd_sstfrz , sn_rcv_ts_ice,   sn_snd_ttilyr, &
          &                  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
+         &                  sn_rcv_icb , sn_rcv_isf, ln_iceshelf_init_atmos,        nn_cats_cpl  
 
       !!---------------------------------------------------------------------
@@ -274,4 +282,9 @@
          WRITE(numout,*)'~~~~~~~~~~~~'
       ENDIF
+
+      !!!!! Getting NEMO4-LIM working at the Met Office: Hardcode number of ice cats to 5 during the initialisation
+      jpl = nn_cats_cpl
+      !!!!!
+
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
          WRITE(numout,*)'  received fields (mutiple ice categogies)'
@@ -300,6 +313,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 ), ')'
@@ -310,4 +325,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  ), ')' 
@@ -318,4 +336,5 @@
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
       ENDIF
 
@@ -511,6 +530,10 @@
       !
       ! non solar sensitivity mandatory for LIM ice model
-      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4 .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' ) &
@@ -557,4 +580,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 = jpl
+      IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl
+      !!!!!
+
       !                                                      ! ------------------------- !
       !                                                      !      Wave breaking        !    
@@ -719,9 +751,10 @@
       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
@@ -746,7 +779,7 @@
       !     2. receiving mixed oce-ice solar radiation 
       IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
-         CALL albedo_oce( zaos, zacs )
+         CALL oce_alb( zaos, zacs )
          ! Due to lack of information on nebulosity : mean clear/overcast sky
-         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
+         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
       ENDIF
 
@@ -757,9 +790,15 @@
       ssnd(jps_ficet)%clname = 'OIceFrcT' 
       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_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = jpl
       ENDIF
       
@@ -779,4 +818,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 = jpl 
+            ssnd(jps_ht_p)%nct = jpl 
+         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 = jpl  
+            ssnd(jps_ht_p)%nct = jpl  
+         ENDIF 
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes; '//sn_snd_mpnd%cldes ) 
+      END SELECT 
+      !!!!!
+ 
       !                                                      ! ------------------------- !
       !                                                      !      Surface current      !
@@ -828,5 +901,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 = jpl 
+         ELSE 
+            IF ( jpl > 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 = jpl 
+      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 = jpl 
+         ELSE 
+            IF ( jpl > 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 = jpl 
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes;'//sn_snd_cond%cldes ) 
+      END SELECT 
+      ! 
+	    
       !                                                      ! ------------------------- ! 
       !                                                      !     Sea surface height    ! 
@@ -1161,4 +1288,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) 
@@ -1248,5 +1396,5 @@
       IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
          ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
-         ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
+         ub (:,:,1) = ssu_m(:,:)                             ! will be used in icestp in the call of lim_sbc_tau
          un (:,:,1) = ssu_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
          CALL iom_put( 'ssu_m', ssu_m )
@@ -1254,5 +1402,5 @@
       IF( srcv(jpr_ocy1)%laction ) THEN
          ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
-         vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
+         vb (:,:,1) = ssv_m(:,:)                             ! will be used in icestp in the call of lim_sbc_tau
          vn (:,:,1) = ssv_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
          CALL iom_put( 'ssv_m', ssv_m )
@@ -1558,5 +1706,5 @@
    
 
-   SUBROUTINE sbc_cpl_ice_flx( p_frld, palbi, psst, pist )
+   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
@@ -1591,8 +1739,8 @@
       !!
       !! ** Details
-      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
+      !!             qns_tot = (1-a) * qns_oce + a * qns_ice               => provided
       !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
       !!
-      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
+      !!             qsr_tot = (1-a) * qsr_oce + a * qsr_ice               => provided
       !!
       !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce).
@@ -1608,26 +1756,34 @@
       !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
+      REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
       ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
-      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      !
-      INTEGER ::   jl         ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, zicefr, zmsk, zsnw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
+      REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
+      !
+      INTEGER ::   ji,jj,jl         ! dummy loop index
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zdevap_ice
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i, zevap_ice
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, zicefr, zmsk, zsnw )
-      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zdevap_ice )
       CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
+      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i )
+
+      IF (sn_rcv_emp%clcat == 'yes') THEN
+         CALL wrk_alloc( jpi, jpj, jpl, zevap_ice)
+      ELSE
+         CALL wrk_alloc( jpi, jpj, 1, zevap_ice)
+      ENDIF
 
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
-      zicefr(:,:) = 1.- p_frld(:,:)
+      ziceld(:,:) = 1. - picefr(:,:)
       zcptn(:,:) = rcp * sst_m(:,:)
       !
@@ -1645,8 +1801,8 @@
          ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
          zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
-         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
+         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:)
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
-         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
-         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
+         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
+         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
          zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
@@ -1654,16 +1810,23 @@
 
 #if defined key_lim3
-      ! zsnw = snow fraction over ice after wind blowing (=zicefr if no blowing)
-      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw )
+      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
+      zsnw(:,:) = 0._wp  ;  CALL ice_thd_snwblow( ziceld, zsnw )
       
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
-      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( picefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
       zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
 
       ! --- evaporation over ocean (used later for qemp) --- !
-      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:)
 
       ! --- evaporation over ice (kg/m2/s) --- !
-      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
+      IF (sn_rcv_emp%clcat == 'yes') THEN
+         DO jl=1,jpl
+            zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl)
+         ENDDO
+      ELSE
+         zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1)
+      ENDIF
+
       ! 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.
@@ -1693,5 +1856,5 @@
          tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
          DO jl=1,jpl
-            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
+            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:,1) * zmsk(:,:)
             devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
          ENDDO
@@ -1703,5 +1866,10 @@
          tprecip(:,:) =         ztprecip(:,:)
          DO jl=1,jpl
-            evap_ice (:,:,jl) = zevap_ice (:,:)
+            IF (sn_rcv_emp%clcat == 'yes') THEN
+               evap_ice (:,:,jl) = zevap_ice (:,:,jl)
+            ELSE
+               evap_ice (:,:,jl) = zevap_ice (:,:,1)
+            ENDIF
+
             devap_ice(:,:,jl) = zdevap_ice(:,:)
          ENDDO
@@ -1709,5 +1877,5 @@
 
 #else
-      zsnw(:,:) = zicefr(:,:)
+      zsnw(:,:) = picefr(:,:)
       ! --- Continental fluxes --- !
       IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
@@ -1748,7 +1916,7 @@
       IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
       IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
       IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
-         &                                                        - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
+         &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
       ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
       !
@@ -1768,5 +1936,5 @@
          ENDIF
       CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
-         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
+         zqns_tot(:,:) =  ziceld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
@@ -1775,7 +1943,7 @@
             ENDDO
          ELSE
-            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
-               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
             ENDDO
@@ -1785,6 +1953,6 @@
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
          zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
-            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
-            &                                           + pist(:,:,1) * zicefr(:,:) ) )
+            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * ziceld(:,:)   &
+            &                                           + pist(:,:,1) * picefr(:,:) ) )
       END SELECT
       !                                     
@@ -1797,7 +1965,7 @@
 #if defined key_lim3      
       ! --- non solar flux over ocean --- !
-      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
+      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
       zqns_oce = 0._wp
-      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
+      WHERE( ziceld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
 
       ! Heat content per unit mass of snow (J/kg)
@@ -1806,5 +1974,5 @@
       ENDWHERE
       ! Heat content per unit mass of rain (J/kg)
-      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * p_frld(:,:) ) 
+      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * ziceld(:,:) ) 
 
       ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
@@ -1821,5 +1989,5 @@
          &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - lfus )   ! solid precip over ocean + snow melting
       zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - lfus )   ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
-!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
+!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
 !!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice
       
@@ -1854,10 +2022,10 @@
       ! clem: this formulation is certainly wrong... but better than it was...
       zqns_tot(:,:) = zqns_tot(:,:)                            &          ! zqns_tot update over free ocean with:
-         &          - (  p_frld(:,:) * zsprecip(:,:) * lfus )  &          ! remove the latent heat flux of solid precip. melting
+         &          - (  ziceld(:,:) * zsprecip(:,:) * lfus )  &          ! remove the latent heat flux of solid precip. melting
          &          - (  zemp_tot(:,:)                         &          ! remove the heat content of mass flux (assumed to be at SST)
          &             - zemp_ice(:,:) ) * zcptn(:,:) 
 
      IF( ln_mixcpl ) THEN
-         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+         qns_tot(:,:) = qns(:,:) * ziceld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
          qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
          DO jl=1,jpl
@@ -1875,5 +2043,5 @@
       IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea',  sprecip(:,:) * ( zcptsnw(:,:) - Lfus )                           ) ! heat flux from snow (cell average)
       IF( iom_use('hflx_rain_cea') ) CALL iom_put('hflx_rain_cea',( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:)                    ) ! heat flux from rain (cell average)
-      IF( iom_use('hflx_evap_cea') ) CALL iom_put('hflx_evap_cea',(frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) & ! heat flux from from evap (cell average)
+      IF( iom_use('hflx_evap_cea') ) CALL iom_put('hflx_evap_cea',(frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) & ! heat flux from from evap (cell average)
          &                                                        ) * zcptn(:,:) * tmask(:,:,1) )
       IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea',sprecip(:,:) * (zcptsnw(:,:) - Lfus) * (1._wp - zsnw(:,:))   ) ! heat flux from snow (over ocean)
@@ -1899,5 +2067,5 @@
          zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
       CASE( 'oce and ice' )
-         zqsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
+         zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
             DO jl=1,jpl
@@ -1906,7 +2074,7 @@
             ENDDO
          ELSE
-            qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
+            qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
             DO jl=1,jpl
-               zqsr_tot(:,:   ) = zqsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
+               zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             ENDDO
@@ -1918,6 +2086,6 @@
 !       ( see OASIS3 user guide, 5th edition, p39 )
          zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
-            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
-            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
+            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
+            &                     + palbi         (:,:,1) * picefr(:,:) ) )
       END SELECT
       IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
@@ -1930,7 +2098,7 @@
 #if defined key_lim3
       ! --- solar flux over ocean --- !
-      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
+      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
       zqsr_oce = 0._wp
-      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
+      WHERE( ziceld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / ziceld(:,:)
 
       IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
@@ -1939,5 +2107,5 @@
 
       IF( ln_mixcpl ) THEN
-         qsr_tot(:,:) = qsr(:,:) * p_frld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
          qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
          DO jl=1,jpl
@@ -1971,23 +2139,67 @@
       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
-
-      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
-      ! Used for LIM2 and LIM3
-      ! Coupled case: since cloud cover is not received from atmosphere 
-      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-
-      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, zicefr, zmsk, zsnw )
-      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
+      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
+
+      ! --- Transmitted shortwave radiation (W/m2) --- !
+      
+      IF ( nice_jules == 0 ) THEN
+               
+         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter
+     
+         ! former coding was    
+         ! Coupled case: since cloud cover is not received from atmosphere 
+         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+         
+         ! to retrieve that coding, we needed to access h_i & h_s from here
+         ! we could even retrieve cloud fraction from the coupler
+               
+         DO jl = 1, jpl
+            DO jj = 1 , jpj
+               DO ji = 1, jpi
+            
+                  !--- surface transmission parameter (Grenfell Maykut 77) --- !
+                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice 
+               
+                  ! --- influence of snow and thin ice --- !
+                  IF ( phs(ji,jj,jl) >= 0.0_wp ) zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque
+                  IF ( phi(ji,jj,jl) <= 0.1_wp ) zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation
+               END DO
+            END DO
+         END DO
+         
+         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation 
+               
+      ENDIF
+      
+      IF ( nice_jules == 2 ) THEN
+      
+         ! here we must receive the qsr_ice_tr array from the coupler
+         ! for now just assume zero
+         
+         qsr_ice_tr(:,:,:) = 0.0_wp
+      
+      ENDIF
+
+
+
+      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zdevap_ice )
       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
+      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i )
+      IF (sn_rcv_emp%clcat == 'yes') THEN
+         CALL wrk_dealloc( jpi,jpj,jpl,zevap_ice)
+      ELSE
+         CALL wrk_dealloc( jpi,jpj,1,zevap_ice)
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
@@ -2033,5 +2245,5 @@
             ! we must send the surface potential temperature 
             IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
-            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
+            ELSE                   ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
             ENDIF
             !
@@ -2061,4 +2273,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(:,:) 
@@ -2073,4 +2296,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          !
@@ -2086,5 +2323,5 @@
                    ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
                 ELSEWHERE
-                   ztmp1(:,:) = albedo_oce_mix(:,:)
+                   ztmp1(:,:) = alb_oce_mix(:,:)
                 END WHERE
              CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
@@ -2114,5 +2351,5 @@
 
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
-         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
+         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
          DO jl=1,jpl
             ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
@@ -2132,4 +2369,13 @@
          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
       
       ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
@@ -2146,11 +2392,11 @@
             SELECT CASE( sn_snd_thick%clcat )
             CASE( 'yes' )   
-               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
-               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
+               ztmp3(:,:,1:jpl) =  h_i(:,:,1:jpl) * a_i(:,:,1:jpl)
+               ztmp4(:,:,1:jpl) =  h_s(:,:,1:jpl) * a_i(:,:,1:jpl)
             CASE( 'no' )
                ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
                DO jl=1,jpl
-                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
-                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
+                  ztmp3(:,:,1) = ztmp3(:,:,1) + h_i(:,:,jl) * a_i(:,:,jl)
+                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_s(:,:,jl) * a_i(:,:,jl)
                ENDDO
             CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
@@ -2159,10 +2405,10 @@
             SELECT CASE( sn_snd_thick%clcat )
             CASE( 'yes' )
-               ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
-               ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
+               ztmp3(:,:,1:jpl) = h_i(:,:,1:jpl)
+               ztmp4(:,:,1:jpl) = h_s(:,:,1:jpl)
             CASE( 'no' )
                WHERE( SUM( a_i, dim=3 ) /= 0. )
-                  ztmp3(:,:,1) = SUM( ht_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
-                  ztmp4(:,:,1) = SUM( ht_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
+                  ztmp3(:,:,1) = SUM( h_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
+                  ztmp4(:,:,1) = SUM( h_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
                ELSEWHERE
                  ztmp3(:,:,1) = 0.
@@ -2176,4 +2422,50 @@
          IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
       ENDIF
+
+      ! 
+      ! Send meltpond fields  
+      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
+         SELECT CASE( sn_snd_mpnd%cldes)  
+         CASE( 'weighted ice' )  
+            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     ! 
@@ -2545,4 +2837,8 @@
       IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
 
+      CALL eos_fzp(tsn(:,:,1,jp_sal), sstfrz)
+      ztmp1(:,:) = sstfrz(:,:)
+      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 )
