Index: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8738)
+++ branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8752)
@@ -17,4 +17,5 @@
    !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
+   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle) 
    !!----------------------------------------------------------------------
 
@@ -42,4 +43,5 @@
    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
    USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icethd_zdf, ONLY:   rn_cnd_s  ! for blk_ice_qcn
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
@@ -64,5 +66,6 @@
    PUBLIC   blk_ice_tau   ! routine called in icestp module
    PUBLIC   blk_ice_flx   ! routine called in icestp module
-#endif
+   PUBLIC   blk_ice_qcn   ! routine called in icestp module
+#endif 
 
 !!Lolo: should ultimately be moved in the module with all physical constants ?
@@ -697,5 +700,5 @@
 
 
-   SUBROUTINE blk_ice_flx( ptsu, palb )
+   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_ice_flx  ***
@@ -710,4 +713,6 @@
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
       !!
@@ -716,4 +721,6 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, z1_lsub           !   -      -
+      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice
+      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice
@@ -722,4 +729,5 @@
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3)
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa
+
       !!---------------------------------------------------------------------
       !
@@ -830,13 +838,36 @@
       CALL wrk_dealloc( jpi,jpj,   zevap, zsnw )
 
-      !--------------------------------------------------------------------
-      ! FRACTIONs of net shortwave radiation which is not absorbed in the
-      ! thin surface layer and penetrates inside the ice cover
-      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
-      !
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-      !
-      !
+      ! --- absorbed and transmitted shortwave radiation (W/m2) --- !
+      !
+      ! former coding was
+      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+
+      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1
+
+      qsr_ice_tr(:,:,:) = 0._wp
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1
+               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
+
+               IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN 
+                  zfrqsr_tr_i  = zfr1            + zfac * zfr2
+               ELSE 
+                  zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque
+               ENDIF
+
+               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation 
+
+            END DO
+         END DO
+      END DO
+                                                                                      
+
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
@@ -854,4 +885,123 @@
 
    END SUBROUTINE blk_ice_flx
+   
+   
+
+   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi )
+   
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_qcn  ***
+      !!
+      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
+      !!                to force sea ice / snow thermodynamics
+      !!                in the case JULES coupler is emulated
+      !!                
+      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
+      !!                following the 0-layer Semtner (1976) approach
+      !!
+      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
+      !!              - qcn_ice : surface inner conduction flux (W/m2)
+      !!
+      !!---------------------------------------------------------------------
+      !!
+      INTEGER, INTENT(in)                        ::   k_monocat                 ! single-category option
+
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   ptsu                      ! sea ice / snow surface temperature
+
+      REAL(wp), DIMENSION(:,:)  , INTENT(in)     ::   ptb                       ! sea ice base temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phs                       ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phi                       ! sea ice thickness
+      
+      INTEGER                                    ::   ji, jj, jl                ! dummy loop indices
+      INTEGER                                    ::   iter                      !
+      REAL(wp)                                   ::   zfac, zfac2, zfac3        ! dummy factors
+      REAL(wp)                                   ::   zkeff_h, ztsu             !
+      REAL(wp)                                   ::   zqc, zqnet                !
+      REAL(wp)                                   ::   zhe, zqa0                 !
+           
+      INTEGER , PARAMETER                        ::   nit = 10                  ! number of iterations
+      REAL(wp), PARAMETER                        ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
+
+      REAL(wp), DIMENSION(:,:,:), POINTER        ::   zgfac                    ! enhanced conduction factor
+      
+      !!---------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn')
+      !
+      CALL wrk_alloc( jpi,jpj,jpl,   zgfac )
+
+      ! -------------------------------------!
+      !      I   Enhanced conduction factor  !
+      ! -------------------------------------!
+      !
+      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3)
+      ! Fichefet and Morales Maqueda, JGR 1997
+      !
+      zgfac(:,:,:) = 1._wp
+      
+      SELECT CASE ( k_monocat ) 
+      
+         CASE ( 1 , 3 )
+      
+            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 ) ) )
+                  END DO
+               END DO
+            END DO
+      
+      END SELECT
+      
+      ! -------------------------------------------------------------!
+      !      II   Surface temperature and conduction flux            !
+      ! -------------------------------------------------------------!
+
+      zfac   =   rcdic * rn_cnd_s 
+                                      ! ========================== !
+      DO jl = 1, jpl                  !  Loop over ice categories  !
+         !                            ! ========================== !
+         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)
+               
+               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 )
+               
+               END DO
+               
+               ptsu(ji,jj,jl)    = MIN( rt0, ztsu )
+               
+               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
+               
+            END DO
+         END DO
+         
+      END DO 
+      
+      CALL wrk_dealloc( jpi,jpj,jpl,   zgfac )
+      
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn')
+
+   END SUBROUTINE blk_ice_qcn
    
 #endif
@@ -1093,8 +1243,6 @@
       zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
       !
-!!      DO jj = 2, jpjm1
-!!         DO ji = fs_2, fs_jpim1
-      DO jj = 1, jpj
-         DO ji = 1, jpi
+      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
+         DO ji = fs_2, fs_jpim1
             ! Virtual potential temperature [K]
             zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
@@ -1140,5 +1288,5 @@
          END DO
       END DO
-!!      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. )
+      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. )
       !
    END SUBROUTINE Cdn10_Lupkes2015
