Index: branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 5975)
+++ branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 5976)
@@ -32,4 +32,16 @@
 !------------------------------------------------------------------------------
    ln_iceini      = .true.         !  activate ice initialization (T) or not (F)
+   ln_iceini_file = .true.         !  ice initialization from a netcdf file (T) or constant (F)
+!
+!         !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!         !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_hti = 'Ice_initialization'   , -12    ,'hti'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_hts = 'Ice_initialization'   , -12    ,'hts'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_ati = 'Ice_initialization'   , -12    ,'ati'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_tsu = 'Ice_initialization'   , -12    ,'tsu'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_tmi = 'Ice_initialization'   , -12    ,'tmi'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_smi = 'Ice_initialization'   , -12    ,'smi'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , ''
+   cn_dir='./'
+   !
    rn_thres_sst   =  2.0           !  maximum water temperature with initial ice (degC)
    rn_hts_ini_n   =  0.3           !  initial real snow thickness (m), North
Index: branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5975)
+++ branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5976)
@@ -28,4 +28,6 @@
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE wrk_nemo         ! work arrays
+   USE fldread          ! read input fields
+   USE iom
 
    IMPLICIT NONE
@@ -47,5 +49,15 @@
    REAL(wp) ::   rn_tmi_ini_s   ! initial temperature
 
-   LOGICAL  ::  ln_iceini    ! initialization or not
+   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
+   INTEGER , PARAMETER ::   jp_hti = 1           ! index of thick (m)    at T-point
+   INTEGER , PARAMETER ::   jp_hts = 2           ! index of thick (m)    at T-point
+   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
+   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
+   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
+   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si    ! structure of input fields (file informations, fields read)
+
+   LOGICAL  ::  ln_iceini        ! initialization or not
+   LOGICAL  ::  ln_iceini_file   ! Ice initialization state from 2D netcdf file
    !!----------------------------------------------------------------------
    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
@@ -88,15 +100,12 @@
       REAL(wp)   :: ztmelts, zdh
       INTEGER    :: i_hemis, i_fill, jl0  
-      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
-      REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
-      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
-      !--------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, zswitch )
-      CALL wrk_alloc( jpi, jpj, zhemis )
-      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
-      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
+      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini               !data by cattegories to fill
+      !--------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
+      CALL wrk_alloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
 
       IF(lwp) WRITE(numout,*)
@@ -117,266 +126,288 @@
 
       ! basal temperature (considered at freezing point)
-      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
-      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
+      t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1) 
 
       IF( ln_iceini ) THEN
 
-      !--------------------------------------------------------------------
-      ! 2) Basal temperature, ice mask and hemispheric index
-      !--------------------------------------------------------------------
-
-      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
-         DO ji = 1, jpi
-            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 
-               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
-            ELSE                                                                                   
-               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
-            ENDIF
+         !--------------------------------------------------------------------
+         ! 2) Basal temperature, ice mask and hemispheric index
+         !--------------------------------------------------------------------
+
+         DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
+            DO ji = 1, jpi
+               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 
+                  zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
+               ELSE                                                                                   
+                  zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
+               ENDIF
+            END DO
          END DO
-      END DO
-
-
-      ! Hemispheric index
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( fcor(ji,jj) >= 0._wp ) THEN    
-               zhemis(ji,jj) = 1 ! Northern hemisphere
-            ELSE
-               zhemis(ji,jj) = 2 ! Southern hemisphere
-            ENDIF
-         END DO
-      END DO
-
-      !--------------------------------------------------------------------
-      ! 3) Initialization of sea ice state variables
-      !--------------------------------------------------------------------
-
-      !-----------------------------
-      ! 3.1) Hemisphere-dependent arrays
-      !-----------------------------
-      ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
-      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness
-      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth
-      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration
-      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity
-      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow)
-
-      zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume
-
-      !---------------------------------------------------------------------
-      ! 3.2) Distribute ice concentration and thickness into the categories
-      !---------------------------------------------------------------------
-      ! a gaussian distribution for ice concentration is used
-      ! then we check whether the distribution fullfills
-      ! volume and area conservation, positivity and ice categories bounds
-      DO i_hemis = 1, 2
-
-      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
-
-      ! note for the great nemo engineers: 
-      ! only very few of the WRITE statements are necessary for the reference version
-      ! they were one day useful, but now i personally doubt of their
-      ! potential for bringing anything useful
-
-      DO i_fill = jpl, 1, -1
-         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
-            !----------------------------
-            ! fill the i_fill categories
-            !----------------------------
-            ! *** 1 category to fill
-            IF ( i_fill .EQ. 1 ) THEN
-               zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis)
-               za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis)
-               zh_i_ini(2:jpl,i_hemis)   = 0._wp
-               za_i_ini(2:jpl,i_hemis)   = 0._wp
-            ELSE
-
-               ! *** >1 categores to fill
-               !--- Ice thicknesses in the i_fill - 1 first categories
-               DO jl = 1, i_fill - 1
-                  zh_i_ini(jl,i_hemis) = hi_mean(jl)
-               END DO
-               
-               !--- jl0: most likely index where cc will be maximum
-               DO jl = 1, jpl
-                  IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. &
-                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN
-                     jl0 = jl
+
+         !--------------------------------------------------------------------
+         ! 3) Initialization of sea ice state variables
+         !--------------------------------------------------------------------
+         IF( ln_iceini_file )THEN
+
+            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
+            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
+            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
+            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
+            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
+            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
+
+         ELSE ! ln_iceini_file = F
+
+            !-----------------------------
+            ! 3.1) Hemisphere-dependent arrays
+            !-----------------------------
+            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( fcor(ji,jj) >= 0._wp ) THEN
+                     zht_i_ini(ji,jj) = rn_hti_ini_n
+                     zht_s_ini(ji,jj) = rn_hts_ini_n
+                     zat_i_ini(ji,jj) = rn_ati_ini_n
+                     zsm_i_ini(ji,jj) = rn_smi_ini_n
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_n
+                  ELSE
+                     zht_i_ini(ji,jj) = rn_hti_ini_s
+                     zht_s_ini(ji,jj) = rn_hts_ini_s
+                     zat_i_ini(ji,jj) = rn_ati_ini_s
+                     zsm_i_ini(ji,jj) = rn_smi_ini_s
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_s
                   ENDIF
                END DO
-               jl0 = MIN(jl0, i_fill)
-               
-               !--- Concentrations
-               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
-               DO jl = 1, i_fill - 1
-                  IF ( jl .NE. jl0 ) THEN
-                     zsigma               = 0.5 * zht_i_ini(i_hemis)
-                     zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma
-                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
-                  ENDIF
-               END DO
-               
-               zA = 0. ! sum of the areas in the jpl categories 
-               DO jl = 1, i_fill - 1
-                 zA = zA + za_i_ini(jl,i_hemis)
-               END DO
-               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category
-               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
-         
-               !--- Ice thickness in the last category
-               zV = 0. ! sum of the volumes of the N-1 categories
-               DO jl = 1, i_fill - 1
-                  zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis)
-               END DO
-               zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
-               IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
-
-               !--- volumes
-               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis)
-               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
-
-            ENDIF ! i_fill
-
-            !---------------------
-            ! Compatibility tests
-            !---------------------
-            ! Test 1: area conservation
-            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
-            IF ( zconv .LT. 1.0e-6 ) THEN
-               ztest_1 = 1
-            ELSE 
-              ! this write is useful
-              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
-               ztest_1 = 0
-            ENDIF
-
-            ! Test 2: volume conservation
-            zV_cons = SUM(zv_i_ini(:,i_hemis))
-            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)
-
-            IF ( zconv .LT. 1.0e-6 ) THEN
-               ztest_2 = 1
-            ELSE
-              ! this write is useful
-              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
-                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
-               ztest_2 = 0
-            ENDIF
-
-            ! Test 3: thickness of the last category is in-bounds ?
-            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
-               ztest_3 = 1
-            ELSE
-               ! this write is useful
-               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
-               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
-               ztest_3 = 0
-            ENDIF
-
-            ! Test 4: positivity of ice concentrations
-            ztest_4 = 1
-            DO jl = 1, jpl
-               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
-                  ! this write is useful
-                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
-                  ztest_4 = 0
-               ENDIF
-            END DO
-
-         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
- 
-         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
-
-      END DO ! i_fill
-
-      IF(lwp) THEN 
-         WRITE(numout,*) ' ztests : ', ztests
-         IF ( ztests .NE. 4 ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) ' !!!! ALERT                  !!! '
-            WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
-            WRITE(numout,*)
-            WRITE(numout,*) ' *** ztests is not equal to 4 '
-            WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
-            WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis)
-            WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis)
-         ENDIF ! ztests .NE. 4
-      ENDIF
-      
-      END DO ! i_hemis
-
-      !---------------------------------------------------------------------
-      ! 3.3) Space-dependent arrays for ice state variables
-      !---------------------------------------------------------------------
-
-      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
-      DO jl = 1, jpl ! loop over categories
+            END DO
+
+         ENDIF ! ln_iceini_file
+
+         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
+
+         !---------------------------------------------------------------------
+         ! 3.2) Distribute ice concentration and thickness into the categories
+         !---------------------------------------------------------------------
+         ! a gaussian distribution for ice concentration is used
+         ! then we check whether the distribution fullfills
+         ! volume and area conservation, positivity and ice categories bounds
+         zh_i_ini(:,:,:) = 0._wp 
+         za_i_ini(:,:,:) = 0._wp
+         zv_i_ini(:,:,:) = 0._wp
+
          DO jj = 1, jpj
             DO ji = 1, jpi
-               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
-               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness
-               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth
-               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
-               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day)
-               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
-
-               ! This case below should not be used if (ht_s/ht_i) is ok in namelist
-               ! In case snow load is in excess that would lead to transformation from snow to ice
-               ! Then, transfer the snow excess into the ice (different from limthd_dh)
-               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
-               ! recompute ht_i, ht_s avoiding out of bounds values
-               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
-               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
-
-               ! ice volume, salt content, age content
-               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
-               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
-               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
-               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
-            END DO
-         END DO
-      END DO
-
-      ! Snow temperature and heat content
-      DO jk = 1, nlay_s
+
+               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
+
+                  ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
+                  ztests  = 0 
+
+                  DO i_fill = jpl, 1, -1
+
+                     IF( ztests .NE. 4 ) THEN
+                        !----------------------------
+                        ! fill the i_fill categories
+                        !----------------------------
+                        ! *** 1 category to fill
+                        IF ( i_fill .EQ. 1 ) THEN
+                           zh_i_ini(ji,jj,    1)   = zht_i_ini(ji,jj)
+                           za_i_ini(ji,jj,    1)   = zat_i_ini(ji,jj)
+                           zh_i_ini(ji,jj,2:jpl)   = 0._wp
+                           za_i_ini(ji,jj,2:jpl)   = 0._wp
+                        ELSE
+
+                           ! *** >1 categores to fill
+                           !--- Ice thicknesses in the i_fill - 1 first categories
+                           DO jl = 1, i_fill - 1
+                              zh_i_ini(ji,jj,jl) = hi_mean(jl)
+                           END DO
+               
+                           !--- jl0: most likely index where cc will be maximum
+                           DO jl = 1, jpl
+                              IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. &
+                                 & ( zht_i_ini(ji,jj) <= hi_max(jl)   ) ) THEN
+                                 jl0 = jl
+                              ENDIF
+                           END DO
+                           jl0 = MIN(jl0, i_fill)
+               
+                           !--- Concentrations
+                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
+                           DO jl = 1, i_fill - 1
+                              IF( jl .NE. jl0 )THEN
+                                 zsigma             = 0.5 * zht_i_ini(ji,jj)
+                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma
+                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
+                              ENDIF
+                           END DO
+               
+                           zA = 0. ! sum of the areas in the jpl categories 
+                           DO jl = 1, i_fill - 1
+                              zA = zA + za_i_ini(ji,jj,jl)
+                           END DO
+                           za_i_ini(ji,jj,i_fill)   = zat_i_ini(ji,jj) - zA ! ice conc in the last category
+                           IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
+         
+                           !--- Ice thickness in the last category
+                           zV = 0. ! sum of the volumes of the N-1 categories
+                           DO jl = 1, i_fill - 1
+                              zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl)
+                           END DO
+                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 
+                           IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
+
+                           !--- volumes
+                           zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:)
+                           IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
+
+                        ENDIF ! i_fill
+
+                        !---------------------
+                        ! Compatibility tests
+                        !---------------------
+                        ! Test 1: area conservation
+                        zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons )
+                        IF ( zconv .LT. 1.0e-6 ) THEN
+                           ztest_1 = 1
+                        ELSE 
+                          !this write is useful
+                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 
+                          ztest_1 = 0
+                        ENDIF
+
+                        ! Test 2: volume conservation
+                        zV_cons = SUM(zv_i_ini(ji,jj,:))
+                        zconv = ABS(zvt_i_ini(ji,jj) - zV_cons)
+
+                        IF( zconv .LT. 1.0e-6 ) THEN
+                           ztest_2 = 1
+                        ELSE
+                           !this write is useful
+                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
+                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj)
+                           ztest_2 = 0
+                        ENDIF
+
+                        ! Test 3: thickness of the last category is in-bounds ?
+                        IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN
+                           ztest_3 = 1
+                        ELSE
+                           ! this write is useful
+                           IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', &
+                           zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
+                           ztest_3 = 0
+                        ENDIF
+
+                        ! Test 4: positivity of ice concentrations
+                        ztest_4 = 1
+                        DO jl = 1, jpl
+                           IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 
+                              ! this write is useful
+                              IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl)
+                              ztest_4 = 0
+                           ENDIF
+                        END DO
+
+                     ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
+ 
+                     ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
+
+                  END DO ! i_fill
+
+                  IF(lwp) THEN 
+                     WRITE(numout,*) ' ztests : ', ztests
+                     IF( ztests .NE. 4 )THEN
+                        WRITE(numout,*)
+                        WRITE(numout,*) ' !!!! ALERT                  !!! '
+                        WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
+                        WRITE(numout,*)
+                        WRITE(numout,*) ' *** ztests is not equal to 4 '
+                        WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
+                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
+                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
+                     ENDIF ! ztests .NE. 4
+                  ENDIF
+      
+               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp
+
+            ENDDO   
+         ENDDO   
+
+         !---------------------------------------------------------------------
+         ! 3.3) Space-dependent arrays for ice state variables
+         !---------------------------------------------------------------------
+
+         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
          DO jl = 1, jpl ! loop over categories
             DO jj = 1, jpj
                DO ji = 1, jpi
-                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
-                   ! Snow energy of melting
-                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
-
-                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
-                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
+                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
+                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
+                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
+                  o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                                    ! age (1 day)
+                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
+
+                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
+                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
+                  ELSE
+                    ht_s(ji,jj,jl)= 0._wp
+                  ENDIF
+
+                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
+                  ! In case snow load is in excess that would lead to transformation from snow to ice
+                  ! Then, transfer the snow excess into the ice (different from limthd_dh)
+                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
+                  ! recompute ht_i, ht_s avoiding out of bounds values
+                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
+                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
+
+                  ! ice volume, salt content, age content
+                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
+                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
+                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
+                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
                END DO
             END DO
          END DO
-      END DO
-
-      ! Ice salinity, temperature and heat content
-      DO jk = 1, nlay_i
-         DO jl = 1, jpl ! loop over categories
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
-                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
-                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
-
-                   ! heat content per unit volume
-                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
-                      -   rcp     * ( ztmelts - rt0 ) )
-
-                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
-                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
+
+         ! Snow temperature and heat content
+         DO jk = 1, nlay_s
+            DO jl = 1, jpl ! loop over categories
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
+                     ! Snow energy of melting
+                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+
+                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
+                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
+                  END DO
                END DO
             END DO
          END DO
-      END DO
-
-      tn_ice (:,:,:) = t_su (:,:,:)
-
-      ELSE 
-         ! if ln_iceini=false
+
+         ! Ice salinity, temperature and heat content
+         DO jk = 1, nlay_i
+            DO jl = 1, jpl ! loop over categories
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
+                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
+                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
+
+                     ! heat content per unit volume
+                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
+                        -   rcp     * ( ztmelts - rt0 ) )
+
+                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
+                  END DO
+               END DO
+            END DO
+         END DO
+
+         tn_ice (:,:,:) = t_su (:,:,:)
+
+      ELSE ! if ln_iceini=false
          a_i  (:,:,:) = 0._wp
          v_i  (:,:,:) = 0._wp
@@ -400,5 +431,5 @@
             END DO
          END DO
-      
+
       ENDIF ! ln_iceini
       
@@ -452,8 +483,6 @@
 
 
-      CALL wrk_dealloc( jpi, jpj, zswitch )
-      CALL wrk_dealloc( jpi, jpj, zhemis )
-      CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
-      CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
+      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
+      CALL wrk_dealloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
 
    END SUBROUTINE lim_istate
@@ -475,7 +504,17 @@
       !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
       !!-----------------------------------------------------------------------------
-      NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  &
-         &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s
-      INTEGER :: ios                 ! Local integer output status for namelist read
+      !
+      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
+      INTEGER :: ji,jj
+      INTEGER :: ifpr, ierror
+      !
+      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files
+      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
+      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
+      !
+      NAMELIST/namiceini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
+         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
+         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
+         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
       !!-----------------------------------------------------------------------------
       !
@@ -489,4 +528,8 @@
       IF(lwm) WRITE ( numoni, namiceini )
 
+      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
+      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
+      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
+
       ! Define the initial parameters
       ! -------------------------
@@ -497,4 +540,5 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
+         WRITE(numout,*) '   ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file
          WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
          WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
@@ -510,4 +554,24 @@
       ENDIF
 
+      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
+         !
+         ! set si structure
+         ALLOCATE( si(jpfldi), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN
+         ENDIF
+
+         DO ifpr = 1, jpfldi
+            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
+            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
+         END DO
+
+         ! fill si with slf_i and control print
+         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
+
+         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
+
+      ENDIF
+
    END SUBROUTINE lim_istate_init
 
