Index: NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/ICB/icbclv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/ICB/icbclv.F90	(revision 13777)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/ICB/icbclv.F90	(revision 13778)
@@ -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,6 @@
       !
       REAL(wp)      ::   zcalving_used, zdist, zfact
+      REAL(wp), DIMENSION(1)      ::   zgreenland_calving_sum, zantarctica_calving_sum 
+      LOGICAL       ::   ll_write
       INTEGER       ::   jn, ji, jj                    ! loop counters
       INTEGER       ::   imx                           ! temporary integer for max berg class
@@ -63,4 +69,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 .AND. ((nn_print>0))
+          ! 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(1) = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
+           IF( lk_mpp ) CALL mpp_sum( 'icbclv', 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) + 1.0e-10_wp )
+
+           ! check
+           IF(ll_write) WRITE(numout, *) 'Greenland iceberg calving climatology (kg/s) : ',zgreenland_calving_sum(1)
+           zgreenland_calving_sum(1) = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
+           IF( lk_mpp ) CALL mpp_sum( 'icbclv', zgreenland_calving_sum )
+           IF(ll_write) WRITE(numout, *) 'Greenland iceberg calving adjusted value (kg/s) : ',zgreenland_calving_sum(1)
+
+           zantarctica_calving_sum(1) = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
+           IF( lk_mpp ) CALL mpp_sum( 'icbclv', 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) + 1.0e-10_wp )
+
+           ! check
+           IF(ll_write) WRITE(numout, *) 'Antarctica iceberg calving climatology (kg/s) : ',zantarctica_calving_sum(1)
+           zantarctica_calving_sum(1) = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
+           IF( lk_mpp ) CALL mpp_sum( 'icbclv', zantarctica_calving_sum )
+           IF(ll_write) WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum(1)
+
+        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.3_icesheet_and_river_coupling/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/IOM/restart.F90	(revision 13777)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/IOM/restart.F90	(revision 13778)
@@ -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.3_icesheet_and_river_coupling/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/cpl_oasis3.F90	(revision 13777)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/cpl_oasis3.F90	(revision 13778)
@@ -29,10 +29,19 @@
 #endif
    USE par_oce                      ! ocean parameters
+   USE cpl_rnf_1d, ONLY: nn_cpl_river   ! Variables used in 1D river outflow 
    USE dom_oce                      ! ocean space and time domain
    USE in_out_manager               ! I/O manager
    USE lbclnk                       ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp
 
    IMPLICIT NONE
    PRIVATE
+
+#if ! defined key_oasis3 
+   ! Dummy interface to oasis_get if not using oasis 
+   INTERFACE oasis_get 
+      MODULE PROCEDURE oasis_get_1d, oasis_get_2d 
+   END INTERFACE 
+#endif 
 
    PUBLIC   cpl_init
@@ -40,7 +49,13 @@
    PUBLIC   cpl_snd
    PUBLIC   cpl_rcv
+   PUBLIC   cpl_rcv_1d
    PUBLIC   cpl_freq
    PUBLIC   cpl_finalize
 
+#if defined key_mpp_mpi
+   INCLUDE 'mpif.h'
+#endif
+
+   INTEGER, PARAMETER         ::   localRoot  = 0
    INTEGER, PUBLIC            ::   OASIS_Rcv  = 1    !: return code if received field
    INTEGER, PUBLIC            ::   OASIS_idle = 0    !: return code if nothing done by oasis
@@ -66,5 +81,5 @@
    INTEGER                    ::   nsnd         ! total number of fields sent 
    INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
-   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=60   ! Maximum number of coupling fields
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=61   ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields
@@ -81,4 +96,5 @@
       INTEGER               ::   nct       ! Number of categories in field
       INTEGER               ::   ncplmodel ! Maximum number of models to/from which this variable may be sent/received
+      INTEGER               ::   dimensions ! Number of dimensions of coupling field 
    END TYPE FLD_CPL
 
@@ -105,4 +121,5 @@
       CHARACTER(len = *), INTENT(in   ) ::   cd_modname   ! model name as set in namcouple file
       INTEGER           , INTENT(  out) ::   kl_comm      ! local communicator of the model
+      INTEGER                           ::   error
       !!--------------------------------------------------------------------
 
@@ -141,6 +158,16 @@
       !
       INTEGER :: id_part
+      INTEGER :: id_part_0d     ! Partition for 0d fields 
+      INTEGER :: id_part_rnf_1d ! Partition for 1d river outflow fields 
+      INTEGER :: id_part_temp   ! Temperary partition used to choose either 0d or 1d partitions 
       INTEGER :: paral(5)       ! OASIS3 box partition
-      INTEGER :: ishape(4)    ! shape of arrays passed to PSMILe
+      INTEGER :: ishape(4)      ! shape of 2D arrays passed to PSMILe
+      INTEGER :: ishape0d1d(2)  ! Shape of 0D or 1D arrays passed to PSMILe.
+      INTEGER :: var_nodims(2)  ! Number of coupling field dimensions.
+                                ! var_nodims(1) is redundant from OASIS3-MCT vn4.0 onwards
+                                ! but retained for backward compatibility. 
+                                ! var_nodims(2) is the number of fields in a bundle 
+                                ! or 1 for unbundled fields (bundles are not yet catered for
+                                ! in NEMO hence we default to 1).  
       INTEGER :: ji,jc,jm       ! local loop indicees
       CHARACTER(LEN=64) :: zclname
@@ -185,4 +212,7 @@
       ishape(3) = 1
       ishape(4) = nlej-nldj+1
+
+      ishape0d1d(1) = 0 
+      ishape0d1d(2) = 0 
       !
       ! ... Allocate memory for data exchange
@@ -211,4 +241,34 @@
    
       CALL oasis_def_partition ( id_part, paral, nerror, jpiglo*jpjglo )
+
+      ! A special partition is needed for 0D fields
+     
+      paral(1) = 0                                       ! serial partitioning
+      paral(2) = 0   
+      IF ( nproc == 0) THEN
+         paral(3) = 1                   ! Size of array to couple (scalar)
+      ELSE
+         paral(3) = 0                   ! Dummy size for PE's not involved
+      END IF
+      paral(4) = 0
+      paral(5) = 0
+	     
+      CALL oasis_def_partition ( id_part_0d, paral, nerror, 1 )
+
+      ! Another partition is needed for 1D river routing fields
+      
+      paral(1) = 0                                       ! serial partitioning
+      paral(2) = 0   
+      IF ( nproc == 0) THEN
+         paral(3) = nn_cpl_river                   ! Size of array to couple (vector)
+      ELSE
+         paral(3) = 0                   ! Dummy size for PE's not involved
+      END IF
+      paral(4) = 0
+      paral(5) = 0
+
+
+      CALL oasis_def_partition ( id_part_rnf_1d, paral, nerror, nn_cpl_river )
+
       !
       ! ... Announce send variables. 
@@ -289,6 +349,41 @@
 #endif
                   IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In
-                  CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
-                     &                OASIS_In           , ishape , OASIS_REAL, nerror )
+                  flush(numout)
+
+                  ! Define 0D (Greenland or Antarctic ice mass) or 1D (river outflow) coupling fields
+                  IF (srcv(ji)%dimensions <= 1) THEN
+                    var_nodims(1) = 1
+                    var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
+                    IF (nproc == 0) THEN
+
+                       IF (srcv(ji)%dimensions == 0) THEN
+
+                          ! If 0D then set temporary variables to 0D components
+                          id_part_temp = id_part_0d
+                          ishape0d1d(2) = 1
+                       ELSE
+
+                          ! If 1D then set temporary variables to river outflow components
+                          id_part_temp = id_part_rnf_1d
+                          ishape0d1d(2)= nn_cpl_river
+
+                       END IF
+
+                       CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_temp   , var_nodims,   &
+	                                OASIS_In           , ishape0d1d(1:2) , OASIS_REAL, nerror )
+                    ELSE
+                       ! Dummy call to keep OASIS3-MCT happy.
+                       CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_0d   , var_nodims,   &
+	                                OASIS_In           , ishape0d1d(1:2) , OASIS_REAL, nerror )
+                    END IF
+                  ELSE
+                    ! It's a "normal" 2D (or pseudo 3D) coupling field. 
+                    ! ... Set the field dimension and bundle count
+                    var_nodims(1) = 2
+                    var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
+
+                    CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , var_nodims,   &
+                       &                OASIS_In           , ishape , OASIS_REAL, nerror )
+                  ENDIF
                   IF ( nerror /= OASIS_Ok ) THEN
                      WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
@@ -471,4 +566,111 @@
 
 
+   SUBROUTINE cpl_rcv_1d( kid, kstep, pdata, nitems, kinfo )
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE cpl_rcv_1d  ***
+      !!
+      !! ** Purpose : - A special version of cpl_rcv to deal exclusively with
+      !! receipt of 0D or 1D fields. 
+      !! The fields are recieved into a 1D array buffer which is simply a 
+      !! dynamically sized sized array (which may be of size 1)
+      !! of 0 dimensional fields. This allows us to pass miltiple 0D 
+      !! fields via a single put/get operation.  
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in   ) ::   nitems      ! Number of 0D items to recieve 
+                                               ! during this get operation. i.e.
+                                               ! The size of the 1D array in which
+                                               ! 0D items are passed.   
+      INTEGER , INTENT(in   ) ::   kid         ! ID index of the incoming
+                                               ! data.  
+      INTEGER , INTENT(in   ) ::   kstep       ! ocean time-step in seconds
+      REAL(wp), INTENT(inout) ::   pdata(1:nitems) ! The original value(s),  
+                                                   ! unchanged if nothing is 
+                                                   ! received
+      INTEGER , INTENT(  out) ::   kinfo       ! OASIS3 info argument
+      !! 
+      REAL(wp) ::   recvfld(1:nitems)          ! Local receive field buffer
+      INTEGER  ::   jc,jm     ! local loop index
+      INTEGER  ::   ierr
+      LOGICAL  ::   llaction
+      INTEGER  ::   MPI_WORKING_PRECISION
+      INTEGER  ::   number_to_print 
+      !!--------------------------------------------------------------------
+      !
+      ! receive local data from OASIS3 on every process
+      !
+      kinfo = OASIS_idle
+      !
+      ! 0D and 1D fields won't have categories or any other form of "pseudo level" 
+      ! so we only cater for a single set of values and thus don't bother 
+      ! with a loop over the jc index 
+      jc = 1
+
+      DO jm = 1, srcv(kid)%ncplmodel
+
+         IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
+
+            IF ( ( srcv(kid)%dimensions <= 1) .AND. (nproc == 0) ) THEN
+               ! Since there is no concept of data decomposition for zero 
+               ! dimension fields, they must only be exchanged through the master PE, 
+               ! unlike "normal" 2D field cases where every PE is involved. 
+
+               CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, recvfld, kinfo )   
+               
+               llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   &
+                           kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
+               
+               IF ( ln_ctl ) WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , &
+                                     llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
+               
+               IF ( llaction ) THEN
+                  
+                  kinfo = OASIS_Rcv
+                  pdata(1:nitems) = recvfld(1:nitems) 
+                  
+                  IF ( ln_ctl ) THEN        
+                     number_to_print = 10
+                     IF ( nitems < number_to_print ) number_to_print = nitems
+                     WRITE(numout,*) '****************'
+                     WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname
+                     WRITE(numout,*) 'oasis_get: ivarid '  , srcv(kid)%nid(jc,jm)
+                     WRITE(numout,*) 'oasis_get:   kstep', kstep
+                     WRITE(numout,*) 'oasis_get:   info ', kinfo
+                     WRITE(numout,*) '     - Minimum Value is ', MINVAL(pdata(:))
+                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(:))
+                     WRITE(numout,*) '     - Start of data is ', pdata(1:number_to_print)
+                     WRITE(numout,*) '****************'
+                  ENDIF
+                  
+               ENDIF
+            ENDIF   
+          ENDIF
+            
+       ENDDO
+
+#if defined key_mpp_mpi
+       ! Set the precision that we want to broadcast using MPI_BCAST
+       SELECT CASE( wp )
+       CASE( sp ) 
+         MPI_WORKING_PRECISION = MPI_REAL                ! Single precision
+       CASE( dp )
+         MPI_WORKING_PRECISION = MPI_DOUBLE_PRECISION    ! Double precision
+       CASE default
+         CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Could not find precision for coupling 0d or 1d field" )
+       END SELECT
+
+       ! We have to broadcast (potentially) received values from PE 0 to all 
+       ! the others. If no new data has been received we're just 
+       ! broadcasting the existing values but there's no more efficient way 
+       ! to deal with that w/o NEMO adopting a UM-style test mechanism
+       ! to determine active put/get timesteps. 
+       CALL mpi_bcast( pdata, nitems, MPI_WORKING_PRECISION, localRoot, mpi_comm_oce, ierr )
+#else
+       CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Unable to use mpi_bcast without key_mpp_mpi present. Please add key_mpp_mpi to your list of NEMO keys." )
+#endif
+
+      !
+   END SUBROUTINE cpl_rcv_1d
+
+
    INTEGER FUNCTION cpl_freq( cdfieldname )  
       !!---------------------------------------------------------------------
@@ -578,5 +780,5 @@
    SUBROUTINE oasis_def_var(k1,cd1,k2,k3,k4,k5,k6,k7)
       CHARACTER(*), INTENT(in   ) ::  cd1
-      INTEGER     , INTENT(in   ) ::  k2,k3(2),k4,k5(2,2),k6
+      INTEGER     , INTENT(in   ) ::  k2,k3(2),k4,k5(*),k6
       INTEGER     , INTENT(  out) ::  k1,k7
       k1 = -1 ; k7 = -1
@@ -598,11 +800,19 @@
    END SUBROUTINE oasis_put
 
-   SUBROUTINE oasis_get(k1,k2,p1,k3)
+   SUBROUTINE oasis_get_1d(k1,k2,p1,k3)
+      REAL(wp), DIMENSION(:)  , INTENT(  out) ::  p1
+      INTEGER                 , INTENT(in   ) ::  k1,k2
+      INTEGER                 , INTENT(  out) ::  k3
+      p1(1) = -1. ; k3 = -1
+      WRITE(numout,*) 'oasis_get_1d: Error you sould not be there...'
+   END SUBROUTINE oasis_get_1d
+
+   SUBROUTINE oasis_get_2d(k1,k2,p1,k3)
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::  p1
       INTEGER                 , INTENT(in   ) ::  k1,k2
       INTEGER                 , INTENT(  out) ::  k3
       p1(1,1) = -1. ; k3 = -1
-      WRITE(numout,*) 'oasis_get: Error you sould not be there...'
-   END SUBROUTINE oasis_get
+      WRITE(numout,*) 'oasis_get_2d: Error you sould not be there...'
+   END SUBROUTINE oasis_get_2d
 
    SUBROUTINE oasis_get_freqs(k1,k5,k2,k3,k4)
Index: NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/cpl_rnf_1d.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/cpl_rnf_1d.F90	(revision 13778)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/cpl_rnf_1d.F90	(revision 13778)
@@ -0,0 +1,205 @@
+MODULE cpl_rnf_1d
+   !!======================================================================
+   !!                       ***  MODULE  cpl_rnf_1d  ***
+   !! Ocean forcing:  River runoff passed from the atmosphere using
+   !!                 1D array. One value per river.
+   !!=====================================================================
+   !! History : ?.?  ! 2018-01 (D. Copsey) Initial setup
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   cpl_rnf_1d_init : runoffs initialisation
+   !!---------------------------------------------------------------------- 
+
+#if defined key_oasis3
+   USE mod_oasis                    ! OASIS3-MCT module
+#endif
+   USE timing          ! Timing
+   USE in_out_manager  ! I/O units
+   USE lib_mpp         ! MPP library
+   USE iom
+   USE dom_oce         ! Domain sizes (for grid box area e1e2t)
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE lib_fortran,    ONLY: DDPDD
+   
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   cpl_rnf_1d_init     ! routine called in nemo_init
+   PUBLIC   cpl_rnf_1d_to_2d      ! routine called in sbccpl.F90
+   
+   TYPE, PUBLIC ::   RIVERS_DATA     !: Storage for river outflow data
+      INTEGER, ALLOCATABLE, DIMENSION(:,:)    ::   river_number       !: River outflow number
+      REAL(wp), ALLOCATABLE, DIMENSION(:)     ::   river_area         ! 1D array listing areas of each river outflow (m2)
+      COMPLEX(wp), ALLOCATABLE, DIMENSION(:)  ::   river_area_c       ! Comlex version of river_area for use in bit reproducible sums (m2)
+   END TYPE RIVERS_DATA
+   
+   TYPE(RIVERS_DATA), PUBLIC, TARGET :: rivers  !: River data
+   
+   INTEGER, PUBLIC            :: nn_cpl_river   ! Maximum number of rivers being passed through the coupler
+   INTEGER, PUBLIC            :: runoff_id      ! OASIS coupling id used in oasis_get command
+   LOGICAL                    :: ln_print_river_info  ! Diagnostic prints of river coupling information
+   
+CONTAINS
+
+   SUBROUTINE cpl_rnf_1d_init
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE cpl_rnf_1d_init  ***
+      !!                     
+      !! ** Purpose : - Read in file for river outflow numbers.
+      !!                Calculate 2D area of river outflow points.
+      !!                Called from nemo_init (nemogcm.F90).
+      !!
+      !!----------------------------------------------------------------------
+      !! namelist variables
+      !!-------------------
+      CHARACTER(len=200)                        ::   file_riv_number             !: Filename for river numbers
+      INTEGER                                   ::   ios                 ! Local integer output status for namelist read
+      INTEGER                                   ::   inum
+      INTEGER                                   ::   ii, jj, rr          !: Loop indices
+      INTEGER                                   ::   max_river
+      REAL(wp), DIMENSION(jpi,jpj)              ::   river_number        ! 2D array containing the river outflow numbers
+      
+      NAMELIST/nam_cpl_rnf_1d/file_riv_number, nn_cpl_river, ln_print_river_info
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing ) CALL timing_start('cpl_rnf_1d_init')
+      
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'cpl_rnf_1d_init : initialization of river runoff coupling'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      IF(lwp) CALL flush(numout)
+      
+      REWIND(numnam_cfg)
+      
+      ! Read the namelist
+      READ  ( numnam_ref, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in reference namelist' )
+      READ  ( numnam_cfg, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in configuration namelist' )
+      IF(lwm .AND. nprint > 2) WRITE ( numond, nam_cpl_rnf_1d )
+
+      !                                               ! Parameter control and print
+      IF(lwp) WRITE(numout,*) '  '
+      IF(lwp) WRITE(numout,*) '          Namelist nam_cpl_rnf_1d : Coupled runoff using 1D array'
+      IF(lwp) WRITE(numout,*) '             Input file that contains river numbers = ',file_riv_number
+      IF(lwp) WRITE(numout,*) '             Maximum number of rivers to couple = ',nn_cpl_river
+      IF(lwp) WRITE(numout,*) '             Print river information = ',ln_print_river_info
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) CALL flush(numout)
+
+      ! Assign space for river numbers
+      ALLOCATE( rivers%river_number( jpi, jpj ) )
+      
+      ! Read the river numbers from netcdf file
+      CALL iom_open (file_riv_number , inum )
+      CALL iom_get  ( inum, jpdom_data, 'river_number', river_number )
+      CALL iom_close( inum )
+      
+      ! Convert from a real array to an integer array
+      max_river=0
+      DO ii = 1, jpi
+        DO jj = 1, jpj
+          rivers%river_number(ii,jj) = INT(river_number(ii,jj))
+          IF ( rivers%river_number(ii,jj) > max_river ) THEN
+            max_river = rivers%river_number(ii,jj)
+          END IF
+        END DO
+      END DO
+      
+      ! Print out the largest river number
+      IF ( ln_print_river_info .AND. lwp) THEN
+         WRITE(numout,*) 'Maximum river number in input file = ',max_river
+         CALL flush(numout)
+      END IF
+      
+      ! Get the area of each river outflow
+      ALLOCATE( rivers%river_area( nn_cpl_river ) )
+      ALLOCATE( rivers%river_area_c( nn_cpl_river ) )
+      rivers%river_area_c(:) = CMPLX( 0.e0, 0.e0, wp )
+      DO ii = nldi, nlei      
+        DO jj = nldj, nlej
+          IF ( tmask_i(ii,jj) > 0.5 ) THEN  ! This makes sure we are not at a duplicated point (at north fold or east-west cyclic point)
+            IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN
+              ! Add the area of each grid box (e1e2t) into river_area_c using DDPDD which should maintain bit reproducibility (needs to be checked)
+              CALL DDPDD( CMPLX( e1e2t(ii,jj), 0.e0, wp ), rivers%river_area_c(rivers%river_number(ii,jj) ) )
+            END IF
+          END IF
+        END DO
+      END DO
+      
+      ! Use mpp_sum to add together river areas on other processors
+      CALL mpp_sum( 'cpl_rnf_1d', rivers%river_area_c )
+      
+      ! Convert from complex number to real
+      rivers%river_area(:) = REAL(rivers%river_area_c(:),wp)
+      
+      IF ( ln_print_river_info .AND. lwp) THEN
+        WRITE(numout,*) 'Area of rivers 1 to 10 are ',rivers%river_area(1:10)
+        WRITE(numout,*) 'Maximum river area = ',MAXVAL(rivers%river_area)
+        WRITE(numout,*) 'Minimum river area = ',MINVAL(rivers%river_area)
+        WRITE(numout,*) 'Sum of river areas = ',SUM(rivers%river_area)
+        CALL flush(numout)
+      END IF
+
+      IF ( MINVAL(rivers%river_area) <= 0 ) THEN
+         WRITE(numout,*) 'ERROR: There is at least one river with a river outflow area of zero. Please check your file_riv_number file'
+         WRITE(numout,*) 'that all the allocated river numbers are at ocean points as defined by the bathymetry file with no river'
+         WRITE(numout,*) 'numbers within the north fold or wraparound points.'
+         DO rr = 1,nn_cpl_river
+           IF ( rivers%river_area(rr) <= 0 ) THEN
+             WRITE(numout,*) 'The first river with this problem is river number ',rr
+             CALL ctl_stop ( 'STOP', 'ERROR: There is at least one river with a river outflow area of zero. See stdout.')
+           END IF
+         END DO
+      END IF
+      
+   END SUBROUTINE cpl_rnf_1d_init
+   
+   SUBROUTINE cpl_rnf_1d_to_2d( runoff_1d )
+   
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE cpl_rnf_1d_to_2d  ***
+      !!                     
+      !! ** Purpose : - Convert river outflow from 1D array (passed from the
+      !!                atmosphere) to the 2D NEMO runoff field.
+      !!                Called from sbc_cpl_ice_flx (sbccpl.F90).
+      !!
+      !!----------------------------------------------------------------------
+      
+      REAL                   , INTENT(in   ) ::   runoff_1d(nn_cpl_river)    ! River runoff. One value per river.
+      
+      INTEGER  ::   ii, jj, rr                 ! Loop indices
+      LOGICAL  ::   found_nan
+            
+      ! Convert the 1D total runoff per river to 2D runoff flux by
+      ! dividing by the area of each runoff zone.
+      DO ii = 1, jpi
+        DO jj = 1, jpj
+          IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river .AND. tmask_i(ii,jj) > 0.5 ) THEN
+            rnf(ii,jj) = runoff_1d(rivers%river_number(ii,jj)) / rivers%river_area(rivers%river_number(ii,jj))
+          ELSE
+            rnf(ii,jj) = 0.0
+          END IF
+            
+        END DO
+      END DO
+
+      IF ( ln_print_river_info ) THEN
+        WRITE(numout,*) 'SUM of river runoff on 1D array = ',SUM(runoff_1d)
+        WRITE(numout,*) 'SUM of river runoff on 2D array = ',SUM(rnf)
+        found_nan = .false.
+        DO rr = 1, nn_cpl_river
+          IF (runoff_1d(rr) /= runoff_1d(rr)) THEN
+            IF (.NOT. found_nan) THEN
+              WRITE(numout,*) 'WARNING: NAN found at river number ',rr
+              found_nan = .true.
+            END IF
+          END IF
+        END DO
+      END IF     
+
+
+   END SUBROUTINE cpl_rnf_1d_to_2d
+
+END MODULE cpl_rnf_1d
Index: NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbc_oce.F90	(revision 13777)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbc_oce.F90	(revision 13778)
@@ -137,4 +137,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra         !: cloud cover (fraction of cloud in a gridcell) [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_icesheet_mask 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mask 
 
    !!----------------------------------------------------------------------
@@ -150,4 +152,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"
@@ -183,4 +212,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_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) ) 
          !
       ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
Index: NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90	(revision 13777)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90	(revision 13778)
@@ -36,4 +36,5 @@
    USE eosbn2         ! 
    USE sbcrnf  , ONLY : l_rnfcpl
+   USE cpl_rnf_1d, ONLY: nn_cpl_river, cpl_rnf_1d_init, cpl_rnf_1d_to_2d   ! Variables used in 1D river outflow 
    USE sbcisf  , ONLY : l_isfcpl
 #if defined key_cice
@@ -120,6 +121,10 @@
    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 ::   jpr_rnf_1d = 60   ! 1D river runoff 
+   INTEGER, PARAMETER ::   jpr_qtr    = 61   ! Transmitted solar
+
+   INTEGER, PARAMETER ::   jprcv      = 61   ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -186,5 +191,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 
@@ -271,5 +277,10 @@
          &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
          &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
-         &                  sn_rcv_ts_ice
+         &                  sn_rcv_ts_ice, sn_rcv_grnm  , sn_rcv_antm  ,                               &
+         &                  nn_coupled_iceshelf_fluxes  , ln_iceshelf_init_atmos ,                     &
+         &                  rn_greenland_total_fw_flux  , rn_greenland_calving_fraction  ,             &
+         &                  rn_antarctica_total_fw_flux , rn_antarctica_calving_fraction ,             &
+         &                  rn_iceshelf_fluxes_tolerance
+
       !!---------------------------------------------------------------------
       !
@@ -310,4 +321,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   ), ')'
@@ -345,4 +358,11 @@
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor 
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
+         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
 
@@ -360,5 +380,9 @@
 
       ! default definitions of srcv
-      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
+      srcv(:)%laction = .FALSE. 
+      srcv(:)%clgrid = 'T' 
+      srcv(:)%nsgn = 1. 
+      srcv(:)%nct = 1 
+      srcv(:)%dimensions = 2 
 
       !                                                      ! ------------------------- !
@@ -479,6 +503,11 @@
       !                                                      ! ------------------------- !
       srcv(jpr_rnf   )%clname = 'O_Runoff'
-      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
-         srcv(jpr_rnf)%laction = .TRUE.
+      srcv(jpr_rnf_1d   )%clname = 'ORunff1D' 
+      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' .OR. TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN  
+         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) srcv(jpr_rnf)%laction = .TRUE. 
+         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN 
+            srcv(jpr_rnf_1d)%laction = .TRUE. 
+            srcv(jpr_rnf_1d)%dimensions = 1 ! 1D field passed through coupler 
+         END IF 
          l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
          ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
@@ -487,5 +516,19 @@
       ENDIF
       !
-      srcv(jpr_cal)%clname = 'OCalving'   ;  IF( TRIM( sn_rcv_cal%cldes) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
+      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' ) THEN 
+         srcv(jpr_grnm)%laction = .TRUE.  
+         srcv(jpr_grnm)%dimensions = 0 ! Scalar field
+      ENDIF 
+       
+      srcv(jpr_antm  )%clname = 'OAntmass' 
+      IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' ) THEN
+         srcv(jpr_antm)%laction = .TRUE.
+         srcv(jpr_antm)%dimensions = 0 ! Scalar field
+      ENDIF
+
       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.
@@ -738,23 +781,4 @@
          ENDIF
       ENDIF
-      
-      ! =================================================== !
-      ! Allocate all parts of frcv used for received fields !
-      ! =================================================== !
-      DO jn = 1, jprcv
-         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
-      END DO
-      ! Allocate taum part of frcv which is used even when not received as coupling field
-      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
-      ! Allocate w10m part of frcv which is used even when not received as coupling field
-      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
-      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
-      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
-      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
-      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
-      IF( k_ice /= 0 ) THEN
-         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
-         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
-      END IF
 
       ! ================================ !
@@ -766,5 +790,9 @@
       
       ! default definitions of nsnd
-      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
+      ssnd(:)%laction = .FALSE. 
+      ssnd(:)%clgrid = 'T' 
+      ssnd(:)%nsgn = 1. 
+      ssnd(:)%nct = 1 
+      ssnd(:)%dimensions = 2 
          
       !                                                      ! ------------------------- !
@@ -1049,4 +1077,42 @@
       ENDIF
 
+      ! Initialise 1D river outflow scheme 
+      nn_cpl_river = 1 
+      IF ( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) CALL cpl_rnf_1d_init   ! Coupled runoff using 1D array
+      
+      ! =================================================== !
+      ! Allocate all parts of frcv used for received fields !
+      ! =================================================== !
+      DO jn = 1, jprcv
+
+         IF ( srcv(jn)%laction ) THEN 
+            SELECT CASE( srcv(jn)%dimensions )
+            !
+            CASE( 0 )   ! Scalar field
+               ALLOCATE( frcv(jn)%z3(1,1,1) )
+               
+            CASE( 1 )   ! 1D field
+               ALLOCATE( frcv(jn)%z3(nn_cpl_river,1,1) )
+               
+            CASE DEFAULT  ! 2D (or pseudo 3D) field.
+               ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
+               
+            END SELECT
+         END IF
+
+      END DO
+      ! Allocate taum part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
+      ! Allocate w10m part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
+      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
+      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
+      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
+      IF( k_ice /= 0 ) THEN
+         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
+         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
+      END IF
+
       !
       ! ================================ !
@@ -1066,4 +1132,23 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
+      !
+      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 
+  
+          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
@@ -1126,9 +1211,15 @@
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
       REAL(wp) ::   zcoef                  ! temporary scalar
+      LOGICAL  ::   ll_wrtstp              ! write diagnostics?
       REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
+      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) ::   zzx, zzy               ! temporary variables
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
       !!----------------------------------------------------------------------
+      !
+      ll_wrtstp  = (( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )) .AND. (nn_print>0)
       !
       IF( kt == nit000 ) THEN
@@ -1147,5 +1238,13 @@
       isec = ( kt - nit000 ) * NINT( rdt )                      ! date of exchanges
       DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
-         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
+        IF( srcv(jn)%laction ) THEN  
+ 
+          IF ( srcv(jn)%dimensions <= 1 ) THEN 
+            CALL cpl_rcv_1d( jn, isec, frcv(jn)%z3, SIZE(frcv(jn)%z3), nrcvinfo(jn) ) 
+          ELSE 
+            CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) ) 
+          END IF 
+
+        END IF 
       END DO
 
@@ -1474,4 +1573,80 @@
          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
+         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
+         ENDIF
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
       ENDIF
       !
@@ -1746,6 +1921,9 @@
       
       ! --- Continental fluxes --- !
-      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
+      IF( srcv(jpr_rnf)%laction ) THEN   ! 2D runoffs (included in emp later on)
          rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      ENDIF
+      IF( srcv(jpr_rnf_1d)%laction ) THEN ! 1D runoff
+         CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:))
       ENDIF
       IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot and emp_oce)
@@ -1786,6 +1964,9 @@
       zsnw(:,:) = picefr(:,:)
       ! --- Continental fluxes --- !
-      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
+      IF( srcv(jpr_rnf)%laction ) THEN   ! 2D runoffs (included in emp later on)
          rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      ENDIF
+      IF( srcv(jpr_rnf_1d)%laction ) THEN  ! 1D runoff
+         CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:)) 
       ENDIF
       IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot)
Index: NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90	(revision 13777)
+++ NEMO/branches/UKMO/NEMO_4.0.3_icesheet_and_river_coupling/src/OCE/SBC/sbcisf.F90	(revision 13778)
@@ -92,9 +92,12 @@
       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
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zfwfisf3d, zqhcisf3d, zqlatisf3d
+      LOGICAL :: ll_wrtstp
       !!---------------------------------------------------------------------
       !
+      ll_wrtstp  = (( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )) .AND. (nn_print>0)
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN    ! compute salt and heat flux
          !
@@ -127,4 +130,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( 'sbcisf', 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( 'sbcisf', 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( 'sbcisf', 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( 'sbcisf', 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 +181,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( 'sbcisf', 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( 'sbcisf', 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( 'sbcisf', 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( 'sbcisf', 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
