Index: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/ASM/asminc.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/ASM/asminc.F90	(revision 15096)
+++ NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/ASM/asminc.F90	(revision 15486)
@@ -35,6 +35,9 @@
    USE sbc_oce         ! Surface boundary condition variables.
    USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step
-#if defined key_si3
-   USE ice      , ONLY : hm_i, at_i, at_i_b
+#if defined key_si3 && defined key_asminc
+   USE phycst         ! physical constants
+   USE ice1D          ! sea-ice: thermodynamics variables
+   USE icetab         ! sea-ice: 3D <==> 2D, 2D <==> 1D
+   USE ice
 #endif
    !
@@ -89,4 +92,9 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
+   REAL(wp) :: zhi_damin                                            ! Ice thickness for new sea ice from DA increment
+#if defined key_si3 && defined key_asminc
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   a_i_bkginc          ! Increment to the background sea ice conc categories
+   LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice    ! Mask .TRUE.=DA positive ice increment to open water
+#endif
 #if defined key_cice && defined key_asminc
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE
@@ -112,9 +120,9 @@
       !! ** Action  : 
       !!----------------------------------------------------------------------
-      INTEGER :: ji, jj, jk, jt  ! dummy loop indices
-      INTEGER :: imid, inum      ! local integers
-      INTEGER :: ios             ! Local integer output status for namelist read
-      INTEGER :: iiauper         ! Number of time steps in the IAU period
-      INTEGER :: icycper         ! Number of time steps in the cycle
+      INTEGER :: ji, jj, jk, jt, jl  ! dummy loop indices
+      INTEGER :: imid, inum          ! local integers
+      INTEGER :: ios                 ! Local integer output status for namelist read
+      INTEGER :: iiauper             ! Number of time steps in the IAU period
+      INTEGER :: icycper             ! Number of time steps in the cycle
       REAL(KIND=dp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step
       REAL(KIND=dp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term
@@ -131,10 +139,12 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
+      REAL(wp) :: zremaining_increment
+
       !!
       NAMELIST/nam_asminc/ ln_bkgwri,                                      &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
-         &                 ln_asmdin, ln_asmiau,                           &
+         &                 ln_seaiceinc, ln_asmdin, ln_asmiau,             &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
-         &                 ln_salfix, salfixmin, nn_divdmp
+         &                 ln_salfix, salfixmin, nn_divdmp, zhi_damin
       !!----------------------------------------------------------------------
 
@@ -173,4 +183,5 @@
          WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
          WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
+         WRITE(numout,*) '      Minimum ice thickness for new ice from DA                zhi_damin = ', zhi_damin
       ENDIF
 
@@ -325,4 +336,8 @@
       ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp
       ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp
+#if defined key_si3 && defined key_asminc
+      ALLOCATE( a_i_bkginc   (jpi,jpj,jpl) )   ;   a_i_bkginc   (:,:,:) = 0._wp
+      ALLOCATE( incr_newice(jpi,jpj,jpl) )     ;   incr_newice(:,:,:) = .FALSE.
+#endif
 #if defined key_asminc
       ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp
@@ -397,6 +412,62 @@
             ! Set missing increments to 0.0 rather than 1e+20
             ! to allow for differences in masks
-            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
-         ENDIF
+            ! very small increments are also set to zero
+            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 .OR. &
+              &    ABS( seaice_bkginc(:,:) ) < 1.0e-03_wp ) seaice_bkginc(:,:) = 0.0_wp
+            
+#if defined key_si3 && defined key_asminc
+            IF (lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'asm_inc_init : Converting single category increment to multi-category'
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            END IF
+
+            !single category increment for sea ice conc
+            !convert single category increment to multi category
+            at_i = SUM( a_i(:,:,:), DIM=3 )
+
+            ! ensure zhi_damin lies within 1st category
+            IF ( zhi_damin > hi_max(1) ) THEN 
+               IF (lwp) THEN
+                  WRITE(numout,*)
+                  WRITE(numout,*) 'minimum DA thickness > hmax(1): setting minimum DA thickness to hmax(1)'
+                  WRITE(numout,*) '~~~~~~~~~~~~'
+               END IF
+               zhi_damin = hi_max(1) - epsi02
+            END IF
+            
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( seaice_bkginc(ji,jj) > 0.0_wp ) THEN
+                     !positive ice concentration increments are always 
+                     !added to the thinnest category of ice
+                     a_i_bkginc(ji,jj,1) = seaice_bkginc(ji,jj)
+                  ELSE
+                     !negative increments are first removed from the thinnest
+                     !available category until it reaches zero concentration
+                     !and then progressively removed from thicker categories
+
+                     !NOTE: might consider the remote possibility that zremaining_increment 
+                     !left over is not zero, particulary if SI3 is being run 
+                     !as a single category 
+                     zremaining_increment = seaice_bkginc(ji,jj)
+                     DO jl = 1, jpl
+                        ! assign as much increment as possible to current category
+                        a_i_bkginc(ji,jj,jl) = -MIN( a_i(ji,jj,jl), -zremaining_increment )
+                        ! update remaining amount of increment
+                        zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl)
+                     END DO
+                  END IF
+               END DO
+            END DO
+            ! find model points where DA new ice should be added to open water
+            ! in any ice category
+            DO jl = 1,jpl
+               WHERE ( at_i(:,:) < epsi02 .AND. seaice_bkginc(:,:) > 0.0_wp )
+                  incr_newice(:,:,jl) = .TRUE.
+               END WHERE
+            END DO
+         ENDIF
+#endif
          !
          CALL iom_close( inum )
@@ -823,9 +894,12 @@
       INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
       !
-      INTEGER  ::   it
-      REAL(wp) ::   zincwgt   ! IAU weight for current time step
-#if defined key_si3
-      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc
-      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
+      INTEGER  ::   it, jk
+      REAL(wp) ::   zincwgt                            ! IAU weight for current time step
+#if defined key_si3 && defined key_asminc
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i         ! change in ice concentration
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: dv_i         ! change in ice volume
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i       ! inverse of ice concentration before current IAU step
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_v_i       ! inverse of ice volume before current IAU step
+      REAL(wp), DIMENSION(jpi,jpj)     :: zhi_damin_2D ! array with DA thickness for incr_newice
 #endif
       !!----------------------------------------------------------------------
@@ -849,27 +923,94 @@
             ! Sea-ice : SI3 case
             !
-#if defined key_si3
-            zofrld (:,:) = 1._wp - at_i(:,:)
-            zohicif(:,:) = hm_i(:,:)
-            !
-            at_i  (:,:) = 1. - MIN( MAX( 1.-at_i  (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
-            at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
-            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
-            !
-            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
-            !
-            ! Nudge sea ice depth to bring it up to a required minimum depth
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt    
+#if defined key_si3 && defined key_asminc
+            ! compute the inverse of key sea ice variables
+            ! to be used later in the code 
+            WHERE( a_i(:,:,:) > epsi10 )
+               z1_a_i(:,:,:) = 1.0_wp / a_i(:,:,:)
+               z1_v_i(:,:,:) = 1.0_wp / v_i(:,:,:)
             ELSEWHERE
-               zhicifinc(:,:) = 0.0_wp
+               z1_a_i(:,:,:) = 0.0_wp
+               z1_v_i(:,:,:) = 0.0_wp
             END WHERE
-            !
-            ! nudge ice depth
-            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
-            !
-            ! seaice salinity balancing (to add)
-#endif
-            !
+
+            ! add positive concentration increments to regions where ice
+            ! is already present and bound them to 1
+            ! ice volume is added based on zhi_damin
+            WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) > 0.0_wp )
+               a_i(:,:,:) = a_i(:,:,:) + MIN( 1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt )
+               v_i(:,:,:) = v_i(:,:,:) + MIN( 1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin
+            END WHERE
+
+            ! add negative concentration increments to regions where ice
+            ! is already present and bound them to 0
+            ! in this case ice volume is changed based on the current thickness
+            WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) < 0.0_wp )
+               a_i(:,:,:) = MAX( a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp )
+               v_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:)
+            END WHERE
+           
+            ! compute changes in ice concentration and volume
+            WHERE ( incr_newice )
+               da_i(:,:,:) = 1.0_wp
+               dv_i(:,:,:) = 1.0_wp
+            ELSEWHERE
+               da_i(:,:,:) = a_i(:,:,:) * z1_a_i(:,:,:)
+               dv_i(:,:,:) = v_i(:,:,:) * z1_v_i(:,:,:)
+            END WHERE
+
+            ! initialise thermodynamics of new ice being added to open water
+            ! just do it once since next IAU steps assume that new ice has 
+            ! already been added in
+            IF ( kt == nitiaustr_r ) THEN
+
+               ! assign zhi_damin to ice forming in open water
+               WHERE ( ANY( incr_newice, DIM=3 ) )
+                  zhi_damin_2D(:,:) = zhi_damin
+               ELSEWHERE
+                  zhi_damin_2D(:,:) = 0.0_wp
+               END WHERE
+
+               ! add ice concentration and volume
+               ! ensure the other prognostic variables are set to zero 
+               WHERE ( incr_newice )
+                  a_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt )
+                  v_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin
+                  v_s (:,:,:) = 0.0_wp
+                  a_ip(:,:,:) = 0.0_wp
+                  v_ip(:,:,:) = 0.0_wp
+                  sv_i(:,:,:) = 0.0_wp
+               END WHERE
+               DO jk = 1, nlay_i
+                  WHERE ( incr_newice )
+                     e_i(:,:,jk,:) = 0.0_wp
+                  END WHERE
+               END DO
+               DO jk = 1, nlay_s
+                  WHERE ( incr_newice )
+                     e_s(:,:,jk,:) = 0.0_wp
+                  END WHERE
+               END DO
+
+               ! Initialisation of the salt content and ice enthalpy
+               ! set flag of new ice to false after this
+               CALL init_new_ice_thd( zhi_damin_2D )
+               incr_newice(:,:,:) = .FALSE.
+            END IF
+
+            ! maintain equivalent values for key prognostic variables
+            v_s(:,:,:) = v_s(:,:,:) * da_i(:,:,:)
+            DO jk = 1, nlay_s
+               e_s(:,:,jk,:) = e_s(:,:,jk,:) * da_i(:,:,:)
+            END DO
+            a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:)
+            v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:)
+ 
+            ! ice volume dependent variables
+            sv_i (:,:,:) = sv_i(:,:,:) * dv_i(:,:,:)
+            DO jk = 1, nlay_i
+               e_i(:,:,jk,:) = e_i(:,:,jk,:) * dv_i(:,:,:)
+            END DO
+#endif
+
 #if defined key_cice && defined key_asminc
             ! Sea-ice : CICE case. Pass ice increment tendency into CICE
@@ -879,4 +1020,8 @@
             IF ( kt == nitiaufin_r ) THEN
                DEALLOCATE( seaice_bkginc )
+#if defined key_si3 && defined key_asminc
+               DEALLOCATE( incr_newice )
+               DEALLOCATE( a_i_bkginc )
+#endif
             ENDIF
             !
@@ -895,30 +1040,10 @@
             !
             neuler = 0                    ! Force Euler forward step
-            !
-            ! Sea-ice : SI3 case
-            !
-#if defined key_si3
-            zofrld (:,:) = 1._wp - at_i(:,:)
-            zohicif(:,:) = hm_i(:,:)
-            ! 
-            ! Initialize the now fields the background + increment
-            at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
-            at_i_b(:,:) = at_i(:,:) 
-            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
-            !
-            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
-            !
-            ! Nudge sea ice depth to bring it up to a required minimum depth
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = zhicifmin - hm_i(:,:)
-            ELSEWHERE
-               zhicifinc(:,:) = 0.0_wp
-            END WHERE
-            !
-            ! nudge ice depth
-            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
-            !
-            ! seaice salinity balancing (to add)
-#endif
+
+            IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'seaice_asm_inc : sea ice direct initialization at time step = ', kt
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
             !
 #if defined key_cice && defined key_asminc
@@ -937,84 +1062,95 @@
             !
          ENDIF
-
-!#if defined defined key_si3 || defined key_cice
-!
-!            IF (ln_seaicebal ) THEN       
-!             !! balancing salinity increments
-!             !! simple case from limflx.F90 (doesn't include a mass flux)
-!             !! assumption is that as ice concentration is reduced or increased
-!             !! the snow and ice depths remain constant
-!             !! note that snow is being created where ice concentration is being increased
-!             !! - could be more sophisticated and
-!             !! not do this (but would need to alter h_snow)
-!
-!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
-!
-!             DO jj = 1, jpj
-!               DO ji = 1, jpi 
-!           ! calculate change in ice and snow mass per unit area
-!           ! positive values imply adding salt to the ocean (results from ice formation)
-!           ! fwf : ice formation and melting
-!
-!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
-!
-!           ! change salinity down to mixed layer depth
-!                 mld=hmld_kara(ji,jj)
-!
-!           ! prevent small mld
-!           ! less than 10m can cause salinity instability 
-!                 IF (mld < 10) mld=10
-!
-!           ! set to bottom of a level 
-!                 DO jk = jpk-1, 2, -1
-!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 
-!                     mld=gdepw(ji,jj,jk+1)
-!                     jkmax=jk
-!                   ENDIF
-!                 ENDDO
-!
-!            ! avoid applying salinity balancing in shallow water or on land
-!            ! 
-!
-!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
-!
-!                 dsal_ocn=0.0_wp
-!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
-!
-!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
-!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
-!
-!           ! put increments in for levels in the mixed layer
-!           ! but prevent salinity below a threshold value 
-!
-!                   DO jk = 1, jkmax              
-!
-!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 
-!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
-!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
-!                     ENDIF
-!
-!                   ENDDO
-!
-!      !            !  salt exchanges at the ice/ocean interface
-!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
-!      !
-!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
-!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt  
-!      !!               
-!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d) 
-!      !!                                                     ! E-P (kg m-2 s-2)
-!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
-!               ENDDO !ji
-!             ENDDO !jj!
-!
-!            ENDIF !ln_seaicebal
-!
-!#endif
          !
       ENDIF
       !
    END SUBROUTINE seaice_asm_inc
-   
+
+
+   SUBROUTINE init_new_ice_thd( hi_new )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE init_new_ice_thd  ***
+      !!
+      !! ** Purpose :   Initialise thermodynamics of new ice 
+      !!                forming at 1st category with thickness hi_new
+      !!
+      !! ** Method  :   Apply SI3  equations to initialise
+      !!                thermodynamics of new ice
+      !!
+      !! ** Action  :   update sea ice thermodynamics
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    ::   hi_new  ! total thickness of new ice
+
+      INTEGER :: jj, ji, jk
+      REAL(wp) ::   ztmelts         ! melting point
+      REAL(wp) ::   Sice_Fz=2.3_wp  ! Salinity of the ice = F(z) [multiyear ice]
+
+      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! 1d version of hi_new
+      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of new ice
+      !!----------------------------------------------------------------------
+
+      ! Identify grid points where new ice forms
+      npti = 0   ;   nptidx(:) = 0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( hi_new(ji,jj) > 0._wp ) THEN
+               npti = npti + 1
+               nptidx( npti ) = (jj - 1) * jpi + ji
+            ENDIF
+         END DO
+      END DO
+
+      ! Move from 2-D to 1-D vectors
+      IF ( npti > 0 ) THEN
+         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , hi_new        )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m         )
+         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo          )
+         DO jk = 1, nlay_i
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) )
+         END DO
+
+         ! --- Salinity of new ice --- ! 
+         SELECT CASE ( nn_icesal )
+         CASE ( 1 )                    ! Sice = constant 
+            zs_newice(1:npti) = rn_icesal
+         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
+            DO ji = 1, npti
+               zs_newice(ji) = MIN(  4.606_wp + 0.91_wp / zh_newice(ji) , rn_simax , 0.5_wp * sss_1d(ji) )
+            END DO
+         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
+            zs_newice(1:npti) = Sice_Fz
+         END SELECT
+
+         ! --- Update ice salt content --- !
+         DO ji = 1, npti
+            sv_i_2d(ji,1) = sv_i_2d(ji,1) + zs_newice(ji) * ( v_i_2d(ji,1) )
+         END DO
+
+         ! --- Heat content of new ice --- !
+         ! We assume that new ice is formed at the seawater freezing point
+         DO ji = 1, npti
+               ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
+               e_i_1d(ji,:)  =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                        &
+                  &                      + rLfus * ( 1.0_wp - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+                  &                      - rcp   *         ztmelts )
+         END DO
+
+         ! Change units for e_i
+         DO jk = 1, nlay_i
+            e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * v_i_2d(1:npti,1) * r1_nlay_i 
+         END DO
+
+         ! Reforming full thermodynamic variables
+         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
+         DO jk = 1, nlay_i
+               CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) )
+         END DO
+      END IF
+
+   END SUBROUTINE init_new_ice_thd
+
+
    !!======================================================================
 END MODULE asminc
