Index: NEMO/branches/2020/dev_14237_KERNEL-01_IMMERSE_SEAMOUNT/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/2020/dev_14237_KERNEL-01_IMMERSE_SEAMOUNT/src/ICE/iceupdate.F90	(revision 14238)
+++ NEMO/branches/2020/dev_14237_KERNEL-01_IMMERSE_SEAMOUNT/src/ICE/iceupdate.F90	(revision 14676)
@@ -104,5 +104,10 @@
       ! Net heat flux on top of the ice-ocean (W.m-2)
       !----------------------------------------------
-      qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:)
+      IF( ln_cndflx ) THEN   ! ice-atm interface = conduction (and melting) fluxes
+         qt_atm_oi(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) + &
+            &             SUM( a_i_b * ( qcn_ice + qml_ice + qtr_ice_top ), dim=3 ) + qemp_ice(:,:)
+      ELSE                   ! ice-atm interface = solar and non-solar fluxes
+         qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 
+      ENDIF
 
       ! --- case we bypass ice thermodynamics --- !
@@ -119,13 +124,19 @@
          ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2)
          !---------------------------------------------------
-         zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
+         IF( ln_cndflx ) THEN   ! ice-atm interface = conduction (and melting) fluxes
+            zqsr = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) + SUM( a_i_b (ji,jj,:) * qtr_ice_bot(ji,jj,:) )
+         ELSE                   ! ice-atm interface = solar and non-solar fluxes
+            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
+         ENDIF
 
          ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2)
          !---------------------------------------------------
-         qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) &
-            &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) &
-            &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) &
-            &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)
-
+         IF( ln_icethd ) THEN
+            qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) &
+               &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) &
+               &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) &
+               &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)
+         ENDIF
+         
          ! New qsr and qns used to compute the oceanic heat flux at the next time step
          !----------------------------------------------------------------------------
@@ -144,5 +155,5 @@
          !
          ! the non-solar is simply derived from the solar flux
-         qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr
+         qns(ji,jj) = qt_oce_ai(ji,jj) - qsr(ji,jj)
 
          ! Mass flux at the atm. surface
@@ -275,7 +286,10 @@
 
       ! other heat fluxes
-      IF( iom_use('hfxsensib'  ) )   CALL iom_put( 'hfxsensib'  ,     -qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux
+      IF( iom_use('hfxsensib'  ) )   CALL iom_put( 'hfxsensib'  ,      qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux
       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
 
       ! controls
@@ -345,5 +359,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.0_wp, tmod_io, 'T', 1.0_wp )
+         CALL lbc_lnk( 'iceupdate', taum, 'T', 1.0_wp, tmod_io, 'T', 1.0_wp )
          !
          utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
@@ -374,5 +388,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.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition
+      CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition
       !
       IF( ln_timing )   CALL timing_stop('ice_update')
