Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/cpl_oasis3.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/cpl_oasis3.F90	(revision 13571)
@@ -165,7 +165,5 @@
       ENDIF
       !
-      ! ... Define the shape for the area that excludes the halo
-      !     For serial configuration (key_mpp_mpi not being active)
-      !     nl* is set to the global values 1 and jp*glo.
+      ! ... Define the shape for the area that excludes the halo as we don't want them to be "seen" by oasis
       !
       ishape(1) = 1
@@ -176,5 +174,5 @@
       ! ... Allocate memory for data exchange
       !
-      ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror)
+      ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror)        ! allocate only inner domain (without halos)
       IF( nerror > 0 ) THEN
          CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN
@@ -182,21 +180,21 @@
       !
       ! -----------------------------------------------------------------
-      ! ... Define the partition 
+      ! ... Define the partition, excluding halos as we don't want them to be "seen" by oasis    
       ! -----------------------------------------------------------------
       
-      paral(1) = 2                                              ! box partitioning
-      paral(2) = jpiglo * (Njs0-1+njmpp-1) + (Nis0-1+nimpp-1)   ! NEMO lower left corner global offset    
-      paral(3) = Ni_0                                           ! local extent in i 
-      paral(4) = Nj_0                                           ! local extent in j
-      paral(5) = jpiglo                                         ! global extent in x
+      paral(1) = 2                                      ! box partitioning
+      paral(2) = Ni0glo * mjg0(nn_hls) + mig0(nn_hls)   ! NEMO lower left corner global offset, without halos 
+      paral(3) = Ni_0                                   ! local extent in i, excluding halos
+      paral(4) = Nj_0                                   ! local extent in j, excluding halos
+      paral(5) = Ni0glo                                 ! global extent in x, excluding halos
       
       IF( sn_cfctl%l_oasout ) THEN
          WRITE(numout,*) ' multiexchg: paral (1:5)', paral
-         WRITE(numout,*) ' multiexchg: jpi, jpj =', jpi, jpj
+         WRITE(numout,*) ' multiexchg: Ni_0, Nj_0 =', Ni_0, Nj_0
          WRITE(numout,*) ' multiexchg: Nis0, Nie0, nimpp =', Nis0, Nie0, nimpp
          WRITE(numout,*) ' multiexchg: Njs0, Nje0, njmpp =', Njs0, Nje0, njmpp
       ENDIF
    
-      CALL oasis_def_partition ( id_part, paral, nerror, jpiglo*jpjglo )
+      CALL oasis_def_partition ( id_part, paral, nerror, Ni0glo*Nj0glo )   ! global number of points, excluding halos
       !
       ! ... Announce send variables. 
@@ -327,5 +325,5 @@
          DO jm = 1, ssnd(kid)%ncplmodel
         
-            IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN
+            IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN   ! exclude halos from data sent to oasis
                CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0:Nie0, Njs0:Nje0,jc), kinfo )
                
@@ -386,7 +384,8 @@
                   &        kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
                
-               IF ( sn_cfctl%l_oasout )   WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
+               IF ( sn_cfctl%l_oasout )   &
+                  &  WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
                
-               IF( llaction ) THEN
+               IF( llaction ) THEN   ! data received from oasis do not include halos
                   
                   kinfo = OASIS_Rcv
@@ -417,6 +416,5 @@
          ENDDO
 
-         !--- Fill the overlap areas and extra hallows (mpp)
-         !--- check periodicity conditions (all cases)
+         !--- we must call lbc_lnk to fill the halos that where not received.
          IF( .NOT. ll_1st ) THEN
             CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )   
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/fldread.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/fldread.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/fldread.F90	(revision 13571)
@@ -216,5 +216,5 @@
                   WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &            
                      & sd(jf)%nrec(1,ibb), sd(jf)%nrec(1,iaa), REAL(sd(jf)%nrec(2,ibb),wp)/rday, REAL(sd(jf)%nrec(2,iaa),wp)/rday
-                  WRITE(numout, *) '      zt_offset is : ',zt_offset
+                  IF( zt_offset /= 0._wp )   WRITE(numout, *) '      zt_offset is : ', zt_offset
                ENDIF
                ! temporal interpolation weights
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbc_ice.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbc_ice.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbc_ice.F90	(revision 13571)
@@ -69,5 +69,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(:,:)   ::   sstfrz         !: sea surface freezing temperature            [degC]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rCdU_ice       !: ice-ocean drag at T-point (<0)               [m/s]
 #endif
 
@@ -89,5 +90,5 @@
    ! variables used in the coupled interface
    INTEGER , PUBLIC, PARAMETER ::   jpl = ncat
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice 
    
    ! already defined in ice.F90 for SI3
@@ -98,5 +99,5 @@
 #endif
 
-   REAL(wp), PUBLIC, SAVE ::   cldf_ice = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-]
+   REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-]
 
    !! arrays relating to embedding ice in the ocean
@@ -131,5 +132,5 @@
          &      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)     , sstfrz   (jpi,jpj)     , STAT= ierr(2) )
+         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , rCdU_ice   (jpi,jpj)     , STAT= ierr(2) )
 #endif
 
@@ -167,5 +168,5 @@
    LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .FALSE.  !: no SI3 ice model
    LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE ice model
-   REAL(wp)        , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-]
+   REAL(wp)        , PUBLIC, PARAMETER ::   pp_cldf    = 0.81     !: cloud fraction over sea ice, summer CLIO value   [-]
    INTEGER         , PUBLIC, PARAMETER ::   jpl = 1 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbc_oce.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbc_oce.F90	(revision 13571)
@@ -136,4 +136,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra         !: cloud cover (fraction of cloud in a gridcell) [-]
 
    !!---------------------------------------------------------------------
@@ -188,5 +189,5 @@
       !
       ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
-         &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) ,                       &
+         &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) , cloud_fra(jpi,jpj),   &
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      &
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk.F90	(revision 13571)
@@ -44,6 +44,6 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_si3
-   USE ice     , ONLY :   jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif
-   USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice
+   USE icevar         ! for CALL ice_var_snwblow
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
@@ -87,7 +87,8 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component)
    !                                                !          seen by the atmospheric forcing (m/s) at T-point
-   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 12   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
-   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 13   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
-   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 13   ! maximum number of files to read
+   INTEGER , PUBLIC, PARAMETER ::   jp_cc    = 12   ! index of cloud cover                     (-)      range:0-1
+   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 13   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 14   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 14   ! maximum number of files to read
 
    ! Warning: keep this structure allocatable for Agrif...
@@ -175,9 +176,9 @@
       TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        "
       TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        "
-      TYPE(FLD_N) ::   sn_hpgi, sn_hpgj                        !       "                        "
+      TYPE(FLD_N) ::   sn_cc, sn_hpgi, sn_hpgj                 !       "                        "
       INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable
       NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     &
-         &                 sn_hpgi, sn_hpgj,                                          &
+         &                 sn_cc, sn_hpgi, sn_hpgj,                                   &
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , rn_zqt, rn_zu,                                    &
@@ -260,5 +261,5 @@
       slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi
       slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow
-      slf_i(jp_slp  ) = sn_slp
+      slf_i(jp_slp  ) = sn_slp     ;   slf_i(jp_cc   ) = sn_cc
       slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm
       slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj
@@ -289,10 +290,16 @@
          !
          IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default)
-            IF(     jfpr == jp_slp  ) THEN
+            IF(     jfpr == jp_slp ) THEN
                sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa
             ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN
                sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents
-            ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN
-               DEALLOCATE( sf(jfpr)%fnow )              ! deallocate as not used in this case
+            ELSEIF( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) THEN
+               IF( .NOT. ln_abl ) THEN
+                  DEALLOCATE( sf(jfpr)%fnow )   ! deallocate as not used in this case
+               ELSE
+                  sf(jfpr)%fnow(:,:,1:ipka) = 0._wp
+               ENDIF
+            ELSEIF( jfpr == jp_cc  ) THEN
+               sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf
             ELSE
                WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr
@@ -303,6 +310,6 @@
             !
             IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   &
-               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
-               &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )
+         &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
+         &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )
          ENDIF
       END DO
@@ -559,4 +566,7 @@
       ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!)
 
+      ! --- cloud cover --- !
+      cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
+
       ! ----------------------------------------------------------------------------- !
       !      0   Wind components and module at T-point relative to the moving ocean   !
@@ -568,5 +578,5 @@
       zwnd_j(:,:) = 0._wp
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
          zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
@@ -576,5 +586,5 @@
 #else
       ! ... scalar wind module at T-point (not masked)
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
       END_2D
@@ -628,5 +638,5 @@
          !     use scalar version of gamma_moist() ...
          IF( ln_tpot ) THEN
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt
             END_2D
@@ -690,5 +700,5 @@
 
       IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zztmp = zU_zu(ji,jj)
             wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod
@@ -710,5 +720,5 @@
          pevp(:,:) = pevp(:,:) * tmask(:,:,1)
 
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF( wndm(ji,jj) > 0._wp ) THEN
                zztmp = taum(ji,jj) / wndm(ji,jj)
@@ -828,5 +838,5 @@
 
       ! use scalar version of L_vap() for AGRIF compatibility
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update
       END_2D
@@ -933,5 +943,5 @@
       ! ------------------------------------------------------------ !
       ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
       END_2D
@@ -978,5 +988,5 @@
          zztmp1 = 11637800.0_wp
          zztmp2 =    -5897.8_wp
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             pcd_dui(ji,jj) = zcd_dui (ji,jj)
             pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
@@ -1019,5 +1029,4 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      -
-      REAL(wp) ::   zfr1, zfr2               ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
@@ -1028,4 +1037,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
       !!---------------------------------------------------------------------
       !
@@ -1112,5 +1122,5 @@
       ! --- evaporation minus precipitation --- !
       zsnw(:,:) = 0._wp
-      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
+      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
       emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
@@ -1139,17 +1149,27 @@
       END DO
 
-      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
-      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
-      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
-      !
-      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
-         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
-      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
-         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
-      ELSEWHERE                                                         ! zero when hs>0
-         qtr_ice_top(:,:,:) = 0._wp
-      END WHERE
-      !
-
+      ! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- !
+      IF( nn_qtrice == 0 ) THEN
+         ! formulation derived from Grenfell and Maykut (1977), where transmission rate
+         !    1) depends on cloudiness
+         !    2) is 0 when there is any snow
+         !    3) tends to 1 for thin ice
+         ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
+         DO jl = 1, jpl
+            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm  
+               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm
+               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
+            ELSEWHERE                                                         ! zero when hs>0
+               qtr_ice_top(:,:,jl) = 0._wp 
+            END WHERE
+         ENDDO
+      ELSEIF( nn_qtrice == 1 ) THEN
+         ! formulation is derived from the thesis of M. Lebrun (2019).
+         !    It represents the best fit using several sets of observations
+         !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
+         qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:)
+      ENDIF
+      !
       IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
          ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
@@ -1233,5 +1253,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
                IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
@@ -1248,5 +1268,5 @@
       !
       DO jl = 1, jpl
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             !
             zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 13571)
@@ -394,5 +394,5 @@
       !!-------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zw = pwnd(ji,jj)   ! wind speed
@@ -430,5 +430,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zta = pzeta(ji,jj)
@@ -481,5 +481,5 @@
       REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zta = pzeta(ji,jj)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 13571)
@@ -430,5 +430,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zta = pzeta(ji,jj)
@@ -481,5 +481,5 @@
       REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zta = pzeta(ji,jj)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 13571)
@@ -410,5 +410,5 @@
       REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab
       !!----------------------------------------------------------------------------------
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
@@ -455,5 +455,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
       !
       zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 13571)
@@ -241,5 +241,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          !
          zw  = pw10(ji,jj)
@@ -277,5 +277,5 @@
       REAL(wp) :: zx2, zx, zstab   ! local scalars
       !!----------------------------------------------------------------------------------
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
          zx2 = MAX( zx2 , 1._wp )
@@ -308,5 +308,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
          zx2 = MAX( zx2 , 1._wp )
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_phy.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_phy.F90	(revision 13571)
@@ -181,5 +181,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
          ztc2 = ztc*ztc
@@ -270,5 +270,5 @@
       INTEGER  ::   ji, jj         ! dummy loop indices
       !!----------------------------------------------------------------------------------
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
       END_2D
@@ -315,5 +315,5 @@
       !!-------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          !
          zqa = (1._wp + rctv0*pqa(ji,jj))
@@ -351,5 +351,5 @@
       !!-------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          !
          zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer...
@@ -448,5 +448,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          !
          ze_sat =  e_sat_sclr( ptak(ji,jj) )
@@ -473,5 +473,5 @@
       !!----------------------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
          q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
@@ -511,5 +511,5 @@
       INTEGER  ::   ji, jj     ! dummy loop indices
       !!----------------------------------------------------------------------------------
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
          zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
@@ -621,5 +621,5 @@
       IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
 
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
          CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_skin_coare.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 13571)
@@ -89,5 +89,5 @@
       REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus
       !!---------------------------------------------------------------------
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
          zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
@@ -156,5 +156,5 @@
       ztime = REAL(nsec_day,wp)/(24._wp*3600._wp) ! time of current time step since 00:00 for current day (UTC) -> ztime = 0 -> 00:00 / ztime = 0.5 -> 12:00 ...
 
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
          l_exit       = .FALSE.
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_skin_ecmwf.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 13571)
@@ -95,5 +95,5 @@
       REAL(wp) :: zQabs, zdlt, zfr, zalfa, zus
       !!---------------------------------------------------------------------
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
          zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
@@ -173,5 +173,5 @@
       IF( PRESENT(pustk) ) l_pustk_known = .TRUE.
 
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
          zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbccpl.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbccpl.F90	(revision 13571)
@@ -41,5 +41,5 @@
 #endif
 #if defined key_si3
-   USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icevar         ! for CALL ice_var_snwblow
 #endif
    !
@@ -48,4 +48,8 @@
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+
+#if defined key_oasis3 
+   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut 
+#endif 
 
    IMPLICIT NONE
@@ -152,5 +156,5 @@
    INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level 
    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_a_p    = 34   ! meltpond area fraction
    INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness
    INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity
@@ -159,4 +163,12 @@
 
    INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent 
+
+#if ! defined key_oasis3 
+   ! Dummy variables to enable compilation when oasis3 is not being used 
+   INTEGER                    ::   OASIS_Sent        = -1 
+   INTEGER                    ::   OASIS_SentOut     = -1 
+   INTEGER                    ::   OASIS_ToRest      = -1 
+   INTEGER                    ::   OASIS_ToRestOut   = -1 
+#endif 
 
    !                                  !!** namelist namsbc_cpl **
@@ -184,4 +196,6 @@
    LOGICAL     ::   ln_usecplmask         !  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)
+   LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration) 
+
    TYPE ::   DYNARR     
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3   
@@ -191,4 +205,7 @@
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+#if defined key_si3 || defined key_cice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time
+#endif
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
@@ -211,5 +228,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(4)
+      INTEGER :: ierr(5)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -221,6 +238,9 @@
 #endif
       ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
-      !
-      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
+#if defined key_si3 || defined key_cice
+      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(4) )
+#endif
+      !
+      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(5) )
 
       sbc_cpl_alloc = MAXVAL( ierr )
@@ -249,14 +269,14 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   & 
+      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             &
+         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  & 
          &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  & 
-         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   & 
-         &                  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_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  & 
          &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  &
-         &                  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_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   &
+         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  &
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
+         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
          &                  sn_rcv_ts_ice
-
       !!---------------------------------------------------------------------
       !
@@ -278,4 +298,8 @@
       ENDIF
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
+         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
+         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux
+         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
          WRITE(numout,*)'  received fields (mutiple ice categogies)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
@@ -326,7 +350,4 @@
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor 
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
-         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
-         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
-         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
       ENDIF
 
@@ -367,5 +388,5 @@
       IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  &
            .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled
-
+      !
       IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
       
@@ -698,5 +719,5 @@
          ! Change first letter to couple with atmosphere if already coupled OPA
          ! this is nedeed as each variable name used in the namcouple must be unique:
-         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
+         ! for example O_Runoff received by OPA from SAS and therefore S_Runoff received by SAS from the Atmosphere
          DO jn = 1, jprcv
             IF( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
@@ -822,4 +843,8 @@
       END SELECT
 
+      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office)
+#if defined key_si3 || defined key_cice
+       a_i_last_couple(:,:,:) = 0._wp
+#endif
       !                                                      ! ------------------------- ! 
       !                                                      !      Ice Meltponds        ! 
@@ -1110,5 +1135,5 @@
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
       REAL(wp) ::   zzx, zzy               ! temporary variables
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
       !!----------------------------------------------------------------------
       !
@@ -1170,5 +1195,5 @@
             !                              
             IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 0, 0 )                                        ! T ==> (U,V)
                   frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
                   frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
@@ -1224,5 +1249,17 @@
          ENDIF
       ENDIF
-
+!!$      !                                                      ! ========================= !
+!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      !
+!!$      !                                                      ! ========================= !
+!!$      cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1)
+!!$      END SELECT
+!!$
+      zcloud_fra(:,:) = pp_cldf   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+      IF( ln_mixcpl ) THEN
+         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:)
+      ELSE
+         cloud_fra(:,:) = zcloud_fra(:,:)
+      ENDIF
+      !                                                      ! ========================= !
       ! u(v)tau and taum will be modified by ice model
       ! -> need to be reset before each call of the ice/fsbc      
@@ -1549,5 +1586,5 @@
             p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
          CASE( 'T' )
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )                    ! T ==> (U,V)
                ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
                zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
@@ -1623,9 +1660,10 @@
       !
       INTEGER  ::   ji, jj, jl   ! dummy loop index
-      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, 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)     ::   zevap_ice_total
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
       !!----------------------------------------------------------------------
       !
@@ -1647,5 +1685,4 @@
          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) ) * picefr(:,:)
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
@@ -1659,6 +1696,45 @@
 
 #if defined key_si3
+
+      ! --- evaporation over ice (kg/m2/s) --- !
+      IF (ln_scale_ice_flux) THEN ! typically met-office requirements
+         IF (sn_rcv_emp%clcat == 'yes') THEN
+            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp
+            END WHERE
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
+            END WHERE
+         ELSE
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp
+            END WHERE
+            zevap_ice_total(:,:) = zevap_ice(:,:,1)
+            DO jl = 2, jpl
+               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
+            ENDDO
+         ENDIF
+      ELSE
+         IF (sn_rcv_emp%clcat == 'yes') THEN
+            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
+            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
+            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
+            END WHERE
+         ELSE
+            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1)
+            zevap_ice_total(:,:) = zevap_ice(:,:,1)
+            DO jl = 2, jpl
+               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
+            ENDDO
+         ENDIF
+      ENDIF
+
+      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
+         ! For conservative case zemp_ice has not been defined yet. Do it now.
+         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
+      ENDIF
+
       ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
-      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw )
+      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw )
       
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
@@ -1667,11 +1743,5 @@
 
       ! --- evaporation over ocean (used later for qemp) --- !
-      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:)
-
-      ! --- evaporation over ice (kg/m2/s) --- !
-      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
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
 
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
@@ -1751,15 +1821,15 @@
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
 !!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
-      IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
-      IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
-      IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
-      IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
-      IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
-      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('rain_ao_cea') )  CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (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) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
+      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
+      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
+      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
+      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
+      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
+      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('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (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) * 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
       !
@@ -1769,4 +1839,7 @@
       CASE( 'oce only' )         ! the required field is directly provided
          zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero
+         ! here so the only flux is the ocean only one.
+         zqns_ice(:,:,:) = 0._wp 
       CASE( 'conservative' )     ! the required fields are directly provided
          zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
@@ -1798,5 +1871,5 @@
                zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    &
                   &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
-                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
             END DO
          ELSE
@@ -1804,5 +1877,5 @@
                zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    &
                   &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
-                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
             END DO
          ENDIF
@@ -1910,4 +1983,7 @@
       CASE( 'oce only' )
          zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
+         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero
+         ! here so the only flux is the ocean only one.
+         zqsr_ice(:,:,:) = 0._wp
       CASE( 'conservative' )
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1995,4 +2071,6 @@
             ENDDO
          ENDIF
+      CASE( 'none' ) 
+         zdqns_ice(:,:,:) = 0._wp
       END SELECT
       
@@ -2010,9 +2088,12 @@
       !                                                      ! ========================= !
       CASE ('coupled')
-         IF( ln_mixcpl ) THEN
-            DO jl=1,jpl
-               qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:)
-               qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:)
-            ENDDO
+         IF (ln_scale_ice_flux) THEN
+            WHERE( a_i(:,:,:) > 1.e-10_wp )
+               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE
+               qml_ice(:,:,:) = 0.0_wp
+               qcn_ice(:,:,:) = 0.0_wp
+            END WHERE
          ELSE
             qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
@@ -2025,19 +2106,32 @@
       IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==!
          !
-         !                    ! ===> 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 when hi>10cm (Grenfell Maykut 77)
-         !
-         WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
-            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) )
-         ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm
-            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri
-         ELSEWHERE                                                         ! zero when hs>0
-            zqtr_ice_top(:,:,:) = 0._wp
-         END WHERE
+         IF( nn_qtrice == 0 ) THEN
+            ! formulation derived from Grenfell and Maykut (1977), where transmission rate
+            !    1) depends on cloudiness
+            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
+            !    2) is 0 when there is any snow
+            !    3) tends to 1 for thin ice
+            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
+            DO jl = 1, jpl
+               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
+               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm
+                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:)
+               ELSEWHERE                                                           ! zero when hs>0
+                  zqtr_ice_top(:,:,jl) = 0._wp 
+               END WHERE
+            ENDDO
+         ELSEIF( nn_qtrice == 1 ) THEN
+            ! formulation is derived from the thesis of M. Lebrun (2019).
+            !    It represents the best fit using several sets of observations
+            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
+            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:)
+         ENDIF
          !     
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
          !
-         !                    ! ===> here we must receive the qtr_ice_top array from the coupler
-         !                           for now just assume zero (fully opaque ice)
+         !          ! ===> here we must receive the qtr_ice_top array from the coupler
+         !                 for now just assume zero (fully opaque ice)
          zqtr_ice_top(:,:,:) = 0._wp
          !
@@ -2096,4 +2190,5 @@
       !
       isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges
+      info = OASIS_idle
 
       zfr_l(:,:) = 1.- fr_i(:,:)
@@ -2234,4 +2329,16 @@
       ENDIF
 
+#if defined key_si3 || defined key_cice
+      ! If this coupling was successful then save ice fraction for use between coupling points. 
+      ! This is needed for some calculations where the ice fraction at the last coupling point 
+      ! is needed. 
+      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. & 
+         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN 
+         IF ( sn_snd_thick%clcat == 'yes' ) THEN 
+           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl)
+         ENDIF
+      ENDIF
+#endif
+
       IF( ssnd(jps_fice1)%laction ) THEN
          SELECT CASE( sn_snd_thick1%clcat )
@@ -2297,5 +2404,5 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl)
+               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl)
                ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)  
             CASE( 'no' )  
@@ -2303,6 +2410,6 @@
                ztmp4(:,:,:) = 0.0  
                DO jl=1,jpl  
-                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)  
-                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)
                ENDDO  
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )  
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcdcy.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcdcy.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcdcy.F90	(revision 13571)
@@ -110,5 +110,5 @@
 
       imask_night(:,:) = 0
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          ztmpm = 0._wp
          IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h
@@ -193,5 +193,5 @@
 
          zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             ztmp = rad * gphit(ji,jj)
             raa(ji,jj) = SIN( ztmp ) * zsin
@@ -202,5 +202,5 @@
          ! rab to test if the day time is equal to 0, less than 24h of full day
          rab(:,:) = -raa(:,:) / rbb(:,:)
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
                ! When is it night?
@@ -226,5 +226,5 @@
          !         Avoid possible infinite scaling factor, associated with very short daylight
          !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040)
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
                rscal(ji,jj) = 0.0_wp
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcflx.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcflx.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcflx.F90	(revision 13571)
@@ -29,5 +29,4 @@
    PUBLIC sbc_flx       ! routine called by step.F90
 
-   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
    INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
@@ -35,4 +34,6 @@
    INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
    INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
+ !!INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux flux
+   INTEGER , PARAMETER ::   jpfld   = 5 !! 6 ! maximum number of files to read 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
 
@@ -59,4 +60,5 @@
       !!                   net downward radiative flux            qsr   (watt/m2)
       !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s)
+      !!                   salt flux                              sfx   (pss*dh*rho/dt => g/m2/s)
       !!
       !!      CAUTION :  - never mask the surface stress fields
@@ -71,5 +73,5 @@
       !!              - emp         upward mass flux (evap. - precip.)
       !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
-      !!                            if ice is present
+      !!                            if ice
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
@@ -85,6 +87,6 @@
       CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
-      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  ! informations about the fields to be read
-      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp
+      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx ! informations about the fields to be read
+      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx
       !!---------------------------------------------------------------------
       !
@@ -105,5 +107,5 @@
          slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
          slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
-         slf_i(jp_emp ) = sn_emp
+         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx
          !
          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
@@ -118,6 +120,4 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
          !
-         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present)
-         !
       ENDIF
 
@@ -126,31 +126,26 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
 
-         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
-         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
+         IF( ln_dm2dc ) THEN   ! modify now Qsr to include the diurnal cycle
+            qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(ji,jj,1)
+         ELSE
+            DO_2D( 0, 0, 0, 0 )
+               qsr(ji,jj) =          sf(jp_qsr)%fnow(ji,jj,1)   * tmask(ji,jj,1)
+            END_2D
          ENDIF
-         DO_2D( 1, 1, 1, 1 )
-            utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
-            vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
-            qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
-            emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
+         DO_2D( 0, 0, 0, 0 )                                      ! set the ocean fluxes from read fields
+            utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
+            vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
+            qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
+            emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
+            !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1) 
          END_2D
          !                                                        ! add to qns the heat due to e-p
-         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
+         !!clem: I do not think it is needed
+         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
          !
-         qns(:,:) = qns(:,:) * tmask(:,:,1)
-         emp(:,:) = emp(:,:) * tmask(:,:,1)
+         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x) 
+         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
+            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
          !
-         !                                                        ! module of wind stress and wind speed at T-point
-         zcoef = 1. / ( zrhoa * zcdrag )
-         DO_2D( 0, 0, 0, 0 )
-            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-            taum(ji,jj) = zmod
-            wndm(ji,jj) = SQRT( zmod * zcoef )
-         END_2D
-         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
-         CALL lbc_lnk( 'sbcflx', taum(:,:), 'T', 1.0_wp )   ;   CALL lbc_lnk( 'sbcflx', wndm(:,:), 'T', 1.0_wp )
-
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
             WRITE(numout,*) 
@@ -166,4 +161,16 @@
          !
       ENDIF
+      !                                                           ! module of wind stress and wind speed at T-point
+      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+      zcoef = 1. / ( zrhoa * zcdrag )
+      DO_2D( 0, 0, 0, 0 )
+         ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
+         zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
+         zmod = 0.5_wp * SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
+         taum(ji,jj) = zmod
+         wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
+      END_2D
+      !
+      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
       !
    END SUBROUTINE sbc_flx
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90	(revision 13571)
@@ -99,6 +99,6 @@
          &             nn_ice   , ln_ice_embd,                                       &
          &             ln_traqsr, ln_dm2dc ,                                         &
-         &             ln_rnf   , nn_fwb     , ln_ssr   , ln_apr_dyn,              &
-         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor  ,     &
+         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_apr_dyn,                  &
+         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc , ln_stcor  ,      &
          &             ln_tauw  , nn_lsm, nn_sdrift
       !!----------------------------------------------------------------------
@@ -119,4 +119,7 @@
 #if defined key_mpp_mpi
       ncom_fsbc = nn_fsbc    ! make nn_fsbc available for lib_mpp
+#endif
+#if ! defined key_si3
+      IF( nn_ice == 2 )    nn_ice = 0  ! without key key_si3 you cannot use si3...
 #endif
       !
@@ -243,5 +246,4 @@
       ENDIF
       !
-
       IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
          IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
@@ -250,4 +252,5 @@
       sfx   (:,:) = 0._wp           !* salt flux due to freezing/melting
       fmmflx(:,:) = 0._wp           !* freezing minus melting flux
+      cloud_fra(:,:) = pp_cldf      !* cloud fraction over sea ice (used in si3)
 
       taum(:,:) = 0._wp             !* wind stress module (needed in GLS in case of reduced restart)
@@ -334,9 +337,9 @@
       IF( l_sbc_clo   )   CALL sbc_clo_init              ! closed sea surface initialisation
       !
-      IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization
-
-      IF( ln_abl      )   CALL sbc_abl_init            ! Atmospheric Boundary Layer (ABL)
-
-      IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization
+      IF( ln_blk      )   CALL sbc_blk_init              ! bulk formulae initialization
+
+      IF( ln_abl      )   CALL sbc_abl_init              ! Atmospheric Boundary Layer (ABL)
+
+      IF( ln_ssr      )   CALL sbc_ssr_init              ! Sea-Surface Restoring initialization
       !
       !
@@ -563,7 +566,4 @@
       ENDIF
       !
-      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice)
-      CALL iom_put( "vtau", vtau )   ! j-wind stress
-      !
       IF(sn_cfctl%l_prtctl) THEN     ! print mean trends (used for debugging)
          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask )
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcrnf.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcrnf.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcrnf.F90	(revision 13571)
@@ -215,7 +215,7 @@
             END_2D
          ELSE                    !* variable volume case
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )              ! update the depth over which runoffs are distributed
                h_rnf(ji,jj) = 0._wp
-               DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
+               DO jk = 1, nk_rnf(ji,jj)                             ! recalculates h_rnf to be the depth in metres
                   h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box
                END DO
@@ -374,5 +374,5 @@
             ENDIF
          END_2D
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                           ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
@@ -404,5 +404,5 @@
          WHERE( zrnfcl(:,:,1) > 0._wp )  h_rnf(:,:) = zacoef * zrnfcl(:,:,1)   ! compute depth for all runoffs
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                ! take in account min depth of ocean rn_hmin
             IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
                jk = mbkt(ji,jj)
@@ -423,5 +423,5 @@
          END_2D
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                          ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcwave.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcwave.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcwave.F90	(revision 13571)
@@ -106,5 +106,5 @@
       !!---------------------------------------------------------------------
       !
-      ALLOCATE( ze3divh(jpi,jpj,jpk) )
+      ALLOCATE( ze3divh(jpi,jpj,jpkm1) )   ! jpkm1 -> avoid lbc_lnk on jpk that is not defined
       ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) )
       !
@@ -121,5 +121,5 @@
             zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
          END_2D
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )          ! exp. wave number & Stokes drift velocity at u- & v-points
             zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
             zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
@@ -164,5 +164,5 @@
          zsqrtpi = SQRT(rpi)
          z_two_thirds = 2.0_wp / 3.0_wp
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! exp. wave number & Stokes drift velocity at u- & v-points
             zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth
             zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth
@@ -204,5 +204,5 @@
       !                       !==  vertical Stokes Drift 3D velocity  ==!
       !
-      DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+      DO_3D( 0, 1, 0, 1, 1, jpkm1 )    ! Horizontal e3*divergence
          ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * usd(ji  ,jj,jk)    &
             &                 - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk)    &
