Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 9157)
@@ -448,4 +448,6 @@
        <field id= "DMS_HALL"   long_name="DMS Surface Concentration, Halloran"       unit="nmol/L"      />
        <field id= "DMS_ANDM"   long_name="DMS Surface Concentration, Anderson modif" unit="nmol/L"      />
+       <field id= "CHL_MLD"    long_name="MLD averaged Chlorophyll"                  unit="mg Chl/m3"   />
+       <field id= "CHL_CPL"    long_name="Scaled Surf or MLD Chlorophyll to UM"      unit="kg Chl/m3"   />
        <field id= "ATM_XCO2"   long_name="Atmospheric xCO2"                          unit="ppm"         />
        <field id= "OCN_FCO2"   long_name="Surface ocean fCO2"                        unit="uatm"        />
@@ -702,5 +704,4 @@
       <field field_ref= "OCN_CO3"    name="OCN_CO3"    />
       <field field_ref= "CO2FLUX"    name="CO2FLUX"    />
-      <field field_ref= "FGCO2"      name="FGCO2"      />
       <field field_ref= "OM_CAL"     name="OM_CAL"     />
       <field field_ref= "OM_ARG"     name="OM_ARG"     />
@@ -784,4 +785,13 @@
       <field field_ref= "CO2STARAIR" name="CO2STARAIR" />
       <field field_ref= "OCN_DPCO2"  name="OCN_DPCO2"  />
+      <field field_ref= "CHL_MLD"    name="CHL_MLD"    />
+    </field_group>
+
+    <field_group id="groupMEDUSA_cpl" >
+      <field field_ref= "CHL_CPL"    name="CHL_CPL"    />
+      <field field_ref= "FGCO2"      name="FGCO2_CPL"  />
+      <field field_ref= "DMS_SURF"   name="DMS_CPL"    />
+      <field field_ref= "ATM_XCO2"   name="AXCO2_CPL"  />
+      <field field_ref= "AEOLIAN"    name="DUST_CPL"   />
     </field_group>
 
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 9157)
@@ -84,4 +84,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
@@ -93,5 +94,4 @@
       REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse
       REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
 
 
@@ -130,5 +130,4 @@
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
-            zvn(:,:,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
@@ -138,5 +137,4 @@
                      zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
                      zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                     zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc
                   ENDDO
                ENDDO
@@ -151,5 +149,5 @@
              tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1)
              sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1)
-             v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) )
+             v_msf(:,:,1) = ptr_sjk( pvtr(:,:,:) )
 
              htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 )
@@ -177,5 +175,5 @@
                     tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
                     sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
-                    v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
+                    v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn) ) 
                     htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 )
                     str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 )
@@ -202,5 +200,5 @@
              WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1)
             
-            vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1))
+            vsum = ptr_sj( pvtr(:,:,:), btmsk(:,:,1))
             tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) )
             tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) )
@@ -224,5 +222,5 @@
                     r1_sjk(:,1,jn) = 0._wp
                     WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
-                    vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn))
+                    vsum = ptr_sj( pvtr(:,:,:), btmsk(:,:,jn))
                     tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
                     tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
@@ -408,5 +406,12 @@
             ENDIF
             IF( iom_use("zomsfeivglo") ) THEN
-               z3d(1,:,:) = ptr_sjk( v_eiv(:,:,:) )  ! zonal cumulative effective transport
+               DO jk=1,jpk
+                  DO jj=1,jpj
+                     DO ji=1,jpi
+                        zvn(ji,jj,jk) = v_eiv(ji,jj,jk) * fse3v(ji,jj,jk) * e1v(ji,jj)
+                     ENDDO
+                  ENDDO
+               ENDDO
+               z3d(1,:,:) = ptr_sjk( zvn(:,:,:) )  ! zonal cumulative effective transport
                DO jk = jpkm1,1,-1
                  z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)   ! effective j-Stream-Function (MSF)
@@ -419,5 +424,5 @@
                IF( ln_subbas ) THEN
                   DO jn = 2, nptr                                    ! by sub-basins
-                     z3d(1,:,:) =  ptr_sjk( v_eiv(:,:,:), btmsk(:,:,jn) ) 
+                     z3d(1,:,:) =  ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
                      DO jk = jpkm1,1,-1
                         z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)    ! effective j-Stream-Function (MSF)
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 9157)
@@ -549,10 +549,34 @@
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
 
-                                        !* sign of grad(H) at u- and v-points
-      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
+      !! AXY (16/08/17): remove the following per George and Andrew bug-hunt
+      !!                                   !* sign of grad(H) at u- and v-points
+      !! mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
+      !! DO jj = 1, jpjm1
+      !!    DO ji = 1, jpim1
+      !!       mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+      !!       mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+      !!    END DO
+      !! END DO
+
+      !! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt
+                                        !* sign of grad(H) at u- and  v-points; zero if grad(H) = 0
+      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
-            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
-            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+#if defined key_bbl_old_nonconserve
+             ! This key allows old (non conservative version) to be used for continuity of results
+             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+#else
+            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
+               mgrhu(ji,jj) = INT(  SIGN( 1.e0, &
+               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            ENDIF
+            !     
+            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
+               mgrhv(ji,jj) = INT(  SIGN( 1.e0, &
+               gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            ENDIF
+#endif
          END DO
       END DO
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 9157)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Add air-sea flux kill switch
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -372,4 +373,9 @@
          ENDDO
       ENDDO
+#  endif
+
+#  if defined key_axy_killco2flux
+      !! AXY (18/08/17): single kill switch on air-sea CO2 flux for budget checking
+      f_co2flux(:,:) = 0.
 #  endif
 
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 9157)
@@ -39,6 +39,5 @@
 # endif
       USE lbclnk,            ONLY: lbc_lnk
-      USE trc,               ONLY: trn
-      USE oce,               ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl
+      USE oce,               ONLY: CO2Flux_out_cpl, DMS_out_cpl
       USE par_oce,           ONLY: jpi, jpj
       USE sbc_oce,           ONLY: lk_oasis, qsr, wndm
@@ -48,5 +47,5 @@
                                    jdms, ocal_ccd, xpar, xze,              &
                                    zb_co2_flx, zb_dms_srf,                 &
-                                   zn_co2_flx, zn_dms_srf, zn_chl_srf
+                                   zn_co2_flx, zn_dms_srf
       USE trc,               ONLY: med_diag
 
@@ -66,11 +65,4 @@
       !!
       IF (jk.eq.1) THEN
-         !! JPALM -- 02-06-2017 --
-         !! add Chl surf coupling
-         !! no need to output, just pass to cpl var
-         IF (lk_oasis) THEN
-            zn_chl_srf(:,:) = (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6  !! surf Chl in Kg-chl/m3 as needed for cpl
-            chloro_out_cpl(:,:) = zn_chl_srf(:,:)        !! Coupling Chl
-         END IF
          IF( med_diag%MED_QSR%dgsave ) THEN
             CALL iom_put( "MED_QSR"  , qsr ) !
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 9157)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Amend bethic reservoir updating
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -38,8 +39,10 @@
 # endif
       USE lbclnk,            ONLY: lbc_lnk
+      USE oce,               ONLY: chloro_out_cpl 
       USE par_medusa,        ONLY: jp_medusa_2d, jp_medusa_3d,          &
-                                   jp_medusa_trd
+                                   jp_medusa_trd, jpchd, jpchn
       USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk
       USE phycst,            ONLY: rsmall
+      USE sbc_oce,           ONLY: lk_oasis
       USE sms_medusa,        ONLY: jinorgben, jorgben,                  &
                                    f3_co3, f3_h2co3, f3_hco3,           &
@@ -50,6 +53,7 @@
                                    zb_sed_n, zb_sed_si,                 &
                                    zn_sed_c, zn_sed_ca, zn_sed_fe,      &
-                                   zn_sed_n, zn_sed_si
-      USE trc,               ONLY: med_diag, nittrc000 
+                                   zn_sed_n, zn_sed_si, zn_chl_srf,     &
+                                   scl_chl, chl_out
+      USE trc,               ONLY: med_diag, nittrc000, trn 
       USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
  
@@ -62,159 +66,52 @@
       REAL(wp) :: fq0,fq1,fq2,fq3
 
+# if defined key_roam                     
+      !!----------------------------------------------------------------------
+      !! AXY (09/08/17): fix benthic submodel
       !!----------------------------------------------------------------------
       !! Process benthic in/out fluxes
       !! These can be handled outside of the 3D calculations since the
-      !! benthic pools (and fluxes) are 2D in nature; this code is
-      !! (shamelessly) borrowed from corresponding code in the LOBSTER
-      !! model
+      !! benthic pools (and fluxes) are 2D in nature; this code was
+      !! developed with help from George Nurser (NOC); it cannot be run
+      !! in a configuration with variable time-stepping with depth
       !!----------------------------------------------------------------------
       !!
-      !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt
+      !! time-step calculation
       if (jorgben.eq.1) then
-         za_sed_n(:,:)  = zn_sed_n(:,:)  +                                 & 
-                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -          &
-                            f_benout_n(:,:)  ) * (rdt / 86400.)
+         za_sed_n(:,:)  = zb_sed_n(:,:)  + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ))
+         za_sed_fe(:,:) = zb_sed_fe(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ))
+         za_sed_c(:,:)  = zb_sed_c(:,:)  + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ))
+      endif
+      if (jinorgben.eq.1) then
+         za_sed_si(:,:) = zb_sed_si(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ))
+         za_sed_ca(:,:) = zb_sed_ca(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ))
+      endif
+      !!
+      !! time-level calculation
+      if (jorgben.eq.1) then
+         zb_sed_n(:,:)  = zn_sed_n(:,:)  + (atfp * &
+                          ( za_sed_n(:,:)  - (2. * zn_sed_n(:,:))  + zb_sed_n(:,:)  ))
          zn_sed_n(:,:)  = za_sed_n(:,:)
-         !!
-         za_sed_fe(:,:) = zn_sed_fe(:,:) +                                 &
-                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -          &
-                            f_benout_fe(:,:) ) * (rdt / 86400.)
+         zb_sed_fe(:,:) = zn_sed_fe(:,:) + (atfp * &
+                          ( za_sed_fe(:,:) - (2. * zn_sed_fe(:,:)) + zb_sed_fe(:,:) ))
          zn_sed_fe(:,:) = za_sed_fe(:,:)
-         !!
-         za_sed_c(:,:)  = zn_sed_c(:,:)  +                                 &
-                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -          &
-                            f_benout_c(:,:)  ) * (rdt / 86400.)
+         zb_sed_c(:,:)  = zn_sed_c(:,:)  + (atfp * &
+                          ( za_sed_c(:,:)  - (2. * zn_sed_c(:,:))  + zb_sed_c(:,:)  ))
          zn_sed_c(:,:)  = za_sed_c(:,:)
       endif
       if (jinorgben.eq.1) then
-         za_sed_si(:,:) = zn_sed_si(:,:) +                                 & 
-                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        &
-                          (rdt / 86400.)
+         zb_sed_si(:,:) = zn_sed_si(:,:) + (atfp * &
+                          ( za_sed_si(:,:) - (2. * zn_sed_si(:,:)) + zb_sed_si(:,:) ))
          zn_sed_si(:,:) = za_sed_si(:,:)
-         !!
-         za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 &
-                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        &
-                          (rdt / 86400.)
+         zb_sed_ca(:,:) = zn_sed_ca(:,:) + (atfp * &
+                          ( za_sed_ca(:,:) - (2. * zn_sed_ca(:,:)) + zb_sed_ca(:,:) ))
          zn_sed_ca(:,:) = za_sed_ca(:,:)
       endif
-      !!
-      if (ibenthic.eq.2) then
-         !! The code below (in this if ... then ... endif loop) is
-         !! effectively commented out because it does not work as 
-         !! anticipated; it can be deleted at a later date
-         if (jorgben.eq.1) then
-            za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -       &
-                               f_benout_n(:,:)  ) * rdt
-            za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -       &
-                               f_benout_fe(:,:) ) * rdt
-            za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -       &
-                               f_benout_c(:,:)  ) * rdt
-         endif
-         if (jinorgben.eq.1) then
-            za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt
-            za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt
-         endif
-         !!
-         !! Leap-frog scheme - only in explicit case, otherwise the 
-         !! time stepping is already being done in trczdf
-         !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN
-         !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc )
-         !!    IF( neuler == 0 .AND. kt == nittrc000 )  zfact = rdttra(jk) * 
-         !!                                             FLOAT(ndttrc)
-         !!    if (jorgben.eq.1) then
-         !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  )
-         !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) )
-         !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  )
-         !!    endif
-         !!    if (jinorgben.eq.1) then
-         !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) )
-         !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) )
-         !!    endif
-         !! ENDIF
-         !! 
-         !! Time filter and swap of arrays
-         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme
-            IF( neuler == 0 .AND. kt == nittrc000 ) THEN
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = zn_sed_n(:,:)
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = zn_sed_fe(:,:)
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = zn_sed_c(:,:)
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = zn_sed_si(:,:)
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = zn_sed_ca(:,:)
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ELSE
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = (atfp  *                                 &
-                                    ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) +  &
-                                      (atfp1 * zn_sed_n(:,:) )
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = (atfp  *                                 &
-                                    ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) +  &
-                                      (atfp1 * zn_sed_fe(:,:))
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = (atfp  *                                 &
-                                    ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) +  &
-                                      (atfp1 * zn_sed_c(:,:) )
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = (atfp  *                                 &
-                                    ( zb_sed_si(:,:) + za_sed_si(:,:) )) +  &
-                                      (atfp1 * zn_sed_si(:,:))
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = (atfp  *                                 &
-                                    ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) +  &
-                                      (atfp1 * zn_sed_ca(:,:))
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ENDIF
-         ELSE                   !  case of smolar scheme or muscl
-            if (jorgben.eq.1) then
-               zb_sed_n(:,:)  = za_sed_n(:,:)
-               zn_sed_n(:,:)  = za_sed_n(:,:)
-               za_sed_n(:,:)  = 0.0
-               !!
-               zb_sed_fe(:,:) = za_sed_fe(:,:)
-               zn_sed_fe(:,:) = za_sed_fe(:,:)
-               za_sed_fe(:,:) = 0.0
-               !!
-               zb_sed_c(:,:)  = za_sed_c(:,:)
-               zn_sed_c(:,:)  = za_sed_c(:,:)
-               za_sed_c(:,:)  = 0.0
-            endif
-            if (jinorgben.eq.1) then
-               zb_sed_si(:,:) = za_sed_si(:,:)
-               zn_sed_si(:,:) = za_sed_si(:,:)
-               za_sed_si(:,:) = 0.0
-               !!
-               zb_sed_ca(:,:) = za_sed_ca(:,:)
-               zn_sed_ca(:,:) = za_sed_ca(:,:)
-               za_sed_ca(:,:) = 0.0
-            endif
-         ENDIF
-      endif
+# endif      
 
 #  if defined key_debug_medusa
@@ -253,7 +150,8 @@
                   fq1 = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj)
                   fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err)=',      &
-                     fq0,fq1,fq2
+                  fq3 = f_benout_n(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      &
+                     fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -266,7 +164,8 @@
                   fq1 = f_fbenin_si(ji,jj)
                   fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err)=',     &
-                     fq0,fq1,fq2
+                  fq3 = f_benout_si(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
+                     fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -278,22 +177,24 @@
                   fq0 = fflx_c(ji,jj)
                   fq1 = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) + f_fbenin_ca(ji,jj)
-                  fq2 = f_co2flux(ji,jj) * e3t_n(ji,jj,1)
+                  fq2 = f_co2flux(ji,jj) * fse3t(ji,jj,1)
                   fq3 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,4f15.10)')               &
-                    'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err)=',  &
-                    fq0,fq1,fq2,fq3
-               ENDIF
-            ENDDO
-         ENDDO   
-         !! alkalinity
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
-                  fq0 = fflx_a(ji,jj)
-                  fq1 = 2.0 * f_fbenin_ca(ji,jj)
-                  fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err)=',     &
-                     fq0,fq1,fq2
+                  fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   &
+                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', &
+                     fq0,fq1,fq2,fq3,fq4
+                ENDIF
+             ENDDO
+          ENDDO   
+          !! alkalinity
+          DO jj = 2,jpjm1
+             DO ji = 2,jpim1
+                if (tmask(ji,jj,1) == 1) then
+                   fq0 = fflx_a(ji,jj)
+                   fq1 = 2.0 * f_fbenin_ca(ji,jj)
+                   fq2 = fq0 + fq1
+                   fq3 = 2.0 * f_benout_ca(ji,jj)
+                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
+                      fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -333,4 +234,24 @@
             ENDDO
          ENDDO
+
+         !!!---------------------------------------------------------------
+         !! Calculates Chl diag for UM coupling 
+         !!!---------------------------------------------------------------
+         !! JPALM -- 02-06-2017 --
+         !! add Chl surf coupling
+         !! no need to output, just pass to cpl var
+         IF (lk_oasis) THEN
+            IF (chl_out.eq.1) THEN
+               !! export and scale surface chl
+               zn_chl_srf(:,:) = MAX( 0.0, (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6 )
+                                 !! surf Chl in Kg-chl/m3 as needed for cpl
+            ELSEIF (chl_out.eq.2) THEN
+               !! export and scale mld chl
+               zn_chl_srf(:,:) = MAX( 0.0, fchl_ml(:,:) * 1.0E-6 )
+                                 !! mld Chl in Kg-chl/m3 as needed for cpl
+            ENDIF
+            chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling Chl
+         END IF
+
          !!----------------------------------------------------------------
          !! Add in XML diagnostics stuff
@@ -360,4 +281,12 @@
             CALL iom_put( "OCAL_LVL"  , fccd )
          ENDIF
+         IF ( med_diag%CHL_MLD%dgsave ) THEN
+            CALL iom_put( "CHL_MLD"  , fchl_ml )
+         ENDIF
+         IF (lk_oasis) THEN
+            IF ( med_diag%CHL_CPL%dgsave ) THEN
+               CALL iom_put( "CHL_CPL"  , chloro_out_cpl )
+            ENDIF
+         ENDIF
          IF ( med_diag%PN_JLIM%dgsave ) THEN
             CALL iom_put( "PN_JLIM"  , fjln2d )
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 9157)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Add slow-sinking detrius variables
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -161,4 +162,6 @@
       fprn_ml(:,:) = 0.0        !! mixed layer PP diagnostics
       fprd_ml(:,:) = 0.0        !! mixed layer PP diagnostics
+      !! AXY (16/08/17)
+      fchl_ml(:,:) = 0.0	!! mixed layer chlorophyll diagnostics
       !! 
       fslownflux(:,:) = 0.0 
@@ -183,5 +186,10 @@
       f_omarg(:,:)    = 0.0
       f_omcal(:,:)    = 0.0
-
+      !!
+      !! AXY (08/08/17): zero slow detritus fluxes
+      fslowsink(:,:)  = 0.0
+# if defined key_roam
+      fslowsinkc(:,:) = 0.0
+# endif      
       !!
       !! allocate and initiate 2D diag
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 9157)
@@ -7,4 +7,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Slow detritus, ML-avg chl variables
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -54,4 +55,6 @@
    !! AXY (01/03/10): add in mixed layer PP diagnostics
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fprn_ml,fprd_ml
+   !! AXY (16/08/17): add in mixed layer chlorophyll diagnostic
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fchl_ml
    !!
    !! nutrient limiting factors
@@ -94,5 +97,11 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregenfastc
 # endif
-
+   !!
+   !! AXY (08/08/17): sinking of detritus moved here
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    fslowsink, fslowgain, fslowloss
+# if defined key_roam
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    fslowsinkc, fslowgainc, fslowlossc
+# endif
+   !!
    !! Particle flux
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdep1
@@ -287,5 +296,5 @@
                fjlim_pn(jpi,jpj),fjlim_pd(jpi,jpj),                   &
                fun_T(jpi,jpj),fun_Q10(jpi,jpj),                       &
-               fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),                     &
+               fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),fchl_ml(jpi,jpj),    &
                fnln(jpi,jpj),ffln2(jpi,jpj),                          &
                fnld(jpi,jpj),ffld(jpi,jpj),fsld(jpi,jpj),             &
@@ -316,4 +325,10 @@
                fregenfastc(jpi,jpj),                                  &
 # endif
+	       fslowsink(jpi,jpj),fslowgain(jpi,jpj),                 &
+               fslowloss(jpi,jpj),                                    &
+# if defined key_roam
+	       fslowsinkc(jpi,jpj),fslowgainc(jpi,jpj),               &
+               fslowlossc(jpi,jpj),                                   &
+# endif
                fdep1(jpi,jpj),                                        &
                ftempn(jpi,jpj),ftempsi(jpi,jpj),ftempfe(jpi,jpj),     &
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 9157)
@@ -6,4 +6,6 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Amend slow-detritus bug
+   !!   -   ! 2017-08 (A. Yool)            Reformatting for clarity
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -60,5 +62,8 @@
                                    fsil_cons, fsil_prod, fsdiss,             &
                                    ftempca, fthetad, fthetan,                &
+                                   fslowsink, fslowgain, fslowloss,          & ! AXY (22/08/17)
+                                   f_sbenin_n, f_sbenin_c,                   &
 # if defined key_roam
+                                   fslowsinkc, fslowgainc, fslowlossc,       & ! AXY (22/08/17)
                                    fcar_cons, fcar_prod, fcomm_resp,         &
                                    fddc, fflx_a, fflx_c, fflx_o2, zoxy,      &
@@ -181,25 +186,21 @@
       ENDDO
 
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            if (tmask(ji,jj,jk) == 1) then
-               !!
-               !!----------------------------------------------------------
-               !! detritus
-               btra(ji,jj,jpdet_lc) = b0 *                                      &
-                                          ! mort. losses 
-                                   (fdpn(ji,jj) + ((1.0 - xfdfrac1) *        &
-                                                   fdpd(ji,jj)) +            &
-                                    fdzmi(ji,jj) +                           &
-                                    ((1.0 - xfdfrac2) * fdzme(ji,jj)) +      &
-                                          ! assim. inefficiency
-                                    ((1.0 - xbetan) * (finmi(ji,jj) +        &
-                                                       finme(ji,jj))) -      &
-                                          ! grazing and remin.
-                                    fgmid(ji,jj) - fgmed(ji,jj) -            &
-                                    fdd(ji,jj) +                             &
-                                          ! seafloor fast->slow
-                                    ffast2slown(ji,jj))
-               !!
+      !!----------------------------------------------------------
+      !! detritus
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,jk) == 1) then
+               !!
+               btra(ji,jj,jpdet_lc) = b0 * (                           &
+                   fdpn(ji,jj)                                         & ! mort. losses 
+                 + ((1.0 - xfdfrac1) * fdpd(ji,jj))                    & ! mort. losses 
+                 + fdzmi(ji,jj)                                        & ! mort. losses
+                 + ((1.0 - xfdfrac2) * fdzme(ji,jj))                   & ! mort. losses
+                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))    & ! assim. inefficiency
+                 - fgmid(ji,jj) - fgmed(ji,jj)                         & ! grazing
+                 - fdd(ji,jj)                                          & ! remin.
+                 + fslowgain(ji,jj) - fslowloss(ji,jj)                 & ! slow-sinking
+                 - (f_sbenin_n(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor
+                 + ffast2slown(ji,jj) )                                  ! seafloor fast->slow 
             ENDIF
          ENDDO
@@ -305,24 +306,26 @@
                                           ffetop(ji,jj) + ffebot(ji,jj) -    &
                                           ffescav(ji,jj) )
+            ENDIF
+         ENDDO
+      ENDDO
+
 # if defined key_roam
-               !!
-               !!----------------------------------------------------------
-               !! AXY (26/11/08): implicit detrital carbon change
-               btra(ji,jj,jpdtc_lc) = b0 * (                                    &
-                                            ! mort. losses
-                                         (xthetapn * fdpn(ji,jj)) +          &
-                                         ((1.0 - xfdfrac1) *                 &
-                                          (xthetapd * fdpd(ji,jj))) +        &
-                                         (xthetazmi * fdzmi(ji,jj)) +        &
-                                         ((1.0 - xfdfrac2) *                 &
-                                          (xthetazme * fdzme(ji,jj))) +      &
-                                             ! assim. inefficiency
-                                         ((1.0 - xbetac) *                   &
-                                          (ficmi(ji,jj) + ficme(ji,jj))) -   &
-                                             ! grazing and remin.
-                                         fgmidc(ji,jj) - fgmedc(ji,jj) -     &
-                                         fddc(ji,jj) +                       &
-                                             ! seafloor fast->slow
-                                         ffast2slowc(ji,jj) )
+      !!----------------------------------------------------------
+      !! AXY (26/11/08): implicit detrital carbon change
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,jk) == 1) then 
+               !!
+               btra(ji,jj,jpdtc_lc) = b0 * (                           &
+                   (xthetapn * fdpn(ji,jj))                            & ! mort. losses 
+                 + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))       & ! mort. losses 
+                 + (xthetazmi * fdzmi(ji,jj))                          & ! mort. losses 
+                 + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))     & ! mort. losses 
+                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))    & ! assim. inefficiency
+                 - fgmidc(ji,jj) - fgmedc(ji,jj)                       & ! grazing
+                 - fddc(ji,jj)                                         & ! remin.
+                 + fslowgainc(ji,jj) - fslowlossc(ji,jj)               & ! slow-sinking
+                 - (f_sbenin_c(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor
+                 + ffast2slowc(ji,jj) )                                  ! seafloor fast->slow
             ENDIF
          ENDDO
@@ -575,8 +578,8 @@
                                                            f_o2flux(ji,jj))
                endif
+            ENDIF
+         ENDDO
+      ENDDO
 # endif
-            ENDIF
-         ENDDO
-      ENDDO
 
 # if defined key_debug_medusa
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 9157)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Revise slow-sinking of detritus
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -35,14 +36,21 @@
                                         f_sbenin_n, fdd,                   &
                                         idf, idfval,                       &   
-# if defined key_roam
+                                        fslowsink,                         &
+                                        fslowgain, fslowloss,              &
+# if defined key_roam
+                                        fslowsinkc,                        &
+                                        fslowgainc, fslowlossc,            &
                                         fddc,                              &
 # endif
                                         fun_T, fun_Q10, zdet, zdtc
       USE detritus_fast_sink_mod, ONLY: detritus_fast_sink
-      USE dom_oce,                ONLY: mbathy, tmask
+      USE dom_oce,                ONLY: mbathy, e3t_0, e3t_n, gphit, tmask
       USE in_out_manager,         ONLY: lwp, numout
       USE par_oce,                ONLY: jpim1, jpjm1
       USE sms_medusa,             ONLY: jmd, jorgben, jsfd, vsed,          &
                                         xrfn, xmd, xmdc, xthetad
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
 
       !! Level
@@ -123,4 +131,39 @@
          DO ji = 2,jpim1
             if (tmask(ji,jj,jk) == 1) then
+               !!----------------------------------------------------------------------
+               !! Detritus sinking (AXY, 08/08/18)
+	       !! Replaces slow-sinking done in trcsed_medusa.F90
+               !!
+               !! Uses the fslowsink variable to carry slow-sinking detritus from one
+               !! grid level to the next, variable fslowgain to "add" detritus sinking
+               !! from above and variable fslowloss to "subtract" detritus sinking out
+               !! to below; these variables appear in the differential equations of
+               !! detrital nitrogen and carbon below
+               !!----------------------------------------------------------------------
+               !!
+               fslowgain(ji,jj)  = fslowsink(ji,jj) / fse3t(ji,jj,jk)                  ! = mmol N / m3 / d
+               if (jk.lt.mbathy(ji,jj)) then
+                  fslowloss(ji,jj)  = (zdet(ji,jj) * vsed * 86400.) / fse3t(ji,jj,jk)  ! = mmol N / m3 / d
+               else
+                  fslowloss(ji,jj)  = 0.                                               ! = mmol N / m3 / d
+               endif
+               fslowsink(ji,jj) = fslowloss(ji,jj) * fse3t(ji,jj,jk)                   ! = mmol N / m2 / d
+               !!
+#  if defined key_roam
+               fslowgainc(ji,jj) = fslowsinkc(ji,jj) / fse3t(ji,jj,jk)                 ! = mmol C / m3 / d
+               if (jk.lt.mbathy(ji,jj)) then
+                  fslowlossc(ji,jj) = (zdtc(ji,jj) * vsed * 86400.) / fse3t(ji,jj,jk)  ! = mmol C / m3 / d
+               else
+                  fslowlossc(ji,jj) = 0.                                               ! = mmol C / m3 / d
+               endif
+               fslowsinkc(ji,jj) = fslowlossc(ji,jj) * fse3t(ji,jj,jk)                 ! = mmol C / m2 / d
+#  endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,jk) == 1) then
                !!---------------------------------------------------------
                !! Detritus addition to benthos
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 9157)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Mean mixed layer chlorophyll
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -42,5 +43,5 @@
                                    zchd, zchn, zdet, zdin, zdtc,         &
                                    zfer, zpds, zphd, zphn, zsil,         &
-                                   zzme, zzmi
+                                   zzme, zzmi, fchl_ml
       USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
       USE in_out_manager,    ONLY: lwp, numout
@@ -55,4 +56,5 @@
                                    xvpd, xvpn, xxi
       USE zdfmxl,            ONLY: hmld
+      USE lbclnk,            ONLY: lbc_lnk
 
    !!* Substitution
@@ -373,7 +375,11 @@
                fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * &
                                                   fse3t(ji,jj,jk) * fq0)
-            ENDIF
-         ENDDO
-      ENDDO
+	       !! AXY (16/08/17)
+	       fchl_ml(ji,jj) = fchl_ml(ji,jj) + ((zchn(ji,jj) + zchd(ji,jj)) * &
+	                                          (fse3t(ji,jj,jk) * fq0) / hmld(ji,jj))
+            ENDIF
+         ENDDO
+      ENDDO
+      CALL lbc_lnk(fchl_ml(:,:),'T',1. )
 
       DO jj = 2,jpjm1
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 9157)
@@ -177,5 +177,8 @@
    INTEGER  ::  jdms_model   !: choice of DMS model passed to atmosphere
 !!                              1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL
-!!
+!! FOR UKESM   
+   REAL(wp) ::  scl_chl      !: scaling factor for tuned Chl passed to the UM 
+   INTEGER  ::  chl_out      !: select Chl field exported and scaled for UM:
+                             !: 1- Surface Chl ; 2- MLD Chl
 !!
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   remdmp   !: depth dependent damping coefficient of passive tracers 
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 9157)
@@ -271,5 +271,5 @@
       zn_dms_srf(:,:)  = 0.0
       za_dms_srf(:,:)  = 0.0
-      zn_chl_srf(:,:)  = 2.0E-8 !! Chl srf
+      zn_chl_srf(:,:)  = 2.0E-8 !! Chl cpl - set first as surf
       !!
       IF(lwp) WRITE(numout,*) ' trc_ini_medusa: DMS and CO2 flux (UKESM) initialised to zero'
@@ -278,5 +278,5 @@
          CO2Flux_out_cpl(:,:) =  zn_co2_flx(:,:)   !! Coupling variable
          DMS_out_cpl(:,:)     =  zn_dms_srf(:,:)   !! Coupling variable
-         chloro_out_cpl(:,:)  =  zn_chl_srf(:,:)   !! Coupling variable
+         chloro_out_cpl(:,:)  =  zn_chl_srf(:,:) * scl_chl   !! Coupling variable
       END IF
       !!
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 9157)
@@ -24,4 +24,5 @@
    USE sms_medusa      ! sms trends
    USE iom             ! I/O manager
+   USE sbc_oce, ONLY: lk_oasis
    !!USE trc_nam_dia     ! JPALM 13-11-2015 -- if iom_use for diag
 
@@ -85,5 +86,6 @@
       &  jriver_n,jriver_si,jriver_c,jriver_alk,jriver_dep,   &
       &  xsdiss,                                              &
-      &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model
+      &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model,   &
+      &  scl_chl, chl_out
 #if defined key_roam
       NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        &
@@ -246,4 +248,6 @@
       jdms_input  = 0
       jdms_model  = 0
+      scl_chl     = 1.
+      chl_out     = 1
             
       !REWIND(numnatm)
@@ -400,4 +404,9 @@
 !!                      1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL, 5 = ANDM
 !!
+!! UKESM1 - exported Chl to UM
+!!       scl_chl     : scaling factor to tune the chl field sent to the UM
+!!       chl_out     : select the chl field to send at the UM:
+!!                     1- Surf Chl ; 2- MLD Chl 
+
       IF(lwp) THEN
 !!
@@ -909,30 +918,46 @@
 !!
 !! UKESM1 - new diagnostics  !! Jpalm; AXY (08/07/15)
-         WRITE(numout,*) '=== UKESM1-related parameters'
-         WRITE(numout,*)     &
-         &   ' include DMS diagnostic?,                                   jdms        = ', jdms
-         if (jdms_input .eq. 0) then
-            WRITE(numout,*)     &
-            &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous'
-         else
-            WRITE(numout,*)     &
-            &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average'
-         endif
-	 if (jdms_model .eq. 1) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)'
-	 elseif (jdms_model .eq. 2) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)'
-	 elseif (jdms_model .eq. 3) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)'
-	 elseif (jdms_model .eq. 4) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)'
-	 elseif (jdms_model .eq. 5) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001; modified)'
-         endif
+         WRITE(numout,*) '=== UKESM1-related parameters ==='
+         WRITE(numout,*) ' ---- --- ---'
+
+         IF (lk_oasis) THEN
+            WRITE(numout,*) '=== UKESM1 --  coupled DMS to the atmosphere'
+            WRITE(numout,*)     &
+            &   ' include DMS diagnostic?,                                   jdms        = ', jdms
+            if (jdms_input .eq. 0) then
+               WRITE(numout,*)     &
+               &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous'
+            else
+               WRITE(numout,*)     &
+               &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average'
+            endif
+   	    if (jdms_model .eq. 1) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)'
+	    elseif (jdms_model .eq. 2) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)'
+	    elseif (jdms_model .eq. 3) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)'
+	    elseif (jdms_model .eq. 4) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)'
+	    elseif (jdms_model .eq. 5) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001; modified)'
+            endif
+         
+            WRITE(numout,*) '=== UKESM1 --  coupled Chl to the atmosphere'
+            WRITE(numout,*)        &
+               &   ' Scaling factor to export tuned Chl to the atmosphere       scl_chl  = ', scl_chl
+            IF (chl_out .eq. 1) THEN
+               WRITE(numout,*)        &
+               &   ' Chl field to be scaled and sent to the atmosphere:         chl_out  = Surface Chl field '
+            ELSEIF (chl_out .eq. 2) THEN
+               WRITE(numout,*)        &
+               &   ' Chl field to be scaled and sent to the atmosphere:         chl_out  = MLD Chl field '
+            ENDIF
+         ENDIF ! IF lk_oasis=true
 !!
       ENDIF
@@ -2053,5 +2078,16 @@
           med_diag%OCN_DPCO2%dgsave = .FALSE.
       ENDIF
-      !!
+      !! UKESM additional
+      IF  (iom_use("CHL_MLD")) THEN 
+          med_diag%CHL_MLD%dgsave = .TRUE.
+      ELSE 
+          med_diag%CHL_MLD%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("CHL_CPL")) THEN 
+          med_diag%CHL_CPL%dgsave = .TRUE.
+      ELSE 
+          med_diag%CHL_CPL%dgsave = .FALSE.
+      ENDIF
+      !! 3D
       IF  (iom_use("TPP3")) THEN 
           med_diag%TPP3%dgsave = .TRUE.
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 9157)
@@ -8,4 +8,5 @@
    !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
+   !!              -   !  2017-08  (A. Yool) amend for slow detritus bug
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -88,9 +89,10 @@
 #  endif
       
-      CALL trc_sed_medusa( kt ) ! sedimentation model
-#  if defined key_debug_medusa
-      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa'
-      CALL flush(numout)
-#  endif
+!! AXY (08/08/2017): remove call to buggy subroutine (now handled by detritus.F90)
+!!       CALL trc_sed_medusa( kt ) ! sedimentation model
+!! #  if defined key_debug_medusa
+!!       IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa'
+!!       CALL flush(numout)
+!! #  endif
 # endif
 
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 9157)
@@ -134,4 +134,6 @@
                   OCN_KWCO2, OCN_K0, CO2STARAIR, OCN_DPCO2,                                          & ! end of regular 2D
                   TPP3, DETFLUX3, REMIN3N, PH3, OM_CAL3,                                             & ! end of regular 3D
+! JPALM (01/09/17): additional UKESM 2D diag
+                  CHL_MLD, CHL_CPL,                                                                  &
 ! AXY (11/11/16): additional CMIP6 2D diagnostics
                   epC100, epCALC100, epN100, epSI100,                                                &
Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9133)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9157)
@@ -276,13 +276,17 @@
       !!                     as proxy of org matter from the ocean
       !!                  -- needed for the coupling with atm
+      !!       07-12-2017 -- To make things cleaner, we want to store an  
+      !!                     unscaled Chl field in the restart and only 
+      !!                     scale it when reading it in.
+
       IF( iom_varid( numrtr, 'N_CHL_srf', ldstop = .FALSE. ) > 0 ) THEN
-         IF(lwp) WRITE(numout,*) 'Chl surf concentration - reading in ...'
+         IF(lwp) WRITE(numout,*) 'Chl cpl concentration - reading in ... - scale by ', scl_chl
          CALL iom_get( numrtr, jpdom_autoglo, 'N_CHL_srf',  zn_chl_srf(:,:)  )
       ELSE
-         IF(lwp) WRITE(numout,*) 'Chl surf concentration - setting to zero ...'
-         zn_chl_srf(:,:)  = (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6
+         IF(lwp) WRITE(numout,*) 'set Chl coupled concentration - scaled by ', scl_chl
+         zn_chl_srf(:,:)  = MAX( 0.0, (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6 )
       ENDIF
       IF (lk_oasis) THEN
-         chloro_out_cpl(:,:) = zn_chl_srf(:,:)        !! Coupling variable
+         chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling variable
       END IF
       !!
@@ -296,5 +300,7 @@
       call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
       call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
-      call trc_rst_dia_stat(zn_chl_srf(:,:), 'CHL surf')
+      IF (lk_oasis) THEN
+         call trc_rst_dia_stat(chloro_out_cpl(:,:), 'CHL  cpl')
+      END IF
       !!  
       !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
@@ -456,4 +462,7 @@
       CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx',  zb_co2_flx(:,:)  )
       CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx',  zn_co2_flx(:,:)  )
+      !! JPALM 07-12-2017 -- To make things cleaner, we want to store an  
+      !!                     unscaled Chl field in the restart and only 
+      !!                     scale it when reading it in.
       CALL iom_rstput( kt, nitrst, numrtw, 'N_CHL_srf',  zn_chl_srf(:,:)  )
       !!
@@ -467,5 +476,5 @@
       call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
       call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
-      call trc_rst_dia_stat(zn_chl_srf(:,:), 'CHL surf')
+      call trc_rst_dia_stat(zn_chl_srf(:,:), 'unscaled CHL cpl')
       !!
       IF(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.'
