Index: /CONFIG/UNIFORM/v6/IPSLCM6/AA_make
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLCM6/AA_make	(revision 2229)
+++ /CONFIG/UNIFORM/v6/IPSLCM6/AA_make	(revision 2230)
@@ -76,8 +76,8 @@
 
 orca2lim3:  ../../modeles/NEMOGCM/CONFIG
-	(cd ../../modeles/NEMOGCM/CONFIG; cp ../../../config/IPSLCM6/SOURCES/NEMO/arch-${FCM_ARCH}.fcm ../ARCH/. ; cp ../../../config/IPSLCM6/SOURCES/NEMO/*.cfg ../TOOLS/COMPILE/. ; ./makenemo -m ${FCM_ARCH} -n ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC TOP_SRC" -e "$(MYSRC_DIR)/SOURCES/NEMO" -j16 add_key "key_lim3 key_coupled key_oasis3 key_oasis_mct key_top key_pisces key_diahth key_diaar5" del_key "key_lim2 key_lim2_vp key_nosignedzero" ; cp ORCA2_LIM3_PISCES/BLD/bin/nemo.exe ../../../bin/opa ; )
+	(cd ../../modeles/NEMOGCM/CONFIG; cp ../../../config/IPSLCM6/SOURCES/NEMO/arch-${FCM_ARCH}.fcm ../ARCH/. ; cp ../../../config/IPSLCM6/SOURCES/NEMO/*.cfg ../TOOLS/COMPILE/. ; ./makenemo -m ${FCM_ARCH} -n ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC TOP_SRC" -e "$(MYSRC_DIR)/SOURCES/NEMO" -j16 add_key "key_lim3 key_coupled key_oasis3 key_oasis_mct key_top key_pisces key_diahth key_diaar5" del_key "key_lim2 key_nosignedzero" ; cp ORCA2_LIM3_PISCES/BLD/bin/nemo.exe ../../../bin/opa ; )
 
 orca1lim3: ../../modeles/NEMOGCM/CONFIG
-	(cd ../../modeles/NEMOGCM/CONFIG; cp ../../../config/IPSLCM6/SOURCES/NEMO/arch-${FCM_ARCH}.fcm ../ARCH/. ; cp ../../../config/IPSLCM6/SOURCES/NEMO/*.cfg ../TOOLS/COMPILE/. ; ./makenemo -m ${FCM_ARCH} -n ORCA1_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC TOP_SRC" -e "$(MYSRC_DIR)/SOURCES/NEMO" -j16 add_key "key_vvl key_lim3 key_oasis3 key_coupled key_oasis_mct  key_top key_pisces key_diahth key_diaar5" del_key "key_lim2 key_lim2_vp key_nosignedzero" ; cp ORCA1_LIM3_PISCES/BLD/bin/nemo.exe ../../../bin/opa ; )
+	(cd ../../modeles/NEMOGCM/CONFIG; cp ../../../config/IPSLCM6/SOURCES/NEMO/arch-${FCM_ARCH}.fcm ../ARCH/. ; cp ../../../config/IPSLCM6/SOURCES/NEMO/*.cfg ../TOOLS/COMPILE/. ; ./makenemo -m ${FCM_ARCH} -n ORCA1_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC TOP_SRC" -e "$(MYSRC_DIR)/SOURCES/NEMO" -j16 add_key "key_lim3 key_oasis3 key_coupled key_oasis_mct  key_top key_pisces key_diahth key_diaar5" del_key "key_lim2 key_lim2_vp key_nosignedzero" ; cp ORCA1_LIM3_PISCES/BLD/bin/nemo.exe ../../../bin/opa ; )
 
 verif: ../../bin
Index: /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90	(revision 2230)
+++ /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90	(revision 2230)
@@ -0,0 +1,538 @@
+MODULE limistate
+   !!======================================================================
+   !!                     ***  MODULE  limistate  ***
+   !!              Initialisation of diagnostics ice variables
+   !!======================================================================
+   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!             -   ! 2012    (C. Rousset) add par_oce (for jp_sal)...bug?
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                    LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_istate      :  Initialisation of diagnostics ice variables
+   !!   lim_istate_init :  initialization of ice state and namelist read
+   !!----------------------------------------------------------------------
+   USE phycst           ! physical constant
+   USE oce              ! dynamics and tracers variables
+   USE dom_oce          ! ocean domain
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice fields
+   USE eosbn2           ! equation of state
+   USE ice              ! sea-ice variables
+   USE par_ice          ! ice parameters
+   USE par_oce          ! ocean parameters
+   USE dom_ice          ! sea-ice domain
+   USE in_out_manager   ! I/O manager
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   USE lib_mpp          ! MPP library
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE wrk_nemo         ! work arrays
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   lim_istate      ! routine called by lim_init.F90
+
+   !! * Module variables
+   !                          !!** init namelist (namiceini) **
+   REAL(wp) ::   ttest   ! threshold water temperature for initial sea ice
+   REAL(wp) ::   hninn   ! initial snow thickness in the north
+   REAL(wp) ::   hnins   ! initial snow thickness in the south
+   REAL(wp) ::   hginn   ! initial ice thickness in the north
+   REAL(wp) ::   hgins   ! initial ice thickness in the south
+   REAL(wp) ::   aginn   ! initial leads area in the north
+   REAL(wp) ::   agins   ! initial leads area in the south
+   REAL(wp) ::   sinn    ! initial salinity 
+   REAL(wp) ::   sins  
+
+   !!----------------------------------------------------------------------
+   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
+   !! $Id: limistate.F90 4335 2013-12-17 14:40:19Z clem $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE lim_istate
+      !!-------------------------------------------------------------------
+      !!                    ***  ROUTINE lim_istate  ***
+      !!
+      !! ** Purpose :   defined the sea-ice initial state
+      !!
+      !! ** Method  :   
+      !!                This routine will put some ice where ocean
+      !!                is at the freezing point, then fill in ice 
+      !!                state variables using prescribed initial 
+      !!                values in the namelist            
+      !!
+      !! ** Steps   :   
+      !!                1) Read namelist
+      !!                2) Basal temperature; ice and hemisphere masks
+      !!                3) Fill in the ice thickness distribution using gaussian
+      !!                4) Fill in space-dependent arrays for state variables
+      !!                5) Diagnostic arrays
+      !!                6) Lateral boundary conditions
+      !!
+      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
+      !!              where there is no ice (clem: I do not know why but it is mandatory) 
+      !!
+      !! History :
+      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
+      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
+      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
+      !!--------------------------------------------------------------------
+
+      !! * Local variables
+      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
+      REAL(wp)   :: epsi20, 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(:)     :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, za_i_ini, zv_i_ini
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zidto    ! ice indicator
+      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
+      !--------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, zidto )
+      CALL wrk_alloc( jpi, jpj, zhemis )
+      CALL wrk_alloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini )
+      CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
+
+      epsi20   = 1.0e-20
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
+
+      !--------------------------------------------------------------------
+      ! 1) Read namelist
+      !--------------------------------------------------------------------
+
+      CALL lim_istate_init     !  reading the initials parameters of the ice
+
+!!gm  in lim2  the initialisation if only done if required in the namelist :
+!!gm      IF( .NOT. ln_limini ) THEN
+!!gm  this should be added in lim3 namelist...
+
+      !--------------------------------------------------------------------
+      ! 2) Basal temperature, ice mask and hemispheric index
+      !--------------------------------------------------------------------
+
+      ! Basal temperature is set to the freezing point of seawater in Celsius
+      t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
+
+      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
+         DO ji = 1, jpi
+            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0._wp      ! no ice
+            ELSE                                                     ;   zidto(ji,jj) = 1._wp      !    ice
+            ENDIF
+         END DO
+      END DO
+
+      t_bo(:,:) = t_bo(:,:) + rt0                          ! conversion to Kelvin
+
+      ! Hemispheric index
+      ! MV 2011 new initialization
+      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
+      ! END MV 2011 new initialization
+
+      !--------------------------------------------------------------------
+      ! 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
+      zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins  ! ice thickness
+      zat_i_ini(1) = aginn ; zat_i_ini(2) = agins  ! ice concentration
+      zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:)   ! ice volume
+      zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins  ! snow depth
+      zsm_i_ini(1) = sinn  ; zsm_i_ini(2) = sins   ! bulk ice salinity
+
+      !---------------------------------------------------------------------
+      ! 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
+               zht_i_ini(1,i_hemis)       = zhm_i_ini(i_hemis)
+               za_i_ini(1,i_hemis)        = zat_i_ini(i_hemis)
+               zht_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
+                  zht_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) )
+               END DO
+
+            !--- jl0: most likely index where cc will be maximum
+               DO jl = 1, jpl
+                  IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. &
+                       ( zhm_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN
+                     jl0 = jl
+                  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 * zhm_i_ini(i_hemis)
+                     zarg                 = ( zht_i_ini(jl,i_hemis) - zhm_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)*zht_i_ini(jl,i_hemis)
+               END DO
+               zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
+               IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
+
+            !--- volumes
+               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_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 ( zht_i_ini(i_fill, i_hemis) .GT. 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 *** zht_i_ini(i_fill,i_hemis) = ', &
+               zht_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,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+            WRITE(numout,*), ' !!!! RED ALERT                  !!! '
+            WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!'
+            WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure '
+            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+            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,*), ' zhm_i_ini : ', zhm_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
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
+               ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness
+               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) )  ! snow depth
+               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity
+               o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age
+               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! 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 / 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 ! ji
+         END DO ! jj
+      END DO ! jl
+
+      ! 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) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt
+                   ! Snow energy of melting
+                   e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
+                   ! Change dimensions
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
+                   ! Multiply by volume, so that heat content in 10^9 Joules
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
+               END DO ! ji
+            END DO ! jj
+         END DO ! jl
+      END DO ! jk
+
+      ! 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) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt 
+                   s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min
+                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
+
+                   ! heat content per unit volume
+                   e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                      +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
+                      -   rcp     * ( ztmelts - rtt ) )
+
+                   ! Correct dimensions to avoid big values
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
+
+                   ! Mutliply by ice volume, and divide by number of layers 
+                   ! to get heat content in 10^9 J
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i
+               END DO ! ji
+            END DO ! jj
+         END DO ! jl
+      END DO ! jk
+
+      !--------------------------------------------------------------------
+      ! 4) Global ice variables for output diagnostics                    | 
+      !--------------------------------------------------------------------
+      fsbbq (:,:)     = 0._wp
+      u_ice (:,:)     = 0._wp
+      v_ice (:,:)     = 0._wp
+      stress1_i(:,:)  = 0._wp
+      stress2_i(:,:)  = 0._wp
+      stress12_i(:,:) = 0._wp
+
+
+      !--------------------------------------------------------------------
+      ! 5) Moments for advection
+      !--------------------------------------------------------------------
+
+      sxopw (:,:) = 0._wp 
+      syopw (:,:) = 0._wp 
+      sxxopw(:,:) = 0._wp 
+      syyopw(:,:) = 0._wp 
+      sxyopw(:,:) = 0._wp
+
+      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
+      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
+      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
+      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
+      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
+
+      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
+      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
+      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
+      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
+      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
+
+      sxsal  (:,:,:)  = 0._wp
+      sysal  (:,:,:)  = 0._wp
+      sxxsal (:,:,:)  = 0._wp
+      syysal (:,:,:)  = 0._wp
+      sxysal (:,:,:)  = 0._wp
+
+      sxage  (:,:,:)  = 0._wp
+      syage  (:,:,:)  = 0._wp
+      sxxage (:,:,:)  = 0._wp
+      syyage (:,:,:)  = 0._wp
+      sxyage (:,:,:)  = 0._wp
+
+      !--------------------------------------------------------------------
+      ! 6) Lateral boundary conditions                                    | 
+      !--------------------------------------------------------------------
+
+      DO jl = 1, jpl
+
+         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( smv_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. )
+         CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. )
+         DO jk = 1, nlay_s
+            CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. )
+            CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. )
+         END DO
+         DO jk = 1, nlay_i
+            CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. )
+            CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
+         END DO
+         !
+         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
+      END DO
+      
+      at_i (:,:) = 0.0_wp
+      DO jl = 1, jpl
+         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
+      END DO
+
+      CALL lbc_lnk( at_i , 'T', 1. )
+      at_i(:,:) = tms(:,:) * at_i(:,:)                       ! put 0 over land
+      !
+      CALL lbc_lnk( fsbbq  , 'T', 1. )
+      !
+      !--------------------------------------------------------------------
+      ! 6) ????                                                           | 
+      !--------------------------------------------------------------------
+      tn_ice (:,:,:) = t_su (:,:,:)
+
+      CALL wrk_dealloc( jpi, jpj, zidto )
+      CALL wrk_dealloc( jpi, jpj, zhemis )
+      CALL wrk_dealloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini )
+      CALL wrk_dealloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
+
+   END SUBROUTINE lim_istate
+
+   SUBROUTINE lim_istate_init
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_istate_init  ***
+      !!        
+      !! ** Purpose : Definition of initial state of the ice 
+      !!
+      !! ** Method : Read the namiceini namelist and check the parameter 
+      !!       values called at the first timestep (nit000)
+      !!
+      !! ** input : 
+      !!        Namelist namiceini
+      !!
+      !! history :
+      !!  8.5  ! 03-08 (C. Ethe) original code 
+      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
+      !!-----------------------------------------------------------------------------
+      NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins
+      !
+      INTEGER :: ios                 ! Local integer output status for namelist read
+      !!-----------------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
+      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
+      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
+      WRITE ( numoni, namiceini )
+
+      ! Define the initial parameters
+      ! -------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest
+         WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn
+         WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins 
+         WRITE(numout,*) '   initial ice thickness  in the north          hginn      = ', hginn
+         WRITE(numout,*) '   initial ice thickness  in the south          hgins      = ', hgins
+         WRITE(numout,*) '   initial ice concentr.  in the north          aginn      = ', aginn
+         WRITE(numout,*) '   initial ice concentr.  in the north          agins      = ', agins
+         WRITE(numout,*) '   initial  ice salinity  in the north          sinn       = ', sinn
+         WRITE(numout,*) '   initial  ice salinity  in the south          sins       = ', sins
+      ENDIF
+
+   END SUBROUTINE lim_istate_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_istate          ! Empty routine
+   END SUBROUTINE lim_istate
+#endif
+
+   !!======================================================================
+END MODULE limistate
