Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 6754)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 6755)
@@ -162,4 +162,5 @@
       INTEGER             ::   iyear, imonth, iday
       REAL (wp)           ::   zsec
+      REAL (wp)           ::   zfjulday
       CHARACTER(len=256)  :: cl_path
       CHARACTER(len=256)  :: cl_filename
@@ -173,5 +174,7 @@
       IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
       IF ( ln_rstdate ) THEN
-         CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )           
+         zfjulday = fjulday + rdttra(1) / rday
+         IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
+         CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )           
          WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
       ELSE
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6754)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6755)
@@ -58,4 +58,5 @@
       INTEGER             ::   iyear, imonth, iday
       REAL (wp)           ::   zsec
+      REAL (wp)           ::   zfjulday
       !!
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
@@ -86,5 +87,7 @@
          IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
             IF ( ln_rstdate ) THEN
-               CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )           
+               zfjulday = fjulday + rdttra(1) / rday
+               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
+               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )           
                WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
             ELSE
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6754)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6755)
@@ -33,5 +33,7 @@
    USE cpl_oasis3      ! OASIS3 coupling
    USE geo2ocean       ! 
-   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
+   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev,            &
+                      CO2Flux_out_cpl, DMS_out_cpl, PCO2a_in_cpl, Dust_in_cpl, &
+                      ln_medusa
    USE albedo          !
    USE in_out_manager  ! I/O manager
@@ -49,4 +51,5 @@
    USE limthd_dh       ! for CALL lim_thd_snwblow
 #endif
+   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -105,5 +108,7 @@
    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 ::   jpr_atm_pco2 = 46          ! Incoming atm CO2 flux
+   INTEGER, PARAMETER ::   jpr_atm_dust = 47          ! Incoming atm aggregate dust 
+   INTEGER, PARAMETER ::   jprcv      = 47            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -140,5 +145,10 @@
    INTEGER, PARAMETER ::   jps_sstfrz = 32            ! sea-surface freezing temperature
    INTEGER, PARAMETER ::   jps_fice1  = 33            ! first-order ice concentration (for time-travelling ice coupling)
-   INTEGER, PARAMETER ::   jpsnd      = 33            ! total number of fields sended
+   INTEGER, PARAMETER ::   jps_bio_co2 = 34           ! MEDUSA air-sea CO2 flux in
+   INTEGER, PARAMETER ::   jps_bio_dms = 35           ! MEDUSA DMS surface concentration in
+   INTEGER, PARAMETER ::   jpsnd      = 35            ! total number of fields sent
+
+   REAL(wp), PARAMETER :: dms_unit_conv = 1.0e+6      ! Coversion factor to get outgong DMS in standard units for coupling
+	                                              ! i.e. specifically nmol/L (= umol/m3)
 
    !                                                         !!** namelist namsbc_cpl **
@@ -152,8 +162,11 @@
    ! Send to the atmosphere                           !
    TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2, sn_snd_cond, sn_snd_mpnd, sn_snd_sstfrz, sn_snd_thick1
+   TYPE(FLD_C) ::   sn_snd_bio_co2, sn_snd_bio_dms                       
 
    ! 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, sn_rcv_grnm, sn_rcv_antm
+   TYPE(FLD_C) ::   sn_rcv_atm_pco2, sn_rcv_atm_dust                         
+
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -230,4 +243,11 @@
          &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
       !!---------------------------------------------------------------------
+
+      ! Add MEDUSA related fields to namelist
+      NAMELIST/namsbc_cpl/  sn_snd_bio_co2, sn_snd_bio_dms,                                           &
+         &                  sn_rcv_atm_pco2, sn_rcv_atm_dust
+
+      !!---------------------------------------------------------------------
+
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
@@ -254,5 +274,5 @@
       ENDIF
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
-         WRITE(numout,*)'  received fields (mutiple ice categogies)'
+         WRITE(numout,*)'  received fields (mutiple ice categories)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
          WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
@@ -271,4 +291,6 @@
          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   ), ')'
+         WRITE(numout,*)'      atm pco2                        = ', TRIM(sn_rcv_atm_pco2%cldes), ' (', TRIM(sn_rcv_atm_pco2%clcat), ')'
+         WRITE(numout,*)'      atm dust                        = ', TRIM(sn_rcv_atm_dust%cldes), ' (', TRIM(sn_rcv_atm_dust%clcat), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
@@ -279,4 +301,6 @@
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
+         WRITE(numout,*)'      bio co2 flux                    = ', TRIM(sn_snd_bio_co2%cldes), ' (', TRIM(sn_snd_bio_co2%clcat), ')'
+         WRITE(numout,*)'      bio dms flux                    = ', TRIM(sn_snd_bio_dms%cldes), ' (', TRIM(sn_snd_bio_dms%clcat), ')'
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
          WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes   ), ' (', TRIM(sn_snd_cond%clcat   ), ')'
@@ -494,4 +518,20 @@
       !                                                      ! ------------------------- !
       srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
+
+
+      !                                                      ! --------------------------------------- !    
+      !                                                      ! Incoming CO2 and DUST fluxes for MEDUSA !
+      !                                                      ! --------------------------------------- !  
+      srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
+
+      IF (TRIM(sn_rcv_atm_pco2%cldes) == 'medusa') THEN
+        srcv(jpr_atm_pco2)%laction = .TRUE.
+      END IF
+               
+      srcv(jpr_atm_dust)%clname = 'OATMDUST'   
+      IF (TRIM(sn_rcv_atm_dust%cldes) == 'medusa')  THEN
+        srcv(jpr_atm_dust)%laction = .TRUE.
+      END IF
+	 
       !                                                      ! ------------------------- !
       !                                                      !   topmelt and botmelt     !   
@@ -733,5 +773,5 @@
       END SELECT
 #else
-      IF( TRIM( sn_snd_mpnd%cldes /= 'none' ) THEN
+      IF( TRIM( sn_snd_mpnd%cldes ) /= 'none' ) THEN
          CALL ctl_stop('Meltponds can only be used with CICEv5')
       ENDIF
@@ -770,4 +810,15 @@
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
       !
+
+      !                                                      ! ------------------------- !
+      !                                                      !   MEDUSA output fields    !
+      !                                                      ! ------------------------- !
+      ! Surface dimethyl sulphide from Medusa
+      ssnd(jps_bio_dms)%clname = 'OBioDMS'   
+      IF( TRIM(sn_snd_bio_dms%cldes) == 'medusa' )    ssnd(jps_bio_dms )%laction = .TRUE.
+
+      ! Surface CO2 flux from Medusa
+      ssnd(jps_bio_co2)%clname = 'OBioCO2'   
+      IF( TRIM(sn_snd_bio_co2%cldes) == 'medusa' )    ssnd(jps_bio_co2 )%laction = .TRUE.
       
       !                                                      ! ------------------------- !
@@ -990,4 +1041,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
       !!----------------------------------------------------------------------
+
+      ! RSRH temporary arrays for testing, just to recieve incoming MEDUSA related fields
+      ! until we know where they need to go.
+      REAL(wp), ALLOCATABLE :: atm_pco2(:,:)
+      REAL(wp), ALLOCATABLE :: atm_dust(:,:)
+
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
@@ -1128,4 +1185,9 @@
       ENDIF
 
+      IF (ln_medusa) THEN
+        IF( srcv(jpr_atm_pco2)%laction) PCO2a_in_cpl(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
+        IF( srcv(jpr_atm_dust)%laction) Dust_in_cpl(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
+      ENDIF
+
 #if defined key_cpl_carbon_cycle
       !                                                      ! ================== !
@@ -1273,11 +1335,17 @@
       zepsilon = rn_iceshelf_fluxes_tolerance
 
+
+      ! See if we need zmask_sum...
+      IF ( srcv(jpr_grnm)%laction .OR. srcv(jpr_antm)%laction ) THEN
+         zmask_sum = glob_sum( tmask(:,:,1) )
+      ENDIF
+
       IF( srcv(jpr_grnm)%laction ) THEN
          greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
          ! take average over ocean points of input array to avoid cumulative error over time
-         zgreenland_icesheet_mass_in = SUM( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zgreenland_icesheet_mass_in )
-         zmask_sum = SUM( tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
+
+         ! The following must be bit reproducible over different PE decompositions
+         zgreenland_icesheet_mass_in = glob_sum( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
+
          zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum
          greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
@@ -1304,8 +1372,7 @@
          antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
          ! take average over ocean points of input array to avoid cumulative error from rounding errors over time
-         zantarctica_icesheet_mass_in = SUM( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zantarctica_icesheet_mass_in )
-         zmask_sum = SUM( tmask(:,:,1) )
-         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
+         ! The following must be bit reproducible over different PE decompositions
+         zantarctica_icesheet_mass_in = glob_sum( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
+
          zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum
          antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
@@ -2184,4 +2251,5 @@
       ENDIF
       !
+#if defined key_cice && ! defined key_cice4
       ! Send meltpond fields 
       IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN
@@ -2226,4 +2294,5 @@
          IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info )
       ENDIF
+#endif
       !
       !
@@ -2235,4 +2304,20 @@
       !
 #endif
+
+
+
+      IF (ln_medusa) THEN
+      !                                                      ! --------------------------------- !
+      !                                                      !  CO2 flux and DMS from MEDUSA     ! 
+      !                                                      ! --------------------------------- !
+         IF ( ssnd(jps_bio_co2)%laction ) THEN
+            CALL cpl_snd( jps_bio_co2, isec, RESHAPE( CO2Flux_out_cpl, (/jpi,jpj,1/) ), info )
+         ENDIF
+
+         IF ( ssnd(jps_bio_dms)%laction )  THEN
+            CALL cpl_snd( jps_bio_dms, isec, RESHAPE( DMS_out_cpl, (/jpi,jpj,1/) ), info )
+         ENDIF
+      ENDIF
+
       !                                                      ! ------------------------- !
       IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
@@ -2407,6 +2492,8 @@
       IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
       
+#if defined key_cice
       ztmp1(:,:) = sstfrz(:,:) + rt0
       IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
+#endif
       !
       CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6754)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6755)
@@ -26,4 +26,5 @@
    USE zdfbfr
    USE fldread         ! read input field at current time step
+   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -257,7 +258,8 @@
               ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
               ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
-
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+ 
+              ! All related global sums must be done bit reproducibly
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
@@ -267,10 +269,11 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
 
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
@@ -280,6 +283,7 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
 
@@ -304,6 +308,7 @@
               ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
 
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+              ! All related global sums must be done bit reproducibly
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
@@ -313,10 +318,11 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
-               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
 
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                ! use ABS function because we need to preserve the sign of fwfisf
                WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
@@ -326,6 +332,7 @@
                ! check
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
-               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
-               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
                IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
 
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 6754)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 6755)
@@ -71,5 +71,19 @@
    !! Energy budget of the leads (open water embedded in sea ice)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraqsr_1lev        !: fraction of solar net radiation absorbed in the first ocean level [-]
+ 
+   !! Arrays used in coupling when MEDUSA is present. These arrays need to be declared
+   !! even if MEDUSA is not active, to allow compilation, in which case they will not be allocated. 
+   !! ---------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE :: CO2Flux_out_cpl(:,:)  ! Output coupling CO2 flux  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE :: DMS_out_cpl(:,:)      ! Output coupling DMS  
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE :: PCO2a_in_cpl(:,:)     ! Input coupling CO2 partial pressure 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE :: Dust_in_cpl(:,:)      ! Input coupling dust 
+
+#if defined key_medusa
+   LOGICAL, PUBLIC, PARAMETER :: ln_medusa=.TRUE. ! Medusa switched on or off. 
+#else
+   LOGICAL, PUBLIC, PARAMETER :: ln_medusa=.FALSE. ! Medusa switched on or off. 
+#endif
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -119,4 +133,14 @@
       ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(4) )
          !
+#if defined key_oasis3
+      IF (ln_medusa) THEN
+         ! We only actually need these arrays to be allocated if coupling and MEDUSA 
+         ! are enabled
+         ALLOCATE( CO2Flux_out_cpl(jpi,jpj), DMS_out_cpl(jpi,jpj),               &
+                   PCO2a_in_cpl(jpi,jpj), Dust_in_cpl(jpi,jpj),     STAT=ierr(5) )
+
+      ENDIF
+#endif
+
       oce_alloc = MAXVAL( ierr )
       IF( oce_alloc /= 0 )   CALL ctl_warn('oce_alloc: failed to allocate arrays')
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6754)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6755)
@@ -33,4 +33,5 @@
    USE step_oce         ! time stepping definition modules
    USE iom
+   USE lbclnk
 
    IMPLICIT NONE
@@ -73,4 +74,5 @@
       !!----------------------------------------------------------------------
       INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   tind     ! tracer loop index
       INTEGER ::   indic    ! error indicator if < 0
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
@@ -109,4 +111,10 @@
                          CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
       ENDIF
+
+      ! We must ensure that tsb halos are up to date on EVERY timestep.
+      DO tind = 1, jpts
+         CALL lbc_lnk( tsb(:,:,:,tind), 'T', 1. )
+      END DO
+
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
