Index: /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/ice.F90
===================================================================
--- /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/ice.F90	(revision 13639)
+++ /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/ice.F90	(revision 13640)
@@ -391,4 +391,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice         !: ice volume variation   [m/s] 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_aice         !: ice conc.  variation   [s-1] 
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_mass     !: advection of mass (kg/m2/s)
@@ -497,7 +498,7 @@
       ! * Ice diagnostics
       ii = ii + 1
-      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),   & 
-         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),   &
-         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj),   &
+      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),                      & 
+         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),                      &
+         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj),  &
          &      diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) )
 
Index: /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icecor.F90
===================================================================
--- /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icecor.F90	(revision 13639)
+++ /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icecor.F90	(revision 13640)
@@ -55,5 +55,4 @@
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) ::   zsal, zzc
-      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag
       !!----------------------------------------------------------------------
       ! controls
@@ -123,47 +122,4 @@
          CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1._wp, v_ice, 'V', -1._wp )
       ENDIF
-
-      !                             !-----------------------------------------------------
-      SELECT CASE( kn )             !  Diagnostics                                       !
-      !                             !-----------------------------------------------------
-      CASE( 1 )                        !--- dyn trend diagnostics
-         !
-         IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
-            diag_heat(:,:) = - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice &      ! W.m-2
-               &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice
-            diag_sice(:,:) =   SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice * rhoi
-            diag_vice(:,:) =   SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhoi
-            diag_vsnw(:,:) =   SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhos
-         ENDIF
-         !                       ! concentration tendency (dynamics)
-         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
-            zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
-            CALL iom_put( 'afxdyn' , zafx )
-         ENDIF
-         !
-      CASE( 2 )                        !--- thermo trend diagnostics & ice aging
-         !
-         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation
-         !
-         IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
-            diag_heat(:,:) = diag_heat(:,:) &
-               &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice &
-               &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice
-            diag_sice(:,:) = diag_sice(:,:) &
-               &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice * rhoi
-            diag_vice(:,:) = diag_vice(:,:) &
-               &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhoi
-            diag_vsnw(:,:) = diag_vsnw(:,:) &
-               &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhos
-            CALL iom_put ( 'hfxdhc' , diag_heat ) 
-         ENDIF
-         !                       ! concentration tendency (total + thermo)
-         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
-            zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
-            CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice )
-            CALL iom_put( 'afxtot' , zafx )
-         ENDIF
-         !
-      END SELECT
       !
       ! controls
Index: /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icestp.F90
===================================================================
--- /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icestp.F90	(revision 13639)
+++ /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icestp.F90	(revision 13640)
@@ -55,5 +55,4 @@
    USE icedyn         ! sea-ice: dynamics
    USE icethd         ! sea-ice: thermodynamics
-   USE icecor         ! sea-ice: corrections
    USE iceupdate      ! sea-ice: sea surface boundary condition update
    USE icedia         ! sea-ice: budget diagnostics
@@ -162,4 +161,6 @@
             &                           CALL ice_dyn( kt )            ! -- Ice dynamics
          !
+                                        CALL diag_trends( 1 )         ! record dyn trends
+         !
          !                          !==  lateral boundary conditions  ==!
          IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
@@ -189,6 +190,5 @@
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
-                                        CALL ice_cor( kt , 2 )        ! -- Corrections
-         !
+                                        CALL diag_trends( 2 )         ! record thermo trends
                                         CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
                                         CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
@@ -391,55 +391,112 @@
       !!               of the time step
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj      ! dummy loop index
-      !!----------------------------------------------------------------------
-      sfx    (:,:) = 0._wp   ;
-      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
-      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
-      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
-      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
-      !
-      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
-      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
-      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
-      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
-      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
-      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
-      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
-      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
-      wfx_snw_sni(:,:) = 0._wp 
-      wfx_pnd(:,:) = 0._wp
-
-      hfx_thd(:,:) = 0._wp   ;
-      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
-      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
-      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
-      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
-      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
-      hfx_err_dif(:,:) = 0._wp
-      wfx_err_sub(:,:) = 0._wp
-      !
-      diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp
-      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
-
-      ! SIMIP diagnostics
-      qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes
-      t_si       (:,:,:) = rt0   ! temp at the ice-snow interface
-
-      tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
-      cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
-      qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T)
-      qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
-      qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F)
-      !
-      ! for control checks (ln_icediachk)
-      diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp
-      diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp
-      diag_trp_sv(:,:) = 0._wp
-      diag_adv_mass(:,:) = 0._wp
-      diag_adv_salt(:,:) = 0._wp
-      diag_adv_heat(:,:) = 0._wp
+      INTEGER  ::   ji, jj, jl      ! dummy loop index
+      !!----------------------------------------------------------------------
+
+      DO jj = 1, jpj 
+         DO ji = 1, jpi
+            sfx    (ji,jj) = 0._wp   ;
+            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
+            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
+            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
+            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
+            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
+            !
+            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
+            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
+            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
+            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
+            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
+            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp  
+            wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
+            wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
+            wfx_snw_sni(ji,jj) = 0._wp 
+            wfx_pnd(ji,jj) = 0._wp
+
+            hfx_thd(ji,jj) = 0._wp   ;
+            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
+            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
+            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
+            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
+            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
+            hfx_err_dif(ji,jj) = 0._wp
+            wfx_err_sub(ji,jj) = 0._wp
+            !
+            diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp
+            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
+            diag_aice(ji,jj) = 0._wp
+
+            tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here)
+            qsb_ice_bot(ji,jj) = 0._wp   ! (needed if ln_icethd=F)
+
+            fhld(ji,jj) = 0._wp   ! needed if ln_icethd=F
+
+            ! for control checks (ln_icediachk)
+            diag_trp_vi(ji,jj) = 0._wp   ;   diag_trp_vs(ji,jj) = 0._wp
+            diag_trp_ei(ji,jj) = 0._wp   ;   diag_trp_es(ji,jj) = 0._wp
+            diag_trp_sv(ji,jj) = 0._wp
+            !
+            diag_adv_mass(ji,jj) = 0._wp
+            diag_adv_salt(ji,jj) = 0._wp
+            diag_adv_heat(ji,jj) = 0._wp
+         END DO
+      END DO
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+               ! SIMIP diagnostics
+               t_si       (ji,jj,jl) = rt0     ! temp at the ice-snow interface
+               qcn_ice_bot(ji,jj,jl) = 0._wp
+               qcn_ice_top(ji,jj,jl) = 0._wp   ! conductive fluxes
+               cnd_ice    (ji,jj,jl) = 0._wp   ! effective conductivity at the top of ice/snow (ln_cndflx=T)
+               qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T)
+               qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs
+            END DO
+         END DO
+      END DO
       
    END SUBROUTINE diag_set0
+
+
+   SUBROUTINE diag_trends( kn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE diag_trends  ***
+      !!
+      !! ** purpose : diagnostics of the trends. Used for conservation purposes
+      !!              and outputs
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
+      !!----------------------------------------------------------------------
+      !
+      ! --- trends of heat, salt, mass (used for conservation controls)
+      IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
+         !
+         diag_heat(:,:) = diag_heat(:,:) &
+            &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice &
+            &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice
+         diag_sice(:,:) = diag_sice(:,:) &
+            &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice * rhoi
+         diag_vice(:,:) = diag_vice(:,:) &
+            &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhoi
+         diag_vsnw(:,:) = diag_vsnw(:,:) &
+            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhos
+         !
+         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend
+         !
+      ENDIF
+      !
+      ! --- trends of concentration (used for simip outputs)
+      IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN
+         !
+         diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
+         !
+         IF( kn == 1 )   CALL iom_put( 'afxdyn' , diag_aice )                                           ! dyn trend
+         IF( kn == 2 )   CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice ) ! thermo trend
+         IF( kn == 2 )   CALL iom_put( 'afxtot' , diag_aice )                                           ! total trend
+         !
+      ENDIF
+      !
+   END SUBROUTINE diag_trends
 
 #else
Index: /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90
===================================================================
--- /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90	(revision 13639)
+++ /NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90	(revision 13640)
@@ -30,4 +30,5 @@
    USE icethd_pnd     ! sea-ice: melt ponds
    USE iceitd         ! sea-ice: remapping thickness distribution
+   USE icecor         ! sea-ice: corrections
    USE icetab         ! sea-ice: 1D <==> 2D transformation
    USE icevar         ! sea-ice: operations
@@ -271,4 +272,8 @@
       IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- !
       !
+                              CALL ice_cor( kt , 2 )                ! --- Corrections --- !
+      !
+      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation     
+      !
       ! convergence tests
       IF( ln_zdf_chkcvg ) THEN
