Index: /branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8422)
+++ /branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8423)
@@ -32,5 +32,5 @@
       !!----------------------------------------------------------------------
       USE bio_medusa_mod
-      USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt
+      USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt, e3t_n, tmask
       USE in_out_manager,    ONLY: lwp, numout
 # if defined key_iomput
@@ -61,4 +61,7 @@
       INTEGER :: ji, jj
       INTEGER :: jn
+
+      REAL(wp) :: fq0,fq1,fq2,fq3
+
       !!----------------------------------------------------------------------
       !! Process benthic in/out fluxes
@@ -235,5 +238,88 @@
          ENDIF
       endif
-      
+     
+#  if defined key_debug_medusa
+         !! AXY (12/07/17)
+         !!-----------------------------------------------------------------
+         !! Check conservation of MEDUSA's sinks-minus-sources using fflx_X 
+         !! diagnostics (i.e. biogeochemical processes only)
+         !!   - fflx_X diagnostics *should* include all transfers between
+         !!     modelled components
+         !!   - they should also include gains / losses due to air-sea 
+         !!     fluxes of C and O2, aeolian and seafloor inputs of Fe, and 
+         !!     inputs from seafloor "benthic buckets" (N, Si, Fe, C and 
+         !!     alkalinity)
+         !!   - however, they do not include the transfer of material to 
+         !!     "benthic buckets" by sedimenting slow- and fast-sinking
+         !!     detritus since these are separate 2D reservoirs
+         !!   - consequently, for a given water column, the integrated
+         !!     fluxes should sum to the "loss" of material to the "benthic
+         !!     buckets"
+         !!   - if they do not, this suggests that MEDUSA contains errors
+         !!     in its accounting (e.g. processes omitted from calculated
+         !!     fluxes)
+         !!   - here, the local integrated fluxes and benthic inputs (plus
+         !!     air-sea fluxes in the case of C) are reported together with
+         !!     the resulting error
+         !!   - only N, Si, C and alkalinity inventories considered; Fe and
+         !!     O2 overlooked because of wholesale loss (and addition, in 
+         !!     the case of O2) of these tracers within the water column
+         !!-----------------------------------------------------------------
+         !!
+         !! nitrogen
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  fq0 = fflx_n(ji,jj)
+                  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
+               ENDIF
+            ENDDO
+         ENDDO
+         !! silicon
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  fq0 = fflx_si(ji,jj)
+                  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
+               ENDIF
+            ENDDO
+         ENDDO
+         !! carbon
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  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)
+                  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
+               ENDIF
+            ENDDO
+         ENDDO
+#  endif
+ 
       IF( ln_diatrc ) THEN
          !!-----------------------------------------------------------------
Index: /branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 8422)
+++ /branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 8423)
@@ -73,5 +73,9 @@
                                    jpalk, jpchd, jpchn, jpdet, jpdic,        &
                                    jpdin, jpdtc, jpfer, jpoxy, jppds,        &
-                                   jpphd, jpphn, jpsil, jpzme, jpzmi
+                                   jpphd, jpphn, jpsil, jpzme, jpzmi,        &
+                                   jpalk_lc, jpchd_lc, jpchn_lc, jpdet_lc,   & 
+                                   jpdic_lc, jpdin_lc, jpdtc_lc, jpfer_lc,   & 
+                                   jpoxy_lc, jppds_lc, jpphd_lc, jpphn_lc,   &
+                                   jpsil_lc, jpzme_lc, jpzmi_lc
       USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
       USE par_trc,           ONLY: jptra
@@ -130,10 +134,10 @@
                !!----------------------------------------------------------
                !! chlorophyll
-               btra(ji,jj,jpchn) = b0 * ( ( (frn(ji,jj) * fprn(ji,jj) *      &
+               btra(ji,jj,jpchn_lc) = b0 * ( ( (frn(ji,jj) * fprn(ji,jj) *      &
                                              zphn(ji,jj) ) -                 &
                                            fgmipn(ji,jj) - fgmepn(ji,jj) -   &
                                            fdpn(ji,jj) - fdpn2(ji,jj) ) *    &
                                           (fthetan(ji,jj) / xxi) )
-               btra(ji,jj,jpchd) = b0 * ( ( (frd(ji,jj) * fprd(ji,jj) *      &
+               btra(ji,jj,jpchd_lc) = b0 * ( ( (frd(ji,jj) * fprd(ji,jj) *      &
                                              zphd(ji,jj) ) -                 &
                                            fgmepd(ji,jj) - fdpd(ji,jj) -     &
@@ -150,11 +154,11 @@
                !!----------------------------------------------------------
                !! phytoplankton
-               btra(ji,jj,jpphn) = b0 * ( (fprn(ji,jj) * zphn(ji,jj)) -      &
+               btra(ji,jj,jpphn_lc) = b0 * ( (fprn(ji,jj) * zphn(ji,jj)) -      &
                                           fgmipn(ji,jj) - fgmepn(ji,jj) -    &
                                           fdpn(ji,jj) - fdpn2(ji,jj) )
-               btra(ji,jj,jpphd) = b0 * ( (fprd(ji,jj) * zphd(ji,jj)) -      &
+               btra(ji,jj,jpphd_lc) = b0 * ( (fprd(ji,jj) * zphd(ji,jj)) -      &
                                           fgmepd(ji,jj) - fdpd(ji,jj) -      &
                                           fdpd2(ji,jj) )
-               btra(ji,jj,jppds) = b0 * ( (fprds(ji,jj) * zpds(ji,jj)) -     &
+               btra(ji,jj,jppds_lc) = b0 * ( (fprds(ji,jj) * zpds(ji,jj)) -     &
                                           fgmepds(ji,jj) - fdpds(ji,jj) -    &
                                           fsdiss(ji,jj) - fdpds2(ji,jj) )
@@ -169,7 +173,7 @@
                !!----------------------------------------------------------
                !! zooplankton
-               btra(ji,jj,jpzmi) = b0 * (fmigrow(ji,jj) - fgmezmi(ji,jj) -   &
+               btra(ji,jj,jpzmi_lc) = b0 * (fmigrow(ji,jj) - fgmezmi(ji,jj) -   &
                                          fdzmi(ji,jj) - fdzmi2(ji,jj))
-               btra(ji,jj,jpzme) = b0 * (fmegrow(ji,jj) - fdzme(ji,jj) -     &
+               btra(ji,jj,jpzme_lc) = b0 * (fmegrow(ji,jj) - fdzme(ji,jj) -     &
                                          fdzme2(ji,jj))
             ENDIF
@@ -183,5 +187,5 @@
                !!----------------------------------------------------------
                !! detritus
-               btra(ji,jj,jpdet) = b0 *                                      &
+               btra(ji,jj,jpdet_lc) = b0 *                                      &
                                           ! mort. losses 
                                    (fdpn(ji,jj) + ((1.0 - xfdfrac1) *        &
@@ -236,5 +240,5 @@
                endif
                !!
-               btra(ji,jj,jpdin) = b0 * ( fn_prod + fn_cons )
+               btra(ji,jj,jpdin_lc) = b0 * ( fn_prod + fn_cons )
                !!
                !! consumption of dissolved nitrogen
@@ -281,5 +285,5 @@
                endif
                !!
-               btra(ji,jj,jpsil) = b0 * ( &
+               btra(ji,jj,jpsil_lc) = b0 * ( &
                  fs_prod + fs_cons )
                !! consumption of dissolved silicon
@@ -298,5 +302,5 @@
                !!----------------------------------------------------------
                !! dissolved "iron" nutrient
-               btra(ji,jj,jpfer) = b0 * ( (xrfn * btra(ji,jj,jpdin)) +       &
+               btra(ji,jj,jpfer_lc) = b0 * ( (xrfn * btra(ji,jj,jpdin_lc)) +       &
                                           ffetop(ji,jj) + ffebot(ji,jj) -    &
                                           ffescav(ji,jj) )
@@ -305,5 +309,5 @@
                !!----------------------------------------------------------
                !! AXY (26/11/08): implicit detrital carbon change
-               btra(ji,jj,jpdtc) = b0 * (                                    &
+               btra(ji,jj,jpdtc_lc) = b0 * (                                    &
                                             ! mort. losses
                                          (xthetapn * fdpn(ji,jj)) +          &
@@ -381,5 +385,5 @@
                endif
                !!
-               btra(ji,jj,jpdic) = b0 * ( fc_prod + fc_cons )
+               btra(ji,jj,jpdic_lc) = b0 * ( fc_prod + fc_cons )
                !! consumption of dissolved carbon
                fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fse3t(ji,jj,jk) *     &
@@ -418,5 +422,5 @@
                endif
                !!
-               btra(ji,jj,jpalk) = b0 * ( fa_prod + fa_cons )
+               btra(ji,jj,jpalk_lc) = b0 * ( fa_prod + fa_cons )
             ENDIF
          ENDDO
@@ -538,5 +542,5 @@
                !! deficient O2; production fluxes only
                if (zoxy(ji,jj).lt.xo2min) then
-                  btra(ji,jj,jpoxy) = b0 * fo2_prod(ji,jj)
+                  btra(ji,jj,jpoxy_lc) = b0 * fo2_prod(ji,jj)
                   foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) *  &
                                                           fo2_prod(ji,jj) )
@@ -545,5 +549,5 @@
                else
                   !! sufficient O2; production + consumption fluxes
-                  btra(ji,jj,jpoxy) = b0 * ( fo2_prod(ji,jj) +               &
+                  btra(ji,jj,jpoxy_lc) = b0 * ( fo2_prod(ji,jj) +               &
                                              fo2_cons(ji,jj) )
                   foxy_prod(ji,jj) = foxy_prod(ji,jj) +                      &
@@ -564,9 +568,9 @@
                   !!
                   !! CO2 flux
-                  btra(ji,jj,jpdic) = btra(ji,jj,jpdic) + (b0 *              &
+                  btra(ji,jj,jpdic_lc) = btra(ji,jj,jpdic_lc) + (b0 *              &
                                                            f_co2flux(ji,jj))
                   !!
                   !! O2 flux (mol/m3/s -> mmol/m3/d)
-                  btra(ji,jj,jpoxy) = btra(ji,jj,jpoxy) + (b0 *              &
+                  btra(ji,jj,jpoxy_lc) = btra(ji,jj,jpoxy_lc) + (b0 *              &
                                                            f_o2flux(ji,jj))
                endif
@@ -586,35 +590,35 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchn)(',jk,')  = ', &
-                                            btra(ji,jj,jpchn)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchd)(',jk,')  = ', &
-                                            btra(ji,jj,jpchd)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphn)(',jk,')  = ', &
-                                            btra(ji,jj,jpphn)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphd)(',jk,')  = ', &
-                                            btra(ji,jj,jpphd)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jppds)(',jk,')  = ', &
-                                            btra(ji,jj,jppds)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzmi)(',jk,')  = ', &
-                                            btra(ji,jj,jpzmi)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzme)(',jk,')  = ', &
-                                            btra(ji,jj,jpzme)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdet)(',jk,')  = ', &
-                                            btra(ji,jj,jpdet)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdin)(',jk,')  = ', &
-                                            btra(ji,jj,jpdin)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpsil)(',jk,')  = ', &
-                                            btra(ji,jj,jpsil)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpfer)(',jk,')  = ', &
-                                            btra(ji,jj,jpfer)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchn_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpchn_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchd_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpchd_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphn_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpphn_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphd_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpphd_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jppds_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jppds_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzmi_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpzmi_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzme_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpzme_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdet_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpdet_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdin_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpdin_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpsil_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpsil_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpfer_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpfer_lc)
 #  if defined key_roam
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdtc)(',jk,')  = ', &
-                                            btra(ji,jj,jpdtc)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdic)(',jk,')  = ', &
-                                            btra(ji,jj,jpdic)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpalk)(',jk,')  = ', &
-                                            btra(ji,jj,jpalk)
-                  IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy)(',jk,')  = ', &
-                                            btra(ji,jj,jpoxy)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdtc_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpdtc_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdic_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpdic_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpalk_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpalk_lc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpoxy_lc)
 #  endif
                endif
@@ -633,30 +637,30 @@
                !! === nitrogen ===
                fflx_n(ji,jj)  = fflx_n(ji,jj) + fse3t(ji,jj,jk) *            &
-                                ( btra(ji,jj,jpphn) + btra(ji,jj,jpphd) +    &
-                                  btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) +    &
-                                  btra(ji,jj,jpdet) + btra(ji,jj,jpdin) )
+                                ( btra(ji,jj,jpphn_lc) + btra(ji,jj,jpphd_lc) +    &
+                                  btra(ji,jj,jpzmi_lc) + btra(ji,jj,jpzme_lc) +    &
+                                  btra(ji,jj,jpdet_lc) + btra(ji,jj,jpdin_lc) )
                !! === silicon ===
                fflx_si(ji,jj) = fflx_si(ji,jj) + fse3t(ji,jj,jk) *           &
-                                ( btra(ji,jj,jppds) + btra(ji,jj,jpsil) )
+                                ( btra(ji,jj,jppds_lc) + btra(ji,jj,jpsil_lc) )
                !! === iron ===
                fflx_fe(ji,jj) = fflx_fe(ji,jj) + fse3t(ji,jj,jk) *           &
                                 ( (xrfn *                                    &
-                                   (btra(ji,jj,jpphn) + btra(ji,jj,jpphd) +  &
-                                    btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) +  &
-                                    btra(ji,jj,jpdet))) + btra(ji,jj,jpfer) )
+                                   (btra(ji,jj,jpphn_lc) + btra(ji,jj,jpphd_lc) +  &
+                                    btra(ji,jj,jpzmi_lc) + btra(ji,jj,jpzme_lc) +  &
+                                    btra(ji,jj,jpdet_lc))) + btra(ji,jj,jpfer_lc) )
 # if defined key_roam
                !! === carbon ===
                fflx_c(ji,jj)  = fflx_c(ji,jj) + fse3t(ji,jj,jk) *            &
-                                ( (xthetapn * btra(ji,jj,jpphn)) +           &
-                                  (xthetapd * btra(ji,jj,jpphd)) +           &
-                                  (xthetazmi * btra(ji,jj,jpzmi)) +          &
-                                  (xthetazme * btra(ji,jj,jpzme)) +          &
-                                  btra(ji,jj,jpdtc) + btra(ji,jj,jpdic) )
+                                ( (xthetapn * btra(ji,jj,jpphn_lc)) +           &
+                                  (xthetapd * btra(ji,jj,jpphd_lc)) +           &
+                                  (xthetazmi * btra(ji,jj,jpzmi_lc)) +          &
+                                  (xthetazme * btra(ji,jj,jpzme_lc)) +          &
+                                  btra(ji,jj,jpdtc_lc) + btra(ji,jj,jpdic_lc) )
                !! === alkalinity ===
                fflx_a(ji,jj)  = fflx_a(ji,jj) + fse3t(ji,jj,jk) *            &
-                                btra(ji,jj,jpalk)
+                                btra(ji,jj,jpalk_lc)
                !! === oxygen ===
                fflx_o2(ji,jj) = fflx_o2(ji,jj) + fse3t(ji,jj,jk) *           &
-                                btra(ji,jj,jpoxy)
+                                btra(ji,jj,jpoxy_lc)
 # endif
             ENDIF
@@ -694,34 +698,34 @@
                if (ibio_switch.eq.1) then
                   tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) +                &
-                                        (btra(ji,jj,jpchn) / 86400.)
+                                        (btra(ji,jj,jpchn_lc) / 86400.)
                   tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) +                &
-                                        (btra(ji,jj,jpchd) / 86400.)
+                                        (btra(ji,jj,jpchd_lc) / 86400.)
                   tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) +                &
-                                        (btra(ji,jj,jpphn) / 86400.)
+                                        (btra(ji,jj,jpphn_lc) / 86400.)
                   tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) +                &
-                                        (btra(ji,jj,jpphd) / 86400.)
+                                        (btra(ji,jj,jpphd_lc) / 86400.)
                   tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) +                &
-                                        (btra(ji,jj,jppds) / 86400.)
+                                        (btra(ji,jj,jppds_lc) / 86400.)
                   tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) +                &
-                                        (btra(ji,jj,jpzmi) / 86400.)
+                                        (btra(ji,jj,jpzmi_lc) / 86400.)
                   tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) +                &
-                                        (btra(ji,jj,jpzme) / 86400.)
+                                        (btra(ji,jj,jpzme_lc) / 86400.)
                   tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) +                &
-                                        (btra(ji,jj,jpdet) / 86400.)
+                                        (btra(ji,jj,jpdet_lc) / 86400.)
                   tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) +                &
-                                        (btra(ji,jj,jpdin) / 86400.)
+                                        (btra(ji,jj,jpdin_lc) / 86400.)
                   tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +                &
-                                        (btra(ji,jj,jpsil) / 86400.)
+                                        (btra(ji,jj,jpsil_lc) / 86400.)
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +                &
-                                        (btra(ji,jj,jpfer) / 86400.)
+                                        (btra(ji,jj,jpfer_lc) / 86400.)
 # if defined key_roam
                   tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) +                &
-                                        (btra(ji,jj,jpdtc) / 86400.)
+                                        (btra(ji,jj,jpdtc_lc) / 86400.)
                   tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +                &
-                                        (btra(ji,jj,jpdic) / 86400.)
+                                        (btra(ji,jj,jpdic_lc) / 86400.)
                   tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) +                &
-                                        (btra(ji,jj,jpalk) / 86400.)
+                                        (btra(ji,jj,jpalk_lc) / 86400.)
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) +                &
-                                        (btra(ji,jj,jpoxy) / 86400.)
+                                        (btra(ji,jj,jpoxy_lc) / 86400.)
 # endif
                endif
@@ -737,37 +741,37 @@
                IF( med_diag%FBDDTALK%dgsave )  THEN
                   fbddtalk(ji,jj)  =  fbddtalk(ji,jj)  +                     &
-                                      (btra(ji,jj,jpalk) * fse3t(ji,jj,jk))
+                                      (btra(ji,jj,jpalk_lc) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDIC%dgsave )  THEN
                   fbddtdic(ji,jj)  =  fbddtdic(ji,jj)  +                     &
-                                      (btra(ji,jj,jpdic) * fse3t(ji,jj,jk))
+                                      (btra(ji,jj,jpdic_lc) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDIFE%dgsave ) THEN
                   fbddtdife(ji,jj) =  fbddtdife(ji,jj) +                     &
-                                      (btra(ji,jj,jpfer) * fse3t(ji,jj,jk))
+                                      (btra(ji,jj,jpfer_lc) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDIN%dgsave )  THEN
                   fbddtdin(ji,jj)  =  fbddtdin(ji,jj)  +                     &
-                                      (btra(ji,jj,jpdin) * fse3t(ji,jj,jk))
+                                      (btra(ji,jj,jpdin_lc) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDISI%dgsave ) THEN
                   fbddtdisi(ji,jj) =  fbddtdisi(ji,jj) +                     &
-                                      (btra(ji,jj,jpsil) * fse3t(ji,jj,jk))
+                                      (btra(ji,jj,jpsil_lc) * fse3t(ji,jj,jk))
                ENDIF
 	       !!
                IF( med_diag%BDDTALK3%dgsave )  THEN
-                  bddtalk3(ji,jj,jk)  =  btra(ji,jj,jpalk)
+                  bddtalk3(ji,jj,jk)  =  btra(ji,jj,jpalk_lc)
                ENDIF
                IF( med_diag%BDDTDIC3%dgsave )  THEN
-                  bddtdic3(ji,jj,jk)  =  btra(ji,jj,jpdic)
+                  bddtdic3(ji,jj,jk)  =  btra(ji,jj,jpdic_lc)
                ENDIF
                IF( med_diag%BDDTDIFE3%dgsave ) THEN
-                  bddtdife3(ji,jj,jk) =  btra(ji,jj,jpfer)
+                  bddtdife3(ji,jj,jk) =  btra(ji,jj,jpfer_lc)
                ENDIF
                IF( med_diag%BDDTDIN3%dgsave )  THEN
-                  bddtdin3(ji,jj,jk)  =  btra(ji,jj,jpdin)
+                  bddtdin3(ji,jj,jk)  =  btra(ji,jj,jpdin_lc)
                ENDIF
                IF( med_diag%BDDTDISI3%dgsave ) THEN
-                  bddtdisi3(ji,jj,jk) =  btra(ji,jj,jpsil)
+                  bddtdisi3(ji,jj,jk) =  btra(ji,jj,jpsil_lc)
                ENDIF
             ENDIF
@@ -790,5 +794,5 @@
             if (tmask(ji,jj,jk) == 1) then
                !!
-               DO jn = 1,jptra
+               DO jn = 1,jp_medusa
                   fq0 = btra(ji,jj,jn)
                   !! AXY (30/01/14): "isnan" problem on HECTOR
@@ -801,5 +805,5 @@
                   endif
                ENDDO
-               DO jn = 1,jptra
+               DO jn = jp_msa0,jp_msa1
                   fq0 = tra(ji,jj,jk,jn)
                   !! AXY (30/01/14): "isnan" problem on HECTOR
@@ -818,46 +822,4 @@
 # endif
 
-      !!----------------------------------------------------------
-      !! Check model conservation
-      !! these terms merely sum up the tendency terms of the relevant
-      !! state variables, which should sum to zero; the iron cycle is
-      !! complicated by fluxes that add (aeolian deposition and 
-      !! seafloor remineralisation) and remove (scavenging) dissolved 
-      !! iron from the model (i.e. the sum of iron fluxes is unlikely 
-      !! to be zero)
-      !!----------------------------------------------------------
-!!      DO jj = 2,jpjm1
-!!         DO ji = 2,jpim1
-!!            if (tmask(ji,jj,jk) == 1) then
-               !!
-               !! fnit0 = btra(ji,jj,jpphn) + btra(ji,jj,jpphd) +             &
-               !!         btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) +             &
-               !!         btra(ji,jj,jpdet) + btra(ji,jj,jpdin)  ! +          &
-               !!         ftempn(ji,jj)
-               !! fsil0 = btra(ji,jj,jppds) + btra(ji,jj,jpsil) ! +           &
-               !!         ftempsi(ji,jj)
-               !! ffer0 = (xrfn * fnit0) + btra(ji,jj,jpfer)
-# if defined key_roam
-               !! fcar0 = 0.
-               !! falk0 = 0.
-               !! foxy0 = 0.
-# endif
-               !!
-               !! if (kt/240*240.eq.kt) then
-               !!    if (ji.eq.2.and.jj.eq.2.and.jk.eq.1) then
-               !!       IF (lwp) write (*,*)                                  &
-               !!          '*******!MEDUSA Conservation!*******',kt
-# if defined key_roam
-               !!       IF (lwp) write (*,*) fnit0,fsil0,ffer0,fcar0,falk0,   &
-               !!                            foxy0
-# else
-               !!       IF (lwp) write (*,*) fnit0,fsil0,ffer0
-# endif
-               !!    endif
-               !! endif     
-
-!!            ENDIF
-!!         ENDDO
-!!      ENDDO
 
    END SUBROUTINE bio_medusa_update
