Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icecor.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icecor.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icecor.F90	(revision 13228)
@@ -114,5 +114,5 @@
             ENDIF
          END_2D
-         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1., v_ice, 'V', -1. )
+         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
       ENDIF
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn.F90	(revision 13228)
@@ -129,6 +129,6 @@
             zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
             zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
-            u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
-            v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
+            u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1.0_wp, zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
+            v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1.0_wp, zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
          END_2D
          ! ---
@@ -159,5 +159,5 @@
                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
             END_2D
-            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
+            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.0_wp )
             ! output
             CALL iom_put( 'icediv' , zdivu_i )
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_adv_pra.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_adv_pra.F90	(revision 13228)
@@ -117,5 +117,5 @@
          END_2D
       END DO
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp )
       !
       ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
@@ -254,5 +254,5 @@
                &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
          END_2D
-         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. )
+         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp )
          !
          ! --- Ensure non-negative fields --- !
@@ -425,7 +425,7 @@
 
       !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
-         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
+         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
+         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
       !
    END SUBROUTINE adv_x
@@ -584,7 +584,7 @@
 
       !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
-         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
+         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
+         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
       !
    END SUBROUTINE adv_y
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_adv_umx.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_adv_umx.F90	(revision 13228)
@@ -122,5 +122,5 @@
          END_2D
       END DO
-      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
+      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp )
       !
       !
@@ -336,5 +336,5 @@
                &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
          END_2D
-         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. )
+         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1.0_wp )
          !
          !
@@ -469,5 +469,5 @@
             END_2D
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. )
+         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1.0_wp )
          !
          IF    ( np_limiter == 1 ) THEN
@@ -500,5 +500,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1.0_wp )
       !
    END SUBROUTINE adv_umx
@@ -552,5 +552,5 @@
                END_2D
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
             !
             DO jl = 1, jpl              !-- flux in y-direction
@@ -576,5 +576,5 @@
                END_2D
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
             !
             DO jl = 1, jpl              !-- flux in x-direction
@@ -598,5 +598,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1.0_wp )
 
    END SUBROUTINE upstream
@@ -660,5 +660,5 @@
                END_2D
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
 
             DO jl = 1, jpl              !-- flux in y-direction
@@ -686,5 +686,5 @@
                END_2D
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
             !
             DO jl = 1, jpl              !-- flux in x-direction
@@ -744,5 +744,5 @@
             END_2D
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
          !
          !                                                        !--  ultimate interpolation of pt at v-point  --!
@@ -771,5 +771,5 @@
             END_2D
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
          !
          !                                                        !--  ultimate interpolation of pt at u-point  --!
@@ -824,5 +824,5 @@
          END DO
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1.0_wp )
       !
       !                                                     !--  BiLaplacian in i-direction  --!
@@ -838,5 +838,5 @@
          END DO
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1.0_wp )
       !
       !
@@ -964,5 +964,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1.0_wp )
       !
       !                                                     !--  BiLaplacian in j-direction  --!
@@ -975,5 +975,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1.0_wp )
       !
       !
@@ -1114,5 +1114,5 @@
             END_2D
          END DO
-         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1.0_wp, ztj_ups, 'T', 1.0_wp )
 
          DO jl = 1, jpl
@@ -1136,5 +1136,5 @@
             END_2D
          END DO
-         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp, pfv_ho, 'V', -1.0_wp )   ! lateral boundary cond.
 
       ENDIF
@@ -1193,5 +1193,5 @@
          END_2D
       END DO
-      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
+      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1.0_wp, zbetdo, 'T', 1.0_wp )   ! lateral boundary cond. (unchanged sign)
 
       
@@ -1248,5 +1248,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
+      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.0_wp)   ! lateral boundary cond.
       
       DO jl = 1, jpl
@@ -1312,5 +1312,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond.
+      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp)   ! lateral boundary cond.
       !
    END SUBROUTINE limiter_x
@@ -1339,5 +1339,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
+      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.0_wp)   ! lateral boundary cond.
 
       DO jl = 1, jpl
@@ -1404,5 +1404,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
+      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.0_wp)   ! lateral boundary cond.
       !
    END SUBROUTINE limiter_y
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_rdgrft.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_rdgrft.F90	(revision 13228)
@@ -300,5 +300,13 @@
 
       !                       ! Ice thickness needed for rafting
+      ! In single precision there were floating point invalids due a sqrt of zhi which happens to have negative values
+      ! To solve that an extra check about the value of pv_i was added.
+      ! Although adding this condition is safe, the double definition (one for single other for double) has been kept to preserve the results of the sette test.
+#if defined key_single
+
+      WHERE( pa_i(1:npti,:) > epsi10 .and. pv_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+#else
       WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+#endif
       ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
@@ -780,5 +788,5 @@
             strength(ji,jj) = zworka(ji,jj)
          END_2D
-         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
          !
       CASE( 2 )               !--- Temporal smoothing
@@ -799,5 +807,5 @@
             ENDIF
          END_2D
-         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
          !
       END SELECT
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_rhg_evp.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icedyn_rhg_evp.F90	(revision 13228)
@@ -300,5 +300,5 @@
 
       END_2D
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
       !
       !                                  !== Landfast ice parameterization ==!
@@ -319,5 +319,5 @@
             tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
          END_2D
-         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp )
          !
       ELSE                               !-- no landfast
@@ -353,5 +353,5 @@
 
          END_2D
-         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp )
 
          DO_2D_01_01
@@ -395,5 +395,5 @@
           
          END_2D
-         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp )
 
          DO_2D_10_10
@@ -484,5 +484,5 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
+            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
             !
 #if defined key_agrif
@@ -533,5 +533,5 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
+            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
             !
 #if defined key_agrif
@@ -584,5 +584,5 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
+            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
             !
 #if defined key_agrif
@@ -633,5 +633,5 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
+            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
             !
 #if defined key_agrif
@@ -694,8 +694,8 @@
 
       END_2D
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp )
       
       ! --- Store the stress tensor for the next time step --- !
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
       pstress1_i (:,:) = zs1 (:,:)
       pstress2_i (:,:) = zs2 (:,:)
@@ -714,6 +714,6 @@
          & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
          !
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
-            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, &
+            &                                  ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
          !
          CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
@@ -752,5 +752,5 @@
             zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
          END_2D
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp )
          !
          CALL iom_put( 'isig1' , zsig1 )
@@ -769,6 +769,6 @@
          & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
          !
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., &
-            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
+            &                                  zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
 
          CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
@@ -802,7 +802,7 @@
          END_2D
 
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
-            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
-            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
+            &                                  zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
+            &                                  zdiag_xatrp    , 'U', -1.0_wp, zdiag_yatrp    , 'V', -1.0_wp )
 
          CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceitd.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceitd.F90	(revision 13228)
@@ -148,6 +148,11 @@
                !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
                !          in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+# if defined key_single
+               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi06 ) )   nptidx(ji) = 0
+               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi06 ) )   nptidx(ji) = 0
+# else
                IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0
                IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0
+# endif
                !
                ! 2) Hn-1 < Hn* < Hn+1  
@@ -170,6 +175,11 @@
             !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
             !    in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+# if defined key_single
+            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi06 ) )   nptidx(ji) = 0
+            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi06 ) )   nptidx(ji) = 0
+# else
             IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0
             IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0
+# endif
          END DO
          !
@@ -538,5 +548,9 @@
       ! 4) Update ice thickness and temperature
       !-------------------------------------------------------------------------------
+# if defined key_single
+      WHERE( a_i_2d(1:npti,:) >= epsi06 )
+# else
       WHERE( a_i_2d(1:npti,:) >= epsi20 )
+# endif
          h_i_2d (1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
          t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icesbc.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icesbc.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icesbc.F90	(revision 13228)
@@ -86,5 +86,5 @@
             vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
          END_2D
-         CALL lbc_lnk_multi( 'icesbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+         CALL lbc_lnk_multi( 'icesbc', utau_ice, 'U', -1.0_wp, vtau_ice, 'V', -1.0_wp )
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd.F90	(revision 13228)
@@ -121,5 +121,5 @@
          END_2D
       ENDIF
-      CALL lbc_lnk( 'icethd', zfric, 'T',  1. )
+      CALL lbc_lnk( 'icethd', zfric, 'T',  1.0_wp )
       !
       !--------------------------------------------------------------------!
@@ -218,5 +218,5 @@
                               CALL ice_thd_dh                           ! Ice-Snow thickness   
                               CALL ice_thd_pnd                          ! Melt ponds formation
-                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
+                              CALL ice_thd_ent( e_i_1d(1:npti,:), .true. )      ! Ice enthalpy remapping
             ENDIF
                               CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !    
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_dh.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_dh.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_dh.F90	(revision 13228)
@@ -186,5 +186,5 @@
       ! Snow precipitation
       !-------------------
-      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
+      CALL ice_thd_snwblow( 1.0_wp - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
 
       zdeltah(1:npti,:) = 0._wp
@@ -442,5 +442,5 @@
                
                zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
-                  &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
+                  &            - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp  * ztmelts
 
                zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_do.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_do.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_do.F90	(revision 13228)
@@ -191,5 +191,5 @@
          END_2D
          ! 
-         CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1., ht_i_new, 'T', 1.  )
+         CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  )
 
       ENDIF
@@ -385,5 +385,5 @@
             END DO
             ! --- Ice enthalpy remapping --- !
-            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 
+            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl), .false. ) 
          END DO
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_ent.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_ent.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icethd_ent.F90	(revision 13228)
@@ -38,5 +38,5 @@
 CONTAINS
  
-   SUBROUTINE ice_thd_ent( qnew )
+   SUBROUTINE ice_thd_ent( qnew, compute_hfx_err )
       !!-------------------------------------------------------------------
       !!               ***   ROUTINE ice_thd_ent  ***
@@ -64,4 +64,6 @@
       !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   qnew             ! new enthlapies (J.m-3, remapped)
+      LOGICAL, INTENT(in)                     ::   compute_hfx_err  ! determines whether to compute diag.
+                                                                    ! error or not
       !
       INTEGER  :: ji         !  dummy loop indices
@@ -128,9 +130,11 @@
       ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do), 
       ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
-      DO ji = 1, npti
-         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice *  &
-            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 
-      END DO
-      
+      IF( compute_hfx_err ) THEN
+         DO ji = 1, npti
+            hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice *  &
+               &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) )
+         END DO
+      END IF
+ 
    END SUBROUTINE ice_thd_ent
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceupdate.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceupdate.F90	(revision 13228)
@@ -342,5 +342,5 @@
             tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
          END_2D
-         CALL lbc_lnk_multi( 'iceupdate', taum, 'T', 1., tmod_io, 'T', 1. )
+         CALL lbc_lnk_multi( 'iceupdate', taum, 'T', 1.0_wp, tmod_io, 'T', 1.0_wp )
          !
          utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
@@ -364,5 +364,5 @@
          vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
       END_2D
-      CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
+      CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition
       !
       IF( ln_timing )   CALL timing_stop('ice_update_tau')
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icevar.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icevar.F90	(revision 13228)
@@ -635,14 +635,15 @@
       !!-------------------------------------------------------------------
       !
-      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
-      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
-      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
-      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
-      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
-      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
-      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
+
+      WHERE( pa_i (1:npti,:)   < 0._wp )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
+      WHERE( pv_i (1:npti,:)   < 0._wp )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
+      WHERE( pv_s (1:npti,:)   < 0._wp )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
+      WHERE( psv_i(1:npti,:)   < 0._wp )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
+      WHERE( poa_i(1:npti,:)   < 0._wp )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
+      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
+      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
       IF( ln_pnd_H12 ) THEN
-         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
-         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
+         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icewri.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icewri.F90	(revision 13219)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/icewri.F90	(revision 13228)
@@ -135,5 +135,5 @@
             z2d(ji,jj) = 0.5_wp * SQRT( z2da * z2da + z2db * z2db )
          END_2D
-         CALL lbc_lnk( 'icewri', z2d, 'T', 1. )
+         CALL lbc_lnk( 'icewri', z2d, 'T', 1.0_wp )
          CALL iom_put( 'icevel', z2d )
 
