Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8909)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8920)
@@ -72,6 +72,6 @@
       SELECT CASE( ksbc )
          CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
-         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! Bulk formulation
-         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation
+         CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! Bulk         formulation
+         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled      formulation
       END SELECT
 
@@ -134,22 +134,25 @@
 !!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
       alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-      
-      SELECT CASE( ksbc )      !==  fluxes over sea ice  ==!
-      !
-      CASE( jp_usr )                !--- user defined formulation
+      !
+      !
+      SELECT CASE( ksbc )   !== fluxes over sea ice ==!
+      !
+      CASE( jp_usr )              !--- user defined formulation
                                   CALL usrdef_sbc_ice_flx( kt )
-         !
-      CASE( jp_blk )                !--- bulk formulation
-                                  CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    ! 
+      CASE( jp_blk )              !--- bulk formulation
+                                  CALL blk_ice_flx    ( t_su, h_s, h_i, alb_ice )    ! 
          IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
-         IF( nn_flxdist /= -1 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
-         !
-      CASE ( jp_purecpl )           !--- coupled formulation
+         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
+         SELECT CASE( nice_jules )
+         CASE( np_jules_ACTIVE )  !    compute conduction flux and surface temperature (as in Jules surface module)
+                                  CALL blk_ice_qcn    ( nn_monocat, t_su, t_bo, h_s, h_i )
+         END SELECT
+      CASE ( jp_purecpl )         !--- coupled formulation
                                   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
-         IF( nn_flxdist /= -1 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
-         !
-      END SELECT
-
-      IF( iom_use('icealb') ) THEN    !--- output ice albedo
+         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
+      END SELECT
+
+      !--- output ice albedo and surface albedo ---!
+      IF( iom_use('icealb') ) THEN
          WHERE( at_i_b <= epsi06 )   ;   zalb(:,:) = rn_alb_oce
          ELSEWHERE                   ;   zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
@@ -157,6 +160,5 @@
          CALL iom_put( "icealb" , zalb(:,:) )
       ENDIF
-
-      IF( iom_use('albedo') ) THEN    !--- output surface albedo
+      IF( iom_use('albedo') ) THEN
          zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
          CALL iom_put( "albedo" , zalb(:,:) )
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8909)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8920)
@@ -68,7 +68,7 @@
    PUBLIC   sbc_blk       ! called in sbcmod
 #if defined key_lim3
-   PUBLIC   blk_ice_tau   ! routine called in icestp module
-   PUBLIC   blk_ice_flx   ! routine called in icestp module
-   PUBLIC   blk_ice_qcn   ! routine called in icestp module
+   PUBLIC   blk_ice_tau   ! routine called in iceforcing
+   PUBLIC   blk_ice_flx   ! routine called in iceforcing
+   PUBLIC   blk_ice_qcn   ! routine called in iceforcing
 #endif 
 
@@ -1041,16 +1041,13 @@
       CASE ( 1 , 3 )
          !
-         zfac    =   1._wp /  ( rn_cnd_s + rcdic )
-         zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon
-         zfac3   =   2._wp / zepsilon
+         zfac  = 1._wp /  ( rn_cnd_s + rcdic )
+         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
+         zfac3 = 2._wp / zepsilon
          !   
          DO jl = 1, jpl                
             DO jj = 1 , jpj
                DO ji = 1, jpi
-                  !                                ! Effective thickness
-                  zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac
-                  !                                ! Enhanced conduction factor
-                  IF( zhe >=  zfac2 )  &
-                     zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) )
+                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * 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
                END DO
             END DO
@@ -1063,24 +1060,19 @@
       ! -------------------------------------------------------------!
       !
-      zfac   =   rcdic * rn_cnd_s 
-      !                             ! ========================== !
-      DO jl = 1, jpl                !  Loop over ice categories  !
-         !                          ! ========================== !
+      zfac = rcdic * rn_cnd_s
+      !
+      DO jl = 1, jpl
          DO jj = 1 , jpj
             DO ji = 1, jpi
-               !                    ! Effective conductivity of the snow-ice system divided by thickness
-               zkeff_h =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) )
-               !                    ! Store initial surface temperature
-               ztsu    =   ptsu(ji,jj,jl)
-               !                    ! Net initial atmospheric heat flux
-               zqa0    =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)
+               !                    
+               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
+                  &      ( rcdic * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
+               ztsu    = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
+               zqa0    = qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)  ! Net initial atmospheric heat flux
                !
-               DO iter = 1, nit     ! --- Iteration loop
-                   !                      ! Conduction heat flux through snow-ice system (>0 downwards)
-                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) )
-                   !                      ! Surface energy budget
-                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc
-                   !                      ! Temperature update
-                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )
+               DO iter = 1, nit     ! --- Iterative loop
+                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
+                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
+                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
                END DO
                !
