Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 9226)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 9227)
@@ -81,5 +81,6 @@
       REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     -
       REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     -
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d0, z2d1   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)       ::   z2d0, z2d1   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwrk         ! 3D workspace
       !!---------------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_hsb')      
@@ -90,15 +91,15 @@
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes
-      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                               ! heat fluxes
-      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                               ! salt fluxes
-      ! Add runoff    heat & salt input
+      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) )   ! volume fluxes
+      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes
+      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes
+      !                    !  Add runoff    heat & salt input
       IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) )
       IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
-      ! Add ice shelf heat & salt input
-      IF( ln_isf    ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) )
-      ! Add penetrative solar radiation
+      !                    ! Add ice shelf heat & salt input
+      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) )
+      !                    ! Add penetrative solar radiation
       IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
-      ! Add geothermal heat flux
+      !                    ! Add geothermal heat flux
       IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * surf(:,:) )
       !
@@ -129,17 +130,14 @@
 
       ! ------------------------ !
-      ! 2 -  Content variations !
+      ! 2 -  Content variations  !
       ! ------------------------ !
       ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
-      zdiff_v2 = 0._wp
-      zdiff_hc = 0._wp
-      zdiff_sc = 0._wp
-
-      ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) )
-
-      ! heat & salt content variation (associated with ssh)
-      IF( ln_linssh ) THEN
-         IF( ln_isfcav ) THEN
+
+      !                    ! volume variation (calculated with ssh)
+      zdiff_v1 = glob_sum_full( surf(:,:)*sshn(:,:) - surf_ini(:,:)*ssh_ini(:,:) )
+
+      !                    ! heat & salt content variation (associated with ssh)
+      IF( ln_linssh ) THEN       ! linear free surface case
+         IF( ln_isfcav ) THEN          ! ISF case
             DO ji = 1, jpi
                DO jj = 1, jpj
@@ -148,5 +146,5 @@
                END DO
             END DO
-         ELSE
+         ELSE                          ! no under ice-shelf seas
             z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 
             z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
@@ -155,16 +153,17 @@
          z_ssh_sc = glob_sum_full( z2d1 ) 
       ENDIF
-
-      DO jk = 1, jpkm1
-         ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 + glob_sum_full( surf(:,:) * tmask(:,:,jk)            &
-            &                           * e3t_n(:,:,jk) - surf_ini(:,:) * e3t_ini(:,:,jk) )
-         ! heat content variation
-         zdiff_hc = zdiff_hc + glob_sum_full(  surf(:,:) * tmask(:,:,jk)                                   & 
-            &                           * e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 
-         ! salt content variation
-         zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk)                           &
-                                        * e3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - surf_ini(:,:) * sc_loc_ini(:,:,jk) )
-      ENDDO
+      !
+      DO jk = 1, jpkm1     ! volume variation (calculated with scale factors)
+         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk) ) * tmask(:,:,jk)
+      END DO
+      zdiff_v2 = glob_sum_full( zwrk(:,:,:) )
+      DO jk = 1, jpkm1           ! heat content variation
+         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_tem) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
+      END DO
+      zdiff_hc = glob_sum_full( zwrk(:,:,:) )
+      DO jk = 1, jpkm1           ! salt content variation
+         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_sal) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
+      END DO
+      zdiff_sc = glob_sum_full( zwrk(:,:,:) )
 
       ! ------------------------ !
@@ -185,8 +184,8 @@
       ! 4 - Diagnostics writing !
       ! ----------------------- !
-      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors)
-      DO jk = 1, jpkm1
-         zvol_tot  = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * e3t_n(:,:,jk) )
+      DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors)
+         zwrk(:,:,jk) = surf(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
+      zvol_tot = glob_sum_full( zwrk(:,:,:) )
 
 !!gm to be added ?
@@ -203,22 +202,33 @@
 
       IF( .NOT. ln_linssh ) THEN
-        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C) 
-        CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (pss)
-        CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content drift    (1.e20 J) 
-        CALL iom_put( 'bgheatfx' , zdiff_hc * rau0 * rcp /  &         ! Heat flux drift       (W/m2) 
-           &                       ( surf_tot * kt * rdt )        )
-        CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)  
-        CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t drift      (km3)  
+         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C) 
+         CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (PSU)
+         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content drift    (1.e20 J) 
+         CALL iom_put( 'bgheatfx' , zdiff_hc * rau0 * rcp /  &         ! Heat flux drift       (W/m2) 
+            &                       ( surf_tot * kt * rdt )        )
+         CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3)
+         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)  
+         CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t drift      (km3)  
+         !
+         IF( kt == nitend .AND. lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+            WRITE(numout,*) '~~~~~~~'
+            WRITE(numout,*) '   Temperature drift = ', zdiff_hc / zvol_tot, ' C'
+            WRITE(numout,*) '   Salinity    drift = ', zdiff_sc / zvol_tot, ' PSU'
+            WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1 * 1.e-9   , ' km^3'
+            WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2 * 1.e-9   , ' km^3'
+         ENDIF
+         !
       ELSE
-        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C) 
-        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (pss)
-        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content drift    (1.e20 J) 
-        CALL iom_put( 'bgheatfx' , zdiff_hc1 * rau0 * rcp /  &        ! Heat flux drift       (W/m2) 
-           &                       ( surf_tot * kt * rdt )         )
-        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)  
-        CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
-        CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
+         CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C) 
+         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (PSU)
+         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content drift    (1.e20 J) 
+         CALL iom_put( 'bgheatfx' , zdiff_hc1 * rau0 * rcp /  &        ! Heat flux drift       (W/m2) 
+            &                       ( surf_tot * kt * rdt )         )
+         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3)
+         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)  
+         CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
+         CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
       ENDIF
       !
@@ -380,6 +390,6 @@
       ! 2 - Time independant variables and file opening !
       ! ----------------------------------------------- !
-      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)      ! masked surface grid cell area
-      surf_tot  = glob_sum( surf(:,:) )                   ! total ocean surface area
+      surf(:,:) = e1e2t(:,:) * tmask_i(:,:)     ! masked surface grid cell area
+      surf_tot  = glob_sum( surf(:,:) )         ! total ocean surface area
 
       IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )         
