Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6696)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6697)
@@ -49,4 +49,5 @@
    USE limthd_dh       ! for CALL lim_thd_snwblow
 #endif
+   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -1273,11 +1274,17 @@
       zepsilon = rn_iceshelf_fluxes_tolerance
 
+
+      ! See if we need zmask_sum...
+      IF ( srcv(jpr_grnm)%laction .OR. srcv(jpr_antm)%laction ) THEN
+         zmask_sum = glob_sum( tmask(:,:,1) )
+      ENDIF
+
       IF( srcv(jpr_grnm)%laction ) THEN
          greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
          ! take average over ocean points of input array to avoid cumulative error over time
-         zgreenland_icesheet_mass_in = SUM( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zgreenland_icesheet_mass_in )
-         zmask_sum = SUM( tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
+
+         ! The following must be bit reproducible over different PE decompositions
+         zgreenland_icesheet_mass_in = glob_sum( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
+
          zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum
          greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
@@ -1304,8 +1311,7 @@
          antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
          ! take average over ocean points of input array to avoid cumulative error from rounding errors over time
-         zantarctica_icesheet_mass_in = SUM( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zantarctica_icesheet_mass_in )
-         zmask_sum = SUM( tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
+         ! The following must be bit reproducible over different PE decompositions
+         zantarctica_icesheet_mass_in = glob_sum( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
+
          zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum
          antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
Index: /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6696)
+++ /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6697)
@@ -26,4 +26,5 @@
    USE zdfbfr
    USE fldread         ! read input field at current time step
+   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -257,7 +258,8 @@
               ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
               ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
-
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+ 
+              ! All related global sums must be done bit reproducibly
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
@@ -267,10 +269,11 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
 
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
@@ -280,6 +283,7 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
 
@@ -304,6 +308,7 @@
               ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
 
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+              ! All related global sums must be done bit reproducibly
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
@@ -313,10 +318,11 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
 
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
@@ -326,6 +332,7 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
 
