Index: NEMO/trunk/src/ICE/ice.F90
===================================================================
--- NEMO/trunk/src/ICE/ice.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/ice.F90	(revision 9912)
@@ -256,6 +256,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping                [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations      [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_atm_oi      !: heat flux available for thermo transformations      [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_oce_ai   !: heat flux remaining at the end of thermo transformations  [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation                   [kg.m-2.s-1]
    
@@ -387,19 +387,19 @@
 
       ii = ii + 1
-      ALLOCATE( t_bo   (jpi,jpj) , wfx_snw_sni(jpi,jpj) ,                                                &
-         &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) ,  &
-         &      wfx_ice(jpi,jpj) , wfx_sub    (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam    (jpi,jpj) ,  &
-         &      wfx_pnd(jpi,jpj) ,                                                                       &
-         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,              &
-         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,              &
-         &      afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj),                                                  &
-         &      fhtur  (jpi,jpj) , qlead  (jpi,jpj) ,                                                    &
-         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
-         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  &
-         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) ,                        & 
-         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld   (jpi,jpj) ,                        &
-         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     &
-         &      hfx_opw(jpi,jpj) , hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,     &
-         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj)        ,  STAT=ierr(ii) )
+      ALLOCATE( t_bo     (jpi,jpj) , wfx_snw_sni(jpi,jpj) ,                                                &
+         &      wfx_snw  (jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) ,  &
+         &      wfx_ice  (jpi,jpj) , wfx_sub    (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam    (jpi,jpj) ,  &
+         &      wfx_pnd  (jpi,jpj) ,                                                                       &
+         &      wfx_bog  (jpi,jpj) , wfx_dyn   (jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,           &
+         &      wfx_res  (jpi,jpj) , wfx_sni   (jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,           &
+         &      afx_tot  (jpi,jpj) , rn_amax_2d(jpi,jpj),                                                  &
+         &      fhtur    (jpi,jpj) , qlead     (jpi,jpj) ,                                                 &
+         &      sfx_res  (jpi,jpj) , sfx_bri   (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
+         &      sfx_bog  (jpi,jpj) , sfx_bom   (jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  &
+         &      hfx_res  (jpi,jpj) , hfx_snw   (jpi,jpj) , hfx_sub(jpi,jpj) ,                        & 
+         &      qt_atm_oi(jpi,jpj) , qt_oce_ai (jpi,jpj) , fhld   (jpi,jpj) ,                        &
+         &      hfx_sum  (jpi,jpj) , hfx_bom   (jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     &
+         &      hfx_opw  (jpi,jpj) , hfx_thd   (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,     &
+         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj)             ,  STAT=ierr(ii) )
 
       ! * Ice global state variables
Index: NEMO/trunk/src/ICE/ice1d.F90
===================================================================
--- NEMO/trunk/src/ICE/ice1d.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/ice1d.F90	(revision 9912)
@@ -53,5 +53,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_out_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qt_oce_ai_1d
 
    ! heat flux associated with ice-atmosphere mass exchange
@@ -190,5 +190,5 @@
          &      hfx_thd_1d(jpij) , hfx_spr_1d    (jpij) ,                      &
          &      hfx_snw_1d(jpij) , hfx_sub_1d    (jpij) ,                      &
-         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , hfx_out_1d(jpij), STAT=ierr(ii) )
+         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) )
       !
       ii = ii + 1
Index: NEMO/trunk/src/ICE/icectl.F90
===================================================================
--- NEMO/trunk/src/ICE/icectl.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/icectl.F90	(revision 9912)
@@ -189,5 +189,5 @@
 
       ! heat flux
-      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es   &
+      zhfx  = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat - diag_trp_ei - diag_trp_es   &
       !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
          &              ) * e1e2t ) * zconv
@@ -572,6 +572,6 @@
                WRITE(numout,*)
                WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj)
-               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj)
-               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj)
+               WRITE(numout,*) ' qt_atm_oi    : ', qt_atm_oi(ji,jj)
+               WRITE(numout,*) ' qt_oce_ai    : ', qt_oce_ai(ji,jj)
                WRITE(numout,*) ' dhc          : ', diag_heat(ji,jj)              
                WRITE(numout,*)
Index: NEMO/trunk/src/ICE/icedia.F90
===================================================================
--- NEMO/trunk/src/ICE/icedia.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/icedia.F90	(revision 9912)
@@ -95,9 +95,9 @@
       ! 2 - Trends due to forcing  !
       ! ---------------------------!
-      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! freshwater flux ice/snow-ocean 
-      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9                     ! freshwater flux ice/snow-atm
-      z_frc_sal    = r1_rau0 * glob_sum(   - sfx(:,:) * e1e2t(:,:) ) * 1.e-9                                          ! salt fluxes ice/snow-ocean
-      z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ocean (and below ice)
-      z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ice-coean
+      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-ocean 
+      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) )                    * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-atm
+      z_frc_sal    = r1_rau0 * glob_sum( -       sfx(:,:)                                     * e1e2t(:,:) ) * 1.e-9   ! salt fluxes ice/snow-ocean
+      z_frc_tembot =           glob_sum(   qt_oce_ai(:,:)                                     * e1e2t(:,:) ) * 1.e-20  ! heat on top of ocean (and below ice)
+      z_frc_temtop =           glob_sum(   qt_atm_oi(:,:)                                     * e1e2t(:,:) ) * 1.e-20  ! heat on top of ice-coean
       !
       frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
@@ -110,7 +110,7 @@
       ! 3 -  Content variations !
       ! ----------------------- !
-      zdiff_vol = r1_rau0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3) 
-      zdiff_sal = r1_rau0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss)
-      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J)
+      zdiff_vol = r1_rau0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3) 
+      zdiff_sal = r1_rau0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 )  - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss)
+      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J)
       !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
 
Index: NEMO/trunk/src/ICE/icethd.F90
===================================================================
--- NEMO/trunk/src/ICE/icethd.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/icethd.F90	(revision 9912)
@@ -177,5 +177,5 @@
             ! Net heat flux on top of the ice-ocean [W.m-2]
             ! ---------------------------------------------
-            hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
+            qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
          END DO
       END DO
@@ -185,7 +185,7 @@
       ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
       IF( .NOT. ln_icedH ) THEN
-         hfx_in(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
-         fhtur (:,:) = 0._wp
-         fhld  (:,:) = 0._wp
+         qt_atm_oi(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
+         fhtur    (:,:) = 0._wp
+         fhld     (:,:) = 0._wp
       ENDIF
 
@@ -196,9 +196,9 @@
       !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema) 
       !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar
-      hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean               
-         &           - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation
-         &           - at_i (:,:) * fhtur(:,:)                               &  ! heat flux taken by turbulence
-         &           - at_i (:,:) *  fhld(:,:)                                  ! heat flux taken during bottom growth/melt 
-                                                                                !    (fhld should be 0 while bott growth)
+      qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean               
+         &             - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation
+         &             - at_i (:,:) * fhtur(:,:)                               &  ! heat flux taken by turbulence
+         &             - at_i (:,:) *  fhld(:,:)                                  ! heat flux taken during bottom growth/melt 
+      !                                                                            !    (fhld should be 0 while bott growth)
       !-------------------------------------------------------------------------------------------!
       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
@@ -429,5 +429,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   )
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem   )
-         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_out_1d    (1:npti), hfx_out       )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     )
          !
          ! SIMIP diagnostics
@@ -507,17 +507,17 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        )
          !
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub        )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d    (1:npti), hfx_spr     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d    (1:npti), hfx_sum     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d    (1:npti), hfx_bom     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d    (1:npti), hfx_bog     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d    (1:npti), hfx_dif     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d    (1:npti), hfx_opw     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d    (1:npti), hfx_snw     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d    (1:npti), hfx_sub     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     )
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) )
Index: NEMO/trunk/src/ICE/icethd_dh.F90
===================================================================
--- NEMO/trunk/src/ICE/icethd_dh.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/icethd_dh.F90	(revision 9912)
@@ -570,5 +570,5 @@
          !    
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
-         hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
+         qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
 
          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
Index: NEMO/trunk/src/ICE/iceupdate.F90
===================================================================
--- NEMO/trunk/src/ICE/iceupdate.F90	(revision 9910)
+++ NEMO/trunk/src/ICE/iceupdate.F90	(revision 9912)
@@ -107,6 +107,6 @@
       ! --- case we bypass ice thermodynamics --- !
       IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
-         hfx_in     (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
-         hfx_out    (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
+         qt_atm_oi  (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
+         qt_oce_ai  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
          qtr_ice_bot(:,:,:) = 0._wp
          emp_ice    (:,:)   = 0._wp
@@ -122,18 +122,18 @@
             zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
 
-            ! Total heat flux reaching the ocean = hfx_out (W.m-2) 
+            ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 
             !---------------------------------------------------
-            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
-            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
+            zqmass           = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
+            qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr
 
             ! Add the residual from heat diffusion equation and sublimation (W.m-2)
             !----------------------------------------------------------------------
-            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
-               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
+            qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) +   &
+               &             ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
 
             ! New qsr and qns used to compute the oceanic heat flux at the next time step
             !----------------------------------------------------------------------------
             qsr(ji,jj) = zqsr                                      
-            qns(ji,jj) = hfx_out(ji,jj) - zqsr              
+            qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr              
 
             ! Mass flux at the atm. surface       
@@ -254,6 +254,6 @@
       IF( iom_use('qt_oce'     ) )   CALL iom_put( "qt_oce"     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
       IF( iom_use('qt_ice'     ) )   CALL iom_put( "qt_ice"     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice )
-      IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( "qt_oce_ai"  , hfx_out * tmask(:,:,1)                                     )   ! total heat flux at the ocean   surface: interface oce-(ice+atm) 
-      IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( "qt_atm_oi"  , hfx_in * tmask(:,:,1)                                      )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce) 
+      IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( "qt_oce_ai"  , qt_oce_ai * tmask(:,:,1)                                   )   ! total heat flux at the ocean   surface: interface oce-(ice+atm) 
+      IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( "qt_atm_oi"  , qt_atm_oi * tmask(:,:,1)                                   )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce) 
       IF( iom_use('qemp_oce'   ) )   CALL iom_put( "qemp_oce"   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean
       IF( iom_use('qemp_ice'   ) )   CALL iom_put( "qemp_ice"   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice
@@ -267,5 +267,5 @@
       IF( iom_use('hfxdif'     ) )   CALL iom_put ("hfxdif"     , hfx_dif             )   ! heat flux used for ice temperature change
       IF( iom_use('hfxsnw'     ) )   CALL iom_put ("hfxsnw"     , hfx_snw             )   ! heat flux used for snow melt 
-      IF( iom_use('hfxerr'     ) )   CALL iom_put ("hfxerr"     , hfx_err_dif         )   ! heat flux error after heat diffusion (included in hfx_out)
+      IF( iom_use('hfxerr'     ) )   CALL iom_put ("hfxerr"     , hfx_err_dif         )   ! heat flux error after heat diffusion (included in qt_oce_ai)
 
       ! heat fluxes associated with mass exchange (freeze/melt/precip...)
