diff --git a/src/MOD_ICE.F90 b/src/MOD_ICE.F90
index 0150e5004..0375d37d1 100644
--- a/src/MOD_ICE.F90
+++ b/src/MOD_ICE.F90
@@ -129,6 +129,8 @@ MODULE MOD_ICE
 
     ! maEVP variables
     real(kind=WP), allocatable, dimension(:)    :: alpha_evp_array, beta_evp_array
+    ! ice/snow thicknesses in the ice-covered area
+    real(kind=WP), allocatable, dimension(:)    :: h_ice, h_snow    
 
     !___________________________________________________________________________
     ! total number of ice tracers (default=3, 1=area, 2=mice, 3=msnow, (4=ice_temp)
@@ -624,6 +626,8 @@ subroutine ice_init(ice, partit, mesh)
     allocate(ice%stress_iceoce_x(      node_size))
     allocate(ice%stress_atmice_y(      node_size))
     allocate(ice%stress_iceoce_y(      node_size))
+    allocate(ice%h_ice          (      node_size))
+    allocate(ice%h_snow         (      node_size))    
     ice%uice            = 0.0_WP
     ice%uice_rhs        = 0.0_WP
     ice%uice_old        = 0.0_WP
@@ -634,6 +638,8 @@ subroutine ice_init(ice, partit, mesh)
     ice%vice_old        = 0.0_WP
     ice%stress_atmice_y = 0.0_WP
     ice%stress_iceoce_y = 0.0_WP
+    ice%h_ice           = 0.0_WP
+    ice%h_snow          = 0.0_WP    
     if (ice%whichEVP /= 0) then
         allocate(ice%uice_aux(         node_size))
         allocate(ice%vice_aux(         node_size))
diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90
index 2b3aa049c..ea6dcf1d9 100755
--- a/src/ice_setup_step.F90
+++ b/src/ice_setup_step.F90
@@ -125,7 +125,7 @@ subroutine ice_timestep(step, ice, partit, mesh)
     !LA 2023-03-08
     real(kind=WP), dimension(:), pointer  :: u_ice_ib, v_ice_ib
 #if defined (__oifs) || defined (__ifsinterface)
-    real(kind=WP), dimension(:), pointer  :: ice_temp, a_ice
+    real(kind=WP), dimension(:), pointer  :: a_ice, ice_temp
     !LA 2023-03-08
     real(kind=WP), dimension(:), pointer  :: a_ice_ib
 #endif
@@ -170,9 +170,8 @@ subroutine ice_timestep(step, ice, partit, mesh)
 !$omp end parallel sections
   end if
 !---------------------------------------------
-
 #if defined (__oifs) || defined (__ifsinterface)
-    a_ice    => ice%data(1)%values(:)
+    a_ice    => ice%data(1)%values(:)    
     ice_temp => ice%data(4)%values(:)
 #endif
     !___________________________________________________________________________
@@ -300,6 +299,8 @@ subroutine ice_timestep(step, ice, partit, mesh)
     !___________________________________________________________________________
 !$OMP PARALLEL DO
     do i=1,myDim_nod2D+eDim_nod2D
+        ice%h_ice(i) =ice%data(2)%values(i)/max(ice%data(1)%values(i), 1.e-3)
+        ice%h_snow(i)=ice%data(3)%values(i)/max(ice%data(1)%values(i), 1.e-3)
         if ( ( U_ice(i)/=0.0_WP .and. mesh%ulevels_nod2d(i)>1) .or. (V_ice(i)/=0.0_WP .and. mesh%ulevels_nod2d(i)>1) ) then
             write(*,*) " --> found cavity velocity /= 0.0_WP , ", mype
             write(*,*) " ulevels_nod2d(n) = ", mesh%ulevels_nod2d(i)
diff --git a/src/io_meandata.F90 b/src/io_meandata.F90
index 655046628..c5f958320 100644
--- a/src/io_meandata.F90
+++ b/src/io_meandata.F90
@@ -241,7 +241,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh)
     end if
 CASE ('m_ice     ')
     if (use_ice) then
-    call def_stream(nod2D, myDim_nod2D, 'm_ice',    'ice height',                     'm',      ice%data(2)%values(1:myDim_nod2D),      io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+    call def_stream(nod2D, myDim_nod2D, 'm_ice',    'ice height per unit area',       'm',      ice%data(2)%values(1:myDim_nod2D),      io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
     end if
 CASE ('thdgr     ')
     if (use_ice) then
@@ -257,9 +257,17 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh)
     end if
 CASE ('m_snow    ')
     if (use_ice) then
-    call def_stream(nod2D, myDim_nod2D, 'm_snow',   'snow height',                     'm',      ice%data(3)%values(1:myDim_nod2D),     io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+    call def_stream(nod2D, myDim_nod2D, 'm_snow',   'snow height  per unit area',      'm',      ice%data(3)%values(1:myDim_nod2D),     io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
     end if
-
+CASE ('h_ice     ')
+    if (use_ice) then
+    call def_stream(nod2D, myDim_nod2D, 'h_ice',    'ice thickness over ice-covered fraction',   'm',     ice%h_ice(1:myDim_nod2D), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+    end if
+CASE ('h_snow    ')
+    if (use_ice) then
+    call def_stream(nod2D, myDim_nod2D, 'h_snow',   'snow thickness over ice-covered fraction',  'm',     ice%h_snow(1:myDim_nod2D), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+    end if    
+    
 ! Debug ice variables    
 CASE ('strength_ice')
     if (use_ice) then