Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icedia.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icedia.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icedia.F90	(revision 14986)
@@ -67,4 +67,5 @@
       REAL(wp)   ::   zbg_ivol, zbg_item, zbg_area, zbg_isal
       REAL(wp)   ::   zbg_svol, zbg_stem
+      REAL(wp)   ::   zbg_ipvol, zbg_ilvol
       REAL(wp)   ::   z_frc_voltop, z_frc_temtop, z_frc_sal
       REAL(wp)   ::   z_frc_volbot, z_frc_tembot
@@ -87,5 +88,6 @@
       ! ----------------------- !
       IF(  iom_use('ibgvol_tot' ) .OR. iom_use('sbgvol_tot' ) .OR. iom_use('ibgarea_tot') .OR. &
-         & iom_use('ibgsalt_tot') .OR. iom_use('ibgheat_tot') .OR. iom_use('sbgheat_tot') ) THEN
+         & iom_use('ibgsalt_tot') .OR. iom_use('ibgheat_tot') .OR. iom_use('sbgheat_tot') .OR. &
+         & iom_use('ipbgvol_tot' ) .OR. iom_use('ilbgvol_tot' ) ) THEN
 
          zbg_ivol = glob_sum( 'icedia', vt_i(:,:) * e1e2t(:,:) ) * 1.e-9  ! ice volume (km3)
@@ -95,4 +97,7 @@
          zbg_item = glob_sum( 'icedia', et_i(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)
          zbg_stem = glob_sum( 'icedia', et_s(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat content (1.e20 J)
+         ! ponds
+         zbg_ipvol = glob_sum( 'icedia', vt_ip(:,:) * e1e2t(:,:) ) * 1.e-9  ! ice pond volume (km3)
+         zbg_ilvol = glob_sum( 'icedia', vt_il(:,:) * e1e2t(:,:) ) * 1.e-9  ! ice pond lid volume (km3)
 
          CALL iom_put( 'ibgvol_tot'  , zbg_ivol )
@@ -102,4 +107,7 @@
          CALL iom_put( 'ibgheat_tot' , zbg_item )
          CALL iom_put( 'sbgheat_tot' , zbg_stem )
+         ! ponds
+         CALL iom_put( 'ipbgvol_tot' , zbg_ipvol )
+         CALL iom_put( 'ilbgvol_tot' , zbg_ilvol )
 
       ENDIF
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icedyn_rhg_vp.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icedyn_rhg_vp.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icedyn_rhg_vp.F90	(revision 14986)
@@ -676,5 +676,5 @@
          CALL lbc_lnk( 'icedyn_rhg_vp', zrhsu, 'U', -1._wp, zrhsv, 'V',  -1._wp)
          CALL lbc_lnk( 'icedyn_rhg_vp', zmU_t, 'U', -1._wp, zmV_t, 'V',  -1._wp)
-         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.)
+         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1._wp, ztauy_oi_rhsv, 'V',  -1._wp)
 
          CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG
@@ -1110,5 +1110,5 @@
       IF ( lwp ) WRITE(numout,*) ' We are out of outer loop '
 
-      CALL lbc_lnk( 'icedyn_rhg_vp', zFU  , 'U',  1., zFV  , 'V',  1. )
+      CALL lbc_lnk( 'icedyn_rhg_vp', zFU  , 'U',  1.0_wp, zFV  , 'V',  1._wp )
       CALL lbc_lnk( 'icedyn_rhg_vp', zBU_prime  , 'U',  1._wp, zBV_prime  , 'V',  1._wp )
       CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime  , 'U',  1._wp, zFV_prime  , 'V',  1._wp )
@@ -1307,5 +1307,5 @@
          
          !
-         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1._wp, ztauy_oi, 'V', -1._wp, ztaux_ai, 'U', -1._wp, ztauy_ai, 'V', -1. ) !, &
+         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1._wp, ztauy_oi, 'V', -1._wp, ztaux_ai, 'U', -1._wp, ztauy_ai, 'V', -1._wp ) !, &
 !            &                           ztaux_bi, 'U', -1.)_wp, ztauy_bi, 'V', -1._wp )
          !
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icestp.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icestp.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icestp.F90	(revision 14986)
@@ -461,4 +461,5 @@
             qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T)
             qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs
+            qml_ice    (ji,jj,jl) = 0._wp   ! surface melt heat flux
             ! Melt pond surface melt diagnostics (mv - more efficient: grouped into one water volume flux)
             dh_i_sum_2d(ji,jj,jl) = 0._wp
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd.F90	(revision 14986)
@@ -536,4 +536,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), qml_ice_1d    (1:npti), qml_ice    (:,:,kl) )
          ! extensive variables
          CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90	(revision 14986)
@@ -224,8 +224,6 @@
       zevap_rema(1:npti) = 0._wp
       DO ji = 1, npti
-         IF( evap_ice_1d(ji) > 0._wp ) THEN
-            zdeltah   (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) )   ! amount of snw that sublimates, < 0
-            zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on)
-         ENDIF
+         zdeltah   (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) )   ! amount of snw that sublimates, < 0
+         zevap_rema(ji) = evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos               ! remaining evap in kg.m-2 (used for ice sublimation later on)
       END DO
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd_ent.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd_ent.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/icethd_ent.F90	(revision 14986)
@@ -121,5 +121,5 @@
          DO ji = 1, npti
             rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
-            qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
+            qnew(ji,jk1) = rswitch * MAX( 0._wp, zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) ! max for roundoff error
          END DO
       END DO
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/iceupdate.F90	(revision 14850)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/iceupdate.F90	(revision 14986)
@@ -289,7 +289,7 @@
       IF( iom_use('hfxcndbot'  ) )   CALL iom_put( 'hfxcndbot'  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux
       IF( iom_use('hfxcndtop'  ) )   CALL iom_put( 'hfxcndtop'  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux
-!!    IF( iom_use('hfxmelt'    ) )   CALL iom_put( 'hfxmelt'    , SUM( qml_ice     * a_i_b, dim=3 ) )   ! Surface melt flux
-!!    IF( iom_use('hfxldmelt'  ) )   CALL iom_put( 'hfxldmelt'  ,      fhld        * at_i_b         )   ! Heat in lead for ice melting 
-!!    IF( iom_use('hfxldgrow'  ) )   CALL iom_put( 'hfxldgrow'  ,      qlead       * r1_Dt_ice      )   ! Heat in lead for ice growth
+      IF( iom_use('hfxmelt'    ) )   CALL iom_put( 'hfxmelt'    , SUM( qml_ice     * a_i_b, dim=3 ) )   ! Surface melt flux
+      IF( iom_use('hfxldmelt'  ) )   CALL iom_put( 'hfxldmelt'  ,      fhld        * at_i_b         )   ! Heat in lead for ice melting 
+      IF( iom_use('hfxldgrow'  ) )   CALL iom_put( 'hfxldgrow'  ,      qlead       * r1_Dt_ice      )   ! Heat in lead for ice growth
 
       ! controls
