Index: /NEMO/trunk/src/ICE/ice.F90
===================================================================
--- /NEMO/trunk/src/ICE/ice.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/ice.F90	(revision 9916)
@@ -360,6 +360,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_si          !: Temperature at Snow-ice interface (K) 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tm_si         !: mean temperature at the snow-ice interface (K) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qcn_ice_bot   !: Bottom  conduction flux (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qcn_ice_top   !: Surface conduction flux (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_bot   !: Bottom  conduction flux (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_top   !: Surface conduction flux (W/m2)
 
    !
@@ -451,5 +451,5 @@
       ! * SIMIP diagnostics
       ii = ii + 1
-      ALLOCATE( t_si (jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj) , qcn_ice_top(jpi,jpj) , STAT = ierr(ii) )
+      ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
Index: /NEMO/trunk/src/ICE/ice1d.F90
===================================================================
--- /NEMO/trunk/src/ICE/ice1d.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/ice1d.F90	(revision 9916)
@@ -109,6 +109,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ib_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_s_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sum      !: Ice surface ablation [m]
@@ -206,5 +204,5 @@
       ii = ii + 1
       ALLOCATE( t_su_1d (jpij) , t_si_1d (jpij) , a_i_1d    (jpij) , a_ib_1d (jpij) ,                   &
-         &      h_i_1d  (jpij) , h_ib_1d (jpij) , h_s_1d    (jpij) , fc_su   (jpij) , fc_bo_i (jpij) ,  &    
+         &      h_i_1d  (jpij) , h_ib_1d (jpij) , h_s_1d    (jpij) ,                                    &    
          &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &    
          &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  &
Index: /NEMO/trunk/src/ICE/icestp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icestp.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/icestp.F90	(revision 9916)
@@ -427,6 +427,6 @@
 
       ! SIMIP diagnostics
-      qcn_ice_bot(:,:) = 0._wp ; qcn_ice_top(:,:) = 0._wp
-      t_si(:,:,:) = rt0       ! temp at the ice-snow interface
+      qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes
+      t_si       (:,:,:) = rt0   ! temp at the ice-snow interface
 
       tau_icebfr(:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
Index: /NEMO/trunk/src/ICE/icethd.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/icethd.F90	(revision 9916)
@@ -431,7 +431,4 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     )
          !
-         ! SIMIP diagnostics
-         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot )
-         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top )
          ! ocean surface fields
          CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
@@ -527,7 +524,7 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
          ! SIMIP diagnostics         
-         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si(:,:,kl) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot  )
-         CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top  )
+         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) )
+         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) )
          ! extensive variables
          CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
Index: /NEMO/trunk/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_dh.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/icethd_dh.F90	(revision 9916)
@@ -137,5 +137,5 @@
          !
          DO ji = 1, npti
-            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - fc_su(ji)
+            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
             qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
             zq_su(ji)      = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
@@ -145,5 +145,5 @@
          !
          DO ji = 1, npti
-            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - fc_su(ji) 
+            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji) 
             qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
             zq_su(ji)      = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
@@ -153,5 +153,5 @@
       !
       DO ji = 1, npti
-         zf_tt(ji)         = fc_bo_i(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
+         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
          zq_bo(ji)         = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO
@@ -416,7 +416,7 @@
       !------------------
       ! Basal growth is driven by heat imbalance at the ice-ocean interface,
-      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux 
+      ! between the inner conductive flux  (qcn_ice_bot), from the open water heat flux 
       ! (fhld) and the sensible ice-ocean flux (qsb_ice_bot). 
-      ! fc_bo_i is positive downwards. qsb_ice_bot and fhld are positive to the ice 
+      ! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice 
 
       ! If salinity varies in time, an iterative procedure is required, because
Index: /NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_zdf_bl99.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/icethd_zdf_bl99.F90	(revision 9916)
@@ -728,12 +728,12 @@
       !-----------------------------
       !
-      ! --- update conduction fluxes
-      !
+      ! --- calculate conduction fluxes (positive downward)
+
       DO ji = 1, npti
          !                                ! 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) )
+         qcn_ice_top_1d(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
-         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji) - t_i_1d(ji,nlay_i) )
+         qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
       END DO
       
@@ -750,5 +750,5 @@
          !
          DO ji = 1, npti
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji) 
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 
          END DO
          !
@@ -770,12 +770,15 @@
                
                IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-                  zhfx_err = (qns_ice_1d(ji) + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice)*a_i_1d(ji)
+                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
+                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
                ELSE                          ! case T_su = 0degC
-                  zhfx_err = (fc_su(ji)      + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice)*a_i_1d(ji)
+                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
+                     &       + zdq * r1_rdtice ) * a_i_1d(ji)
                ENDIF
                
             ELSEIF( k_jules == np_jules_ACTIVE ) THEN
             
-               zhfx_err = ( fc_su(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
+                  &          + zdq * r1_rdtice ) * a_i_1d(ji)
             
             ENDIF
@@ -787,21 +790,4 @@
             hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
             !
-         END DO
-         !
-         ! --- SIMIP diagnostics
-         !
-         DO ji = 1, npti
-            !--- Conduction fluxes (positive downward)
-            qcn_ice_bot_1d(ji) = qcn_ice_bot_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
-            qcn_ice_top_1d(ji) = qcn_ice_top_1d(ji) + fc_su  (ji) * a_i_1d(ji) / at_i_1d(ji)
-   
-            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-            IF( h_s_1d(ji) >= zhs_min ) THEN
-               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
-            ELSE
-               t_si_1d(ji) = t_su_1d(ji)
-            ENDIF
          END DO
          !
@@ -827,11 +813,23 @@
       IF( k_jules == np_jules_EMULE ) THEN
          ! 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)
+         t_s_1d    (1:npti,:) = ztsold        (1:npti,:)
+         t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
+         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
       ENDIF
       !
+      ! --- SIMIP diagnostics
+      !
+      DO ji = 1, npti         
+         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
+         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
+         IF( h_s_1d(ji) >= zhs_min ) THEN
+            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
+               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
+         ELSE
+            t_si_1d(ji) = t_su_1d(ji)
+         ENDIF
+      END DO
+      !
    END SUBROUTINE ice_thd_zdf_BL99
-
 
 #else
Index: /NEMO/trunk/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceupdate.F90	(revision 9915)
+++ /NEMO/trunk/src/ICE/iceupdate.F90	(revision 9916)
@@ -277,7 +277,7 @@
 
       ! other heat fluxes
-      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"  , qcn_ice_bot  * at_i_b )   ! Bottom conduction flux
-      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( "hfxcndtop"  , qcn_ice_top  * at_i_b )   ! Surface conduction 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
 
       ! diags
