Index: /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6250)
@@ -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/nemo_v3_6_STABLE_pkg/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6250)
@@ -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/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 6250)
@@ -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
+         IF(lwp) 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(lwp) 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(lwp) 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(lwp) 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/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6250)
@@ -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/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 6250)
@@ -125,4 +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_icesheet_mass_array, greenland_icesheet_mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mass_array, antarctica_icesheet_mask
 
    !!----------------------------------------------------------------------
@@ -137,4 +139,23 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m]
    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.
+   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   ! Absolute tolerance for detecting differences in icesheet masses. 
 
    !! * Substitutions
@@ -172,4 +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_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/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6250)
@@ -103,5 +103,7 @@
    INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
    INTEGER, PARAMETER ::   jpr_ts_ice = 43            ! skin temperature of sea-ice (used for melt-ponds)
-   INTEGER, PARAMETER ::   jprcv      = 43            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_grnm   = 44            ! Greenland ice mass
+   INTEGER, PARAMETER ::   jpr_antm   = 45            ! Antarctic ice mass
+   INTEGER, PARAMETER ::   jprcv      = 45            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -153,5 +155,5 @@
    ! Received from the atmosphere                     !
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice 
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -224,5 +226,7 @@
          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     &
          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
-         &                  sn_rcv_co2 , sn_rcv_ts_ice, nn_cplmodel  , ln_usecplmask
+         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
+         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, &
+         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
       !!---------------------------------------------------------------------
       !
@@ -263,4 +267,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,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
@@ -280,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
 
@@ -421,4 +431,7 @@
       !
       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.
+
 
       !                                                      ! ------------------------- !
@@ -884,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 )
       !
@@ -945,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
@@ -1228,4 +1267,63 @@
 
       ENDIF
+      
+      !                                                        ! 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 ) > 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) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used 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 ) > 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) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used 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/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 6250)
@@ -172,5 +172,4 @@
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
       REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
       INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices
       !!---------------------------------------------------------------------
@@ -185,29 +184,8 @@
       jj_off = INT ( (jpjglo - ny_global) / 2 )
 
-      ! Set freezing temperatures and ensure consistencey between NEMO and CICE
-      CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d ) 
-      DO jk=1,jpk
-         ztfrz3d(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept_n(:,:,jk) )
-      ENDDO
-
-      !Ensure that no temperature points are below freezing if not a NEMO restart
-      IF( .NOT. ln_rstart ) THEN
-         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),ztfrz3d)
-         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
-      ENDIF
-
-#if defined key_nemocice_decomp
-      ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
-      ! there is no restart file.
-      ! Values from a CICE restart file would overwrite this
-      IF ( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
-      ENDIF  
-#endif
-
-! Initialize CICE
+      ! Initialize CICE
       CALL CICE_Initialize
 
-! Do some CICE consistency checks
+      ! Do some CICE consistency checks
       IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          IF ( calc_strair .OR. calc_Tsfc ) THEN
@@ -221,13 +199,31 @@
 
 
-! allocate sbc_ice and sbc_cice arrays
-      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
+      ! allocate sbc_ice and sbc_cice arrays
+      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' )
       IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
 
-      ! Populate the surface freezing temperature array
-      sstfrz(:,:)=ztfrz3d(:,:,1)
-
-      fr_iu(:,:)=0.0
-      fr_iv(:,:)=0.0
+      ! Ensure that no temperature points are below freezing if not a NEMO restart
+      IF( .NOT. ln_rstart ) THEN
+
+         CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d ) 
+         DO jk=1,jpk
+            ztfrz3d(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept_n(:,:,jk) )
+         ENDDO
+         tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
+         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
+         CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d ) 
+
+#if defined key_nemocice_decomp
+         ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
+         ! there is no restart file.
+         ! Values from a CICE restart file would overwrite this
+         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
+#endif
+
+      ENDIF  
+
+      ! calculate surface freezing temperature and send to CICE
+      sstfrz(:,:) = eos_fzp(sss_m(:,:), fsdept_n(:,:,1)) 
+      CALL nemo2cice(sstfrz,Tf, 'T', 1. )
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
@@ -240,4 +236,6 @@
 ! T point to U point
 ! T point to V point
+      fr_iu(:,:)=0.0
+      fr_iv(:,:)=0.0
       DO jj=1,jpjm1
          DO ji=1,jpim1
@@ -303,5 +301,4 @@
  
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
-      CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
@@ -377,9 +374,13 @@
 
 ! Surface downward latent heat flux (CI_5)
-         IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
+         IF (ksbc == jp_flx) THEN
             DO jl=1,ncat
                ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
             ENDDO
-         ELSE
+         ELSE IF (ksbc == jp_purecpl) THEN
+            DO jl=1,ncat
+               ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
+            ENDDO
+	 ELSE
            !In coupled mode - qla_ice calculated in sbc_cpl for each category
            ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
Index: /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6249)
+++ /branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6250)
@@ -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
+               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 )
+               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 )
+               ! 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) 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 )
+               IF(lwp) 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
+               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 )
+               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 )
+               ! 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) 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 )
+               IF(lwp) 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
