Index: NEMO/branches/2021/ticket2607_r14608_halo1_halo2_compatibility/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2021/ticket2607_r14608_halo1_halo2_compatibility/src/OCE/SBC/sbcblk.F90	(revision 14617)
+++ NEMO/branches/2021/ticket2607_r14608_halo1_halo2_compatibility/src/OCE/SBC/sbcblk.F90	(revision 14807)
@@ -892,8 +892,12 @@
       REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable
       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! net long wave radiative heat flux
-      !!---------------------------------------------------------------------
-      !
-      ! local scalars ( place there for vector optimisation purposes)
-
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
+      !!---------------------------------------------------------------------
+      !
+      ! Heat content per unit mass (J/kg)
+      zcptrain(:,:) = (      ptair        - rt0 ) * rcp  * tmask(:,:,1)
+      zcptsnw (:,:) = ( MIN( ptair, rt0 ) - rt0 ) * rcpi * tmask(:,:,1)
+      zcptn   (:,:) =        ptsk                 * rcp  * tmask(:,:,1)
+      !
       ! ----------------------------------------------------------------------------- !
       !     III    Net longwave radiative FLUX                                        !
@@ -907,14 +911,11 @@
       ! ----------------------------------------------------------------------------- !
       !
-      emp (:,:) = (  pevp(:,:)                                       &   ! mass flux (evap. - precip.)
-         &         - pprec(:,:) * rn_pfac  ) * tmask(:,:,1)
-      !
-      qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:)                   &   ! Downward Non Solar
-         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
-         &     - pevp(:,:) * ptsk(:,:) * rcp                         &   ! remove evap heat content at SST
-         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair
-         &     * ( ptair(:,:) - rt0 ) * rcp                          &
-         &     + psnow(:,:) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
+      emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac ) * tmask(:,:,1)      ! mass flux (evap. - precip.)
+      !
+      qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:)                     &   ! Downward Non Solar
+         &     - psnow(:,:) * rn_pfac * rLfus                          &   ! remove latent melting heat for solid precip
+         &     - pevp(:,:) * zcptn(:,:)                                &   ! remove evap heat content at SST
+         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac * zcptrain(:,:) &   ! add liquid precip heat content at Tair
+         &     + psnow(:,:) * rn_pfac * zcptsnw(:,:)                       ! add solid  precip heat content at min(Tair,Tsnow)
       qns(:,:) = qns(:,:) * tmask(:,:,1)
       !
@@ -1000,5 +1001,5 @@
       ! C-grid ice dynamics :   U & V-points (same as ocean)
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-      wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
+         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
       END_2D
       !
@@ -1120,6 +1121,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
-      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)
       !!---------------------------------------------------------------------
       !
@@ -1130,4 +1131,9 @@
       dqla_ice(:,:,:) = 0._wp
 
+      ! Heat content per unit mass (J/kg)
+      zcptrain(:,:) = (      ptair        - rt0 ) * rcp  * tmask(:,:,1)
+      zcptsnw (:,:) = ( MIN( ptair, rt0 ) - rt0 ) * rcpi * tmask(:,:,1)
+      zcptn   (:,:) =        sst_m                * rcp  * tmask(:,:,1)
+      !
       !                                     ! ========================== !
       DO jl = 1, jpl                        !  Loop over ice categories  !
@@ -1205,10 +1211,8 @@
 
       ! --- heat flux associated with emp --- !
-      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
-         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp               & ! liquid precip at Tair
-         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
-         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
-      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
-         &              ( ( MIN( ptair(:,:), 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 --- !
@@ -1218,5 +1222,5 @@
 
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), 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) ---
@@ -1250,22 +1254,27 @@
       !
       IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
-         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
-         IF( iom_use('evap_ao_cea'  ) )  CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
-         IF( iom_use('hflx_evap_cea') )  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(:,:) ) )
-         IF( iom_use('hflx_rain_cea') )  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 )
-         IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
-         IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
-         IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
+         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') .OR. iom_use('hflx_subl_cea') ) THEN
+         CALL iom_put( 'subl_ai_cea'  , SUM( a_i_b(:,:,:) *  evap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Sublimation over sea-ice (cell average)
+         CALL iom_put( 'hflx_subl_cea', SUM( a_i_b(:,:,:) * qevap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Heat flux from sublimation (cell average)
       ENDIF
       !
