Index: NEMO/trunk/src/OCE/DIA/diahsb.F90
===================================================================
--- NEMO/trunk/src/OCE/DIA/diahsb.F90	(revision 15017)
+++ NEMO/trunk/src/OCE/DIA/diahsb.F90	(revision 15048)
@@ -83,44 +83,56 @@
       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,jpkm1) ::   zwrk         ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,13)      ::   ztmp
+      REAL(wp), DIMENSION(jpi,jpj,jpkm1,4) ::   ztmpk
+      REAL(wp), DIMENSION(17)              ::   zbg          
       !!---------------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_hsb')
       !
+      ztmp (:,:,:)   = 0._wp ! should be better coded
+      ztmpk(:,:,:,:) = 0._wp ! should be better coded
+      !
       ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
       ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
+      !
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes
-      z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes
-      z_frc_trd_s =           glob_sum( 'diahsb', 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( 'diahsb', rnf_tsc(:,:,jp_tem) * surf(:,:) )
-      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) )
-      !                    ! Add ice shelf heat & salt input
-      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t &
-         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) )
-      !                    ! Add penetrative solar radiation
-      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rho0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) )
-      !                    ! Add geothermal heat flux
-      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
-      !
-      IF( ln_linssh ) THEN
+      ! prepare trends
+      ztmp(:,:,1)  = - r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) ) * surf(:,:)    ! volume
+      ztmp(:,:,2)  =   sbc_tsc(:,:,jp_tem) * surf(:,:)                                                      ! heat
+      ztmp(:,:,3)  =   sbc_tsc(:,:,jp_sal) * surf(:,:)                                                      ! salt
+      IF( ln_rnf     )    ztmp(:,:,4) =   rnf_tsc(:,:,jp_tem) * surf(:,:)                                   ! runoff temp
+      IF( ln_rnf_sal )    ztmp(:,:,5) =   rnf_tsc(:,:,jp_sal) * surf(:,:)                                   ! runoff salt
+      IF( ln_isf     )    ztmp(:,:,6) = ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ! isf temp
+      IF( ln_traqsr  )    ztmp(:,:,7) =   r1_rho0_rcp * qsr(:,:) * surf(:,:)                                ! penetrative solar radiation
+      IF( ln_trabbc  )    ztmp(:,:,8) =   qgh_trd0(:,:) * surf(:,:)                                         ! geothermal heat
+      !
+      IF( ln_linssh ) THEN   ! Advection flux through fixed surface (z=0)
          IF( ln_isfcav ) THEN
             DO ji=1,jpi
                DO jj=1,jpj
-                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
-                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
+                  ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
                END DO
             END DO
          ELSE
-            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
-            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
+            ztmp(:,:,9 ) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            ztmp(:,:,10) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
          END IF
-         z_wn_trd_t = - glob_sum( 'diahsb', z2d0 )
-         z_wn_trd_s = - glob_sum( 'diahsb', z2d1 )
-      ENDIF
-
+      ENDIF
+      
+      ! global sum
+      zbg(1:10) = glob_sum_vec( 'dia_hsb', ztmp(:,:,1:10) )
+
+      ! adding up
+      z_frc_trd_v = zbg(1)  ! volume fluxes
+      z_frc_trd_t = zbg(2)  ! heat fluxes
+      z_frc_trd_s = zbg(3)  ! salt fluxes
+      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + zbg(4) ! runoff heat
+      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + zbg(5) ! runoff salt
+      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t + zbg(6) ! isf heat
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + zbg(7) ! penetrative solar flux
+      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + zbg(8) ! geothermal heat
+      !
       frc_v = frc_v + z_frc_trd_v * rn_Dt
       frc_t = frc_t + z_frc_trd_t * rn_Dt
@@ -128,15 +140,18 @@
       !                                          ! Advection flux through fixed surface (z=0)
       IF( ln_linssh ) THEN
+         z_wn_trd_t = zbg(9)
+         z_wn_trd_s = zbg(10)
+         !
          frc_wn_t = frc_wn_t + z_wn_trd_t * rn_Dt
          frc_wn_s = frc_wn_s + z_wn_trd_s * rn_Dt
       ENDIF
 
-      ! ------------------------ !
-      ! 2 -  Content variations  !
-      ! ------------------------ !
+      ! --------------------------------- !
+      ! 2 -  Content variations with ssh  !
+      ! --------------------------------- !
       ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
-
+      !
       !                    ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
+      ztmp(:,:,11) = surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:)
 
       !                    ! heat & salt content variation (associated with ssh)
@@ -145,34 +160,55 @@
             DO ji = 1, jpi
                DO jj = 1, jpj
-                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
-                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
+                  ztmp(ji,jj,12) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
+                  ztmp(ji,jj,13) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
                END DO
             END DO
          ELSE                          ! no under ice-shelf seas
-            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
-            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
+            ztmp(:,:,12) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
+            ztmp(:,:,13) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
          END IF
-         z_ssh_hc = glob_sum_full( 'diahsb', z2d0 )
-         z_ssh_sc = glob_sum_full( 'diahsb', z2d1 )
-      ENDIF
-      !
-      DO jk = 1, jpkm1           ! volume variation (calculated with scale factors)
-         zwrk(:,:,jk) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   &
-            &           - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk)
+      ENDIF
+
+      ! global sum
+      zbg(11:13) = glob_sum_full_vec( 'dia_hsb', ztmp(:,:,11:13) )
+      
+      zdiff_v1 = zbg(11)
+      !                    ! heat & salt content variation (associated with ssh)
+      IF( ln_linssh ) THEN       ! linear free surface case
+         z_ssh_hc = zbg(12)
+         z_ssh_sc = zbg(13)
+      ENDIF
+      !
+      ! --------------------------------- !
+      ! 3 -  Content variations with e3t  !
+      ! --------------------------------- !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
+      !
+      DO jk = 1, jpkm1           ! volume
+         ztmpk(:,:,jk,1) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   &
+            &              - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk)
       END DO
-      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different
-      DO jk = 1, jpkm1           ! heat content variation
-         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
-            &           - surf_ini(:,:) *         hc_loc_ini(:,:,jk) )
+      DO jk = 1, jpkm1           ! heat
+         ztmpk(:,:,jk,2) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
+            &              - surf_ini(:,:) *         hc_loc_ini(:,:,jk) )
       END DO
-      zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
-      DO jk = 1, jpkm1           ! salt content variation
-         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
-            &           - surf_ini(:,:) *         sc_loc_ini(:,:,jk) )
+      DO jk = 1, jpkm1           ! salt
+         ztmpk(:,:,jk,3) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
+            &              - surf_ini(:,:) *         sc_loc_ini(:,:,jk) )
       END DO
-      zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+      DO jk = 1, jpkm1           ! total ocean volume
+         ztmpk(:,:,jk,4) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+      END DO
+      
+      ! global sum
+      zbg(14:17) = glob_sum_full_vec( 'dia_hsb', ztmpk(:,:,:,1:4) )
+      
+      zdiff_v2 = zbg(14)     ! glob_sum_full needed as tmask and tmask_ini could be different
+      zdiff_hc = zbg(15)
+      zdiff_sc = zbg(16)
+      zvol_tot = zbg(17)
 
       ! ------------------------ !
-      ! 3 -  Drifts              !
+      ! 4 -  Drifts              !
       ! ------------------------ !
       zdiff_v1 = zdiff_v1 - frc_v
@@ -186,12 +222,4 @@
          zerr_sc1  = z_ssh_sc - frc_wn_s
       ENDIF
-
-      ! ----------------------- !
-      ! 4 - Diagnostics writing !
-      ! ----------------------- !
-      DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors)
-         zwrk(:,:,jk) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
-      END DO
-      zvol_tot = glob_sum( 'diahsb', zwrk(:,:,:) )
 
 !!gm to be added ?
