Index: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12576)
@@ -264,4 +264,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="fwfisf3d"     long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" />
Index: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/cfgs/SHARED/namelist_ref	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/cfgs/SHARED/namelist_ref	(revision 12576)
@@ -290,4 +290,17 @@
    !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
    nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
+   nn_coupled_iceshelf_fluxes = 0 ! =0 : total freshwater input from iceberg calving and ice shelf basal melting 
+                                  ! taken from climatologies used (no action in coupling routines).
+                                  ! =1 :  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.
+                                  ! =2 :  specify constant freshwater inputs in this namelist to set the combined
+                                  ! magnitude of iceberg calving and iceshelf melting freshwater fluxes.
+   ln_iceshelf_init_atmos     = .true.  ! If true force ocean to initialise icesheet masses from atmospheric values rather than
+                                        ! from values in ocean restart file. 
+   rn_greenland_total_fw_flux   = 0.0  ! Constant total rate of freshwater input (kg/s) for Greenland (if nn_coupled_iceshelf_fluxes=2) 
+   rn_greenland_calving_fraction = 0.5  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_antarctica_total_fw_flux  = 0.0  ! Constant total rate of freshwater input (kg/s) for Antarctica (if nn_coupled_iceshelf_fluxes=2)
+   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: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/ICB/icbclv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/ICB/icbclv.F90	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/ICB/icbclv.F90	(revision 12576)
@@ -26,4 +26,8 @@
    USE icb_oce        ! iceberg parameters 
 
+   USE sbc_oce        ! for icesheet freshwater input variables 
+   USE in_out_manager 
+   USE iom 
+
    IMPLICIT NONE
    PRIVATE
@@ -49,4 +53,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
@@ -63,4 +68,42 @@
       ! Heat in units of W/m2, and mask (just in case)
       berg_grid%calving_hflx(:,:) = src_calving_hflx(:,:) * tmask_i(:,:) * tmask(:,:,1)
+
+      IF( lk_oasis) THEN
+        ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+        IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
+          ll_write = ((MOD( kt, sn_cfctl%ptimincr ) == 0) .OR. ( kt == nitend )) .AND. lwp
+          ! 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
+           IF(ll_write) 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 )
+           IF(ll_write) 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
+           IF(ll_write) 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 )
+           IF(ll_write) WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum
+
+        ENDIF
+      ENDIF
+   
+      CALL iom_put( 'berg_calve', berg_grid%calving(:,:) )
+
 
       IF( ll_first_call .AND. .NOT. l_restarted_bergs ) THEN      ! This is a hack to simplify initialization
Index: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/IOM/restart.F90	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/IOM/restart.F90	(revision 12576)
@@ -29,4 +29,5 @@
    USE diurnal_bulk
    USE lib_mpp         ! distribued memory computing library
+   USE sbc_oce         ! for icesheet freshwater input variables 
 
    IMPLICIT NONE
@@ -161,4 +162,16 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn, ldxios = lwxios      )
                      CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      )
+
+                     IF( lk_oasis) THEN
+                     ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+                       IF( nn_coupled_iceshelf_fluxes .eq. 1 ) 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
                   ! extra variable needed for the ice sheet coupling
                   IF ( ln_iscpl ) THEN 
@@ -295,4 +308,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: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbc_oce.F90	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbc_oce.F90	(revision 12576)
@@ -136,4 +136,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+   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 
 
    !!----------------------------------------------------------------------
@@ -149,4 +151,31 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
 
+   !!----------------------------------------------------------------------
+   !!  Surface scalars of total ice sheet mass for Greenland and Antarctica, 
+   !! passed from atmosphere to be converted to dvol and hence a freshwater 
+   !! flux  by using old values. New values are saved in the dump, to become
+   !! old values next coupling timestep. Freshwater fluxes split between 
+   !! sub iceshelf melting and iceberg calving, scalled to flux per second
+   !!----------------------------------------------------------------------
+   
+   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.
+   INTEGER, PUBLIC     ::   nn_coupled_iceshelf_fluxes     ! =0 : total freshwater input from iceberg calving and ice shelf basal melting 
+                                                           ! taken from climatologies used (no action in coupling routines).
+                                                           ! =1 :  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.
+                                                           ! =2 :  specify constant freshwater inputs in this namelist to set the combined
+                                                           ! magnitude of iceberg calving and iceshelf melting freshwater fluxes.
+   LOGICAL, PUBLIC     ::   ln_iceshelf_init_atmos         ! If true force ocean to initialise iceshelf masses from atmospheric values rather
+                                                           ! than values in ocean restart (applicable if nn_coupled_iceshelf_fluxes=1).
+   REAL(wp), PUBLIC    ::   rn_greenland_total_fw_flux    ! Constant total rate of freshwater input (kg/s) for Greenland (if nn_coupled_iceshelf_fluxes=2) 
+   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_total_fw_flux   ! Constant total rate of freshwater input (kg/s) for Antarctica (if nn_coupled_iceshelf_fluxes=2) 
+   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   ! Absolute tolerance for detecting differences in icesheet masses. 
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -182,4 +211,7 @@
          &      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_icesheet_mass_array(jpi,jpj) , antarctica_icesheet_mass_array(jpi,jpj) ) 
+      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) ) 
          !
       ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
Index: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90	(revision 12576)
@@ -116,6 +116,8 @@
    INTEGER, PARAMETER ::   jpr_tauwy  = 56   ! y component of the ocean stress from waves
    INTEGER, PARAMETER ::   jpr_ts_ice = 57   ! Sea ice surface temp
-
-   INTEGER, PARAMETER ::   jprcv      = 57   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_grnm   = 58   ! Greenland ice mass 
+   INTEGER, PARAMETER ::   jpr_antm   = 59   ! Antarctic ice mass 
+
+   INTEGER, PARAMETER ::   jprcv      = 59   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -174,5 +176,6 @@
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr,  &
       &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf,      &
+                    sn_rcv_grnm, sn_rcv_antm
    ! Send to waves 
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
@@ -256,5 +259,7 @@
          &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   &
          &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   &
-         &                  sn_rcv_ts_ice
+         &                  sn_rcv_ts_ice, sn_rcv_grnm  , sn_rcv_antm  , nn_coupled_iceshelf_fluxes,   &
+         &                  rn_greenland_calving_fraction, rn_antarctica_calving_fraction,             &
+         &                  rn_iceshelf_fluxes_tolerance
 
       !!---------------------------------------------------------------------
@@ -292,4 +297,6 @@
          WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
          WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
+         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')' 
+         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')' 
          WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')'
          WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
@@ -330,4 +337,11 @@
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
          WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
+         WRITE(numout,*)'  nn_coupled_iceshelf_fluxes          = ', nn_coupled_iceshelf_fluxes 
+         WRITE(numout,*)'  ln_iceshelf_init_atmos              = ', ln_iceshelf_init_atmos
+         WRITE(numout,*)'  rn_greenland_total_fw_flux         = ', rn_greenland_total_fw_flux
+         WRITE(numout,*)'  rn_antarctica_total_fw_flux        = ', rn_antarctica_total_fw_flux
+         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
 
@@ -471,4 +485,6 @@
       !
       srcv(jpr_cal)%clname = 'OCalving'   ;  IF( TRIM( sn_rcv_cal%cldes) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
+      srcv(jpr_grnm  )%clname = 'OGrnmass'   ;   IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' )   srcv(jpr_grnm)%laction = .TRUE. 
+      srcv(jpr_antm  )%clname = 'OAntmass'   ;   IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' )   srcv(jpr_antm)%laction = .TRUE. 
       srcv(jpr_isf)%clname = 'OIcshelf'   ;  IF( TRIM( sn_rcv_isf%cldes) == 'coupled' )   srcv(jpr_isf)%laction = .TRUE.
       srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE.
@@ -1045,4 +1061,27 @@
          &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
       IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
+
+      IF( nn_coupled_iceshelf_fluxes .gt. 0 ) 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 
       !
    END SUBROUTINE sbc_cpl_init
@@ -1103,5 +1142,8 @@
       INTEGER  ::   ji, jj, jn             ! dummy loop indices
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdt did not change since nit000)
-      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
+      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
@@ -1445,4 +1487,82 @@
          IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
          !
+      ENDIF
+
+      !                                                        ! land ice masses : Greenland
+      zepsilon = rn_iceshelf_fluxes_tolerance
+
+      IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
+      
+         ! This is a zero dimensional, single value field. 
+         zgreenland_icesheet_mass_in =  frcv(jpr_grnm)%z3(1,1,1)
+           
+         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
+
+         IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
+            ! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
+            ! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
+            zgreenland_icesheet_mass_b = zgreenland_icesheet_mass_in
+            greenland_icesheet_mass = zgreenland_icesheet_mass_in
+         ENDIF
+
+         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
+            zgreenland_icesheet_mass_b = greenland_icesheet_mass
+            
+            ! Only update the mass if it has increased.
+            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
+               greenland_icesheet_mass = zgreenland_icesheet_mass_in
+            ENDIF
+            
+            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 .AND. ll_wrtstp) THEN
+            WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
+            WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass
+            WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
+            WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
+            IF(lflush) CALL flush(numout)
+         ENDIF
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         greenland_icesheet_mass_rate_of_change = rn_greenland_total_fw_flux
+      ENDIF
+
+      !                                                        ! land ice masses : Antarctica
+      IF( srcv(jpr_antm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
+         
+         ! This is a zero dimensional, single value field. 
+         zantarctica_icesheet_mass_in = frcv(jpr_antm)%z3(1,1,1)
+           
+         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
+
+         IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
+            ! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
+            ! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
+            zantarctica_icesheet_mass_b = zantarctica_icesheet_mass_in
+            antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+         ENDIF
+
+         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
+            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
+            
+            ! Only update the mass if it has increased.
+            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
+               antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+            END IF
+            
+            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 .AND. ll_wrtstp) THEN
+            WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
+            WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass
+            WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
+            WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
+            IF(lflush) CALL flush(numout)
+         ENDIF
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
       ENDIF
       !
Index: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90	(revision 12575)
+++ NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90	(revision 12576)
@@ -92,4 +92,5 @@
       INTEGER ::   ji, jj, jk   ! loop index
       INTEGER ::   ikt, ikb     ! local integers
+      REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum 
       REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep) 
       REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d
@@ -127,4 +128,45 @@
                fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
             ENDIF
+
+            IF( lk_oasis) THEN
+              ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+              IF( nn_coupled_iceshelf_fluxes .gt. 0 ) 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.
+
+                ! 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 )                                                                  &
+                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+                &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+                 ! check
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
+
+                 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_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 ) &
+                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+                &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+
+                 ! check
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
+
+                 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+              ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * rLfusisf             ! heat flux
             stbl(:,:)   = soce
@@ -137,4 +179,45 @@
                fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf
             ENDIF
+
+            IF( lk_oasis) THEN
+              ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+              IF( nn_coupled_iceshelf_fluxes .gt. 0 ) 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.
+
+                ! 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 )                                                                  &
+                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+                &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+                 ! check
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
+
+                 zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_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 ) &
+                &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+                &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+
+                 ! check
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
+
+                 zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
+                 IF(lwp .AND. ll_wrtstp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+              ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * rLfusisf               ! heat flux
             stbl(:,:)   = soce
