Index: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12854)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12890)
@@ -94,12 +94,11 @@
       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
       REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
-      REAL(wp) ::   zhs_min   =  0.1_wp       ! minimum snow thickness for conductivity calculation 
+      REAL(wp) ::   zhs_ssl   =  0.03_wp      ! surface scattering layer in the snow 
+      REAL(wp) ::   zhi_ssl   =  0.10_wp      ! surface scattering layer in the ice
       REAL(wp) ::   ztmelts                   ! ice melting temperature
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
       REAL(wp) ::   zcpi                      ! Ice specific heat
       REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
-      REAL(wp) ::   zfac                      ! dummy factor
-      !
-      REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only)
+      !
       REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
@@ -150,28 +149,14 @@
       !------------------
       DO ji = 1, npti
-
-         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead
-         !!IF( h_s_1d(ji) < zhs_min ) THEN
-         !!  isnow(ji) = h_s_1d(ji) / zhs_min
-         !!  zh_s(ji) = zhs_min * r1_nlay_s
-         !!ELSE
-         !!  isnow(ji) = 1.0_wp
-         !!  zh_s(ji) = h_s_1d(ji) * r1_nlay_s
-         !!END IF
-         isnow(ji) = za_s_fra(ji) !!clem: 2 variables for the same thing
-         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
-         
-         
-         ! layer thickness
-         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
-
+         zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction
+         zh_i(ji) = MAX( zhi_ssl , h_i_1d(ji) ) * r1_nlay_i ! set a minimum ice thickness for conduction
       END DO
       !
-      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
-      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
+      WHERE( h_i_1d(1:npti) >= zhi_ssl )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
+      ELSEWHERE                            ;   z1_h_i(1:npti) = 0._wp ! put 0 if ice thick < scattering layer
       END WHERE
       !
-      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
-      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
+      WHERE( h_s_1d(1:npti) >= zhs_ssl )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
+      ELSEWHERE                            ;   z1_h_s(1:npti) = 0._wp ! put 0 if snw thick < scattering layer
       END WHERE
       !
@@ -198,5 +183,5 @@
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) )
             !                             ! radiation absorbed by the layer-th snow layer
             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
@@ -204,9 +189,9 @@
       END DO
       !
-      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
+      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * (1._wp - za_s_fra(1:npti))
       DO jk = 1, nlay_i 
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
             !                             ! radiation absorbed by the layer-th ice layer
             zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
@@ -216,5 +201,5 @@
       qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
       !
-      iconv    = 0          ! number of iterations
+      iconv = 0          ! number of iterations
       !
       l_T_converged(:) = .FALSE.
@@ -243,5 +228,5 @@
                DO ji = 1, npti
                   ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
-                     &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+                     &                    MIN( -epsi10, 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 )
                END DO
             END DO
@@ -251,13 +236,13 @@
             DO ji = 1, npti
                ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
-                  &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
+                  &                            - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
                ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
-                  &                           - 0.011_wp * ( t_bo_1d(ji) - rt0 )
+                  &                            - 0.011_wp * ( t_bo_1d(ji) - rt0 )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1, npti
-                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
-                     &                        MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
-                     &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
+                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /       &
+                     &                         MIN( -epsi10, 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 ) &
+                     &                        - 0.011_wp * ( 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 )
                END DO
             END DO
@@ -304,7 +289,7 @@
          DO ji = 1, npti   ! Snow-ice interface
             IF ( .NOT. l_T_converged(ji) ) THEN
-               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
-               IF( zfac > epsi10 ) THEN
-                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
+               IF( h_s_1d(ji) >= zhs_ssl ) THEN
+                  zkappa_s(ji,nlay_s) =   zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / &
+                     &                  ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) )
                ELSE
                   zkappa_s(ji,nlay_s) = 0._wp
@@ -324,5 +309,5 @@
          DO ji = 1, npti   ! Snow-ice interface
             IF ( .NOT. l_T_converged(ji) ) &
-               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
+               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji))
          END DO
          !
@@ -558,5 +543,5 @@
                   IF( t_su_1d(ji) < rt0 ) THEN
                      t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
-                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
+                        &          ( za_s_fra(ji) * t_s_1d(ji,1) + (1._wp - za_s_fra(ji)) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
                   ENDIF
                ENDIF
@@ -790,6 +775,6 @@
          !
          DO ji = 1, npti
-            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
-               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
+            qcn_ice_top_1d(ji) =  -           za_s_fra(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
+               &                  - ( 1._wp - za_s_fra(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
          END DO
          !
@@ -806,7 +791,7 @@
          DO ji = 1, npti
             t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature
-               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &
-               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) &
-               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 )
+               &                +          za_s_fra(ji)  * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &
+               &                + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) &
+               &          ) / MAX( epsi10, za_s_fra(ji)  * zkappa_s(ji,0) * zg1s + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 )
             t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su
          END DO
@@ -866,11 +851,12 @@
       !--------------------------------------------------------------------
       ! effective conductivity and 1st layer temperature (needed by Met Office)
+      ! this is a conductivity at mid-layer, hence the factor 2
       DO ji = 1, npti
-         IF( h_i_1d(ji) > 0.1_wp ) THEN
+         IF( h_i_1d(ji) >= zhi_ssl ) THEN   ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0)
             cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
          ELSE
-            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp ! cnd_ice is capped by: cond_i/0.1m
+            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl
          ENDIF
-         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
+         t1_ice_1d(ji) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1)
       END DO
       !
@@ -886,8 +872,7 @@
       DO ji = 1, npti         
          !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-         IF( h_s_1d(ji) >= zhs_min ) THEN !!clem change that
-            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
+         IF( h_s_1d(ji) >= zhs_ssl ) THEN
+            t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / &
+               &          ( rn_cnd_s * zh_i(ji)                + ztcond_i(ji,1) * zh_s(ji)                )
          ELSE
             t_si_1d(ji) = t_su_1d(ji)
