Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90	(revision 14715)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90	(revision 14717)
@@ -390,5 +390,11 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
-      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
+      !!---------------------------------------------------------------------
+      !
+      ! Heat content per unit mass (J/kg)
+      zcptrain(:,:) = (      sf(jp_tair)%fnow(:,:,1)        - rt0 ) * rcp  * tmask(:,:,1)
+      zcptsnw (:,:) = ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1)
+      zcptn   (:,:) =        pst(:,:)                               * rcp  * tmask(:,:,1)
       !
       ! local scalars ( place there for vector optimisation purposes)
@@ -541,9 +547,9 @@
       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
-         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
+         &     - zevap(:,:) * zcptn(:,:)                                          &   ! remove evap heat content at SST
          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
-         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
+         &     * zcptrain(:,:)                                                    &
          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
+         &     * zcptsnw(:,:)
       qns(:,:) = qns(:,:) * tmask(:,:,1)
       !
@@ -819,6 +825,6 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
-      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
       !!---------------------------------------------------------------------
       !
@@ -827,4 +833,9 @@
       !
       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      !
+      ! Heat content per unit mass (J/kg)
+      zcptrain(:,:) = (      sf(jp_tair)%fnow(:,:,1)        - rt0 ) * rcp  * tmask(:,:,1)
+      zcptsnw (:,:) = ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1)
+      zcptn   (:,:) =        sst_m(:,:)                             * rcp  * tmask(:,:,1)
       !
       zztmp = 1. / ( 1. - albo )
@@ -901,10 +912,8 @@
 
       ! --- heat flux associated with emp --- !
-      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
-         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
-         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
-      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * zcptn(:,:)         & ! evap at sst
+         &          + ( tprecip(:,:) - sprecip(:,:) )   *   zcptrain(:,:)         & ! liquid precip at Tair
+         &          +   sprecip(:,:) * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus ) ! solid precip at min(Tair,Tsnow)
+      qemp_ice(:,:) =   sprecip(:,:) *           zsnw   * ( zcptsnw (:,:) - rLfus ) ! solid precip (only)
 
       ! --- total solar and non solar fluxes --- !
@@ -914,5 +923,5 @@
 
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
+      qprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus )
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
@@ -947,22 +956,28 @@
 
       IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
-         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
-         CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
-         CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average)
-      ENDIF
-      IF( iom_use('hflx_rain_cea') ) THEN
-         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
-         CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average)
-      ENDIF
-      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
-          WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
-          ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)    
-          ENDWHERE
-          ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
-          CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
-          CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
-          CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
-      ENDIF
-      !
+         CALL iom_put( 'evap_ao_cea'  , zevap(:,:) * ( 1._wp - at_i_b(:,:) ) * tmask(:,:,1)              )   ! ice-free oce evap (cell average)
+         CALL iom_put( 'hflx_evap_cea', zevap(:,:) * ( 1._wp - at_i_b(:,:) ) * tmask(:,:,1) * zcptn(:,:) )   ! heat flux from evap (cell average)
+      ENDIF
+      IF( iom_use('rain') .OR. iom_use('rain_ao_cea') .OR. iom_use('hflx_rain_cea') ) THEN
+         CALL iom_put( 'rain'         ,   tprecip(:,:) - sprecip(:,:)                             )          ! liquid precipitation 
+         CALL iom_put( 'rain_ao_cea'  , ( tprecip(:,:) - sprecip(:,:) ) * ( 1._wp - at_i_b(:,:) ) )          ! liquid precipitation over ocean (cell average)
+         CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )                    ! heat flux from rain (cell average)
+      ENDIF
+      IF(  iom_use('snow_ao_cea')   .OR. iom_use('snow_ai_cea')      .OR. &
+         & iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
+         CALL iom_put( 'snow_ao_cea'     , sprecip(:,:)                            * ( 1._wp - zsnw(:,:) ) ) ! Snow over ice-free ocean  (cell average)
+         CALL iom_put( 'snow_ai_cea'     , sprecip(:,:)                            *           zsnw(:,:)   ) ! Snow over sea-ice         (cell average)
+         CALL iom_put( 'hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) )                         ! heat flux from snow (cell average)
+         CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
+         CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
+      ENDIF
+      IF( iom_use('hflx_prec_cea') ) THEN                                                                    ! heat flux from precip (cell average)
+         CALL iom_put('hflx_prec_cea' ,    sprecip(:,:)                  * ( zcptsnw (:,:) - rLfus )  &
+            &                          + ( tprecip(:,:) - sprecip(:,:) ) *   zcptrain(:,:) )
+      ENDIF
+      !
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Sublimation over sea-ice (cell average)
+!!clem IF( iom_use('hflx_subl_cea') ) CALL iom_put( 'hflx_subl_cea', SUM( a_i_b(:,:,:) * qevap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Heat flux from sublimation (cell average)
+
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90	(revision 14715)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90	(revision 14717)
@@ -1855,5 +1855,5 @@
       IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
       IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
-      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average)
+      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * ziceld(:,:)         )  ! liquid precipitation over ocean (cell average)
       IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , zevap_ice_total(:,:) * picefr(:,:) * tmask(:,:,1)     )  ! Sublimation over sea-ice (cell average)
       IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
@@ -2030,4 +2030,6 @@
       IF (        iom_use('hflx_snow_ai_cea') ) &                                                    ! heat flux from snow (over ice)
          &   CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *  zsnw(:,:) )
+      IF(         iom_use('hflx_subl_cea') )    &                                                    ! heat flux from sublimation
+         &   CALL iom_put('hflx_subl_cea' ,   SUM( qevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) * tmask(:,:,1) )
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
       !
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcfwb.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcfwb.F90	(revision 14715)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcfwb.F90	(revision 14717)
@@ -112,4 +112,7 @@
             emp(:,:) = emp(:,:) - z_fwfprv(1)        * tmask(:,:,1)
             qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+            ! outputs
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', zcoef * sst_m(:,:) * tmask(:,:,1) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1)        * tmask(:,:,1) )
          ENDIF
          !
@@ -156,4 +159,7 @@
             emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1)
             qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+            ! outputs
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zcoef * sst_m(:,:) * tmask(:,:,1) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -fwfold             * tmask(:,:,1) )
          ENDIF
          ! Output restart information
@@ -210,4 +216,7 @@
             qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction
             erp(:,:) = erp(:,:) + zerp_cor(:,:)
+            ! outputs
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zerp_cor(:,:) * rcp * sst_m(:,:) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -zerp_cor(:,:) )
             !
             IF( nprint == 1 .AND. lwp ) THEN                   ! control print
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcrnf.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcrnf.F90	(revision 14715)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcrnf.F90	(revision 14717)
@@ -128,6 +128,10 @@
              IF( ln_rnf_icb ) THEN
                 fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt
-                CALL iom_put( 'iceberg_cea'  , fwficb(:,:)  )         ! output iceberg flux
-                CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
+                rnf(:,:) = rnf(:,:) + fwficb(:,:)
+                qns(:,:) = qns(:,:) - fwficb(:,:) * rLfus
+                !!qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus                
+                !!qns_oce(:,:) = qns_oce(:,:) - fwficb(:,:) * rLfus                
+                CALL iom_put( 'iceberg_cea'  ,  fwficb(:,:)  )          ! output iceberg flux
+                CALL iom_put( 'hflx_icb_cea' , -fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
              ENDIF
          ENDIF
@@ -149,4 +153,5 @@
                                          CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
          IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rau0 * rcp )   ! output runoff sensible heat (W/m2)
+         IF( iom_use('sflx_rnf_cea') )   CALL iom_put( 'sflx_rnf_cea', rnf_tsc(:,:,jp_sal) * rau0       )   ! output runoff salt flux (g/m2/s)
       ENDIF
       !
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcssr.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcssr.F90	(revision 14715)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcssr.F90	(revision 14717)
@@ -91,4 +91,6 @@
          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    Add restoring term     !
             !                                      ! ========================= !
+            qrp(:,:) = 0._wp ! necessary init
+            erp(:,:) = 0._wp
             !
             IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
@@ -140,7 +142,13 @@
                      qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj)
                      erp(ji,jj) = zerp
+                     qrp(ji,jj) = qrp(ji,jj) - zerp * rcp * sst_m(ji,jj)
                   END DO
                END DO
             ENDIF
+            !
+            ! outputs
+            CALL iom_put( 'hflx_ssr_cea', qrp(:,:) )
+            IF( nn_sssr == 1 )   CALL iom_put( 'sflx_ssr_cea',  erp(:,:) * sss_m(:,:) )
+            IF( nn_sssr == 2 )   CALL iom_put( 'vflx_ssr_cea', -erp(:,:) )
             !
          ENDIF
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/TRA/tradmp.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/TRA/tradmp.F90	(revision 14715)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/TRA/tradmp.F90	(revision 14717)
@@ -99,5 +99,5 @@
       IF( ln_timing )   CALL timing_start('tra_dmp')
       !
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+      IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN   !* Save ta and sa trends
          ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
          ztrdts(:,:,:,:) = tsa(:,:,:,:) 
@@ -149,4 +149,10 @@
       END SELECT
       !
+      ! outputs
+      IF( iom_use('hflx_dmp_cea') ) &
+         & CALL iom_put('hflx_dmp_cea', SUM( ( tsa(:,:,:,jp_tem) - ztrdts(:,:,:,jp_tem) ) * e3t_n(:,:,:), dim=3 ) * rcp * rau0 ) ! W/m2
+      IF( iom_use('sflx_dmp_cea') ) &
+         & CALL iom_put('sflx_dmp_cea', SUM( ( tsa(:,:,:,jp_sal) - ztrdts(:,:,:,jp_sal) ) * e3t_n(:,:,:), dim=3 ) * rau0 )       ! g/m2/s
+      !
       IF( l_trdtra )   THEN       ! trend diagnostic
          ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)
