Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8938)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8939)
@@ -44,9 +44,9 @@
    INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation
    !                                           ! associated indices:
-   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999)
-
-   INTEGER, PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns
-   INTEGER, PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns
-   INTEGER, PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd
+   INTEGER, PARAMETER ::   np_BL99 = 1         ! Bitz and Lipscomb (1999)
+
+   INTEGER, PARAMETER ::   np_zdf_jules_OFF = 0   !   compute all temperatures from qsr and qns
+   INTEGER, PARAMETER ::   np_zdf_jules_SND = 1   !   compute conductive heat flux and surface temperature from qsr and qns
+   INTEGER, PARAMETER ::   np_zdf_jules_RCV = 2   !   compute snow and inner ice temperatures from qcnd
    
    !!----------------------------------------------------------------------
@@ -143,11 +143,11 @@
       REAL(wp) ::   zfac                      ! dummy factor
       !
-      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
-      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
-      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
-      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
-      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
+      REAL(wp), DIMENSION(jpij) ::   isnow        ! switch for presence (1) or absence (0) of snow
+      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
+      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
+      REAL(wp), DIMENSION(jpij) ::   zh_s, z1_h_s ! snow layer thickness
+      REAL(wp), DIMENSION(jpij) ::   zqns_ice_b   ! solar radiation absorbed at the surface
+      REAL(wp), DIMENSION(jpij) ::   zfnet        ! surface flux function
+      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function
       !
       REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
@@ -265,4 +265,5 @@
       iconv    =  0          ! number of iterations
       zdti_max =  1000._wp   ! maximal value of error on all points
+      !
       !                                                          !============================!
       DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
@@ -292,12 +293,12 @@
             !
             DO ji = 1, npti
-               ztcond_i(ji,0)      = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
+               ztcond_i(ji,0)      = rcdic + 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 )
-               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  &
+               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, 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(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /          &
+                  ztcond_i(ji,jk) = rcdic + 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 )
@@ -321,5 +322,5 @@
             DO ji = 1, npti
                zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity
-               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )   ! Effective thickness he (zhe)
+               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )     ! Effective thickness he (zhe)
                IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
                   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he)
@@ -394,6 +395,6 @@
             !----------------------------
 
-            DO ji = 1, npti
             ! update of the non solar flux according to the update in T_su
+            DO ji = 1, npti
                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
             END DO
@@ -483,6 +484,5 @@
                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
                   ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
-                     &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
+                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
                ENDIF
                !                            !---------------------!
@@ -515,6 +515,5 @@
                      ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
                      ztrid(ji,jm_min(ji)+1,3)  = 0.0
-                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  &
-                        &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
+                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
                   ENDIF
 
@@ -561,35 +560,34 @@
 
             DO jk = jm_mint+1, jm_maxt
-            DO ji = 1, npti
-               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
-               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
-               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
-            END DO
-            END DO
-
-            DO ji = 1, npti
+               DO ji = 1, npti
+                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
+                  zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
+                  zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
+               END DO
+            END DO
+
             ! ice temperatures
-            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+            DO ji = 1, npti
+               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
             END DO
 
             DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-            DO ji = 1, npti
-               jk    =  jm - nlay_s - 1
-               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
-            END DO
-            END DO
-
-            DO ji = 1, npti
-            ! snow temperatures      
-            IF( h_s_1d(ji) > 0._wp ) THEN
-               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
-               &                / zdiagbis(ji,nlay_s+1)
-            ENDIF
-            ! surface temperature
-            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) *  &
-                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
-            ENDIF
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+               END DO
+            END DO
+
+            DO ji = 1, npti
+               ! snow temperatures      
+               IF( h_s_1d(ji) > 0._wp ) THEN
+                  t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
+               ENDIF
+               ! surface temperature
+               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) *  &
+                     &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
+               ENDIF
             END DO
             !
@@ -601,21 +599,21 @@
             zdti_max = 0._wp
             DO ji = 1, npti
-            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) ) )     
-            END DO
-
-            DO jk  =  1, nlay_s
-            DO ji = 1, npti
-               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
-            END DO
-
-            DO jk  =  1, nlay_i
-            DO ji = 1, npti
-               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
-               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
-               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
-            END DO
+               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) ) )     
+            END DO
+
+            DO jk = 1, nlay_s
+               DO ji = 1, npti
+                  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
+            END DO
+
+            DO jk = 1, nlay_i
+               DO ji = 1, npti
+                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
+                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+               END DO
             END DO
 
@@ -652,82 +650,78 @@
 
             DO jm = nlay_s + 2, nlay_s + nlay_i 
-            DO ji = 1, npti
-               jk = jm - nlay_s - 1
-               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
-               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
-               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-            END DO
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
+                  ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+                  ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
+                  zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+               END DO
             ENDDO
 
             jm =  nlay_s + nlay_i + 1
             DO ji = 1, npti
-            !!ice bottom term
-            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
-            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
-            ztrid(ji,jm,3)  = 0.0
-            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
+               !!ice bottom term
+               ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
+               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
+               ztrid(ji,jm,3)  = 0.0
+               zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+                  &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
             ENDDO
 
 
             DO ji = 1, npti
-            !                              !---------------------!
-            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
-               !                           !---------------------!
-               ! snow interior terms (bottom equation has the same form as the others)
-               DO jm = 3, nlay_s + 1
-               jk = jm - 1
-               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
-               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
-               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
-               END DO
-
-               ! case of only one layer in the ice (ice equation is altered)
-               IF ( nlay_i == 1 ) THEN
-               ztrid(ji,nlay_s+2,3)  = 0.0
-               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+               !                              !---------------------!
+               IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
+                  !                           !---------------------!
+                  ! snow interior terms (bottom equation has the same form as the others)
+                  DO jm = 3, nlay_s + 1
+                     jk = jm - 1
+                     ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
+                     ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+                     ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
+                     zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+                  END DO
+                  
+                  ! case of only one layer in the ice (ice equation is altered)
+                  IF ( nlay_i == 1 ) THEN
+                     ztrid(ji,nlay_s+2,3)  = 0.0
+                     zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+                  ENDIF
+                  
+                  jm_min(ji) = 2
+                  jm_max(ji) = nlay_i + nlay_s + 1
+                  
+                  ! first layer of snow equation
+                  ztrid(ji,2,1)  = 0.0
+                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)
+                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
+                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
+                  
+                  !                            !---------------------!
+               ELSE                            ! cells without snow  !
+                  !                            !---------------------!
+                  jm_min(ji)    =  nlay_s + 2
+                  jm_max(ji)    =  nlay_i + nlay_s + 1
+                  
+                  ! first layer of ice equation
+                  ztrid(ji,jm_min(ji),1)  = 0.0
+                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)  
+                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
+                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
+                  
+                  ! case of only one layer in the ice (surface & ice equations are altered)
+                  IF ( nlay_i == 1 ) THEN
+                     ztrid(ji,jm_min(ji),1)  = 0.0
+                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)
+                     ztrid(ji,jm_min(ji),3)  = 0.0
+                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)  &
+                        &                    + qcn_ice_1d(ji) )
+                  ENDIF
+                  
                ENDIF
-
-               jm_min(ji) = 2
-               jm_max(ji) = nlay_i + nlay_s + 1
-
-               ! first layer of snow equation
-               ztrid(ji,2,1)  = 0.0
-               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)
-               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
-               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
-               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
-
-            !                               !---------------------!
-            ELSE                            ! cells without snow  !
-            !                               !---------------------!
-
-               jm_min(ji)    =  nlay_s + 2
-               jm_max(ji)    =  nlay_i + nlay_s + 1
-
-               ! first layer of ice equation
-               ztrid(ji,jm_min(ji),1)  = 0.0
-               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)  
-               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
-               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
-               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) )
-
-               ! case of only one layer in the ice (surface & ice equations are altered)
-               IF ( nlay_i == 1 ) THEN
-               ztrid(ji,jm_min(ji),1)  = 0.0
-               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)
-               ztrid(ji,jm_min(ji),3)  = 0.0
-               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    &
-               &                    + qcn_ice_1d(ji) )
-
-               ENDIF
-
-            ENDIF
-            !
-            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
-            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
-            !
+               !
+               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
+               zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
+               !
             END DO
             !
@@ -740,34 +734,33 @@
             jm_mint = nlay_i+5
             DO ji = 1, npti
-            jm_mint = MIN(jm_min(ji),jm_mint)
-            jm_maxt = MAX(jm_max(ji),jm_maxt)
-            END DO
-
+               jm_mint = MIN(jm_min(ji),jm_mint)
+               jm_maxt = MAX(jm_max(ji),jm_maxt)
+            END DO
+            
             DO jk = jm_mint+1, jm_maxt
-            DO ji = 1, npti
-               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
-               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
-               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
-            END DO
-            END DO
-
-            DO ji = 1, npti
+               DO ji = 1, npti
+                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
+                  zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
+                  zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
+               END DO
+            END DO
+            
             ! ice temperatures
-            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+           DO ji = 1, npti
+               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
             END DO
 
             DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-            DO ji = 1, npti
-               jk    =  jm - nlay_s - 1
-               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
-            END DO
-            END DO
-
-            DO ji = 1, npti
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+               END DO
+            END DO
+            
             ! snow temperatures      
-            IF( h_s_1d(ji) > 0._wp ) THEN
-               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
-               &                / zdiagbis(ji,nlay_s+1)
-            ENDIF
+            DO ji = 1, npti
+               IF( h_s_1d(ji) > 0._wp ) THEN
+                  t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
+               ENDIF
             END DO
             !
@@ -779,17 +772,17 @@
             zdti_max = 0._wp
 
-            DO jk  =  1, nlay_s
-            DO ji = 1, npti
-               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
-            END DO
-
-            DO jk  =  1, nlay_i
-            DO ji = 1, npti
-               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
-               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
-               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
-            END DO
+            DO jk = 1, nlay_s
+               DO ji = 1, npti
+                  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
+            END DO
+            
+            DO jk = 1, nlay_i
+               DO ji = 1, npti
+                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
+                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+               END DO
             END DO
 
@@ -815,8 +808,8 @@
       !
       DO ji = 1, npti
-      !                                ! surface ice conduction flux
+         !                                ! surface ice conduction flux
          fc_su(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))
-      !                                ! bottom ice conduction flux
+         !                                ! bottom ice conduction flux
          fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
       END DO
@@ -892,25 +885,23 @@
       ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
       !---------------------------------------------------------------------------------------
-      ! effective conductivity and 1st layer temperature (Jules coupling)
+      ! effective conductivity and 1st layer temperature (needed by Met Office)
       DO ji = 1, npti
-         IF (h_s_1d(ji) > 0.1 ) THEN 
-             cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
+         IF( h_s_1d(ji) > 0.1_wp ) THEN 
+            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
          ELSE
-             IF (h_i_1d(ji) > 0.1 ) THEN
-                 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
-             ELSE
-                 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1
-             ENDIF
+            IF( h_i_1d(ji) > 0.1_wp ) THEN
+               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
+            ELSE
+               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp
+            ENDIF
          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) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
       END DO
       !
       IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode
-         !
          ! Restore temperatures to their initial values
          t_s_1d    (1:npti,:) = ztsold(1:npti,:)
          t_i_1d    (1:npti,:) = ztiold(1:npti,:)
          qcn_ice_1d(1:npti)   = fc_su (1:npti)
-         !  
       ENDIF
       !
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8938)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8939)
@@ -32,9 +32,9 @@
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      ! 
-   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
+   USE oce     , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
    USE ocealb         ! 
    USE eosbn2         ! 
-   USE sbcrnf, ONLY : l_rnfcpl
-   USE sbcisf   , ONLY : l_isfcpl
+   USE sbcrnf  , ONLY : l_rnfcpl
+   USE sbcisf  , ONLY : l_isfcpl
 #if defined key_cice
    USE ice_domain_size, only: ncat
@@ -1243,17 +1243,8 @@
       ! needed by Met Office
       IF( srcv(jpr_ts_ice)%laction ) THEN 
-        DO jn = 1, jpl 
-          DO jj = 1, jpj 
-            DO ji = 1, jpi 
-              IF (frcv(jpr_ts_ice)%z3(ji,jj,jn) > 0.0) THEN 
-                tsfc_ice(ji,jj,jn) = 0.0 
-              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jn) < -60.0) THEN 
-                tsfc_ice(ji,jj,jn) = -60.0 
-              ELSE 
-                tsfc_ice(ji,jj,jn) = frcv(jpr_ts_ice)%z3(ji,jj,jn) 
-              ENDIF 
-            END DO 
-          END DO 
-        END DO 
+         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   tsfc_ice(:,:,:) = 0.0 
+         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   tsfc_ice(:,:,:) = -60.
+         ELSEWHERE                                        ;   tsfc_ice(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:)
+         END WHERE
       ENDIF 
 #endif
@@ -2066,7 +2057,7 @@
 #if defined key_lim3      
       IF( ln_meto_cpl ) THEN
-         !                                                      ! ========================= !
-         SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
-         !                                                      ! ========================= !
+         !                                                   ! ========================= !
+         SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )          !  ice topmelt and botmelt  !
+         !                                                   ! ========================= !
          CASE ('coupled')
             qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i(:,:,:)
@@ -2186,8 +2177,10 @@
          IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
       ENDIF
-
+      !
+      !                                                      ! ------------------------- !
+      !                                                      ! 1st layer ice/snow temp.  !
+      !                                                      ! ------------------------- !
 #if defined key_lim3
-      !!!!! Getting NEMO4-LIM working at Met Office
-      ! Top layer ice temperature
+      ! needed by  Met Office
       IF( ssnd(jps_ttilyr)%laction) THEN
          SELECT CASE( sn_snd_ttilyr%cldes)
@@ -2198,7 +2191,5 @@
          IF( ssnd(jps_ttilyr)%laction )   CALL cpl_snd( jps_ttilyr, isec, ztmp3, info )
       ENDIF
-      !!!!!
 #endif
-
       !                                                      ! ------------------------- !
       !                                                      !           Albedo          !
@@ -2315,6 +2306,8 @@
 
 #if defined key_lim3
-      ! NEMO4 - Jules coupling - Met Office
-      ! Send meltpond fields  
+      !                                                      ! ------------------------- !
+      !                                                      !      Ice melt ponds       ! 
+      !                                                      ! ------------------------- !
+      ! needed by Met Office
       IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
          SELECT CASE( sn_snd_mpnd%cldes)  
@@ -2331,17 +2324,22 @@
                  ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl) 
                ENDDO  
-            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' )  
+            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )  
             END SELECT  
-         CASE( 'default' )    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%cldes' )     
+         CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%cldes' )     
          END SELECT  
-         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )     
+         IF( ssnd(jps_a_p)%laction  )   CALL cpl_snd( jps_a_p , isec, ztmp3, info )     
          IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )     
-         ! 
-         ! Send ice effective conductivity 
+      ENDIF 
+      ! 
+      !                                                      ! ------------------------- !
+      !                                                      !     Ice conductivity      ! 
+      !                                                      ! ------------------------- !
+      ! needed by Met Office
+      IF( ssnd(jps_kice)%laction ) THEN 
          SELECT CASE( sn_snd_cond%cldes) 
          CASE( 'weighted ice' )    
             SELECT CASE( sn_snd_cond%clcat ) 
             CASE( 'yes' )    
-	               ztmp3(:,:,1:jpl) =  cnd_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
+	       ztmp3(:,:,1:jpl) =  cnd_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
             CASE( 'no' ) 
                ztmp3(:,:,:) = 0.0 
@@ -2349,13 +2347,12 @@
                  ztmp3(:,:,1) = ztmp3(:,:,1) + cnd_ice(:,:,jl) * a_i(:,:,jl) 
                ENDDO 
-            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' ) 
+            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' ) 
             END SELECT 
          CASE( 'ice only' )    
            ztmp3(:,:,1:jpl) = cnd_ice(:,:,1:jpl) 
+         CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%cldes' )     
          END SELECT 
          IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info ) 
       ENDIF 
-      !   
-      !!!!!
 #endif
 
@@ -2682,5 +2679,8 @@
 
 #if defined key_lim3
-      ! NEMO4 - Jules coupling - Met Office
+      !                                                      ! ------------------------- !
+      !                                                      ! Sea surface freezing temp ! 
+      !                                                      ! ------------------------- !
+      ! needed by Met Office
       CALL eos_fzp(tsn(:,:,1,jp_sal), sstfrz)
       ztmp1(:,:) = sstfrz(:,:) + rt0
