Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8011)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8012)
@@ -50,4 +50,5 @@
                                    fccd, fcomm_resp,                     &
                                    fd_cal3, fd_car3, fd_nit3, fd_sil3,   &
+                                   fdep1,                                &
                                    fdd, fdd2d, fddc,                     &
                                    fdpd, fdpd2, fdpd22d, fdpd2d,         &
@@ -74,7 +75,7 @@
                                    fprds, fprds2d,                       & 
                                    fprn, fprn_ml, fprn2d,                &
-                                   fregen2d,                     &
+                                   fregen, fregen2d,                     &
                                    fregenfast, fregenfastsi,             &
-                                   fregensi2d,                 &
+                                   fregensi, fregensi2d,                 &
                                    freminc, freminc2d,                   &
                                    freminca, freminca2d,                 & 
@@ -95,8 +96,11 @@
                                    iben_c2d, iben_ca2d, iben_fe2d,       &
                                    iben_n2d, iben_si2d,                  &
+                                   intdissic, intdissin,                 &
+                                   intdissisi, inttalk,                  &
                                    iters, lyso_ca2d,                     &
                                    mdetc2d, megrazd3, megrazp3,          &
                                    megrazz3,                             & 
                                    migrazd3, migrazp3,                   &
+                                   o2min, o2sat3,                        &
                                    oben_c2d, oben_ca2d, oben_fe2d,       &
                                    oben_n2d, oben_si2d,                  &
@@ -114,10 +118,14 @@
                                    zimesd2d, zimesdc2d, zimesn2d,        &
                                    ziresp2d,                             &
-                                   zpds, zphd, zphn
-      USE dom_oce,           ONLY: e3t_0, e3t_n, mbathy, tmask
+                                   zo2min,                               &
+                                   zalk, zdet, zdic, zdin, zdtc,         &
+                                   zfer, zoxy, zpds, zphd, zphn,         &
+                                   zsal, zsil, ztmp, zzme, zzmi
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,       &
+                                   mbathy, tmask
       USE in_out_manager,    ONLY: lwp, numout
       USE iom,               ONLY: lk_iomput
       USE par_kind,          ONLY: wp
-      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
+      USE par_oce,           ONLY: jpim1, jpjm1
       USE phycst,            ONLY: rsmall
       USE sbc_oce,           ONLY: qsr, wndm
@@ -126,7 +134,8 @@
                                    i0100, i0150, i0200, i0500, i1000,    &
                                    jdms, ocal_ccd,                       &
-                                   xbetac, xbetan, xpar, xphi,           &
-                                   xthetapd, xthetapn, xthetazmi, xze
+                                   xbetac, xbetan, xpar, xphi, xrfn,     &
+                                   xthetapd, xthetapn, xthetazme, xthetazmi, xze
       USE trc,               ONLY: ln_diatrc, med_diag, trc2d, trc3d 
+      USE trcoxy_medusa,     ONLY: oxy_sato
 
    !!* Substitution
@@ -140,6 +149,4 @@
       !! Loop avariables
       INTEGER :: ji, jj, jn
-
-      REAL(wp), DIMENSION(jpi,jpj) :: fregen,fregensi
 
 # if defined key_trc_diabio
@@ -158,4 +165,124 @@
       !!
       IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio'
+# endif
+
+      !! The section below, down to calculation of zo2min, was moved 
+      !! from before the call to AIR_SEA in trcbio_medusa.F90 - marc 9/5/17 
+      IF( lk_iomput ) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  !! sum tracers for inventory checks
+                  IF ( med_diag%INVTN%dgsave )   THEN
+                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        &
+                        (fse3t(ji,jj,jk) * (zphn(ji,jj) + zphd(ji,jj) +      &
+                                            zzmi(ji,jj) + zzme(ji,jj) +      &
+                                            zdet(ji,jj) + zdin(ji,jj)))
+                  ENDIF
+                  IF ( med_diag%INVTSI%dgsave )  THEN
+                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       & 
+                       (fse3t(ji,jj,jk) * (zpds(ji,jj) + zsil(ji,jj)))
+                  ENDIF
+                  IF ( med_diag%INVTFE%dgsave )  THEN
+                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       & 
+                        (fse3t(ji,jj,jk) * (xrfn *                           &
+                                            (zphn(ji,jj) + zphd(ji,jj) +     &
+                                             zzmi(ji,jj) + zzme(ji,jj) +     &
+                                             zdet(ji,jj)) +                  &
+                                            zfer(ji,jj)))
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+
+# if defined key_roam
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  IF ( med_diag%INVTC%dgsave )  THEN
+                     ftot_c(ji,jj)  = ftot_c(ji,jj) +                        & 
+                        (fse3t(ji,jj,jk) * ((xthetapn * zphn(ji,jj)) +       &
+                                            (xthetapd * zphd(ji,jj)) +       &
+                                            (xthetazmi * zzmi(ji,jj)) +      &
+                                            (xthetazme * zzme(ji,jj)) +      &
+                                            zdtc(ji,jj) + zdic(ji,jj)))
+                  ENDIF
+                  IF ( med_diag%INVTALK%dgsave ) THEN
+                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     &
+                                                       zalk(ji,jj))
+                  ENDIF
+                  IF ( med_diag%INVTO2%dgsave )  THEN
+                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    &
+                                                        zoxy(ji,jj))
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  IF ( med_diag%INVTC%dgsave )  THEN
+                     !!
+                     !! AXY (10/11/16): CMIP6 diagnostics
+                     IF ( med_diag%INTDISSIC%dgsave ) THEN
+                        intdissic(ji,jj) = intdissic(ji,jj) +                &
+                                           (fse3t(ji,jj,jk) * zdic(ji,jj))
+                     ENDIF
+                     IF ( med_diag%INTDISSIN%dgsave ) THEN
+                        intdissin(ji,jj) = intdissin(ji,jj) +                &
+                                           (fse3t(ji,jj,jk) * zdin(ji,jj))
+                     ENDIF
+                     IF ( med_diag%INTDISSISI%dgsave ) THEN
+                        intdissisi(ji,jj) = intdissisi(ji,jj) +              &
+                                            (fse3t(ji,jj,jk) * zsil(ji,jj))
+                     ENDIF
+                     IF ( med_diag%INTTALK%dgsave ) THEN
+                        inttalk(ji,jj) = inttalk(ji,jj) +                    &
+                                         (fse3t(ji,jj,jk) * zalk(ji,jj))
+                     ENDIF
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  IF ( med_diag%O2min%dgsave ) THEN
+                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then
+                        o2min(ji,jj)  = zoxy(ji,jj)
+                        IF ( med_diag%ZO2min%dgsave ) THEN
+                           !! layer midpoint
+                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               &
+                                            fdep1(ji,jj)) / 2.0
+                        ENDIF
+                     endif
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+# endif
+      ENDIF
+
+# if defined key_roam
+      !! This section is moved from just below CALL to AIR_SEA in
+      !! trcbio_medusa.F90 - marc 9/5/17
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
+
+               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
+               IF ( med_diag%O2SAT3%dgsave ) THEN
+! Remove f_o2sat3 - marc 9/5/17
+!                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 )
+!                  o2sat3(ji, jj, jk) = f_o2sat3
+                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj),                   &
+                                 o2sat3(ji,jj,jk) )
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
 # endif
 
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 8011)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 8012)
@@ -85,4 +85,5 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fslown, fslowc
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fslownflux, fslowcflux
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregen,fregensi
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregenfast,fregenfastsi
 # if defined key_roam
@@ -309,6 +310,7 @@
                fdpn2(jpi,jpj),fdpd2(jpi,jpj),fdpds2(jpi,jpj),         &
                fdzmi2(jpi,jpj),fdzme2(jpi,jpj),                       &
-               fslown(jpi,jpj), fslowc(jpi,jpj),                      &
+               fslown(jpi,jpj),fslowc(jpi,jpj),                       &
                fslownflux(jpi,jpj),fslowcflux(jpi,jpj),               &
+               fregen(jpi,jpj),fregensi(jpi,jpj),                     &
                fregenfast(jpi,jpj),fregenfastsi(jpi,jpj),             &
 # if defined key_roam
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 8011)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 8012)
@@ -37,17 +37,28 @@
                                    f_fbenin_n, f_fbenin_si,                &
                                    f_omcal, fcaco3,                        &
-                                   fccd, fccd_dep, fdep1,                  &
-                                   fdpd, fdpds, fdzme, fgmepds,            &
+                                   fccd, fccd_dep, fdep1, fdd,             &
+                                   fdpd, fdpd2, fdpds, fdpds2,             &
+                                   fdpn, fdpn2,                            &
+                                   fdzme, fdzme2, fdzmi, fdzmi2,           &
                                    ffast2slowc, ffast2slown,               &
                                    ffastc, ffastca, ffastfe, ffastn,       &
                                    ffastsi,                                &
+                                   fgmed, fgmepd, fgmepds, fgmepn,         &
+                                   fgmezmi,                                &
+                                   fgmid, fgmipn,                          &
+                                   ficme, ficmi,                           &
                                    fifd_c, fifd_fe, fifd_n, fifd_si,       &
+                                   finme, finmi,                           &
+                                   fmeexcr, fmiexcr,                       &
                                    fofd_fe, fofd_n, fofd_si,               &
                                    fprotf,                                 &
-                                   fregenfast, fregenfastsi,               &
+                                   fregen, fregenfast, fregenfastsi,       &
+                                   fregensi,                               &
                                    freminc, freminca, freminfe,            &
                                    freminn, freminsi,                      &
+                                   fsdiss,                                 &
                                    fsedc, fsedca, fsedn, fsedfe, fsedsi,   &
-                                   fslowc, fslown,                         &
+                                   fslowc, fslowcflux,                     &
+                                   fslown, fslownflux,                     &
                                    ftempc, ftempca, ftempfe, ftempn,       &
                                    ftempsi,                                &
@@ -55,5 +66,6 @@
                                    fofd_c, fregenfastc,                    &
 # endif
-                                   idf, idfval
+                                   idf, idfval,                            &
+                                   zdet, zdtc
       USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,         &
                                    gphit, mbathy, tmask
@@ -65,11 +77,15 @@
                                    jexport, jfdfate, jinorgben, jocalccd,  &
                                    jorgben, jp_tem, jrratio,               &
-                                   ocal_ccd, xcaco3a, xcaco3b,             &
+                                   ocal_ccd, vsed,                         &
+                                   xbetac, xbetan,                         &
+                                   xcaco3a, xcaco3b,                       &
                                    xfastc, xfastca, xfastsi,               &
                                    xfdfrac1, xfdfrac2, xfdfrac3,           &
                                    xmassc, xmassca, xmasssi,               &
-                                   xprotca, xprotsi, xrfn, xridg_r0,       &
+                                   xphi, xprotca, xprotsi,                 &
+                                   xrfn, xridg_r0,                         &
                                    xsedc, xsedca, xsedfe,xsedn, xsedsi,    &
-                                   xthetapd, xthetazme,                    &
+                                   xthetapd, xthetapn,                     &
+                                   xthetazme, xthetazmi,                   &
                                    zn_sed_c, zn_sed_ca, zn_sed_fe,         &
                                    zn_sed_n, zn_sed_si
@@ -89,4 +105,93 @@
       !! temporary variables
       REAL(wp) :: fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8
+
+      !! The two sections below, slow detritus creation and Nutrient 
+      !! regeneration are moved from just above the CALL to DETRITUS
+      !! in trcbio_medusa.F90.
+      !!---------------------------------------------------------
+      !! Slow detritus creation
+      !!---------------------------------------------------------
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,1) == 1) THEN
+               !!
+               !! this variable integrates the creation of slow sinking 
+               !! detritus to allow the split between fast and slow 
+               !! detritus to be diagnosed
+               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) +                 &
+                                ((1.0 - xfdfrac1) * fdpd(ji,jj)) +           &
+                                ((1.0 - xfdfrac2) * fdzme(ji,jj)) +          &
+                                ((1.0 - xbetan) *                            &
+                                 (finmi(ji,jj) + finme(ji,jj)))
+               !!
+               !! this variable records the slow detrital sinking flux at 
+               !! this particular depth; it is used in the output of this 
+               !! flux at standard depths in the diagnostic outputs; 
+               !! needs to be adjusted from per second to per day because 
+               !! of parameter vsed
+               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400.
+# if defined key_roam
+               !!
+               !! and the same for detrital carbon
+               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) +                   &
+                                (xthetazmi * fdzmi(ji,jj)) +                 &
+                                (xthetapd * (1.0 - xfdfrac1) *               &
+                                 fdpd(ji,jj)) +                              &
+                                (xthetazme * (1.0 - xfdfrac2) *              &
+                                 fdzme(ji,jj)) +                             &
+                                ((1.0 - xbetac) * (ficmi(ji,jj) +            &
+                                                   ficme(ji,jj)))
+               !!
+               !! this variable records the slow detrital sinking flux 
+               !! at this particular depth; it is used in the output of 
+               !! this flux at standard depths in the diagnostic 
+               !! outputs; needs to be adjusted from per second to per 
+               !! day because of parameter vsed
+               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400.
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      !!---------------------------------------------------------
+      !! Nutrient regeneration
+      !! this variable integrates total nitrogen regeneration down the
+      !! watercolumn; its value is stored and output as a 2D diagnostic;
+      !! the corresponding dissolution flux of silicon (from sources
+      !! other than fast detritus) is also integrated; note that,
+      !! confusingly, the linear loss terms from plankton compartments
+      !! are labelled as fdX2 when one might have expected fdX or fdX1
+      !!---------------------------------------------------------
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,1) == 1) THEN
+               !!
+               !! nitrogen
+               fregen(ji,jj) =                                             &
+                                     ! messy feeding
+                               (((xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + &
+                                 (xphi *                                   &
+                                  (fgmepn(ji,jj) + fgmepd(ji,jj) +         &
+                                   fgmezmi(ji,jj) + fgmed(ji,jj))) +       &
+                                     ! excretion + D remin.
+                                 fmiexcr(ji,jj) + fmeexcr(ji,jj) +         &
+                                 fdd(ji,jj) +                              &
+                                     ! linear mortality
+                                 fdpn2(ji,jj) + fdpd2(ji,jj) +             &
+                                 fdzmi2(ji,jj) + fdzme2(ji,jj)) *          &
+                                fse3t(ji,jj,jk))
+               !!
+               !! silicon
+               fregensi(ji,jj) =                                           &
+                                     ! dissolution + non-lin. mortality
+                                 ((fsdiss(ji,jj) +                         &
+                                   ((1.0 - xfdfrac1) * fdpds(ji,jj)) +     &
+                                     ! egestion by zooplankton
+                                   ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +   &
+                                     ! linear mortality
+                                   fdpds2(ji,jj)) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
 
       !!-------------------------------------------------------------------
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 8011)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 8012)
@@ -35,5 +35,5 @@
       USE bio_medusa_mod,    ONLY: fdpd, fdpd2, fdpds, fdpds2,             &
                                    fdpn, fdpn2, fdzme, fdzme2,             &
-                                   fdzmi, fdzmi2, fsin,                    &
+                                   fdzmi, fdzmi2, fsdiss, fsin,            &
                                    zphd, zphn, zpds, zzme, zzmi
       USE dom_oce,           ONLY: tmask
@@ -43,5 +43,5 @@
                                    xkphd, xkphn, xkzme, xkzmi,             &
                                    xmetapd, xmetapn, xmetazme, xmetazmi,   &
-                                   xmpd, xmpn, xmzme, xmzmi
+                                   xmpd, xmpn, xmzme, xmzmi, xsdiss
       USE zooplankton_mod,   ONLY: zooplankton
 
@@ -178,4 +178,14 @@
       ENDDO
 
+      !! diatom frustule dissolution. This section is moved from just
+      !! below CALL to iron_chem_scav in trcbio_medusa.F90 - marc 9/5/17
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               fsdiss(ji,jj)  = xsdiss * zpds(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
    END SUBROUTINE plankton
 
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8011)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8012)
@@ -243,5 +243,5 @@
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fslown, fslowc
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux
-      REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi
 # if defined key_roam
@@ -397,8 +397,8 @@
       ENDIF
 
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! b0 is present for debugging purposes; using b0 = 0 sets the tendency
       !! terms of all biological equations to 0.
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !!
       !! AXY (03/09/14): probably not the smartest move ever, but it'll fit
@@ -410,18 +410,18 @@
       b0 = 1.
 # endif
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! fast detritus ballast scheme (0 = no; 1 = yes)
       !! alternative to ballast scheme is same scheme but with no ballast
       !! protection (not dissimilar to Martin et al., 1987)
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !!
       iball = 1
 
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
       !! these should *only* be used in 1D since they give comprehensive
       !! output for ecological functions in the model; primarily used in
       !! debugging
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !!
       idf    = 0
@@ -434,17 +434,9 @@
       endif
 
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! Initialise arrays to zero and set up arrays for diagnostics
-      !!----------------------------------------------------------------------
-! tmp - marc
-      write(numout,*) 'bbb13. before call to bio_medusa_init,kt=',kt
-      flush(numout)
-!
+      !!------------------------------------------------------------------
       CALL bio_medusa_init( kt )
-! tmp - marc
-      write(numout,*) 'bbb14. after call to bio_medusa_init,kt=',kt
-      flush(numout)
-!
-       !!
+
 # if defined key_axy_nancheck
        DO jn = 1,jptra
@@ -452,12 +444,14 @@
           !! fq1 = MAXVAL(trn(:,:,:,jn))
           fq2 = SUM(trn(:,:,:,jn))
-          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', &
-          !! &        kt, jn, fq0, fq1, fq2
-          !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR
-          !!                 and has been replaced here with a specialist routine
+          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK',     &
+          !!                kt, jn, fq0, fq1, fq2
+          !! AXY (30/01/14): much to our surprise, the next line doesn't 
+          !!                 work on HECTOR and has been replaced here with 
+          !!                 a specialist routine
           !! if (fq2 /= fq2 ) then
           if ( ieee_is_nan( fq2 ) ) then
              !! there's a NaN here
-             if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:'
+             if (lwp) write(numout,*) 'NAN detected in field', jn,           &
+                                      'at time', kt, 'at position:'
              DO jk = 1,jpk
                 DO jj = 1,jpj
@@ -466,6 +460,6 @@
                       !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then
                       if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
-                         if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', &
-                         &        tmask(ji,jj,jk), ji, jj, jk, jn
+                         if (lwp) write (numout,'(a,1pe12.2,4i6)')           &
+                            'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn
                       endif
                    enddo
@@ -479,12 +473,13 @@
 
 # if defined key_debug_medusa
-      IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked'
+      IF (lwp) write (numout,*)                                              &
+                     'trc_bio_medusa: variables initialised and checked'
       CALL flush(numout)
 # endif 
 
 # if defined key_roam
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! calculate atmospheric pCO2
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !!
       !! what's atmospheric pCO2 doing? (data start in 1859)
@@ -504,7 +499,10 @@
          !! AXY (14/06/12): tweaked to make more sense (and be correct)
 #  if defined key_bs_axy_yrlen
-         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  !! bugfix: for 360d year with HadGEM2-ES forcing
+         !! bugfix: for 360d year with HadGEM2-ES forcing
+         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  
 #  else
-         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0  !! original use of 365 days (not accounting for leap year or 360d year)
+         !! original use of 365 days (not accounting for leap year or 
+         !! 360d year)
+         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0
 #  endif
          fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
@@ -512,5 +510,6 @@
       endif
 #  if defined key_axy_pi_co2
-      f_xco2a(:,:) = 284.725          !! OCMIP pre-industrial pCO2
+      !! OCMIP pre-industrial pCO2
+      f_xco2a(:,:) = 284.725
 #  endif
       !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
@@ -540,17 +539,19 @@
       !!=============================
       !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call :
-      !!          we don't want to call on the first time-step of all run submission, 
-      !!          but only on the very first time-step, and then every month
-      !!          So we call on nittrc000 if not restarted run, 
-      !!          else if one month after last call.
-      !!          assume one month is 30d --> 3600*24*30 : 2592000s
-      !!          try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt)   
+      !!          we don't want to call on the first time-step of all run 
+      !!          submission, but only on the very first time-step, and 
+      !!          then every month. So we call on nittrc000 if not 
+      !!          restarted run, else if one month after last call.
+      !!          Assume one month is 30d --> 3600*24*30 : 2592000s
+      !!          try to call carb-chem at 1st month's tm-stp : 
+      !!          x * 30d + 1*rdt(i.e: mod = rdt)   
       !!          ++ need to pass carb-chem output var through restarts
-      If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN
-         !!----------------------------------------------------------------------
+      If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR.                      &
+           ( mod(kt*rdt,2592000.) == rdt ) ) THEN
+         !!---------------------------------------------------------------
          !! Calculate the carbonate chemistry for the whole ocean on the first
          !! simulation timestep and every month subsequently; the resulting 3D
          !! field of omega calcite is used to determine the depth of the CCD
-         !!----------------------------------------------------------------------
+         !!---------------------------------------------------------------
          CALL carb_chem( kt )
 
@@ -563,9 +564,9 @@
 # endif 
 
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! MEDUSA has unified equation through the water column
       !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers) 
       !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1          
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !!
       !! NOTE: the ordering of the loops below differs from that of some other
@@ -583,391 +584,259 @@
          !! OPEN horizontal loops
          DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,jk) == 1) then               
-               !!===========================================================
-               !! SETUP LOCAL GRID CELL
-               !!===========================================================
-               !!
-               !!-----------------------------------------------------------
-               !! Some notes on grid vertical structure
-               !! - fsdepw(ji,jj,jk) is the depth of the upper surface of 
-               !!   level jk
-               !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of 
-               !!   level jk
-               !! - fse3t(ji,jj,jk)  is the thickness of level jk
-               !!-----------------------------------------------------------
-               !!
-               !! AXY (01/03/10): set up level depth (bottom of level)
-               fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
-               !! AXY (28/11/16): local seafloor depth
-               !!                 previously mbathy(ji,jj) - 1, now 
-               !!                 mbathy(ji,jj)
-! I should be able to remove this - marc 5/5/17
-!               mbathy(ji,jj) = mbathy(ji,jj)
-               !!
-               !! set up model tracers
-               !! negative values of state variables are not allowed to
-               !! contribute to the calculated fluxes
-               !! non-diatom chlorophyll
-               zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
-               !! diatom chlorophyll
-               zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
-               !! non-diatoms
-               zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
-               !! diatoms
-               zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
-               !! diatom silicon
-               zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
-               !! AXY (28/01/10): probably need to take account of 
-               !! chl/biomass connection
-               if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
-               if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
-               if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
-               if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
-	       !! AXY (23/01/14): duh - why did I forget diatom silicon?
-	       if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
-	       if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
-               !! microzooplankton
-               zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
-               !! mesozooplankton
-               zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
-               !! detrital nitrogen
-               zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
-               !! dissolved inorganic nitrogen
-               zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
-               !! dissolved silicic acid
-               zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
-               !! dissolved "iron"
-               zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
+            DO ji = 2,jpim1
+               !! OPEN wet point IF..THEN loop
+               if (tmask(ji,jj,jk) == 1) then               
+                  !!=========================================================
+                  !! SETUP LOCAL GRID CELL
+                  !!=========================================================
+                  !!
+                  !!---------------------------------------------------------
+                  !! Some notes on grid vertical structure
+                  !! - fsdepw(ji,jj,jk) is the depth of the upper surface of 
+                  !!   level jk
+                  !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of 
+                  !!   level jk
+                  !! - fse3t(ji,jj,jk)  is the thickness of level jk
+                  !!---------------------------------------------------------
+                  !!
+                  !! AXY (01/03/10): set up level depth (bottom of level)
+                  fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
+                  !!
+                  !! set up model tracers
+                  !! negative values of state variables are not allowed to
+                  !! contribute to the calculated fluxes
+                  !! non-diatom chlorophyll
+                  zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
+                  !! diatom chlorophyll
+                  zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
+                  !! non-diatoms
+                  zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
+                  !! diatoms
+                  zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
+                  !! diatom silicon
+                  zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
+                  !! AXY (28/01/10): probably need to take account of 
+                  !! chl/biomass connection
+                  if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
+                  if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
+                  if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
+                  if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
+	          !! AXY (23/01/14): duh - why did I forget diatom silicon?
+	          if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
+	          if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
+               ENDIF
+            ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  !! microzooplankton
+                  zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
+                  !! mesozooplankton
+                  zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
+                  !! detrital nitrogen
+                  zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
+                  !! dissolved inorganic nitrogen
+                  zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
+                  !! dissolved silicic acid
+                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
+                  !! dissolved "iron"
+                  zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
+               ENDIF
+            ENDDO
+         ENDDO
+
 # if defined key_roam
-               !! detrital carbon
-               zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
-               !! dissolved inorganic carbon
-               zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
-               !! alkalinity
-               zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
-               !! oxygen
-               zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  !! detrital carbon
+                  zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
+                  !! dissolved inorganic carbon
+                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
+                  !! alkalinity
+                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
+                  !! oxygen
+                  zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
 #  if defined key_axy_carbchem && defined key_mocsy
-               !! phosphate via DIN and Redfield
-               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
+                  !! phosphate via DIN and Redfield
+                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
 #  endif
-               !!
-               !! also need physical parameters for gas exchange calculations
-               ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
-               zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
-               !!
-	       !! AXY (28/02/14): check input fields
-               if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
-                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ',  &
-                     tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',     &
-                     ji, ',', jj, ',', jk, ') at time', kt
-		  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',&
-                     tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
-                  !! temperatur
-                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
-               endif
-               if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
-                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', &
-                     tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    &
-                     ji, ',', jj, ',', jk, ') at time', kt
-               endif
+                  !!
+                  !! also need physical parameters for gas exchange 
+                  !! calculations
+                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
+                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
+                  !!
+	          !! AXY (28/02/14): check input fields
+                  if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
+                     IF(lwp) WRITE(numout,*)                                 &
+                        ' trc_bio_medusa: T WARNING 2D, ',                   &
+                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
+                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
+		     IF(lwp) WRITE(numout,*)                                 &
+                        ' trc_bio_medusa: T SWITCHING 2D, ',                 &
+                        tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
+                     !! temperatur
+                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
+                  endif
+                  if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
+                     IF(lwp) WRITE(numout,*)                                 &
+                        ' trc_bio_medusa: S WARNING 2D, ',                   &
+                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
+                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
+                  endif
+               ENDIF
+            ENDDO
+         ENDDO
 # else
-               !! implicit detrital carbon
-               zdtc(ji,jj) = zdet(ji,jj) * xthetad
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  !! implicit detrital carbon
+                  zdtc(ji,jj) = zdet(ji,jj) * xthetad
+               ENDIF
+            ENDDO
+         ENDDO
 # endif
 # if defined key_debug_medusa
-               if (idf.eq.1) then
-                  !! AXY (15/01/10)
-                  if (trn(ji,jj,jk,jpdin).lt.0.) then
-                     IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',        &
-                        trn(ji,jj,jk,jpdin)
-                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',        &
-                        ji, jj, jk, kt
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  if (idf.eq.1) then
+                     !! AXY (15/01/10)
+                     if (trn(ji,jj,jk,jpdin).lt.0.) then
+                        IF (lwp) write (numout,*)                            &
+                           '------------------------------'
+                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',    &
+                           trn(ji,jj,jk,jpdin)
+                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',    &
+                           ji, jj, jk, kt
+                     endif
+                     if (trn(ji,jj,jk,jpsil).lt.0.) then
+                        IF (lwp) write (numout,*)                            &
+                           '------------------------------'
+                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',    &
+                           trn(ji,jj,jk,jpsil)
+                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',    &
+                           ji, jj, jk, kt
+                     endif
+#  if defined key_roam
+                     if (trn(ji,jj,jk,jpdic).lt.0.) then
+                        IF (lwp) write (numout,*)                            &
+                           '------------------------------'
+                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',    &
+                           trn(ji,jj,jk,jpdic)
+                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',    &
+                           ji, jj, jk, kt
+                     endif
+                     if (trn(ji,jj,jk,jpalk).lt.0.) then
+                        IF (lwp) write (numout,*)                            &
+                           '------------------------------'
+                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',    &
+                           trn(ji,jj,jk,jpalk)
+                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',    &
+                           ji, jj, jk, kt
+                     endif
+                     if (trn(ji,jj,jk,jpoxy).lt.0.) then
+                        IF (lwp) write (numout,*)                            &
+                           '------------------------------'
+                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',    &
+                           trn(ji,jj,jk,jpoxy)
+                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',    &
+                           ji, jj, jk, kt
+                     endif
+#  endif
                   endif
-                  if (trn(ji,jj,jk,jpsil).lt.0.) then
-                     IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',        &
-                        trn(ji,jj,jk,jpsil)
-                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',        &
-                        ji, jj, jk, kt
-                  endif
+               ENDIF
+            ENDDO
+         ENDDO
+# endif
+# if defined key_debug_medusa
+! I'M NOT SURE THIS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
+!         if (idf.eq.1.AND.idfval.eq.1) then
+!            DO jj = 2,jpjm1
+!               DO ji = 2,jpim1
+!                  if (tmask(ji,jj,1) == 1) then
+!                     !! report state variable values
+!                     IF (lwp) write (numout,*)                               &
+!                        '------------------------------'
+!                     IF (lwp) write (numout,*) 'fthk(',jk,') = ',            &
+!                        fse3t(ji,jj,jk)
+!                     IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
+!                     IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
+!                     IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
+!                     IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
+!                     IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
+!                     IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
+!                     IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
+!                     IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
+!                     IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
 #  if defined key_roam
-                  if (trn(ji,jj,jk,jpdic).lt.0.) then
-                     IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',        &
-                        trn(ji,jj,jk,jpdic)
-                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',        &
-                        ji, jj, jk, kt
-                  endif
-                  if (trn(ji,jj,jk,jpalk).lt.0.) then
-                     IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',        &
-                        trn(ji,jj,jk,jpalk)
-                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',        &
-                        ji, jj, jk, kt
-                  endif
-                  if (trn(ji,jj,jk,jpoxy).lt.0.) then
-                     IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',        &
-                        trn(ji,jj,jk,jpoxy)
-                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',        &
-                        ji, jj, jk, kt
-                  endif
+!                     IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
+!                     IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
+!                     IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
+!                     IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)
 #  endif
-               endif
-# endif
+!                  ENDIF
+!               ENDDO
+!            ENDDO
+!         endif
+# endif
+
 # if defined key_debug_medusa
-               !! report state variable values
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fthk(',jk,') = ', fse3t(ji,jj,jk)
-                  IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
-                  IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
-                  IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
-                  IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
-                  IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
-                  IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
-                  IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
-                  IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
-                  IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
-#  if defined key_roam
-                  IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
-                  IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
-                  IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
-                  IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)                  
-#  endif
-               endif
-# endif
-
-# if defined key_debug_medusa
-               if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
-               endif
-# endif
-
-               !! sum tracers for inventory checks
-               IF( lk_iomput ) THEN
-                  IF ( med_diag%INVTN%dgsave )   THEN
-                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        &
-                        (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) +     &
-                                             zzmi(ji,jj) + zzme(ji,jj) +     &
-                                             zdet(ji,jj) + zdin(ji,jj) ) )
-                  ENDIF
-                  IF ( med_diag%INVTSI%dgsave )  THEN
-                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       & 
-                        (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) )
-                  ENDIF
-                  IF ( med_diag%INVTFE%dgsave )  THEN
-                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       & 
-                        (fse3t(ji,jj,jk) * ( xrfn *                          &
-                                            ( zphn(ji,jj) + zphd(ji,jj) +    &
-                                              zzmi(ji,jj) + zzme(ji,jj) +    &
-                                              zdet(ji,jj) ) + zfer(ji,jj) ) )
-                  ENDIF
-# if defined key_roam
-                  IF ( med_diag%INVTC%dgsave )  THEN
-                     ftot_c(ji,jj)  = ftot_c(ji,jj) + & 
-                        (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) +      &
-                                             (xthetapd * zphd(ji,jj)) +      &
-                                             (xthetazmi * zzmi(ji,jj)) +     &
-                                             (xthetazme * zzme(ji,jj)) +     &
-                                             zdtc(ji,jj) + zdic(ji,jj) ) )
-                  ENDIF
-                  IF ( med_diag%INVTALK%dgsave ) THEN
-                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     &
-                                                       zalk(ji,jj))
-                  ENDIF
-                  IF ( med_diag%INVTO2%dgsave )  THEN
-                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    &
-                                                        zoxy(ji,jj))
-                  ENDIF
-                  !!
-                  !! AXY (10/11/16): CMIP6 diagnostics
-                  IF ( med_diag%INTDISSIC%dgsave ) THEN
-                     intdissic(ji,jj) = intdissic(ji,jj) +                   &
-                                        (fse3t(ji,jj,jk) * zdic(ji,jj))
-                  ENDIF
-                  IF ( med_diag%INTDISSIN%dgsave ) THEN
-                     intdissin(ji,jj) = intdissin(ji,jj) +                   &
-                                        (fse3t(ji,jj,jk) * zdin(ji,jj))
-                  ENDIF
-                  IF ( med_diag%INTDISSISI%dgsave ) THEN
-                     intdissisi(ji,jj) = intdissisi(ji,jj) +                 &
-                                         (fse3t(ji,jj,jk) * zsil(ji,jj))
-                  ENDIF
-                  IF ( med_diag%INTTALK%dgsave ) THEN
-                     inttalk(ji,jj) = inttalk(ji,jj) +                       &
-                                      (fse3t(ji,jj,jk) * zalk(ji,jj))
-                  ENDIF
-                  IF ( med_diag%O2min%dgsave ) THEN
-                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then
-                        o2min(ji,jj)  = zoxy(ji,jj)
-                        IF ( med_diag%ZO2min%dgsave ) THEN
-                           !! layer midpoint
-                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               &
-                                            fdep1(ji,jj)) / 2.0
-                        ENDIF
-                     endif
-                  ENDIF
-# endif
-               ENDIF
-
-               CALL flush(numout)
-
-            ENDIF
-         ENDDO
-         ENDDO
-
-         !!----------------------------------------------------------------
+! I'M NOT SURE THIS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
+!         if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
+!            DO jj = 2,jpjm1
+!               DO ji = 2,jpim1
+!                  if (tmask(ji,jj,1) == 1) then
+!                     IF (lwp) write (numout,*)                               &
+!                       '------------------------------'
+!                     IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
+!                  ENDIF
+!               ENDDO
+!            ENDDO
+!         endif
+# endif
+
+         !!---------------------------------------------------------------
          !! Calculate air-sea gas exchange and river inputs
-         !!----------------------------------------------------------------
+         !!---------------------------------------------------------------
          IF ( jk == 1 ) THEN
-            call air_sea( kt )
-         END IF
-
-# if defined key_roam
-
-! Moved from above - marc 21/4/17
-! I think this should be moved into diagnostics at bottom - it
-! doesn't like it is used anyway else - marc 21/4/17
-         DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,jk) == 1) then
-
-               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
-               IF ( med_diag%O2SAT3%dgsave ) THEN
-                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 )
-                  o2sat3(ji, jj, jk) = f_o2sat3
-               ENDIF
-            ENDIF
-         ENDDO
-         ENDDO
-# endif
-
-         !!------------------------------------------------------------------
+            CALL air_sea( kt )
+         ENDIF
+
+         !!---------------------------------------------------------------
          !! Phytoplankton growth, zooplankton grazing and miscellaneous
          !! plankton losses. 
-         !!------------------------------------------------------------------
+         !!---------------------------------------------------------------
          CALL plankton( jk )
 
-         !!----------------------------------------------------------------
+         !!---------------------------------------------------------------
          !! Iron chemistry and scavenging
-         !!----------------------------------------------------------------
+         !!---------------------------------------------------------------
          CALL iron_chem_scav( jk )
 
-! Miscellaneous processes - marc
-
-         DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,jk) == 1) then
-
-               !!----------------------------------------------------------------------
-               !! Miscellaneous
-               !!----------------------------------------------------------------------
-               !!
-               !! diatom frustule dissolution
-               fsdiss(ji,jj)  = xsdiss * zpds(ji,jj)
-
-# if defined key_debug_medusa
-               !! report miscellaneous calculations
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss(ji,jj)
-               endif
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Slow detritus creation
-               !!----------------------------------------------------------------------
-               !! this variable integrates the creation of slow sinking detritus
-               !! to allow the split between fast and slow detritus to be 
-               !! diagnosed
-               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj)) + &
-               ((1.0 - xfdfrac2) * fdzme(ji,jj)) + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))
-               !!
-               !! this variable records the slow detrital sinking flux at this
-               !! particular depth; it is used in the output of this flux at
-               !! standard depths in the diagnostic outputs; needs to be
-               !! adjusted from per second to per day because of parameter vsed
-               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400.
-# if defined key_roam
-               !!
-               !! and the same for detrital carbon
-               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) + (xthetazmi * fdzmi(ji,jj)) + &
-               (xthetapd * (1.0 - xfdfrac1) * fdpd(ji,jj)) + &
-               (xthetazme * (1.0 - xfdfrac2) * fdzme(ji,jj)) + &
-               ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))
-               !!
-               !! this variable records the slow detrital sinking flux at this
-               !! particular depth; it is used in the output of this flux at
-               !! standard depths in the diagnostic outputs; needs to be
-               !! adjusted from per second to per day because of parameter vsed
-               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400.
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Nutrient regeneration
-               !! this variable integrates total nitrogen regeneration down the
-               !! watercolumn; its value is stored and output as a 2D diagnostic;
-               !! the corresponding dissolution flux of silicon (from sources
-               !! other than fast detritus) is also integrated; note that,
-               !! confusingly, the linear loss terms from plankton compartments
-               !! are labelled as fdX2 when one might have expected fdX or fdX1
-               !!----------------------------------------------------------------------
-               !!
-               !! nitrogen
-               fregen(ji,jj)   = (( (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) +                &  ! messy feeding
-               (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) +           &  ! messy feeding
-               fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) +                                &  ! excretion + D remin.
-               fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)) * fse3t(ji,jj,jk))                    ! linear mortality
-               !!
-               !! silicon
-               fregensi(ji,jj) = (( fsdiss(ji,jj) + ((1.0 - xfdfrac1) * fdpds(ji,jj)) +      &  ! dissolution + non-lin. mortality
-               ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +                           &  ! egestion by zooplankton
-               fdpds2(ji,jj)) * fse3t(ji,jj,jk))                                             ! linear mortality
-# if defined key_roam
-               !!
-               !! carbon
-! Doesn't look this is used - marc 10/4/17
-!               fregenc(ji,jj)  = (( (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))) +  &  ! messy feeding
-!               (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +     &  ! messy feeding
-!               (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))) +                       &  ! messy feeding
-!               fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) +                               &  ! respiration + D remin.
-!               (xthetapn * fdpn2(ji,jj)) + (xthetapd * fdpd2(ji,jj)) +                &  ! linear mortality
-!               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality
-# endif
-
-            ENDIF
-         ENDDO
-         ENDDO
-
-         !!------------------------------------------------------------------
-         !! Detritus process
-         !!------------------------------------------------------------------
+         !!---------------------------------------------------------------
+         !! Detritus processes
+         !!---------------------------------------------------------------
          CALL detritus( jk, iball )
 
-         !!------------------------------------------------------------------
+         !!---------------------------------------------------------------
          !! Updating tracers
-         !!------------------------------------------------------------------
+         !!---------------------------------------------------------------
          CALL bio_medusa_update( kt, jk )
 
-         !!------------------------------------------------------------------
+         !!---------------------------------------------------------------
          !! Diagnostics
-         !!------------------------------------------------------------------
+         !!---------------------------------------------------------------
          CALL bio_medusa_diag( kt, jk )
 
+         !!-------------------------------------------------------
+         !! 2d specific k level diags
+         !!-------------------------------------------------------
          IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
-
-            !!-------------------------------------------------------
-            !! 2d specific k level diags
-            !!-------------------------------------------------------
             CALL bio_medusa_diag_slice( jk )
-
          ENDIF
 
@@ -975,7 +844,7 @@
       ENDDO
 
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       !! Final calculations for diagnostics
-      !!----------------------------------------------------------------------
+      !!------------------------------------------------------------------
       CALL bio_medusa_fin( kt )
 
