Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5708)
@@ -193,4 +193,5 @@
 
          <!-- * variable related to ice shelf forcing * -->
+         <field id="berg_calve"   long_name="Iceberg calving"                               unit="kg/m2/s"  />
          <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  />
          <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     />
Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5708)
@@ -378,4 +378,9 @@
    ln_usecplmask = .false.   !  use a coupling mask file to merge data received from several models
                              !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_coupled_iceshelf_fluxes = .false. ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the 
+                                        ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+   rn_greenland_calving_fraction = 0.5  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_antarctica_calving_fraction = 0.5 ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_iceshelf_fluxes_tolerance = 1e-6  ! Fractional threshold for detecting differences in icesheet masses (must be positive definite).
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 5708)
@@ -25,4 +25,8 @@
    USE icbutl         ! iceberg utility routines
 
+   USE sbc_oce        ! for icesheet freshwater input variables
+   USE in_out_manager
+   USE iom
+
    IMPLICIT NONE
    PRIVATE
@@ -48,4 +52,5 @@
       !
       REAL(wp)                        :: zcalving_used, zdist, zfact
+      REAL(wp)                        :: zgreenland_calving_sum, zantarctica_calving_sum
       INTEGER                         :: jn, ji, jj                    ! loop counters
       INTEGER                         :: imx                           ! temporary integer for max berg class
@@ -59,4 +64,41 @@
       zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 850._wp
       berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact
+
+      IF( lk_oasis) THEN
+      ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+      IF( ln_coupled_iceshelf_fluxes ) THEN
+
+        ! Adjust total calving rates so that sum of iceberg calving and iceshelf melting in the northern
+        ! 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_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum )
+         WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                                 &
+        &    berg_grid%calving(:,:) = berg_grid%calving(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction &
+        &                                     / ( zgreenland_calving_sum + 1.0e-10_wp )
+
+         ! check
+         WRITE(numout, *) 'Greenland iceberg calving climatology (kg/s) : ',zgreenland_calving_sum
+         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum )
+         WRITE(numout, *) 'Greenland iceberg calving adjusted value (kg/s) : ',zgreenland_calving_sum
+
+         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum )
+         WHERE( antarctica_icesheet_mask(:,:) == 1.0 )                                                                              &
+         berg_grid%calving(:,:) = berg_grid%calving(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction &
+        &                           / ( zantarctica_calving_sum + 1.0e-10_wp )
+ 
+         ! check
+         WRITE(numout, *) 'Antarctica iceberg calving climatology (kg/s) : ',zantarctica_calving_sum
+         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum )
+         WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum
+
+      ENDIF
+      ENDIF
+   
+      CALL iom_put( 'berg_calve', berg_grid%calving(:,:) )
 
       ! Heat in units of W/m2, and mask (just in case)
Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 5708)
@@ -24,4 +24,5 @@
    USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
    USE divcur          ! hor. divergence and curl      (div & cur routines)
+   USE sbc_oce         ! for icesheet freshwater input variables
 
    IMPLICIT NONE
@@ -145,4 +146,16 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       )
 #endif
+                     IF( lk_oasis) THEN
+                     ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+                     IF( ln_coupled_iceshelf_fluxes ) THEN
+                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass', greenland_icesheet_mass )
+                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed )
+                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change )
+                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass', antarctica_icesheet_mass )
+                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed )
+                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change )
+                     ENDIF
+                     ENDIF
+
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
@@ -258,4 +271,23 @@
 #endif
       !
+      IF( iom_varid( numror, 'greenland_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN
+         CALL iom_get( numror, 'greenland_icesheet_mass', greenland_icesheet_mass )
+         CALL iom_get( numror, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed )
+         CALL iom_get( numror, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change )
+      ELSE
+         greenland_icesheet_mass = 0.0 
+         greenland_icesheet_mass_rate_of_change = 0.0 
+         greenland_icesheet_timelapsed = 0.0
+      ENDIF
+      IF( iom_varid( numror, 'antarctica_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN
+         CALL iom_get( numror, 'antarctica_icesheet_mass', antarctica_icesheet_mass )
+         CALL iom_get( numror, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed )
+         CALL iom_get( numror, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change )
+      ELSE
+         antarctica_icesheet_mass = 0.0 
+         antarctica_icesheet_mass_rate_of_change = 0.0 
+         antarctica_icesheet_timelapsed = 0.0
+      ENDIF
+
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
          tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5708)
@@ -125,6 +125,6 @@
 #endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_mass_array
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_mass_array
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_icesheet_mass_array, greenland_icesheet_mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mass_array, antarctica_icesheet_mask
 
    !!----------------------------------------------------------------------
@@ -148,6 +148,14 @@
    !!----------------------------------------------------------------------
    
-   REAL, PUBLIC  :: greenland_mass
-   REAL, PUBLIC  :: antarctica_mass
+   REAL(wp), PUBLIC  :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed 
+   REAL(wp), PUBLIC  :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed
+
+   ! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to 
+   ! avoid circular dependencies.
+   LOGICAL, PUBLIC     ::   ln_coupled_iceshelf_fluxes     ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the 
+                                                           ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+   REAL(wp), PUBLIC    ::   rn_greenland_calving_fraction  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_iceshelf_fluxes_tolerance   ! Fractional threshold for detecting differences in icesheet masses. 
 
    !! * Substitutions
@@ -185,5 +193,6 @@
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      &
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
-      ALLOCATE( greenland_mass_array(jpi,jpj) , antarctica_mass_array(jpi,jpj) )
+      ALLOCATE( greenland_icesheet_mass_array(jpi,jpj) , antarctica_icesheet_mass_array(jpi,jpj) )
+      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) )
          !
 #if defined key_vvl
Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5708)
@@ -227,5 +227,6 @@
          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
          &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
-         &                  ln_usecplmask
+         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, &
+         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
       !!---------------------------------------------------------------------
       !
@@ -285,4 +286,8 @@
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes
+         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction
+         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction
+         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance
       ENDIF
 
@@ -892,4 +897,27 @@
       ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
+      IF( ln_coupled_iceshelf_fluxes ) THEN
+          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
+          ! more complicated could be done if required.
+          greenland_icesheet_mask = 0.0
+          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0
+          antarctica_icesheet_mask = 0.0
+          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0
+
+          ! initialise other variables
+          greenland_icesheet_mass_array(:,:) = 0.0
+          antarctica_icesheet_mass_array(:,:) = 0.0
+
+          IF( .not. ln_rstart ) THEN
+             greenland_icesheet_mass = 0.0 
+             greenland_icesheet_mass_rate_of_change = 0.0 
+             greenland_icesheet_timelapsed = 0.0
+             antarctica_icesheet_mass = 0.0 
+             antarctica_icesheet_mass_rate_of_change = 0.0 
+             antarctica_icesheet_timelapsed = 0.0
+          ENDIF
+
+      ENDIF
+
       CALL wrk_dealloc( jpi,jpj, zacs, zaos )
       !
@@ -953,4 +981,7 @@
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
+      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in
+      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b
+      REAL(wp) ::   zmask_sum, zepsilon      
       REAL(wp) ::   zcoef                  ! temporary scalar
       REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
@@ -1237,12 +1268,50 @@
       ENDIF
       
-      !                                                        ! land ice masses
-      IF( srcv(jpr_grnm)%laction )   greenland_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
-      IF( srcv(jpr_antm)%laction )   antarctica_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
-      greenland_mass = greenland_mass_array(1,1)
-      antarctica_mass = antarctica_mass_array(1,1)
-         
-      WRITE(numout,*) 'Greenland mass is ', greenland_mass
-      WRITE(numout,*) 'Antarctica mass is ', antarctica_mass
+      !                                                        ! land ice masses : Greenland
+      zepsilon = rn_iceshelf_fluxes_tolerance
+
+      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 ) 
+         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum
+         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
+         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass )/( greenland_icesheet_mass + zepsilon ) > zepsilon ) THEN
+            zgreenland_icesheet_mass_b = greenland_icesheet_mass
+            greenland_icesheet_mass = zgreenland_icesheet_mass_in
+            IF( zgreenland_icesheet_mass_b /= 0.0 ) &
+           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed 
+            greenland_icesheet_timelapsed = 0.0_wp       
+         ENDIF
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) is ', greenland_icesheet_mass
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
+      ENDIF
+
+      !                                                        ! land ice masses : Antarctica
+      IF( srcv(jpr_antm)%laction ) THEN
+         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 ) 
+         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum
+         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
+         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass )/( antarctica_icesheet_mass + zepsilon ) > zepsilon ) THEN
+            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
+            antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+            IF( zantarctica_icesheet_mass_b /= 0.0 ) &
+          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed 
+            antarctica_icesheet_timelapsed = 0.0_wp       
+         ENDIF
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) is ', antarctica_icesheet_mass
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
+      ENDIF
+
       !
       CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
Index: /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5707)
+++ /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5708)
@@ -27,6 +27,4 @@
    USE fldread         ! read input field at current time step
 
-
-
    IMPLICIT NONE
    PRIVATE
@@ -90,4 +88,5 @@
     INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
     INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
     REAL(wp)                     ::   rmin
     REAL(wp)                     ::   zhk
@@ -256,4 +255,42 @@
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+
+            IF( lk_oasis) THEN
+            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( ln_coupled_iceshelf_fluxes ) THEN
+
+              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
+              ! 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 )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+               ! check
+               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 )
+               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 )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+      
+               ! check
+               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 )
+               WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+            ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             stbl(:,:)   = soce
@@ -264,4 +301,42 @@
             !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
             fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+
+            IF( lk_oasis) THEN
+            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( ln_coupled_iceshelf_fluxes ) THEN
+
+              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
+              ! 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 )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+               ! check
+               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 )
+               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 )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+      
+               ! check
+               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 )
+               WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+            ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
