Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 8521)
@@ -448,4 +448,5 @@
        <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= "ATM_XCO2"   long_name="Atmospheric xCO2"                          unit="ppm"         />
        <field id= "OCN_FCO2"   long_name="Surface ocean fCO2"                        unit="uatm"        />
@@ -784,4 +785,5 @@
       <field field_ref= "CO2STARAIR" name="CO2STARAIR" />
       <field field_ref= "OCN_DPCO2"  name="OCN_DPCO2"  />
+      <field field_ref= "CHL_MLD"    name="CHL_MLD"    />
     </field_group>
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 8521)
@@ -549,10 +549,28 @@
       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( 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
          END DO
       END DO
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8521)
@@ -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
@@ -277,4 +278,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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8521)
@@ -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
@@ -62,159 +63,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 +147,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 +161,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 +174,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
@@ -360,4 +258,7 @@
             CALL iom_put( "OCAL_LVL"  , fccd )
          ENDIF
+         IF ( med_diag%CHL_MLD%dgsave ) THEN
+            CALL iom_put( "CHL_MLD"  , fchl_ml )
+         ENDIF
          IF ( med_diag%PN_JLIM%dgsave ) THEN
             CALL iom_put( "PN_JLIM"  , fjln2d )
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 8521)
@@ -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,7 +162,15 @@
       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 
       fslowcflux(:,:) = 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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 8521)
@@ -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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 8521)
@@ -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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 8521)
@@ -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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 8521)
@@ -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
@@ -373,4 +374,7 @@
                fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * &
                                                   fse3t(ji,jj,jk) * fq0)
+	       !! 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
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 8521)
@@ -2053,5 +2053,11 @@
           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
+      !! 3D
       IF  (iom_use("TPP3")) THEN 
           med_diag%TPP3%dgsave = .TRUE.
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 8521)
@@ -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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 8447)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 8521)
@@ -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,                                                                           &
 ! AXY (11/11/16): additional CMIP6 2D diagnostics
                   epC100, epCALC100, epN100, epSI100,                                                &
