Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4566)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4626)
@@ -34,6 +34,6 @@
    !!PUBLIC   lim_diahsb_rst   ! routine called by ice_init.F90
 
-   REAL(dp) ::   frc_sal, frc_vol   ! global forcing trends
-   REAL(dp) ::   bg_grme            ! global ice growth+melt trends
+   REAL(wp) ::   frc_sal, frc_vol   ! global forcing trends
+   REAL(wp) ::   bg_grme            ! global ice growth+melt trends
    REAL(wp) ::   epsi06 = 1.e-6_wp  ! small number
    REAL(wp) ::   epsi03 = 1.e-3_wp  ! small number
@@ -44,9 +44,8 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.4 , NEMO Consortium (2012)
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
    !! $Id: limdiahsb.F90 3294 2012-10-18 16:44:18Z rblod $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -58,11 +57,10 @@
       !!	
       !!---------------------------------------------------------------------------
-      !!
-      REAL(dp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
-      REAL(dp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_thd, zbg_sfx_res, zbg_sfx_mec 
-      REAL(dp)   ::   zbg_emp, zbg_emp_bog, zbg_emp_lag, zbg_emp_sig, zbg_emp_dyg, zbg_emp_bom, zbg_emp_sum, zbg_emp_res 
-      REAL(dp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
-      REAL(dp)   ::   z1_area                     !    -     -
-      REAL(dp)   ::   zinda, zindb
+      REAL(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
+      REAL(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_thd, zbg_sfx_res, zbg_sfx_mec 
+      REAL(wp)   ::   zbg_emp, zbg_emp_bog, zbg_emp_lag, zbg_emp_sig, zbg_emp_dyg, zbg_emp_bom, zbg_emp_sum, zbg_emp_res 
+      REAL(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
+      REAL(wp)   ::   z1_area, zcoef                    
+      REAL(wp)   ::   zinda, zindb
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('lim_diahsb')
@@ -70,8 +68,7 @@
       IF( numit == nstart ) CALL lim_diahsb_init 
 
-      ! 1/area
-      z1_area = 1.d0 / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )
-
-      zinda = MAX( 0.d0 , SIGN( 1.d0 , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )
+      z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )
+
+      zinda = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )
       ! ----------------------- !
       ! 1 -  Content variations !
@@ -80,5 +77,5 @@
       zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow
       zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area
-      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )       ! mean salt content
+      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 )      * area(:,:) * tms(:,:) )  ! mean salt content
       zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content
 
@@ -89,18 +86,22 @@
       zbg_shc = glob_sum( et_s(:,:) * 1.e-11 ) ! snow heat content [10^9*1.e-11 J]
 
-      zbg_emp     = zinda * glob_sum(         emp(:,:) * area(:,:) * tms(:,:) ) * z1_area         * r1_rau0 * rday
-      zbg_emp_bog = zinda * glob_sum( diag_bot_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-      zbg_emp_lag = zinda * glob_sum( diag_lat_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-      zbg_emp_sig = zinda * glob_sum( diag_sni_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-      zbg_emp_dyg = zinda * glob_sum( diag_dyn_gr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-      zbg_emp_bom = zinda * glob_sum( diag_bot_me(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-      zbg_emp_sum = zinda * glob_sum( diag_sur_me(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-      zbg_emp_res = zinda * glob_sum( diag_res_pr(:,:) * area(:,:) * tms(:,:) ) * z1_area * rhoic * r1_rau0 * rday
-
-      zbg_sfx     = zinda * glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
-      zbg_sfx_bri = zinda * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
-      zbg_sfx_thd = zinda * glob_sum( sfx_thd(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
-      zbg_sfx_res = zinda * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
-      zbg_sfx_mec = zinda * glob_sum( sfx_mec(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zcoef = zinda * z1_area * r1_rau0 * rday
+      zbg_emp     =  glob_sum(         emp(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      
+      zcoef = zinda * z1_area * rhoic * r1_rau0 * rday
+      zbg_emp_bog = glob_sum( diag_bot_gr(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_emp_lag = glob_sum( diag_lat_gr(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_emp_sig = glob_sum( diag_sni_gr(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_emp_dyg = glob_sum( diag_dyn_gr(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_emp_bom = glob_sum( diag_bot_me(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_emp_sum = glob_sum( diag_sur_me(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_emp_res = glob_sum( diag_res_pr(:,:) * area(:,:) * tms(:,:) ) * zcoef
+
+      zcoef = zinda * z1_area * r1_rau0 * rday
+      zbg_sfx     = glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_sfx_bri = glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_sfx_thd = glob_sum( sfx_thd(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_sfx_res = glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_sfx_mec = glob_sum( sfx_mec(:,:) * area(:,:) * tms(:,:) ) * zcoef
       
       ! --------------------------------------------- !
@@ -109,10 +110,11 @@
       z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes
       z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes
-      z_bg_grme = glob_sum( ( diag_bot_gr(:,:) + diag_lat_gr(:,:) + diag_sni_gr(:,:) + diag_dyn_gr(:,:) + &
-                          &    diag_bot_me(:,:) + diag_sur_me(:,:) + diag_res_pr(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes
-      !
-      frc_vol  = frc_vol  + z_frc_vol  * rdt_ice
-      frc_sal  = frc_sal  + z_frc_sal  * rdt_ice
-      bg_grme  = bg_grme  + z_bg_grme  * rdt_ice
+      z_bg_grme = glob_sum(  (  diag_bot_gr(:,:) + diag_lat_gr(:,:) + diag_sni_gr(:,:) + diag_dyn_gr(:,:)     &
+         &                    + diag_bot_me(:,:) + diag_sur_me(:,:) + diag_res_pr(:,:)                    )   &
+         &                 * area(:,:) * tms(:,:) ) ! volume fluxes
+      !
+      frc_vol  = frc_vol + z_frc_vol * rdt_ice
+      frc_sal  = frc_sal + z_frc_sal * rdt_ice
+      bg_grme  = bg_grme + z_bg_grme * rdt_ice
       
       ! difference
@@ -123,5 +125,5 @@
       ! 3 - Diagnostics writing !
       ! ----------------------- !
-      zindb = MAX( 0.d0 , SIGN( 1.d0 , zbg_ivo - epsi06 ) )
+      zindb = MAX( 0.d0 , SIGN( 1._wp , zbg_ivo - epsi06 ) )
       !
       CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)         
@@ -156,5 +158,5 @@
       !
       IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb')
-!
+      !
    END SUBROUTINE lim_diahsb
 
@@ -190,11 +192,8 @@
       ! 2 - initial conservation variables !
       ! ---------------------------------- !
-      !frc_vol = 0.d0                                           ! volume       trend due to forcing
-      !frc_sal = 0.d0                                           ! salt content   -    -   -    -         
-      !bg_grme = 0.d0                                           ! ice growth + melt volume trend
-      !
       CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE lim_diahsb_init
+
 
    SUBROUTINE lim_diahsb_rst( kt, cdrw )
@@ -226,7 +225,7 @@
            IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state '
            IF(lwp) WRITE(numout,*) '~~~~~~~'
-           frc_vol  = 0.d0                                           
-           frc_sal  = 0.d0                                                  
-           bg_grme  = 0.d0                                        
+           frc_vol  = 0._wp                                           
+           frc_sal  = 0._wp                                                  
+           bg_grme  = 0._wp                                        
        ENDIF	
 
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4566)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4626)
@@ -145,4 +145,5 @@
       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis
       REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms
+      REAL(wp) ::   ztemp   ! local scalar 
       !!------------------------------------------------------------------     
       ! 
@@ -328,13 +329,15 @@
             DO layer = 1, nlay_i-1
                DO ji = kideb , kiut
-                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
-                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   &
-                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )  
+                  ztemp = t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2._wp * rtt
+                  ztcond_i(ji,layer) = rcdic + 0.0900_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
+                     &                                   / MIN( -2.0_wp * epsi10, ztemp )   &
+                     &                       - 0.0055_wp * ztemp
                   ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
                END DO
             END DO
             DO ji = kideb , kiut
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   &
-                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )  
+               ztemp = t_bo_b(ji) - rtt
+               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN( -epsi10, ztemp )   &
+                  &                        - 0.011_wp * ztemp  
                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
             END DO
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 4566)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 4626)
@@ -53,5 +53,5 @@
       !
       INTEGER  ::   ji, jk     ! dummy loop indices 
-      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
+      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, ice_switch,  ztmelts   ! local scalars
       REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
       REAL(wp), POINTER, DIMENSION(:) ::   ze_init, zhiold, zsiold
@@ -99,18 +99,17 @@
             ! Switches 
             !----------
-            iflush       =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )    ! =1 if summer 
-            igravdr      =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
-            iaccrbo      =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )    ! =1 if bottom accretion
-            i_ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) )
-            isnowic      = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * i_ice_switch   ! =1 if snow ice formation
-
+            iflush     =         MAX( 0._wp , SIGN( 1.0 , t_su_b(ji) - rtt )        )           ! =1 if summer 
+            igravdr    =         MAX( 0._wp , SIGN( 1.0 , t_bo_b(ji) - t_su_b(ji) ) )           ! =1 if t_su < t_bo
+            iaccrbo    =         MAX( 0._wp , SIGN( 1.0 , dh_i_bott(ji) )           )           ! =1 if bottom accretion
+            ice_switch = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + 1.e-2 ) )           ! =1 ice thickness > 1cm
+            isnowic    = 1._wp - MAX ( 0._wp , SIGN( 1._wp , - dh_snowice(ji) ) ) * ice_switch  ! =1 if snow ice formation
+            !
             !---------------------
             ! Salinity tendencies
             !---------------------
-            !                                   ! drainage by gravity drainage
-            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
-            !                                   ! drainage by flushing  
-            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
-
+            !                               
+            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  ! drainage by gravity
+            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice  ! drainage by flushing 
+            !
             !-----------------
             ! Update salinity   
@@ -119,19 +118,19 @@
             ! snow ice / bottom sources are added in lim_thd_ent to conserve energy
             sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
-
+            !
             ! if no ice, salinity = 0.1
-            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
-            sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch )
-
+            ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
+            sm_i_b(ji) = ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - ice_switch )
+            !
             !----------------------------
             ! Heat flux - brine drainage
             !----------------------------
             fhbri_1d(ji) = 0._wp
-
+            !
             !----------------------------
             ! Salt flux - brine drainage
             !----------------------------
-            sfx_bri_1d(ji) = sfx_bri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) * r1_rdtice
-
+            sfx_bri_1d(ji) = sfx_bri_1d(ji) - ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) * r1_rdtice
+            !
          END DO
 
Index: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4566)
+++ branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4626)
@@ -331,6 +331,4 @@
       !!----------------------------------------------------------------------
     
-
-    
       CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
@@ -338,49 +336,49 @@
       CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicestru", "i-Wind stress over ice"  , "Pa"     , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicestrv", "j-Wind stress over ice"  , "Pa"     , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicesflx", "Solar flux over ocean"   , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean","w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
       CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
       CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
       CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxthd", "Salt flux from thermo"   , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxmec", "Salt flux from dynmics"  , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iisfxthd", "Salt flux from thermo"   , ""       , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iisfxmec", "Salt flux from dynmics"  , ""       , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
 
       CALL histend( kid, snc4set )   ! end of the file definition
 
-      CALL histwrite( kid, "iicethic", kt, icethi        , jpi*jpj, (/1/) )    
-      CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicetemp", kt, tm_i - rtt    , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
-
-      CALL histwrite( kid, "iicebopr", kt, diag_bot_gr        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedypr", kt, diag_dyn_gr        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicelapr", kt, diag_lat_gr        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesipr", kt, diag_sni_gr        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicerepr", kt, diag_res_pr        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicebome", kt, diag_bot_me        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesume", kt, diag_sur_me        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxthd", kt, sfx_thd        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxmec", kt, sfx_mec        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicethic", kt, icethi      , jpi*jpj, (/1/) )    
+      CALL histwrite( kid, "iiceconc", kt, at_i        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicetemp", kt, tm_i - rtt  , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicevelu", kt, u_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicevelv", kt, v_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicestru", kt, utau_ice    , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicestrv", kt, vtau_ice    , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesflx", kt, qsr         , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicenflx", kt, qns         , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "isnowpre", kt, sprecip     , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesali", kt, smt_i       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicevolu", kt, vt_i        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8, jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "iicebopr", kt, diag_bot_gr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicedypr", kt, diag_dyn_gr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicelapr", kt, diag_lat_gr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesipr", kt, diag_sni_gr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicerepr", kt, diag_res_pr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicebome", kt, diag_bot_me , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesume", kt, diag_sur_me , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxthd", kt, sfx_thd     , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxmec", kt, sfx_mec     , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxres", kt, sfx_res     , jpi*jpj, (/1/) )
 
       ! Close the file
