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 12919)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 13148)
@@ -95,4 +95,5 @@
       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) ::   zh_min    =  1.e-3_wp     ! minimum ice/snow thickness for conduction
       REAL(wp) ::   ztmelts                   ! ice melting temperature
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
@@ -123,4 +124,5 @@
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), DIMENSION(jpij)            ::   zkappa_comb ! Combined snow and ice surface conductivity
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
@@ -130,4 +132,5 @@
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
       REAL(wp), DIMENSION(jpij)            ::   za_s_fra    ! ice fraction covered by snow 
+      REAL(wp), DIMENSION(jpij)            ::   isnow       ! snow presence (1) or not (0) 
       !
       ! Mono-category
@@ -143,5 +146,6 @@
       END DO
 
-      CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) )   ! calculate ice fraction covered by snow
+      ! calculate ice fraction covered by snow for radiation
+      CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) )
       
       !------------------
@@ -158,21 +162,30 @@
       ENDIF
       !
-      ! layers thicknesses
+      ! thicknesses
       DO ji = 1, npti
-         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
+         ! ice thickness
+         IF( h_i_1d(ji) > 0._wp ) THEN 
+            zh_i  (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction
+            z1_h_i(ji) = 1._wp / zh_i(ji)                       !       it must be very small
+         ELSE
+            zh_i  (ji) = 0._wp
+            z1_h_i(ji) = 0._wp
+         ENDIF
+         ! snow thickness
+         IF( h_s_1d(ji) > 0._wp ) THEN
+            zh_s  (ji) = MAX( zh_min , h_s_1d(ji) ) * r1_nlay_s ! set a minimum thickness for conduction
+            z1_h_s(ji) = 1._wp / zh_s(ji)                       !       it must be very small
+            isnow (ji) = 1._wp
+         ELSE
+            zh_s  (ji) = 0._wp
+            z1_h_s(ji) = 0._wp
+            isnow (ji) = 0._wp
+         ENDIF
       END DO
-      !
-      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( 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
+      ! clem: we should apply correction on snow thickness to take into account snow fraction
+      !       it must be a distribution, so it is a bit complicated
       !
       ! Store initial temperatures and non solar heat fluxes
       IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
-         !
          ztsub      (1:npti) = t_su_1d(1:npti)                          ! surface temperature at iteration n-1
          ztsuold    (1:npti) = t_su_1d(1:npti)                          ! surface temperature initial value
@@ -200,14 +213,12 @@
       END DO
       !
-      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))
+      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) =           za_s_fra(ji)   * zradtr_s(ji,nlay_s)                       &   ! part covered by snow
-               &                                       * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) &
+               &                                       * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min  ) ) &
                &            + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji)                        &   ! part snow free
-               &                                       * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
-            
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
+               &                                       * 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)
@@ -304,12 +315,7 @@
          END DO
          DO ji = 1, npti   ! Snow-ice interface
-            IF ( .NOT. l_T_converged(ji) ) THEN
-               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
-               ENDIF
-            ENDIF
+            IF ( .NOT. l_T_converged(ji) ) &
+               zkappa_s(ji,nlay_s) = isnow(ji) * 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) ) )
          END DO
 
@@ -324,6 +330,10 @@
          END DO
          DO ji = 1, npti   ! Snow-ice interface
-            IF ( .NOT. l_T_converged(ji) ) &
-               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji))
+            IF ( .NOT. l_T_converged(ji) ) THEN
+               ! Calculate combined surface snow and ice conductivity to pass through the coupler
+               zkappa_comb(ji) = isnow(ji) * zkappa_s(ji,0) + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)
+               ! If there is snow then use the same snow-ice interface conductivity for the top layer of ice
+               IF( h_s_1d(ji) > 0._wp )   zkappa_i(ji,0) = zkappa_s(ji,nlay_s)
+           ENDIF
          END DO
          !
@@ -334,5 +344,5 @@
             DO ji = 1, npti
                zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
+               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / zcpi
             END DO
          END DO
@@ -558,9 +568,10 @@
                   ztsub(ji) = t_su_1d(ji)
                   IF( t_su_1d(ji) < rt0 ) THEN
-                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
-                        &          ( 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))
+                     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))
                   ENDIF
                ENDIF
             END DO
+            !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1)
             !
             !--------------------------------------------------------------
@@ -575,9 +586,14 @@
 
                IF ( .NOT. l_T_converged(ji) ) THEN
+
                   t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
                   zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
 
-                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
-                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
+                  IF( h_s_1d(ji) > 0._wp ) THEN
+                     DO jk = 1, nlay_s
+                        t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+                        zdti_max      = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+                     END DO
+                  ENDIF
 
                   DO jk = 1, nlay_i
@@ -586,4 +602,5 @@
                      zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
                   END DO
+                  
                   ! convergence test
                   IF( ln_zdf_chkcvg ) THEN
@@ -745,4 +762,5 @@
                ENDIF
             END DO
+            !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1)
             !
             !--------------------------------------------------------------
@@ -757,8 +775,12 @@
 
                IF ( .NOT. l_T_converged(ji) ) THEN
-                  ! t_s
-                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp )
-                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
-                  ! t_i
+
+                  IF( h_s_1d(ji) > 0._wp ) THEN
+                     DO jk = 1, nlay_s
+                        t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+                        zdti_max      = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+                     END DO
+                  ENDIF
+
                   DO jk = 1, nlay_i
                      ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
@@ -766,4 +788,5 @@
                      zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
                   END DO
+
                   ! convergence test
                   IF( ln_zdf_chkcvg ) THEN
@@ -789,5 +812,5 @@
       !     bottom ice conduction flux
       DO ji = 1, npti
-         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
+         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
       END DO
       !     surface ice conduction flux
@@ -795,6 +818,6 @@
          !
          DO ji = 1, npti
-            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) )
+            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) )
          END DO
          !
@@ -810,8 +833,7 @@
          !
          DO ji = 1, npti
-            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature
-               &                +          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) = ( qcn_ice_top_1d(ji) +          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 )
             t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su
          END DO
@@ -873,10 +895,11 @@
       ! this is a conductivity at mid-layer, hence the factor 2
       DO ji = 1, npti
-         IF( h_i_1d(ji) >= zhi_ssl ) THEN   ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0)
+         IF( h_i_1d(ji) >= zhi_ssl ) THEN
+!!clem metO            cnd_ice_1d(ji) = 2._wp * zkappa_comb(ji)
             cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
          ELSE
             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) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1)
+         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
       END DO
       !
@@ -893,6 +916,8 @@
          !--- Snow-ice interfacial temperature (diagnostic SIMIP)
          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)                )
+            t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,1)   &
+               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) &
+               &          / ( rn_cnd_s       * h_i_1d(ji) * r1_nlay_i &
+               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s )
          ELSE
             t_si_1d(ji) = t_su_1d(ji)
