Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4737)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5837)
@@ -9,5 +9,4 @@
    !!----------------------------------------------------------------------
    !!   lim_wri      : write of the diagnostics variables in ouput file 
-   !!   lim_wri_init : initialization and namelist read
    !!   lim_wri_state : write for initial state or/and abandon
    !!----------------------------------------------------------------------
@@ -25,5 +24,4 @@
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! work arrays
-   USE par_ice
    USE iom
    USE timing          ! Timing
@@ -36,23 +34,4 @@
    PUBLIC lim_wri_state  ! called by dia_wri_state 
 
-   INTEGER, PARAMETER ::   jpnoumax = 43   !: maximum number of variable for ice output
-   
-   INTEGER  ::   noumef             ! number of fields
-   INTEGER  ::   noumefa            ! number of additional fields
-   INTEGER  ::   add_diag_swi       ! additional diagnostics
-   INTEGER  ::   nz                                         ! dimension for the itd field
-
-   REAL(wp) , DIMENSION(jpnoumax) ::   cmulti         ! multiplicative constant
-   REAL(wp) , DIMENSION(jpnoumax) ::   cadd           ! additive constant
-   REAL(wp) , DIMENSION(jpnoumax) ::   cmultia        ! multiplicative constant
-   REAL(wp) , DIMENSION(jpnoumax) ::   cadda          ! additive constant
-   CHARACTER(len = 35), DIMENSION(jpnoumax) ::   titn, titna   ! title of the field
-   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   nam , nama    ! name of the field
-   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   uni , unia    ! unit of the field
-   INTEGER            , DIMENSION(jpnoumax) ::   nc  , nca     ! switch for saving field ( = 1 ) or not ( = 0 )
-
-   REAL(wp)  ::   epsi06 = 1.e-6_wp
-   REAL(wp)  ::   zzero  = 0._wp
-   REAL(wp)  ::   zone   = 1._wp      
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
@@ -78,300 +57,219 @@
       INTEGER, INTENT(in) ::   kindic   ! if kindic < 0 there has been an error somewhere
       !
-      INTEGER ::  ji, jj, jk, jl, jf, ipl ! dummy loop indices
-      INTEGER ::  ierr
-      REAL(wp),DIMENSION(1) ::   zdept
-      REAL(wp) ::  zsto, zjulian, zout, zindh, zinda, zindb, zindc
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zcmo, zcmoa
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zfield
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zmaskitd, zoi, zei
-
-      CHARACTER(len = 60) ::   clhstnam, clop, clhstnama
-
-      INTEGER , SAVE ::   nice, nhorid, ndim, niter, ndepid
-      INTEGER , SAVE ::   nicea, nhorida, ndimitd
-      INTEGER , ALLOCATABLE, DIMENSION(:), SAVE ::   ndex51
-      INTEGER , ALLOCATABLE, DIMENSION(:), SAVE ::   ndexitd
+      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
+      REAL(wp) ::  z1_365
+      REAL(wp) ::  ztmp
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace
       !!-------------------------------------------------------------------
 
       IF( nn_timing == 1 )  CALL timing_start('limwri')
 
-      CALL wrk_alloc( jpi, jpj, zfield )
-      CALL wrk_alloc( jpi, jpj, jpnoumax, zcmo, zcmoa )
-      CALL wrk_alloc( jpi, jpj, jpl, zmaskitd, zoi, zei )
-
-      ipl = jpl
-
-      IF( numit == nstart ) THEN 
-
-         ALLOCATE( ndex51(jpij), ndexitd(jpij*jpl), STAT=ierr )
-         IF( lk_mpp    )   CALL mpp_sum ( ierr )
-         IF( ierr /= 0 ) THEN
-            CALL ctl_stop( 'lim_wri : unable to allocate standard arrays' )   ;   RETURN
-         ENDIF
-
-         CALL lim_wri_init 
-
-         IF(lwp) WRITE(numout,*) ' lim_wri, first time step '
-         IF(lwp) WRITE(numout,*) ' add_diag_swi ', add_diag_swi
-
-         !--------------------
-         !  1) Initialization
-         !--------------------
-
-         !-------------
-         ! Normal file
-         !-------------
-         niter    = ( nit000 - 1 ) / nn_fsbc
-         CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian )
-         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
-!clem
-!         zsto     = rdt_ice
-!         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
-!         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
-!         ENDIF
-!         zout     = nwrite * rdt_ice / nn_fsbc
-!         zdept(1) = 0.
-!
-!         CALL dia_nam ( clhstnam, nwrite, 'icemod_old' )
-!         CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice,   &
-!            &           nhorid, nice, domain_id=nidom, snc4chunks=snc4set )
-!         CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down")
-!         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
-!
-!         DO jf = 1 , noumef
-!            IF(lwp) WRITE(numout,*) 'jf', jf
-!            IF ( nc(jf) == 1 ) THEN
-!               CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj &
-!                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
-!               IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout'
-!               IF(lwp) WRITE(numout,*)  nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout 
-!            ENDIF
-!         END DO
-!
-!         CALL histend(nice, snc4set)
-!clem
-         !
-         !-----------------
-         ! ITD file output
-         !-----------------
-         zsto     = rdt_ice
-         clop     = "ave(x)"
-         zout     = nwrite * rdt_ice / nn_fsbc
-         zdept(1) = 0.
-
-         CALL dia_nam ( clhstnama, nwrite, 'icemoa' )
-         CALL histbeg ( clhstnama, jpi, glamt, jpj, gphit,         &
-            1, jpi, 1, jpj,            & ! zoom
-            niter, zjulian, rdt_ice,   & ! time
-            nhorida,                   & ! ? linked with horizontal ...
-            nicea , domain_id=nidom, snc4chunks=snc4set)                  ! file 
-         CALL histvert( nicea, "icethi", "L levels","m", ipl , hi_mean , nz )
+      CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
+      CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi )
+
+      !-----------------------------
+      ! Mean category values
+      !-----------------------------
+      z1_365 = 1._wp / 365._wp
+
+      CALL lim_var_icetm      ! mean sea ice temperature
+
+      CALL lim_var_bv         ! brine volume
+
+      DO jj = 1, jpj          ! presence indicator of ice
+         DO ji = 1, jpi
+            zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) )
+         END DO
+      END DO
+      !
+      !
+      !                                             
+      IF ( iom_use( "icethic_cea" ) ) THEN                       ! mean ice thickness
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)
+            END DO
+         END DO
+         CALL iom_put( "icethic_cea"  , z2d              )
+      ENDIF
+
+      IF ( iom_use( "snowthic_cea" ) ) THEN                      ! snow thickness = mean snow thickness over the cell 
+         DO jj = 1, jpj                                            
+            DO ji = 1, jpi
+               z2d(ji,jj)  = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)
+            END DO
+         END DO
+         CALL iom_put( "snowthic_cea" , z2d              )       
+      ENDIF
+      !
+      IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 
+         DO jj = 2 , jpjm1
+            DO ji = 2 , jpim1
+               z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
+               z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+           END DO
+         END DO
+         CALL lbc_lnk( z2da, 'T', -1. )
+         CALL lbc_lnk( z2db, 'T', -1. )
+         CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component
+         CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component
+         DO jj = 1, jpj                                 
+            DO ji = 1, jpi
+               z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) ) 
+            END DO
+         END DO
+         CALL iom_put( "icevel"       , z2d              )       ! ice velocity module
+      ENDIF
+      !
+      IF ( iom_use( "miceage" ) ) THEN 
+         z2d(:,:) = 0.e0
          DO jl = 1, jpl
-            zmaskitd(:,:,jl) = tmask(:,:,1)
-         END DO
-         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
-         CALL wheneq( jpi*jpj*jpl, zmaskitd, 1, 1., ndexitd, ndimitd  )  
-         CALL histdef( nicea, "iice_itd", "Ice area in categories"         , "-"    ,   &  
-            jpi, jpj, nhorida, jpl, 1, jpl, nz, 15, clop, zsto, zout )
-         CALL histdef( nicea, "iice_hid", "Ice thickness in categories"    , "m"    ,   &  
-            jpi, jpj, nhorida, jpl, 1, jpl, nz, 15, clop, zsto, zout )
-         CALL histdef( nicea, "iice_hsd", "Snow depth in in categories"    , "m"    ,   &  
-            jpi, jpj, nhorida, jpl, 1, jpl, nz, 15, clop, zsto, zout )
-         CALL histdef( nicea, "iice_std", "Ice salinity distribution"      , "ppt"  ,   &  
-            jpi, jpj, nhorida, jpl, 1, jpl, nz, 15, clop, zsto, zout )
-         CALL histdef( nicea, "iice_otd", "Ice age distribution"               , "days",   &  
-            jpi, jpj, nhorida, jpl, 1, jpl, nz, 15, clop, zsto, zout )
-         CALL histdef( nicea, "iice_etd", "Brine volume distr. "               , "%"    ,   &  
-            jpi, jpj, nhorida, jpl, 1, jpl, nz, 15, clop, zsto, zout )
-         CALL histend(nicea, snc4set)
-      ENDIF
-
-      !     !-----------------------------------------------------------------------!
-      !     !--2. Computation of instantaneous values                               ! 
-      !     !-----------------------------------------------------------------------!
-
-      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      !IF( ln_nicep ) THEN
-      !   WRITE(numout,*)
-      !   WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit
-      !   WRITE(numout,*) '~~~~~~~ '
-      !   WRITE(numout,*) ' kindic = ', kindic
-      !ENDIF
-      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-      !-- calculs des valeurs instantanees
-      zcmo ( 1:jpi, 1:jpj, 1:jpnoumax ) = 0._wp
-      zcmoa( 1:jpi, 1:jpj, 1:jpnoumax ) = 0._wp
-
-      ! Ice surface temperature and some fluxes
-      DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) )
+                  z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age
+      ENDIF
+
+      IF ( iom_use( "micet" ) ) THEN 
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zinda  = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi06 ) )
-               zcmo(ji,jj,17) = zcmo(ji,jj,17) + a_i(ji,jj,jl)*qsr_ice (ji,jj,jl) 
-               zcmo(ji,jj,18) = zcmo(ji,jj,18) + a_i(ji,jj,jl)*qns_ice(ji,jj,jl) 
-               zcmo(ji,jj,27) = zcmo(ji,jj,27) + zinda*(t_su(ji,jj,jl)-rtt)*a_i(ji,jj,jl)/MAX(at_i(ji,jj),epsi06)
-               zcmo(ji,jj,21) = zcmo(ji,jj,21) + zinda*oa_i(ji,jj,jl)/MAX(at_i(ji,jj),epsi06) 
-            END DO
-         END DO
-      END DO
-
-      ! Mean sea ice temperature
-      CALL lim_var_icetm
-
-      ! Brine volume
-      CALL lim_var_bv
-
-      DO jj = 2 , jpjm1
-         DO ji = 2 , jpim1
-            zinda  = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi06 ) )
-            zindb  = MAX( zzero , SIGN( zone , at_i(ji,jj) ) )
-
-            zcmo(ji,jj,1)  = at_i(ji,jj)
-            zcmo(ji,jj,2)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zinda
-            zcmo(ji,jj,3)  = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zinda
-            zcmo(ji,jj,4)  = diag_bot_gr(ji,jj) * rday     ! Bottom thermodynamic ice production
-            zcmo(ji,jj,5)  = diag_dyn_gr(ji,jj) * rday     ! Dynamic ice production (rid/raft)
-            zcmo(ji,jj,22) = diag_lat_gr(ji,jj) * rday     ! Lateral thermodynamic ice production
-            zcmo(ji,jj,23) = diag_sni_gr(ji,jj) * rday     ! Snow ice production ice production
-            zcmo(ji,jj,24) = (tm_i(ji,jj) - rtt) * zinda
-
-            zcmo(ji,jj,6)  = fbif(ji,jj)*at_i(ji,jj)
-            zcmo(ji,jj,7)  = (  u_ice(ji,jj) * tmu(ji,jj) + u_ice(ji-1,jj) * tmu(ji-1,jj) ) * 0.5_wp
-            zcmo(ji,jj,8)  = (  v_ice(ji,jj) * tmv(ji,jj) + v_ice(ji,jj-1) * tmv(ji,jj-1) ) * 0.5_wp
-            zcmo(ji,jj,9)  = sst_m(ji,jj)
-            zcmo(ji,jj,10) = sss_m(ji,jj)
-
-            zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
-            zcmo(ji,jj,12) = qsr(ji,jj)
-            zcmo(ji,jj,13) = qns(ji,jj)
-            zcmo(ji,jj,14) = fhbri(ji,jj)
-            zcmo(ji,jj,15) = utau_ice(ji,jj)
-            zcmo(ji,jj,16) = vtau_ice(ji,jj)
-            zcmo(ji,jj,17) = zcmo(ji,jj,17) + ( 1._wp - at_i(ji,jj) ) * qsr(ji,jj)
-            zcmo(ji,jj,18) = zcmo(ji,jj,18) + ( 1._wp - at_i(ji,jj) ) * qns(ji,jj)
-            zcmo(ji,jj,19) = sprecip(ji,jj)
-            zcmo(ji,jj,20) = smt_i(ji,jj)
-            zcmo(ji,jj,25) = et_i(ji,jj)
-            zcmo(ji,jj,26) = et_s(ji,jj)
-            zcmo(ji,jj,28) = sfx_bri(ji,jj)
-            zcmo(ji,jj,29) = sfx_thd(ji,jj)
-
-            zcmo(ji,jj,30) = bv_i(ji,jj)
-            zcmo(ji,jj,31) = hicol(ji,jj) * zindb
-            zcmo(ji,jj,32) = strength(ji,jj)
-            zcmo(ji,jj,33) = SQRT(  zcmo(ji,jj,7)*zcmo(ji,jj,7) + zcmo(ji,jj,8)*zcmo(ji,jj,8)  )
-            zcmo(ji,jj,34) = diag_sur_me(ji,jj) * rday     ! Surface melt
-            zcmo(ji,jj,35) = diag_bot_me(ji,jj) * rday     ! Bottom melt
-            zcmo(ji,jj,36) = divu_i(ji,jj)
-            zcmo(ji,jj,37) = shear_i(ji,jj)
-            zcmo(ji,jj,38) = diag_res_pr(ji,jj) * rday     ! Bottom melt
-            zcmo(ji,jj,39) = vt_i(ji,jj)  ! ice volume
-            zcmo(ji,jj,40) = vt_s(ji,jj)  ! snow volume
-
-            zcmo(ji,jj,41) = sfx_mec(ji,jj)
-            zcmo(ji,jj,42) = sfx_res(ji,jj)
-
-            zcmo(ji,jj,43) = diag_trp_vi(ji,jj) * rday     ! transport of ice volume
-
-        END DO
-      END DO
-
-      !
-      ! ecriture d'un fichier netcdf
-      !
-      niter = niter + 1
-!clem
-!      DO jf = 1 , noumef
-!         !
-!         zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf)
-!         !
-!         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN   ;   CALL lbc_lnk( zfield, 'T', -1. )
-!         ELSE                                                            ;   CALL lbc_lnk( zfield, 'T',  1. )
-!         ENDIF
-!         !
-!         IF( ln_nicep ) THEN 
-!            WRITE(numout,*)
-!            WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim'
-!            WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim
-!         ENDIF
-!         IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
-!         !
-!      END DO
-!
-!      IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN
-!         IF( lwp) WRITE(numout,*) ' Closing the icemod file '
-!         CALL histclo( nice )
-!      ENDIF
-!clem
-      !
-       CALL iom_put ('iceconc', zcmo(:,:,1) )          ! field1: ice concentration
-       CALL iom_put ('icethic_cea', zcmo(:,:,2) )      ! field2: ice thickness (i.e. icethi(:,:))
-       CALL iom_put ('snowthic_cea', zcmo(:,:,3))      ! field3: snow thickness
-       CALL iom_put ('icebopr', zcmo(:,:,4) )   ! field4: daily bottom thermo ice production
-       CALL iom_put ('icedypr', zcmo(:,:,5) )   ! field5: daily dynamic ice production
-       CALL iom_put ('ioceflxb', zcmo(:,:,6) )         ! field6: Oceanic flux at the ice base
-       CALL iom_put ('uice_ipa', zcmo(:,:,7) )         ! field7: ice velocity u component
-       CALL iom_put ('vice_ipa', zcmo(:,:,8) )         ! field8: ice velocity v component
-       CALL iom_put ('isst', zcmo(:,:,9) )             ! field 9: sea surface temperature
-       CALL iom_put ('isss', zcmo(:,:,10) )            ! field 10: sea surface salinity
-       CALL iom_put ('qt_oce', zcmo(:,:,11) )           ! field 11: total flux at ocean surface
-       CALL iom_put ('qsr_oce', zcmo(:,:,12) )          ! field 12: solar flux at ocean surface
-       CALL iom_put ('qns_oce', zcmo(:,:,13) )          ! field 13: non-solar flux at ocean surface
-       !CALL iom_put ('hfbri', fhbri )                  ! field 14: heat flux due to brine release
-       CALL iom_put( 'utau_ice', zcmo(:,:,15)  )     ! Wind stress over ice along i-axis at I-point
-       CALL iom_put( 'vtau_ice', zcmo(:,:,16) )     ! Wind stress over ice along j-axis at I-point
-       CALL iom_put ('qsr_io', zcmo(:,:,17) )          ! field 17: solar flux at ice/ocean surface
-       CALL iom_put ('qns_io', zcmo(:,:,18) )          ! field 18: non-solar flux at ice/ocean surface
-       !CALL iom_put ('snowpre', zcmo(:,:,19) * rday ! field 19 :snow precip          
-       CALL iom_put ('micesalt', zcmo(:,:,20) )        ! field 20 :mean ice salinity
-       CALL iom_put ('miceage', zcmo(:,:,21) / 365)    ! field 21: mean ice age
-       CALL iom_put ('icelapr',zcmo(:,:,22) )   ! field 22: daily lateral thermo ice prod.
-       CALL iom_put ('icesipr',zcmo(:,:,23) )   ! field 23: daily snowice ice prod.
-       CALL iom_put ('micet', zcmo(:,:,24) )           ! field 24: mean ice temperature
-       CALL iom_put ('icehc', zcmo(:,:,25) )           ! field 25: ice total heat content
-       CALL iom_put ('isnowhc', zcmo(:,:,26) )         ! field 26: snow total heat content
-       CALL iom_put ('icest', zcmo(:,:,27) )           ! field 27: ice surface temperature
-       CALL iom_put ('sfxbri', zcmo(:,:,28) * rday )           ! field 28: brine salt flux
-       CALL iom_put ('sfxthd', zcmo(:,:,29) * rday )           ! field 29: equivalent FW salt flux
-       CALL iom_put ('ibrinv', zcmo(:,:,30) *100 )     ! field 30: brine volume
-       CALL iom_put ('icecolf', zcmo(:,:,31) )         ! field 31: frazil ice collection thickness
-       CALL iom_put ('icestr', zcmo(:,:,32) * 0.001 )  ! field 32: ice strength
-       CALL iom_put ('icevel', zcmo(:,:,33) )          ! field 33: ice velocity
-       CALL iom_put ('isume', zcmo(:,:,34) )    ! field 34: surface melt
-       CALL iom_put ('ibome', zcmo(:,:,35) )     ! field 35: bottom melt
-       CALL iom_put ('idive', zcmo(:,:,36) * 1.0e8)    ! field 36: divergence
-       CALL iom_put ('ishear', zcmo(:,:,37) * 1.0e8 )  ! field 37: shear
-       CALL iom_put ('icerepr', zcmo(:,:,38) ) ! field 38: daily prod./melting due to limupdate
-       CALL iom_put ('icevolu', zcmo(:,:,39) ) ! field 39: ice volume
-       CALL iom_put ('snowvol', zcmo(:,:,40) ) ! field 40: snow volume
-       CALL iom_put ('sfxmec', zcmo(:,:,41) * rday )           ! field 41: salt flux from ridging rafting
-       CALL iom_put ('sfxres', zcmo(:,:,42) * rday )           ! field 42: salt flux from limupdate (resultant)
-       CALL iom_put ('icetrp', zcmo(:,:,43) )    ! field 43: ice volume transport
-
-      !-----------------------------
-      ! Thickness distribution file
-      !-----------------------------
-      IF( add_diag_swi == 1 ) THEN
-
-         DO jl = 1, jpl 
-            CALL lbc_lnk( a_i(:,:,jl)  , 'T' ,  1. )
-            CALL lbc_lnk( sm_i(:,:,jl) , 'T' ,  1. )
-            CALL lbc_lnk( oa_i(:,:,jl) , 'T' ,  1. )
-            CALL lbc_lnk( ht_i(:,:,jl) , 'T' ,  1. )
-            CALL lbc_lnk( ht_s(:,:,jl) , 'T' ,  1. )
-         END DO
-
-         ! Compute ice age
+               z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj)
+            END DO
+         END DO
+         CALL iom_put( "micet"       , z2d               )        ! mean ice temperature
+      ENDIF
+      !
+      IF ( iom_use( "icest" ) ) THEN 
+         z2d(:,:) = 0.e0
+         DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "icest"       , z2d              )        ! ice surface temperature
+      ENDIF
+
+      IF ( iom_use( "icecolf" ) ) THEN 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) ) )
+               z2d(ji,jj) = hicol(ji,jj) * rswitch
+            END DO
+         END DO
+         CALL iom_put( "icecolf"     , z2d              )        ! frazil ice collection thickness
+      ENDIF
+
+      CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature
+      CALL iom_put( "isss"        , sss_m               )        ! sea surface salinity
+      CALL iom_put( "iceconc"     , at_i                )        ! ice concentration
+      CALL iom_put( "icevolu"     , vt_i                )        ! ice volume = mean ice thickness over the cell
+      CALL iom_put( "icehc"       , et_i                )        ! ice total heat content
+      CALL iom_put( "isnowhc"     , et_s                )        ! snow total heat content
+      CALL iom_put( "ibrinv"      , bv_i * 100._wp      )        ! brine volume
+      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point
+      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point
+      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation 
+      CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity
+
+      CALL iom_put( "icestr"      , strength * 0.001    )        ! ice strength
+      CALL iom_put( "idive"       , divu_i * 1.0e8      )        ! divergence
+      CALL iom_put( "ishear"      , shear_i * 1.0e8     )        ! shear
+      CALL iom_put( "snowvol"     , vt_s                )        ! snow volume
+      
+      CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport
+      CALL iom_put( "snwtrp"      , diag_trp_vs * rday  )        ! snw volume transport
+      CALL iom_put( "saltrp"      , diag_trp_smv * rday * rhoic ) ! salt content transport
+      CALL iom_put( "deitrp"      , diag_trp_ei         )        ! advected ice enthalpy (W/m2)
+      CALL iom_put( "destrp"      , diag_trp_es         )        ! advected snw enthalpy (W/m2)
+
+      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting
+      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant)
+      CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines
+      CALL iom_put( "sfx"         , sfx     * rday      )        ! total salt flux
+
+      ztmp = rday / rhoic
+      CALL iom_put( "vfxres"     , wfx_res * ztmp       )        ! daily prod./melting due to limupdate 
+      CALL iom_put( "vfxopw"     , wfx_opw * ztmp       )        ! daily lateral thermodynamic ice production
+      CALL iom_put( "vfxsni"     , wfx_sni * ztmp       )        ! daily snowice ice production
+      CALL iom_put( "vfxbog"     , wfx_bog * ztmp       )        ! daily bottom thermodynamic ice production
+      CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp       )        ! daily dynamic ice production (rid/raft)
+      CALL iom_put( "vfxsum"     , wfx_sum * ztmp       )        ! surface melt 
+      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt 
+      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt 
+      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
+      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow) 
+      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
+      
+      CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total)
+      CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics)
+      CALL iom_put( "afxthd"     , afx_thd * rday       )        ! concentration tendency (thermo)
+
+      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !  
+      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !  
+      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   !  
+      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   !  
+      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   !  
+      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !  
+      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !  
+      CALL iom_put ('hfxerr'     , hfx_err(:,:)         )   !  
+      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !  
+      
+      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   !  
+      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   !  
+      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   !  
+      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !  
+      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !  
+      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base 
+      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
+      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
+      
+      !--------------------------------
+      ! Output values for each category
+      !--------------------------------
+      CALL iom_put( "iceconc_cat"      , a_i         )        ! area for categories
+      CALL iom_put( "icethic_cat"      , ht_i        )        ! thickness for categories
+      CALL iom_put( "snowthic_cat"     , ht_s        )        ! snow depth for categories
+      CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories
+
+      ! ice temperature
+      IF ( iom_use( "icetemp_cat" ) ) THEN 
+         zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i
+         CALL iom_put( "icetemp_cat"   , zt_i - rt0  )
+      ENDIF
+      
+      ! snow temperature
+      IF ( iom_use( "snwtemp_cat" ) ) THEN 
+         zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s
+         CALL iom_put( "snwtemp_cat"   , zt_s - rt0  )
+      ENDIF
+
+      ! Compute ice age
+      IF ( iom_use( "iceage_cat" ) ) THEN 
          DO jl = 1, jpl 
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - epsi06 ) )
-                  zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , epsi06 ) * zinda
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) )
+                  rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) )
+                  zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch
                END DO
             END DO
          END DO
-
-         ! Compute brine volume
+         CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories
+      ENDIF
+
+      ! Compute brine volume
+      IF ( iom_use( "brinevol_cat" ) ) THEN 
          zei(:,:,:) = 0._wp
          DO jl = 1, jpl 
@@ -379,38 +277,21 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zinda = MAX( zzero , SIGN( zone , a_i(ji,jj,jl) - epsi06 ) )
-                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* &
-                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * &
-                        zinda / nlay_i
+                     rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
+                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 *  &
+                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * &
+                        rswitch * r1_nlay_i
                   END DO
                END DO
             END DO
          END DO
-
-         DO jl = 1, jpl 
-            CALL lbc_lnk( zei(:,:,jl) , 'T' ,  1. )
-         END DO
-
-         CALL histwrite( nicea, "iice_itd", niter, a_i  , ndimitd , ndexitd  )   ! area
-         CALL histwrite( nicea, "iice_hid", niter, ht_i , ndimitd , ndexitd  )   ! thickness
-         CALL histwrite( nicea, "iice_hsd", niter, ht_s , ndimitd , ndexitd  )   ! snow depth
-         CALL histwrite( nicea, "iice_std", niter, sm_i , ndimitd , ndexitd  )   ! salinity
-         CALL histwrite( nicea, "iice_otd", niter, zoi  , ndimitd , ndexitd  )   ! age
-         CALL histwrite( nicea, "iice_etd", niter, zei  , ndimitd , ndexitd  )   ! brine volume
-
-         !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
-         !     IF( kindic < 0 )   CALL lim_wri_state( 'output.abort' )
-         !     not yet implemented
-
-         IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN
-            IF(lwp) WRITE(numout,*) ' Closing the icemod file '
-            CALL histclo( nicea ) 
-         ENDIF
-         !
-      ENDIF
-
-      CALL wrk_dealloc( jpi, jpj, zfield )
-      CALL wrk_dealloc( jpi, jpj, jpnoumax, zcmo, zcmoa )
-      CALL wrk_dealloc( jpi, jpj, jpl, zmaskitd, zoi, zei )
+         CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories
+      ENDIF
+
+      !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
+      !     IF( kindic < 0 )   CALL lim_wri_state( 'output.abort' )
+      !     not yet implemented
+      
+      CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
+      CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')
@@ -419,124 +300,4 @@
 #endif
 
-   SUBROUTINE lim_wri_init
-      !!-------------------------------------------------------------------
-      !!                    ***   ROUTINE lim_wri_init  ***
-      !!                
-      !! ** Purpose :   ???
-      !!
-      !! ** Method  : Read the namicewri namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namicewri
-      !!-------------------------------------------------------------------
-      INTEGER ::   nf      ! ???
-      INTEGER ::   ios     ! Local integer output status for namelist read
-
-      TYPE FIELD 
-         CHARACTER(len = 35) :: ztitle 
-         CHARACTER(len = 8 ) :: zname          
-         CHARACTER(len = 8 ) :: zunit
-         INTEGER             :: znc   
-         REAL                :: zcmulti 
-         REAL                :: zcadd        
-      END TYPE FIELD
-
-      TYPE(FIELD) ::  &
-         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
-         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
-         field_13, field_14, field_15, field_16, field_17, field_18,   &
-         field_19, field_20, field_21, field_22, field_23, field_24,   &
-         field_25, field_26, field_27, field_28, field_29, field_30,   &
-         field_31, field_32, field_33, field_34, field_35, field_36,   &
-         field_37, field_38, field_39, field_40, field_41, field_42, field_43
-
-      TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield
-      !
-      NAMELIST/namiceout/ noumef, &
-         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
-         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
-         field_13, field_14, field_15, field_16, field_17, field_18,   &
-         field_19, field_20, field_21, field_22, field_23, field_24,   &
-         field_25, field_26, field_27, field_28, field_29, field_30,   &
-         field_31, field_32, field_33, field_34, field_35, field_36,   &
-         field_37, field_38, field_39, field_40, field_41, field_42, field_43, add_diag_swi
-      !!-------------------------------------------------------------------
-      REWIND( numnam_ice_ref )              ! Namelist namiceout in reference namelist : Ice outputs
-      READ  ( numnam_ice_ref, namiceout, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceout in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namiceout in configuration namelist : Ice outputs
-      READ  ( numnam_ice_cfg, namiceout, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceout in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namiceout )
-
-      zfield(1)  = field_1
-      zfield(2)  = field_2
-      zfield(3)  = field_3
-      zfield(4)  = field_4
-      zfield(5)  = field_5
-      zfield(6)  = field_6
-      zfield(7)  = field_7
-      zfield(8)  = field_8
-      zfield(9)  = field_9
-      zfield(10) = field_10
-      zfield(11) = field_11
-      zfield(12) = field_12
-      zfield(13) = field_13
-      zfield(14) = field_14
-      zfield(15) = field_15
-      zfield(16) = field_16
-      zfield(17) = field_17
-      zfield(18) = field_18
-      zfield(19) = field_19
-      zfield(20) = field_20
-      zfield(21) = field_21
-      zfield(22) = field_22
-      zfield(23) = field_23
-      zfield(24) = field_24
-      zfield(25) = field_25
-      zfield(26) = field_26
-      zfield(27) = field_27
-      zfield(28) = field_28
-      zfield(29) = field_29
-      zfield(30) = field_30
-      zfield(31) = field_31
-      zfield(32) = field_32
-      zfield(33) = field_33
-      zfield(34) = field_34
-      zfield(35) = field_35
-      zfield(36) = field_36
-      zfield(37) = field_37
-      zfield(38) = field_38
-      zfield(39) = field_39
-      zfield(40) = field_40
-      zfield(41) = field_41
-      zfield(42) = field_42
-      zfield(43) = field_43
-
-      DO nf = 1, noumef
-         titn  (nf) = zfield(nf)%ztitle
-         nam   (nf) = zfield(nf)%zname
-         uni   (nf) = zfield(nf)%zunit
-         nc    (nf) = zfield(nf)%znc
-         cmulti(nf) = zfield(nf)%zcmulti
-         cadd  (nf) = zfield(nf)%zcadd
-      END DO
-
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_wri_init : Ice parameters for outputs'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '    number of fields to be stored         noumef = ', noumef
-         WRITE(numout,*) '           title                            name     unit      Saving (1/0) ',   &
-            &            '    multiplicative constant       additive constant '
-         DO nf = 1 , noumef         
-            WRITE(numout,*) '   ', titn(nf), '   '    , nam   (nf), '      '  , uni (nf),   &
-               &            '  ' , nc  (nf),'        ', cmulti(nf), '        ', cadd(nf)
-         END DO
-         WRITE(numout,*) ' add_diag_swi ', add_diag_swi
-      ENDIF
-      !
-   END SUBROUTINE lim_wri_init
  
    SUBROUTINE lim_wri_state( kt, kid, kh_i )
@@ -555,35 +316,49 @@
       INTEGER, INTENT( in ) ::   kid , kh_i       
       !!----------------------------------------------------------------------
-      !CALL histvert( kid, "icethi", "L levels","m", jpl , hi_mean , nz )
-
-      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 )
-      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      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, "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, "iice_itd", "Ice concentration by cat", "%"      , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 
-      !CALL histdef( kid, "iice_hid", "Ice thickness by cat"    , "m"      , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 
-      !CALL histdef( kid, "iice_hsd", "Snow thickness by cat"   , "m"      , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 
-      !CALL histdef( kid, "iice_std", "Ice salinity by cat"     , "PSU"    , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 
+
+      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 )
+      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      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, "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, "iisfxdyn", "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
@@ -591,5 +366,5 @@
       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, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
       CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
       CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
@@ -603,19 +378,17 @@
       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, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) )
       CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
 
-      !CALL histwrite( kid, "iice_itd", kt, a_i  , jpi*jpj*jpl, (/1/)  )   ! area
-      !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/)  )   ! thickness
-      !CALL histwrite( kid, "iice_hsd", kt, ht_s , jpi*jpj*jpl, (/1/)  )   ! snow depth
-      !CALL histwrite( kid, "iice_std", kt, sm_i , jpi*jpj*jpl, (/1/)  )   ! salinity
+      ! Close the file
+      ! -----------------
+      !CALL histclo( kid )
 
     END SUBROUTINE lim_wri_state
