Index: NEMO/releases/release-4.0/src/ICE/ice1d.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/ice1d.F90	(revision 10783)
+++ NEMO/releases/release-4.0/src/ICE/ice1d.F90	(revision 10785)
@@ -123,4 +123,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sv_i_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   o_i_1d        !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !:
@@ -209,5 +210,5 @@
          &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) ,                   &
          &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   &
-         &      sv_i_1d (jpij) , oa_i_1d (jpij) , STAT=ierr(ii) )
+         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
       !
       ii = ii + 1
Index: NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90	(revision 10783)
+++ NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90	(revision 10785)
@@ -23,4 +23,5 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
@@ -160,30 +161,43 @@
          
          zamsk = 1._wp
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area
          zamsk = 0._wp
          !
          zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                ! Ice volume
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                !-- Ice volume
          !
          zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                ! Snw volume
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                !-- Snw volume
          !
          zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               ! Salt content
-         !
-         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )               ! Age content
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               !-- Salt content
          !
          DO jk = 1, nlay_i
             zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) )   ! Ice heat content
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) )   !-- Ice heat content
          END DO
          !
          DO jk = 1, nlay_s
             zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) )   ! Snw heat content
-         END DO
-         !
-         IF ( ln_pnd_H12 ) THEN
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) )   !-- Snw heat content
+         END DO
+         !
+         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN                                                              !-- Ice Age
+            ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with
+            !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that
+            !       spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration
+            !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
+            !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )
+            ! set u*a=u for advection of ice age
+            DO jl = 1, jpl
+               zua_ho(:,:,jl) = zudy(:,:)
+               zva_ho(:,:,jl) = zvdx(:,:)
+            END DO
+            zamsk = 1._wp
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i )
+            zamsk = 0._wp
+         ENDIF
+         !
+         IF ( ln_pnd_H12 ) THEN                                                                                               !-- melt ponds
             ! set u*a=u for advection of Ap only 
             DO jl = 1, jpl
@@ -191,11 +205,11 @@
                zva_ho(:,:,jl) = zvdx(:,:)
             END DO
-            
+            !
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction
             zamsk = 0._wp
             !
             zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! mp volume
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! volume
          ENDIF
          !
@@ -203,5 +217,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                   ! Open water area
+               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                              !-- Open water area
                   &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
             END DO
Index: NEMO/releases/release-4.0/src/ICE/icethd.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icethd.F90	(revision 10783)
+++ NEMO/releases/release-4.0/src/ICE/icethd.F90	(revision 10785)
@@ -244,6 +244,4 @@
          !
       END DO
-      ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting)
-      oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
 
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
@@ -418,5 +416,9 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
          CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
-
+         !
+         ! to update ice age
+         CALL tab_2d_1d( npti, nptidx(1:npti), o_i_1d (1:npti), o_i (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
+         !
          ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
          DO jk = 1, nlay_i
@@ -443,4 +445,5 @@
          sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
          v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
+         oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
          
          CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
@@ -516,4 +519,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
          !
       END SELECT
Index: NEMO/releases/release-4.0/src/ICE/icethd_dh.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icethd_dh.F90	(revision 10783)
+++ NEMO/releases/release-4.0/src/ICE/icethd_dh.F90	(revision 10785)
@@ -614,6 +614,12 @@
       DO jk = 1, nlay_s
          DO ji = 1,npti
-            ! mask enthalpy
-            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - h_s_1d(ji) )  )
+            ! where there is no ice or no snow
+            rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) )
+            ! mass & energy loss to the ocean
+            hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * &
+               &                              ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! heat flux to the ocean [W.m-2], < 0
+            wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * &
+               &                              ( rhos          * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! mass flux
+            ! update energy (mass is updated in the next loop)
             e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
             ! recalculate t_s_1d from e_s_1d
@@ -622,6 +628,9 @@
       END DO
 
-      ! --- ensure that a_i = 0 where h_i = 0 ---
-      WHERE( h_i_1d(1:npti) == 0._wp )   a_i_1d(1:npti) = 0._wp
+      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
+      WHERE( h_i_1d(1:npti) == 0._wp )   
+         a_i_1d(1:npti) = 0._wp
+         h_s_1d(1:npti) = 0._wp
+      END WHERE
       !
    END SUBROUTINE ice_thd_dh
Index: NEMO/releases/release-4.0/src/ICE/icewri.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icewri.F90	(revision 10783)
+++ NEMO/releases/release-4.0/src/ICE/icewri.F90	(revision 10785)
@@ -100,5 +100,5 @@
       IF( iom_use('snwthic'  ) )   CALL iom_put( "snwthic", hm_s  * zmsk00        )   ! snw thickness
       IF( iom_use('icebrv'   ) )   CALL iom_put( "icebrv" , bvm_i * zmsk00 * 100. )   ! brine volume
-      IF( iom_use('iceage'   ) )   CALL iom_put( "iceage" , om_i  * zmsk00 / rday )   ! ice age
+      IF( iom_use('iceage'   ) )   CALL iom_put( "iceage" , om_i  * zmsk15 / rday )   ! ice age
       IF( iom_use('icehnew'  ) )   CALL iom_put( "icehnew", ht_i_new              )   ! new ice thickness formed in the leads
       IF( iom_use('snwvolu'  ) )   CALL iom_put( "snwvolu", vt_s  * zmsksn        )   ! snow volume
