Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 8025)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 8025)
@@ -0,0 +1,723 @@
+MODULE bio_med_diag_iomput_mod
+   !!======================================================================
+   !!                         ***  MODULE bio_med_diag_iomput_mod  ***
+   !! Calculates diagnostics
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   bio_med_diag_iomput        ! Called in bio_medusa_diag.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bio_med_diag_iomput( jk )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE bio_med_diag_iomput  ***
+      !! Calculates the diagnostics used with iom_put
+      !!-------------------------------------------------------------------
+      USE bio_medusa_mod
+      USE dom_oce,           ONLY: e3t_0, e3t_n, mbathy, tmask
+      USE in_out_manager,    ONLY: lwp, numout
+      USE par_oce,           ONLY: jpim1, jpjm1
+      USE phycst,            ONLY: rsmall
+      USE sms_medusa,        ONLY: f3_omarg, f3_omcal,                   &
+                                   i0100, i0500, i1000,                  &
+                                   xbetac, xbetan, xphi,                 &
+                                   xthetapd, xthetapn, xthetazmi
+      USE trc,               ONLY: med_diag
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+
+      !! level
+      INTEGER, INTENT( in ) :: jk
+
+      !! Loop avariables
+      INTEGER :: ji, jj, jn
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !!-------------------------------------------------------
+               !! Add in XML diagnostics stuff
+               !!-------------------------------------------------------
+               !!
+               !! ** 2D diagnostics
+#   if defined key_debug_medusa
+               IF (lwp) write (numout,*)                                     &
+                  'trc_bio_medusa: diag in ij-jj-jk loop'
+               CALL flush(numout)
+#   endif
+               IF ( med_diag%PRN%dgsave ) THEN
+                   fprn2d(ji,jj) = fprn2d(ji,jj) +                           &
+                                   (fprn(ji,jj)  * zphn(ji,jj) *             &
+                                    fse3t(ji,jj,jk)) 
+               ENDIF
+               IF ( med_diag%MPN%dgsave ) THEN
+                   fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn(ji,jj) *            &
+                                                    fse3t(ji,jj,jk))
+               ENDIF
+               IF ( med_diag%PRD%dgsave ) THEN
+                   fprd2d(ji,jj) = fprd2d(ji,jj) +                           &
+                                   (fprd(ji,jj)  * zphd(ji,jj) *             &
+                                    fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%MPD%dgsave ) THEN
+                   fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd(ji,jj) *            &
+                                                    fse3t(ji,jj,jk)) 
+               ENDIF
+               !  IF( med_diag%DSED%dgsave ) THEN
+               !      CALL iom_put( "DSED"  , ftot_n )
+               !  ENDIF
+               IF( med_diag%OPAL%dgsave ) THEN
+                   fprds2d(ji,jj) = fprds2d(ji,jj) +                         &
+                                    (fprds(ji,jj) * zpds(ji,jj) *            &
+                                     fse3t(ji,jj,jk)) 
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%OPALDISS%dgsave ) THEN
+                   fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss(ji,jj) *      &
+                                                        fse3t(ji,jj,jk))  
+               ENDIF
+               IF( med_diag%GMIPn%dgsave ) THEN
+                   fgmipn2d(ji,jj) = fgmipn2d(ji,jj) +                       &
+                                     (fgmipn(ji,jj)  * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%GMID%dgsave ) THEN
+                   fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid(ji,jj) *         &
+                                                      fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%MZMI%dgsave ) THEN
+                   fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi(ji,jj) *         &
+                                                      fse3t(ji,jj,jk)) 
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%GMEPN%dgsave ) THEN
+                   fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn(ji,jj) *      &
+                                                        fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%GMEPD%dgsave ) THEN
+                   fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd(ji,jj) *      &
+                                                        fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%GMEZMI%dgsave ) THEN
+                   fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) +                     &
+                                      (fgmezmi(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%GMED%dgsave ) THEN
+                   fgmed2d(ji,jj) = fgmed2d(ji,jj) +                         &
+                                       (fgmed(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%MZME%dgsave ) THEN
+                  fdzme2d(ji,jj) = fdzme2d(ji,jj) +                          &
+                                    (fdzme(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               !  IF( med_diag%DEXP%dgsave ) THEN
+               !      CALL iom_put( "DEXP"  , ftot_n )
+               !  ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%DETN%dgsave ) THEN
+                  fslown2d(ji,jj) = fslown2d(ji,jj) +                        &
+                                    (fslown(ji,jj) * fse3t(ji,jj,jk))  
+               ENDIF
+               IF( med_diag%MDET%dgsave ) THEN
+                  fdd2d(ji,jj) = fdd2d(ji,jj) +                              &
+                                 (fdd(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%AEOLIAN%dgsave ) THEN
+                  ffetop2d(ji,jj) = ffetop2d(ji,jj) +                        &
+                                    (ffetop(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%BENTHIC%dgsave ) THEN
+                  ffebot2d(ji,jj) = ffebot2d(ji,jj) +                        &
+                                    (ffebot(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%SCAVENGE%dgsave ) THEN
+                  ffescav2d(ji,jj) = ffescav2d(ji,jj) +                      &
+                                     (ffescav(ji,jj) * fse3t(ji,jj,jk))  
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%PN_JLIM%dgsave ) THEN
+                  ! fjln2d(ji,jj) = fjln2d(ji,jj) +                          &
+                  !                 (fjln(ji,jj)  * zphn(ji,jj) *            &
+                  !                  fse3t(ji,jj,jk)) 
+                  fjln2d(ji,jj) = fjln2d(ji,jj) +                            &
+                                  (fjlim_pn(ji,jj) * zphn(ji,jj) *           &
+                                   fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%PN_NLIM%dgsave ) THEN
+                  fnln2d(ji,jj) = fnln2d(ji,jj) +                            &
+                                  (fnln(ji,jj) * zphn(ji,jj) *               &
+                                   fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%PN_FELIM%dgsave ) THEN
+                  ffln2d(ji,jj) = ffln2d(ji,jj) +                            &
+                                  (ffln2(ji,jj) * zphn(ji,jj) *              &
+                                   fse3t(ji,jj,jk)) 
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%PD_JLIM%dgsave ) THEN
+                   ! fjld2d(ji,jj) = fjld2d(ji,jj) +                          &
+                   !                 (fjld(ji,jj)  * zphd(ji,jj) *            &
+                   !                  fse3t(ji,jj,jk)) 
+                   fjld2d(ji,jj) = fjld2d(ji,jj) +                           &
+                                   (fjlim_pd(ji,jj) * zphd(ji,jj) *          &
+                                    fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%PD_NLIM%dgsave ) THEN
+                   fnld2d(ji,jj) = fnld2d(ji,jj) +                           &
+                                   (fnld(ji,jj) * zphd(ji,jj) *              &
+                                    fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%PD_FELIM%dgsave ) THEN
+                   ffld2d(ji,jj) = ffld2d(ji,jj) +                           &
+                                   (ffld(ji,jj) * zphd(ji,jj) *              &
+                                    fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%PD_SILIM%dgsave ) THEN
+                   fsld2d2(ji,jj) = fsld2d2(ji,jj) +                         &
+                                    (fsld2(ji,jj) * zphd(ji,jj) *            &
+                                     fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%PDSILIM2%dgsave ) THEN
+                   fsld2d(ji,jj) = fsld2d(ji,jj) +                           &
+                                   (fsld(ji,jj) * zphd(ji,jj) *              &
+                                    fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! 
+               IF( med_diag%TOTREG_N%dgsave ) THEN
+                  fregen2d(ji,jj) = fregen2d(ji,jj) + fregen(ji,jj)
+               ENDIF
+               IF( med_diag%TOTRG_SI%dgsave ) THEN
+                  fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi(ji,jj)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! 
+               IF( med_diag%FASTN%dgsave ) THEN
+                   ftempn2d(ji,jj) = ftempn2d(ji,jj) +                       &
+                                     (ftempn(ji,jj)  * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%FASTSI%dgsave ) THEN
+                   ftempsi2d(ji,jj) = ftempsi2d(ji,jj) +                     &
+                                      (ftempsi(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%FASTFE%dgsave ) THEN
+                   ftempfe2d(ji,jj) = ftempfe2d(ji,jj) +                     &
+                                      (ftempfe(ji,jj) * fse3t(ji,jj,jk))  
+               ENDIF
+               IF( med_diag%FASTC%dgsave ) THEN
+                   ftempc2d(ji,jj) = ftempc2d(ji,jj) +                       &
+                                     (ftempc(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%FASTCA%dgsave ) THEN
+                   ftempca2d(ji,jj) = ftempca2d(ji,jj) +                     &
+                                      (ftempca(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! 
+               IF( med_diag%REMINN%dgsave ) THEN
+                   freminn2d(ji,jj) = freminn2d(ji,jj) +                     &
+                                      (freminn(ji,jj)  * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%REMINSI%dgsave ) THEN
+                   freminsi2d(ji,jj) = freminsi2d(ji,jj) +                   &
+                                       (freminsi(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%REMINFE%dgsave ) THEN
+                   freminfe2d(ji,jj) = freminfe2d(ji,jj) +                   &
+                                       (freminfe(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%REMINC%dgsave ) THEN
+                   freminc2d(ji,jj) = freminc2d(ji,jj) +                     &
+                                      (freminc(ji,jj)  * fse3t(ji,jj,jk)) 
+               ENDIF
+               IF( med_diag%REMINCA%dgsave ) THEN
+                   freminca2d(ji,jj) = freminca2d(ji,jj) +                   &
+                                       (freminca(ji,jj) * fse3t(ji,jj,jk)) 
+               ENDIF
+               !!
+            ENDIF
+         ENDDO
+      ENDDO
+
+# if defined key_roam
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !!
+               !! AXY (09/11/16): CMIP6 diagnostics
+               IF( med_diag%FD_NIT3%dgsave ) THEN
+                  fd_nit3(ji,jj,jk) = ffastn(ji,jj)
+               ENDIF
+               IF( med_diag%FD_SIL3%dgsave ) THEN
+                  fd_sil3(ji,jj,jk) = ffastsi(ji,jj)
+               ENDIF
+               IF( med_diag%FD_CAR3%dgsave ) THEN
+                  fd_car3(ji,jj,jk) = ffastc(ji,jj)
+               ENDIF
+               IF( med_diag%FD_CAL3%dgsave ) THEN
+                  fd_cal3(ji,jj,jk) = ffastca(ji,jj)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      IF (jk.eq.i0100) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  IF( med_diag%RR_0100%dgsave ) THEN
+                     ffastca2d(ji,jj) =                                      &
+                                ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+      ELSE IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  IF( med_diag%RR_0500%dgsave ) THEN
+                     ffastca2d(ji,jj) =                                      &
+                                ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+      ELSE IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  IF( med_diag%RR_1000%dgsave ) THEN
+                     ffastca2d(ji,jj) =                                      &
+                                ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+      ELSE IF (jk.eq.mbathy(ji,jj)) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  IF( med_diag%IBEN_N%dgsave ) THEN
+                     iben_n2d(ji,jj) = f_sbenin_n(ji,jj) +                   &
+                                       f_fbenin_n(ji,jj)
+                  ENDIF
+                  IF( med_diag%IBEN_FE%dgsave ) THEN
+                     iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) +                 &
+                                        f_fbenin_fe(ji,jj)
+                  ENDIF
+                  IF( med_diag%IBEN_C%dgsave ) THEN
+                     iben_c2d(ji,jj) = f_sbenin_c(ji,jj) +                   &
+                                       f_fbenin_c(ji,jj)
+                  ENDIF
+                  IF( med_diag%IBEN_SI%dgsave ) THEN
+                     iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
+                  ENDIF
+                  IF( med_diag%IBEN_CA%dgsave ) THEN
+                     iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
+                  ENDIF
+                  IF( med_diag%OBEN_N%dgsave ) THEN
+                     oben_n2d(ji,jj) = f_benout_n(ji,jj)
+                  ENDIF
+                  IF( med_diag%OBEN_FE%dgsave ) THEN
+                     oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
+                  ENDIF
+                  IF( med_diag%OBEN_C%dgsave ) THEN
+                     oben_c2d(ji,jj) = f_benout_c(ji,jj)
+                  ENDIF
+                  IF( med_diag%OBEN_SI%dgsave ) THEN
+                     oben_si2d(ji,jj) = f_benout_si(ji,jj)
+                  ENDIF
+                  IF( med_diag%OBEN_CA%dgsave ) THEN
+                     oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
+                  ENDIF
+                  IF( med_diag%SFR_OCAL%dgsave ) THEN
+                     sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
+                  ENDIF
+                  IF( med_diag%SFR_OARG%dgsave ) THEN
+                     sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
+                  ENDIF
+                  IF( med_diag%LYSO_CA%dgsave ) THEN
+                     lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+      !! end bathy-1 diags
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !!
+               IF( med_diag%RIV_N%dgsave ) THEN
+                  rivn2d(ji,jj) = rivn2d(ji,jj) +                            &
+                                  (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%RIV_SI%dgsave ) THEN
+                  rivsi2d(ji,jj) = rivsi2d(ji,jj) +                          &
+                                   (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%RIV_C%dgsave ) THEN
+                  rivc2d(ji,jj) = rivc2d(ji,jj) +                            &
+                                  (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%RIV_ALK%dgsave ) THEN
+                  rivalk2d(ji,jj) = rivalk2d(ji,jj) +                        &
+                                    (f_riv_loc_alk(ji,jj) *                  &
+                                     fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%DETC%dgsave ) THEN
+                  fslowc2d(ji,jj) = fslowc2d(ji,jj) +                        &
+                                    (fslowc(ji,jj)  * fse3t(ji,jj,jk))   
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !!
+               IF( med_diag%PN_LLOSS%dgsave ) THEN
+                  fdpn22d(ji,jj) = fdpn22d(ji,jj) +                          &
+                                   (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%PD_LLOSS%dgsave ) THEN
+                  fdpd22d(ji,jj) = fdpd22d(ji,jj) +                          &
+                                   (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%ZI_LLOSS%dgsave ) THEN
+                  fdzmi22d(ji,jj) = fdzmi22d(ji,jj) +                        &
+                                    (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_LLOSS%dgsave ) THEN
+                  fdzme22d(ji,jj) = fdzme22d(ji,jj) +                        &
+                                    (fdzme2(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%ZI_MES_N%dgsave ) THEN
+                  zimesn2d(ji,jj) = zimesn2d(ji,jj) +                        &
+                                    (xphi * (fgmipn(ji,jj) +                 &
+                                             fgmid(ji,jj)) *                 &
+                                     fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZI_MES_D%dgsave ) THEN
+                  zimesd2d(ji,jj) = zimesd2d(ji,jj) +                        & 
+                                    ((1. - xbetan) * finmi(ji,jj) *          &
+                                     fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZI_MES_C%dgsave ) THEN
+                  zimesc2d(ji,jj) = zimesc2d(ji,jj) +                        &
+                                    (xphi * ((xthetapn * fgmipn(ji,jj)) +    &
+                                             fgmidc(ji,jj)) *                &
+                                             fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZI_MESDC%dgsave ) THEN
+                  zimesdc2d(ji,jj) = zimesdc2d(ji,jj) +                      &
+                                     ((1. - xbetac) * ficmi(ji,jj) *         &
+                                      fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZI_EXCR%dgsave ) THEN
+                  ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +                        &
+                                    (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZI_RESP%dgsave ) THEN
+                  ziresp2d(ji,jj) = ziresp2d(ji,jj) +                        &
+                                    (fmiresp(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZI_GROW%dgsave ) THEN
+                  zigrow2d(ji,jj) = zigrow2d(ji,jj) +                        &
+                                    (fmigrow(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               IF( med_diag%ZE_MES_N%dgsave ) THEN
+                  zemesn2d(ji,jj) = zemesn2d(ji,jj) +                        &
+                                    (xphi *                                  &
+                                     (fgmepn(ji,jj) + fgmepd(ji,jj) +        &
+                                      fgmezmi(ji,jj) + fgmed(ji,jj)) *       &
+                                     fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_MES_D%dgsave ) THEN
+                  zemesd2d(ji,jj) = zemesd2d(ji,jj) +                        &
+                                    ((1. - xbetan) * finme(ji,jj) *          &
+                                     fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_MES_C%dgsave ) THEN
+                  zemesc2d(ji,jj) = zemesc2d(ji,jj) +                        & 
+                                    (xphi *                                  &
+                                     ((xthetapn * fgmepn(ji,jj)) +           &
+                                      (xthetapd * fgmepd(ji,jj)) +           &
+                                      (xthetazmi * fgmezmi(ji,jj)) +         &
+                                      fgmedc(ji,jj)) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_MESDC%dgsave ) THEN
+                  zemesdc2d(ji,jj) = zemesdc2d(ji,jj) +                      &
+                                     ((1. - xbetac) * ficme(ji,jj) *         &
+                                      fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_EXCR%dgsave ) THEN
+                  zeexcr2d(ji,jj) = zeexcr2d(ji,jj) +                        &
+                                    (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_RESP%dgsave ) THEN
+                  zeresp2d(ji,jj) = zeresp2d(ji,jj) +                        &
+                                    (fmeresp(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%ZE_GROW%dgsave ) THEN
+                  zegrow2d(ji,jj) = zegrow2d(ji,jj) +                        &
+                                    (fmegrow(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+              IF( med_diag%MDETC%dgsave ) THEN
+                  mdetc2d(ji,jj) = mdetc2d(ji,jj) +                          &
+                                   (fddc(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%GMIDC%dgsave ) THEN
+                  gmidc2d(ji,jj) = gmidc2d(ji,jj) +                          &
+                                   (fgmidc(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+               IF( med_diag%GMEDC%dgsave ) THEN
+                  gmedc2d(ji,jj) = gmedc2d(ji,jj) +                          &
+                                   (fgmedc(ji,jj) * fse3t(ji,jj,jk))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+# endif                   
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !!
+               !! ** 3D diagnostics
+               IF( med_diag%TPP3%dgsave ) THEN
+                  tpp3d(ji,jj,jk) = (fprn(ji,jj) * zphn(ji,jj)) +            &
+                                    (fprd(ji,jj) * zphd(ji,jj))
+                  !CALL iom_put( "TPP3"  , tpp3d )
+               ENDIF
+               IF( med_diag%TPPD3%dgsave ) THEN
+                  tppd3(ji,jj,jk) = (fprd(ji,jj) * zphd(ji,jj))
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+                  
+               IF( med_diag%REMIN3N%dgsave ) THEN
+                  !! remineralisation
+                  remin3dn(ji,jj,jk) = fregen(ji,jj) +                       &
+                                       (freminn(ji,jj) * fse3t(ji,jj,jk))
+                  !CALL iom_put( "REMIN3N"  , remin3dn )
+               ENDIF
+               !! IF( med_diag%PH3%dgsave ) THEN
+               !!     CALL iom_put( "PH3"  , f3_pH )
+               !! ENDIF
+               !! IF( med_diag%OM_CAL3%dgsave ) THEN
+               !!     CALL iom_put( "OM_CAL3"  , f3_omcal )
+               !! ENDIF
+	       !! 
+	       !! AXY (09/11/16): CMIP6 diagnostics
+	       IF ( med_diag%DCALC3%dgsave   ) THEN
+                  dcalc3(ji,jj,jk) = freminca(ji,jj)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+	       IF ( med_diag%FEDISS3%dgsave  ) THEN
+                  fediss3(ji,jj,jk) = ffetop(ji,jj)
+               ENDIF
+	       IF ( med_diag%FESCAV3%dgsave  ) THEN
+                  fescav3(ji,jj,jk) = ffescav(ji,jj)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+	       IF ( med_diag%MIGRAZP3%dgsave ) THEN
+                  migrazp3(ji,jj,jk) = fgmipn(ji,jj) * xthetapn
+               ENDIF
+	       IF ( med_diag%MIGRAZD3%dgsave ) THEN
+                  migrazd3(ji,jj,jk) = fgmidc(ji,jj)
+               ENDIF
+	       IF ( med_diag%MEGRAZP3%dgsave ) THEN
+                  megrazp3(ji,jj,jk) = (fgmepn(ji,jj) * xthetapn) +          &
+                                       (fgmepd(ji,jj) * xthetapd)
+               ENDIF
+	       IF ( med_diag%MEGRAZD3%dgsave ) THEN
+                  megrazd3(ji,jj,jk) = fgmedc(ji,jj)
+               ENDIF
+	       IF ( med_diag%MEGRAZZ3%dgsave ) THEN
+                   megrazz3(ji,jj,jk) = (fgmezmi(ji,jj) * xthetazmi)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+	       IF ( med_diag%PBSI3%dgsave    ) THEN
+                  pbsi3(ji,jj,jk)    = (fprds(ji,jj) * zpds(ji,jj))
+               ENDIF
+	       IF ( med_diag%PCAL3%dgsave    ) THEN
+                  pcal3(ji,jj,jk)    = ftempca(ji,jj)
+               ENDIF
+	       IF ( med_diag%REMOC3%dgsave   ) THEN
+                  remoc3(ji,jj,jk)   = freminc(ji,jj)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+	       IF ( med_diag%PNLIMJ3%dgsave  ) THEN
+                  ! pnlimj3(ji,jj,jk)  = fjln(ji,jj)
+                  pnlimj3(ji,jj,jk)  = fjlim_pn(ji,jj)
+               ENDIF
+	       IF ( med_diag%PNLIMN3%dgsave  ) THEN
+                  pnlimn3(ji,jj,jk)  = fnln(ji,jj)
+               ENDIF
+	       IF ( med_diag%PNLIMFE3%dgsave ) THEN
+                  pnlimfe3(ji,jj,jk) = ffln2(ji,jj)
+               ENDIF
+	       IF ( med_diag%PDLIMJ3%dgsave  ) THEN
+                  ! pdlimj3(ji,jj,jk)  = fjld(ji,jj)
+                  pdlimj3(ji,jj,jk)  = fjlim_pd(ji,jj)
+               ENDIF
+	       IF ( med_diag%PDLIMN3%dgsave  ) THEN
+                  pdlimn3(ji,jj,jk)  = fnld(ji,jj)
+               ENDIF
+	       IF ( med_diag%PDLIMFE3%dgsave ) THEN
+                  pdlimfe3(ji,jj,jk) = ffld(ji,jj)
+               ENDIF
+	       IF ( med_diag%PDLIMSI3%dgsave ) THEN
+                  pdlimsi3(ji,jj,jk) = fsld2(ji,jj)
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE bio_med_diag_iomput
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE bio_med_diag_iomput( )                    ! Empty routine
+      WRITE(*,*) 'bio_med_diag_iomput: You should not have seen this print! error?'
+   END SUBROUTINE bio_med_diag_iomput
+#endif 
+
+   !!======================================================================
+END MODULE bio_med_diag_iomput_mod
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90	(revision 8025)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90	(revision 8025)
@@ -0,0 +1,1215 @@
+MODULE bio_med_diag_trc_mod
+   !!======================================================================
+   !!                         ***  MODULE bio_med_diag_trc_mod  ***
+   !! Calculates diagnostics
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   bio_med_diag_trc        ! Called in bio_medusa_diag.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bio_med_diag_trc( jk )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE bio_med_diag_trc  ***
+      !! Calculates diagnostics without using iom_use
+      !!-------------------------------------------------------------------
+      USE bio_medusa_mod
+      USE dom_oce,           ONLY: e3t_0, e3t_n, mbathy, tmask
+      USE in_out_manager,    ONLY: lwp, numout
+      USE par_oce,           ONLY: jpim1, jpjm1
+      USE phycst,            ONLY: rsmall
+      USE sbc_oce,           ONLY: qsr, wndm
+      USE sms_medusa,        ONLY: f2_ccd_arg, f2_ccd_cal,               &
+                                   f3_omarg, f3_omcal, f3_pH,            &
+                                   i0100, i0150, i0200, i0500, i1000,    &
+                                   jdms, ocal_ccd,                       &
+                                   xbetac, xbetan, xpar, xphi,           &
+                                   xthetapd, xthetapn, xthetazme,        &
+                                   xthetazmi, xze
+      USE trc,               ONLY: med_diag, trc2d, trc3d 
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+
+      !! level
+      INTEGER, INTENT( in ) :: jk
+
+      !! Loop avariables
+      INTEGER :: ji, jj, jn
+
+      !!
+      !! ** Without using iom_use
+#   if defined key_debug_medusa
+      IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk ln_diatrc'
+      CALL flush(numout)
+#   endif
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) then
+               !!-------------------------------------------------------
+               !! Prepare 2D diagnostics
+               !!-------------------------------------------------------
+               !!
+               !! if ((kt / 240*240).eq.kt) then
+               !!    IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt
+               !! endif     
+               !! nitrogen inventory
+               trc2d(ji,jj,1)  =  ftot_n(ji,jj)
+               !! silicon  inventory
+               trc2d(ji,jj,2)  =  ftot_si(ji,jj)
+               !! iron     inventory
+               trc2d(ji,jj,3)  =  ftot_fe(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! non-diatom production
+               trc2d(ji,jj,4)  = trc2d(ji,jj,4)  +                           &
+                                 (fprn(ji,jj)  * zphn(ji,jj) *               &
+                                  fse3t(ji,jj,jk))
+               !! non-diatom non-grazing losses
+               trc2d(ji,jj,5)  = trc2d(ji,jj,5)  +                           &
+                                 (fdpn(ji,jj) * fse3t(ji,jj,jk))
+               !! diatom production
+               trc2d(ji,jj,6)  = trc2d(ji,jj,6)  +                           &
+                                 (fprd(ji,jj) * zphd(ji,jj) *                &
+                                  fse3t(ji,jj,jk))
+               !! diatom non-grazing losses
+               !! diagnostic field  8 is (ostensibly) supplied by trcsed.F
+               trc2d(ji,jj,7)  = trc2d(ji,jj,7)  +                           &
+                                 (fdpd(ji,jj) * fse3t(ji,jj,jk))
+               !! diatom silicon production
+               trc2d(ji,jj,9)  = trc2d(ji,jj,9)  +                           &
+                                 (fprds(ji,jj) * zpds(ji,jj) *               &
+                                  fse3t(ji,jj,jk))
+               !! diatom silicon dissolution
+               trc2d(ji,jj,10) = trc2d(ji,jj,10) +                           &
+                                 (fsdiss(ji,jj)  * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! microzoo grazing on non-diatoms
+               trc2d(ji,jj,11) = trc2d(ji,jj,11) +                           &
+                                 (fgmipn(ji,jj)  * fse3t(ji,jj,jk))
+               !! microzoo grazing on detrital nitrogen
+               trc2d(ji,jj,12) = trc2d(ji,jj,12) +                           &
+                                 (fgmid(ji,jj) * fse3t(ji,jj,jk))
+               !! microzoo non-grazing losses
+               trc2d(ji,jj,13) = trc2d(ji,jj,13) +                           &
+                                 (fdzmi(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! mesozoo  grazing on non-diatoms
+               trc2d(ji,jj,14) = trc2d(ji,jj,14) +                           &
+                                 (fgmepn(ji,jj)  * fse3t(ji,jj,jk))
+               !! mesozoo  grazing on diatoms
+               trc2d(ji,jj,15) = trc2d(ji,jj,15) +                           &
+                                 (fgmepd(ji,jj)  * fse3t(ji,jj,jk))
+               !! mesozoo  grazing on microzoo
+               trc2d(ji,jj,16) = trc2d(ji,jj,16) +                           &
+                                 (fgmezmi(ji,jj) * fse3t(ji,jj,jk))
+               !! mesozoo  grazing on detrital nitrogen
+               trc2d(ji,jj,17) = trc2d(ji,jj,17) +                           &
+                                 (fgmed(ji,jj) * fse3t(ji,jj,jk))
+               !! mesozoo  non-grazing losses
+               trc2d(ji,jj,18) = trc2d(ji,jj,18) +                           &
+                                 (fdzme(ji,jj)   * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! diagnostic field 19 is (ostensibly) supplied by trcexp.F
+               !! slow sinking detritus N production
+               trc2d(ji,jj,20) = trc2d(ji,jj,20) +                           &
+                                 (fslown(ji,jj) * fse3t(ji,jj,jk))
+               !! detrital remineralisation
+               trc2d(ji,jj,21) = trc2d(ji,jj,21) +                           &
+                                 (fdd(ji,jj) * fse3t(ji,jj,jk))
+               !! aeolian  iron addition
+               trc2d(ji,jj,22) = trc2d(ji,jj,22) +                           &
+                                 (ffetop(ji,jj) * fse3t(ji,jj,jk))
+               !! seafloor iron addition
+               trc2d(ji,jj,23) = trc2d(ji,jj,23) +                           &
+                                 (ffebot(ji,jj) * fse3t(ji,jj,jk))
+               !! "free" iron scavenging
+               trc2d(ji,jj,24) = trc2d(ji,jj,24) +                           &
+                                 (ffescav(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! non-diatom J  limitation term 
+               trc2d(ji,jj,25) = trc2d(ji,jj,25) +                           &
+                                 (fjlim_pn(ji,jj) * zphn(ji,jj) *            &
+                                  fse3t(ji,jj,jk))
+               !! non-diatom N  limitation term 
+               trc2d(ji,jj,26) = trc2d(ji,jj,26) +                           &
+                                 (fnln(ji,jj) * zphn(ji,jj) *                &
+                                  fse3t(ji,jj,jk))
+               !! non-diatom Fe limitation term 
+               trc2d(ji,jj,27) = trc2d(ji,jj,27) +                           &
+                                 (ffln2(ji,jj) * zphn(ji,jj) *               &
+                                  fse3t(ji,jj,jk))
+               !! diatom     J  limitation term 
+               trc2d(ji,jj,28) = trc2d(ji,jj,28) +                           &
+                                 (fjlim_pd(ji,jj) * zphd(ji,jj) *            &
+                                  fse3t(ji,jj,jk))
+               !! diatom     N  limitation term 
+               trc2d(ji,jj,29) = trc2d(ji,jj,29) +                           &
+                                 (fnld(ji,jj) * zphd(ji,jj) *                &
+                                  fse3t(ji,jj,jk))
+               !! diatom     Fe limitation term 
+               trc2d(ji,jj,30) = trc2d(ji,jj,30) +                           &
+                                 (ffld(ji,jj) * zphd(ji,jj) *                &
+                                  fse3t(ji,jj,jk))
+               !! diatom     Si limitation term 
+               trc2d(ji,jj,31) = trc2d(ji,jj,31) +                           &
+                                 (fsld2(ji,jj) * zphd(ji,jj) *               &
+                                  fse3t(ji,jj,jk))
+               !! diatom     Si uptake limitation term
+               trc2d(ji,jj,32) = trc2d(ji,jj,32) +                           &
+                                 (fsld(ji,jj) * zphd(ji,jj) *                &
+                                  fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      IF (jk.eq.i0100) THEN 
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at  100 m
+                  trc2d(ji,jj,33) = fslownflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at  200 m
+                  trc2d(ji,jj,34) = fslownflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at  500 m
+                  trc2d(ji,jj,35) = fslownflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at 1000 m
+                  trc2d(ji,jj,36) = fslownflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! non-fast N  full column regeneration
+               trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen(ji,jj)
+               !! non-fast Si full column regeneration
+               trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi(ji,jj)
+               !! non-fast N  regeneration to  100 m
+            ENDIF
+         ENDDO
+      ENDDO
+
+      IF (jk.eq.i0100) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  trc2d(ji,jj,39) = trc2d(ji,jj,37)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! non-fast N  regeneration to  200 m
+                  trc2d(ji,jj,40) = trc2d(ji,jj,37)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! non-fast N  regeneration to  500 m
+                  trc2d(ji,jj,41) = trc2d(ji,jj,37)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! non-fast N  regeneration to 1000 m
+                  trc2d(ji,jj,42) = trc2d(ji,jj,37)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! fast sinking detritus N production
+               trc2d(ji,jj,43) = trc2d(ji,jj,43) +                           &
+                                 (ftempn(ji,jj) * fse3t(ji,jj,jk))
+               !! fast sinking detritus Si production
+               trc2d(ji,jj,44) = trc2d(ji,jj,44) +                           &
+                                 (ftempsi(ji,jj) * fse3t(ji,jj,jk))
+               !! fast sinking detritus Fe production
+               trc2d(ji,jj,45) = trc2d(ji,jj,45) +                           &
+                                 (ftempfe(ji,jj) * fse3t(ji,jj,jk))
+               !! fast sinking detritus C production
+               trc2d(ji,jj,46) = trc2d(ji,jj,46) +                           &
+                                 (ftempc(ji,jj)  * fse3t(ji,jj,jk))
+               !! fast sinking detritus CaCO3 production
+               trc2d(ji,jj,47) = trc2d(ji,jj,47) +                           &
+                                 (ftempca(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      IF (jk.eq.i0100) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus N  flux at  100 m
+                  trc2d(ji,jj,48) = ffastn(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus N  flux at  200 m
+                  trc2d(ji,jj,49) = ffastn(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus N  flux at  500 m
+                  trc2d(ji,jj,50) = ffastn(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus N  flux at 1000 m
+                  trc2d(ji,jj,51) = ffastn(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0100) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! N  regeneration to  100 m
+                  trc2d(ji,jj,52) = fregenfast(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! N  regeneration to  200 m
+                  trc2d(ji,jj,53) = fregenfast(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! N  regeneration to  500 m
+                  trc2d(ji,jj,54) = fregenfast(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! N  regeneration to 1000 m
+                  trc2d(ji,jj,55) = fregenfast(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0100) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus Si flux at  100 m
+                  trc2d(ji,jj,56) = ffastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus Si flux at  200 m
+                  trc2d(ji,jj,57) = ffastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus Si flux at  500 m
+                  trc2d(ji,jj,58) = ffastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! fast detritus Si flux at 1000 m
+                  trc2d(ji,jj,59) = ffastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0100) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! Si regeneration to  100 m
+                  trc2d(ji,jj,60) = fregenfastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! Si regeneration to  200 m
+                  trc2d(ji,jj,61) = fregenfastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! Si regeneration to  500 m
+                  trc2d(ji,jj,62) = fregenfastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! Si regeneration to 1000 m
+                  trc2d(ji,jj,63) = fregenfastsi(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! sum of fast-sinking N  fluxes
+               trc2d(ji,jj,64) = trc2d(ji,jj,64) +                           &
+                                 (freminn(ji,jj) * fse3t(ji,jj,jk))
+               !! sum of fast-sinking Si fluxes
+               trc2d(ji,jj,65) = trc2d(ji,jj,65) +                           &
+                                 (freminsi(ji,jj) * fse3t(ji,jj,jk))
+               !! sum of fast-sinking Fe fluxes
+               trc2d(ji,jj,66) = trc2d(ji,jj,66) +                           &
+                                 (freminfe(ji,jj) * fse3t(ji,jj,jk))
+               !! sum of fast-sinking C  fluxes
+               trc2d(ji,jj,67) = trc2d(ji,jj,67) +                           &
+                                 (freminc(ji,jj) * fse3t(ji,jj,jk))
+               !! sum of fast-sinking Ca fluxes
+               trc2d(ji,jj,68) = trc2d(ji,jj,68) +                           &
+                                 (freminca(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+
+      if (jk.eq.mbathy(ji,jj)) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! N  sedimentation flux
+                  trc2d(ji,jj,69) = fsedn(ji,jj)
+                  !! Si sedimentation flux
+                  trc2d(ji,jj,70) = fsedsi(ji,jj)
+                  !! Fe sedimentation flux
+                  trc2d(ji,jj,71) = fsedfe(ji,jj)
+                  !! C  sedimentation flux
+                  trc2d(ji,jj,72) = fsedc(ji,jj)
+                  !! Ca sedimentation flux
+                  trc2d(ji,jj,73) = fsedca(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk.eq.1) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  trc2d(ji,jj,74) = qsr(ji,jj)
+                  trc2d(ji,jj,75) = xpar(ji,jj,jk)
+                  !! trc2d(ji,jj,75) = real(iters(ji,jj))
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! diagnostic fields 76 to 80 calculated below
+               !! mixed layer non-diatom production
+               trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj)
+               !! mixed layer     diatom production
+               trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+# if defined key_gulf_finland
+      if (jk.eq.1) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! Gulf of Finland check
+                  trc2d(ji,jj,83) = real(ibio_switch)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+# else
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! calcite CCD depth
+               trc2d(ji,jj,83) = ocal_ccd(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+# endif
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! last model level above calcite CCD depth
+               trc2d(ji,jj,84) = fccd(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      IF (jk.eq.1) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! surface "free" iron
+                  trc2d(ji,jj,85) = xFree(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+! I'm keeping this the same as before, but it looks like it should 
+! be i0100 and not i0200 - marc 8/5/17
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! "free" iron at  100 m
+                  trc2d(ji,jj,86) = xFree(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+
+      IF (jk.eq.i0200) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! "free" iron at  200 m
+                  trc2d(ji,jj,87) = xFree(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+
+      IF (jk.eq.i0500) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! "free" iron at  500 m
+                  trc2d(ji,jj,88) = xFree(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+
+      IF (jk.eq.i1000) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! "free" iron at 1000 m
+                  trc2d(ji,jj,89) = xFree(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+
+      IF (jk.eq.1) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! AXY (27/06/12): extract "euphotic depth"
+                  trc2d(ji,jj,90) = xze(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
+# if defined key_roam
+      if (jk .eq. 1) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! ROAM provisionally has access to a further 20 2D 
+                  !! diagnostics
+                  !! surface wind
+                  trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + wndm(ji,jj)
+                  !! atmospheric pCO2
+                  trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm(ji,jj)
+                  !! ocean pH
+                  trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph(ji,jj)
+                  !! ocean pCO2
+                  trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w(ji,jj)
+                  !! ocean H2CO3 conc.
+                  trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3(ji,jj)
+                  !! ocean HCO3 conc.
+                  trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3(ji,jj)
+                  !! ocean CO3 conc.
+                  trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3(ji,jj)
+                  !! air-sea CO2 flux
+                  trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux(ji,jj)
+               ENDIF
+            ENDDO
+        ENDDO
+
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! ocean omega calcite 
+                  trc2d(ji,jj,99)  = trc2d(ji,jj,99)  + f_omcal(ji,jj)
+                  !! ocean omega aragonite
+                  trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj)
+                  !! ocean TDIC
+                  trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC(ji,jj)
+                  !! ocean TALK
+                  trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK(ji,jj)
+                  !! surface kw660
+                  trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660(ji,jj)
+                  !! surface pressure
+                  trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0(ji,jj)
+                  !! air-sea O2 flux
+                  trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux(ji,jj)
+                  !! ocean O2 saturation
+                  trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat(ji,jj)
+                  !! depth calcite CCD
+                  trc2d(ji,jj,107) = f2_ccd_cal(ji,jj)
+                  !! depth aragonite CCD
+                  trc2d(ji,jj,108) = f2_ccd_arg(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk .eq. mbathy(ji,jj)) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! seafloor omega calcite
+                  trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)
+                  !! seafloor omega aragonite
+                  trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk.eq.i0100) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! diagnostic fields 111 to 117 calculated below
+                  !! rain ratio at  100 m
+                  trc2d(ji,jj,118) =                                         &
+                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk.eq.i0500) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! rain ratio at  500 m
+                  trc2d(ji,jj,119) =                                         &
+                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk.eq.i1000) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! rain ratio at 1000 m
+                  trc2d(ji,jj,120) =                                         &
+                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk.eq.mbathy(ji,jj)) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! AXY (18/01/12): benthic flux diagnostics
+                  trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
+                  trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
+                  trc2d(ji,jj,123) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
+                  trc2d(ji,jj,124) = f_fbenin_si(ji,jj)
+                  trc2d(ji,jj,125) = f_fbenin_ca(ji,jj)
+                  trc2d(ji,jj,126) = f_benout_n(ji,jj)
+                  trc2d(ji,jj,127) = f_benout_fe(ji,jj)
+                  trc2d(ji,jj,128) = f_benout_c(ji,jj)
+                  trc2d(ji,jj,129) = f_benout_si(ji,jj)
+                  trc2d(ji,jj,130) = f_benout_ca(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! diagnostics fields 131 to 135 calculated below
+               trc2d(ji,jj,136) = f_runoff(ji,jj)
+               !! AXY (19/07/12): amended to allow for riverine 
+               !! nutrient addition below surface
+               trc2d(ji,jj,137) = trc2d(ji,jj,137) +                         &
+                                     (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
+               trc2d(ji,jj,138) = trc2d(ji,jj,138) +                         &
+                                     (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
+               trc2d(ji,jj,139) = trc2d(ji,jj,139) +                         &
+                                     (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
+               trc2d(ji,jj,140) = trc2d(ji,jj,140) +                         &
+                                     (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
+               !! slow sinking detritus C production
+               trc2d(ji,jj,141) = trc2d(ji,jj,141) +                         &
+                                     (fslowc(ji,jj)  * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      if (jk.eq.i0100) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at  100 m
+                  trc2d(ji,jj,142) = fslowcflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      if (jk.eq.i0200) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at  200 m
+                  trc2d(ji,jj,143) = fslowcflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+
+      if (jk.eq.i0500) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at  500 m
+                  trc2d(ji,jj,144) = fslowcflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+
+      if (jk.eq.i1000) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! slow detritus flux at 1000 m
+                  trc2d(ji,jj,145) = fslowcflux(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! carbon     inventory
+               trc2d(ji,jj,146)  = trc2d(ji,jj,146)  + ftot_c(ji,jj)
+               !! alkalinity inventory
+               trc2d(ji,jj,147)  = trc2d(ji,jj,147)  + ftot_a(ji,jj)
+               !! oxygen     inventory
+               trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      if (jk.eq.mbathy(ji,jj)) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! community respiration
+               trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+	       !!
+	       !! AXY (14/02/14): a Valentines Day gift to BASIN - a 
+               !!                 shedload of new diagnostics that 
+               !!                 they'll most likely never need!
+               !!                 (actually, as with all such gifts, 
+               !!                 I'm giving them some things I'd like 
+               !!                 myself!)
+               !! 
+               !! ------------------------------------------------------
+               !! linear losses
+               !! non-diatom
+               trc2d(ji,jj,151) = trc2d(ji,jj,151) +                         &
+                                  (fdpn2(ji,jj) * fse3t(ji,jj,jk))
+               !! diatom
+               trc2d(ji,jj,152) = trc2d(ji,jj,152) +                         &
+                                  (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
+               !! microzooplankton
+               trc2d(ji,jj,153) = trc2d(ji,jj,153) +                         &
+                                  (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
+               !! mesozooplankton
+               trc2d(ji,jj,154) = trc2d(ji,jj,154) +                         &
+                                  (fdzme2(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! ------------------------------------------------------
+               !! microzooplankton grazing
+               !! microzooplankton messy -> N
+               trc2d(ji,jj,155) = trc2d(ji,jj,155) +                         &
+                                  (xphi * (fgmipn(ji,jj) +                   &
+                                           fgmid(ji,jj)) * fse3t(ji,jj,jk))
+               !! microzooplankton messy -> D
+               trc2d(ji,jj,156) = trc2d(ji,jj,156) +                         &
+                                  ((1. - xbetan) * finmi(ji,jj) *            &
+                                   fse3t(ji,jj,jk))
+               !! microzooplankton messy -> DIC
+               trc2d(ji,jj,157) = trc2d(ji,jj,157) +                         &
+                                  (xphi * ((xthetapn * fgmipn(ji,jj)) +      &
+                                           fgmidc(ji,jj)) *                  &
+                                   fse3t(ji,jj,jk))
+               !! microzooplankton messy -> Dc
+               trc2d(ji,jj,158) = trc2d(ji,jj,158) +                         &
+                                  ((1. - xbetac) * ficmi(ji,jj) *            &
+                                   fse3t(ji,jj,jk))
+               !! microzooplankton excretion
+               trc2d(ji,jj,159) = trc2d(ji,jj,159) +                         &
+                                  (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
+               !! microzooplankton respiration
+               trc2d(ji,jj,160) = trc2d(ji,jj,160) +                         &
+                                  (fmiresp(ji,jj) * fse3t(ji,jj,jk))
+               !! microzooplankton growth
+               trc2d(ji,jj,161) = trc2d(ji,jj,161) +                         &
+                                  (fmigrow(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! ------------------------------------------------------
+               !! mesozooplankton grazing
+               !! mesozooplankton messy -> N
+               trc2d(ji,jj,162) = trc2d(ji,jj,162) +                         &
+                                  (xphi *                                    &
+                                   (fgmepn(ji,jj) + fgmepd(ji,jj) +          &
+                                    fgmezmi(ji,jj) + fgmed(ji,jj)) *         &
+                                   fse3t(ji,jj,jk))
+               !! mesozooplankton messy -> D
+               trc2d(ji,jj,163) = trc2d(ji,jj,163) +                         &
+                                  ((1. - xbetan) * finme(ji,jj) *            &
+                                   fse3t(ji,jj,jk))
+               !! mesozooplankton messy -> DIC
+               trc2d(ji,jj,164) = trc2d(ji,jj,164) +                         &
+                                  (xphi *                                    &
+                                   ((xthetapn * fgmepn(ji,jj)) +             &
+                                    (xthetapd * fgmepd(ji,jj)) +             &
+                                    (xthetazmi * fgmezmi(ji,jj)) +           &
+                                   fgmedc(ji,jj)) * fse3t(ji,jj,jk))
+               !! mesozooplankton messy -> Dc
+               trc2d(ji,jj,165) = trc2d(ji,jj,165) +                         &
+                                  ((1. - xbetac) * ficme(ji,jj) *            &
+                                   fse3t(ji,jj,jk))
+               !! mesozooplankton excretion
+               trc2d(ji,jj,166) = trc2d(ji,jj,166) +                         &
+                                  (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
+               !! mesozooplankton respiration
+               trc2d(ji,jj,167) = trc2d(ji,jj,167) +                         &
+                                  (fmeresp(ji,jj) * fse3t(ji,jj,jk))
+               !! mesozooplankton growth
+               trc2d(ji,jj,168) = trc2d(ji,jj,168) +                         &
+                                  (fmegrow(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! ------------------------------------------------------
+               !! miscellaneous
+               !! detrital C remineralisation
+               trc2d(ji,jj,169) = trc2d(ji,jj,169) +                         &
+                                  (fddc(ji,jj) * fse3t(ji,jj,jk))
+               !! microzoo grazing on detrital carbon
+               trc2d(ji,jj,170) = trc2d(ji,jj,170) +                         &
+                                  (fgmidc(ji,jj)  * fse3t(ji,jj,jk))
+               !! mesozoo  grazing on detrital carbon
+               trc2d(ji,jj,171) = trc2d(ji,jj,171) +                         &
+                                  (fgmedc(ji,jj)  * fse3t(ji,jj,jk))
+               !!
+            ENDIF
+         ENDDO
+      ENDDO
+
+      !! ------------------------------------------------------
+      !!
+      !! AXY (23/10/14): extract primary production related 
+      !!                 surface fields to deal with diel 
+      !!                 cycle issues; hijacking BASIN 150m
+      !!                 diagnostics to do so (see commented 
+      !!                 out diagnostics below this section)
+      !!
+      !! extract relevant BASIN fields at 150m
+      if (jk .eq. i0150) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! Pn PP
+                  trc2d(ji,jj,172) = trc2d(ji,jj,4)
+                  !! Pn linear loss
+                  trc2d(ji,jj,173) = trc2d(ji,jj,151)
+                  !! Pn non-linear loss
+                  trc2d(ji,jj,174) = trc2d(ji,jj,5)
+                  !! Pn grazing to Zmi
+                  trc2d(ji,jj,175) = trc2d(ji,jj,11)
+                  !! Pn grazing to Zme
+                  trc2d(ji,jj,176) = trc2d(ji,jj,14)
+                  !! Pd PP
+                  trc2d(ji,jj,177) = trc2d(ji,jj,6)
+                  !! Pd linear loss
+                  trc2d(ji,jj,178) = trc2d(ji,jj,152)
+                  !! Pd non-linear loss
+                  trc2d(ji,jj,179) = trc2d(ji,jj,7)
+                  !! Pd grazing to Zme
+                  trc2d(ji,jj,180) = trc2d(ji,jj,15)
+                  !! Zmi grazing on D
+                  trc2d(ji,jj,181) = trc2d(ji,jj,12)
+                  !! Zmi grazing on Dc
+                  trc2d(ji,jj,182) = trc2d(ji,jj,170)
+                  !! Zmi messy feeding loss to N
+                  trc2d(ji,jj,183) = trc2d(ji,jj,155)
+                  !! Zmi messy feeding loss to D
+                  trc2d(ji,jj,184) = trc2d(ji,jj,156)
+                  !! Zmi messy feeding loss to DIC
+                  trc2d(ji,jj,185) = trc2d(ji,jj,157)
+                  !! Zmi messy feeding loss to Dc
+                  trc2d(ji,jj,186) = trc2d(ji,jj,158)
+                  !! Zmi excretion
+                  trc2d(ji,jj,187) = trc2d(ji,jj,159)
+                  !! Zmi respiration
+                  trc2d(ji,jj,188) = trc2d(ji,jj,160)
+                  !! Zmi growth
+                  trc2d(ji,jj,189) = trc2d(ji,jj,161)
+                  !! Zmi linear loss
+                  trc2d(ji,jj,190) = trc2d(ji,jj,153)
+                  !! Zmi non-linear loss
+                  trc2d(ji,jj,191) = trc2d(ji,jj,13)
+                  !! Zmi grazing to Zme
+                  trc2d(ji,jj,192) = trc2d(ji,jj,16)
+                  !! Zme grazing on D
+                  trc2d(ji,jj,193) = trc2d(ji,jj,17)
+                  !! Zme grazing on Dc
+                  trc2d(ji,jj,194) = trc2d(ji,jj,171)
+                  !! Zme messy feeding loss to N
+                  trc2d(ji,jj,195) = trc2d(ji,jj,162)
+                  !! Zme messy feeding loss to D
+                  trc2d(ji,jj,196) = trc2d(ji,jj,163)
+                  !! Zme messy feeding loss to DIC
+                  trc2d(ji,jj,197) = trc2d(ji,jj,164)
+                  !! Zme messy feeding loss to Dc
+                  trc2d(ji,jj,198) = trc2d(ji,jj,165)
+                  !! Zme excretion
+                  trc2d(ji,jj,199) = trc2d(ji,jj,166)
+                  !! Zme respiration
+                  trc2d(ji,jj,200) = trc2d(ji,jj,167)
+                  !! Zme growth
+                  trc2d(ji,jj,201) = trc2d(ji,jj,168)
+                  !! Zme linear loss
+                  trc2d(ji,jj,202) = trc2d(ji,jj,154)
+                  !! Zme non-linear loss
+                  trc2d(ji,jj,203) = trc2d(ji,jj,18)
+                  !! Slow detritus production, N
+                  trc2d(ji,jj,204) = trc2d(ji,jj,20)
+                  !! Slow detritus remineralisation, N
+                  trc2d(ji,jj,205) = trc2d(ji,jj,21)
+                  !! Slow detritus production, C
+                  trc2d(ji,jj,206) = trc2d(ji,jj,141)
+                  !! Slow detritus remineralisation, C
+                  trc2d(ji,jj,207) = trc2d(ji,jj,169)
+                  !! Fast detritus production, N
+                  trc2d(ji,jj,208) = trc2d(ji,jj,43)
+                  !! Fast detritus remineralisation, N
+                  trc2d(ji,jj,209) = trc2d(ji,jj,21)
+                  !! Fast detritus production, C
+                  trc2d(ji,jj,210) = trc2d(ji,jj,64)
+                  !! Fast detritus remineralisation, C
+                  trc2d(ji,jj,211) = trc2d(ji,jj,67)
+                  !! Community respiration
+                  trc2d(ji,jj,212) = trc2d(ji,jj,150)
+                  !! Slow detritus N flux at 150 m
+                  trc2d(ji,jj,213) = fslownflux(ji,jj)
+                  !! Slow detritus C flux at 150 m
+                  trc2d(ji,jj,214) = fslowcflux(ji,jj)
+                  !! Fast detritus N flux at 150 m
+                  trc2d(ji,jj,215) = ffastn(ji,jj)
+                  !! Fast detritus C flux at 150 m
+                  trc2d(ji,jj,216) = ffastc(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+
+      !! 
+      !! Jpalm (11-08-2014)
+      !! Add UKESM1 diagnoatics 
+      !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+      if ((jk .eq. 1) .and.( jdms.eq.1)) then
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! DMS surface concentration 
+                  trc2d(ji,jj,221) = dms_surf(ji,jj)
+                  !! AXY (13/03/15): add in other DMS estimates
+                  !! DMS surface concentration 
+                  trc2d(ji,jj,222) = dms_andr(ji,jj)
+                  !! DMS surface concentration 
+                  trc2d(ji,jj,223) = dms_simo(ji,jj)
+                  !! DMS surface concentration 
+                  trc2d(ji,jj,224) = dms_aran(ji,jj)
+                  !! DMS surface concentration 
+                  trc2d(ji,jj,225) = dms_hall(ji,jj)
+               ENDIF
+            ENDDO
+         ENDDO
+      endif
+# endif
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! other possible future diagnostics include:
+               !!   - integrated tracer values (esp. biological)
+               !!   - mixed layer tracer values
+               !!   - sub-surface chlorophyll maxima (plus depth)
+               !!   - different mixed layer depth criteria (T, sigma, 
+               !!     var. sigma)
+               !!-------------------------------------------------------
+               !! Prepare 3D diagnostics
+               !!-------------------------------------------------------
+               !!
+               !! primary production  
+               trc3d(ji,jj,jk,1)  = ((fprn(ji,jj) + fprd(ji,jj)) *           &
+                                     zphn(ji,jj))
+               !! detrital flux
+               trc3d(ji,jj,jk,2)  = fslownflux(ji,jj) + ffastn(ji,jj)
+               !! remineralisation
+               trc3d(ji,jj,jk,3)  = fregen(ji,jj) +                          &
+                                    (freminn(ji,jj) * fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+# if defined key_roam
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! pH
+               trc3d(ji,jj,jk,4)  = f3_pH(ji,jj,jk)
+               !! omega calcite
+               trc3d(ji,jj,jk,5)  = f3_omcal(ji,jj,jk)
+            ENDIF
+         ENDDO
+      ENDDO
+# else
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               !! fast Si flux
+               trc3d(ji,jj,jk,4)  = ffastsi(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+# endif
+
+   END SUBROUTINE bio_med_diag_trc
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE bio_med_diag_trc( )                    ! Empty routine
+      WRITE(*,*) 'bio_med_diag_trc: You should not have seen this print! error?'
+   END SUBROUTINE bio_med_diag_trc
+#endif 
+
+   !!======================================================================
+END MODULE bio_med_diag_trc_mod
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 8023)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8025)
@@ -25,29 +25,24 @@
 CONTAINS
 
-   SUBROUTINE bio_medusa_diag( kt, jk )
+   SUBROUTINE bio_medusa_diag( jk )
       !!-------------------------------------------------------------------
       !!                     ***  ROUTINE bio_medusa_diag  ***
       !! This called from TRC_BIO_MEDUSA and calculates diagnostics
       !!-------------------------------------------------------------------
+      USE bio_med_diag_iomput_mod,  ONLY: bio_med_diag_iomput
+      USE bio_med_diag_trc_mod,     ONLY: bio_med_diag_trc
       USE bio_medusa_mod
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,       &
-                                   mbathy, tmask
-      USE in_out_manager,    ONLY: lwp, numout
+      USE dom_oce,                  ONLY: e3t_0, e3t_n,                  &
+                                          gdepw_0, gdepw_n, tmask
+      USE in_out_manager,           ONLY: lwp, numout
 # if defined key_iomput
-      USE iom,               ONLY: lk_iomput
-# endif
-      USE par_kind,          ONLY: wp
-      USE par_oce,           ONLY: jpim1, jpjm1
-      USE phycst,            ONLY: rsmall
-      USE sbc_oce,           ONLY: qsr, wndm
-      USE sms_medusa,        ONLY: f2_ccd_arg, f2_ccd_cal,               &
-                                   f3_omarg, f3_omcal, f3_pH,            &
-                                   i0100, i0150, i0200, i0500, i1000,    &
-                                   jdms, ocal_ccd,                       &
-                                   xbetac, xbetan, xpar, xphi, xrfn,     &
-                                   xthetapd, xthetapn, xthetazme, xthetazmi, xze
-      USE trc,               ONLY: ln_diatrc, med_diag, trc2d, trc3d 
+      USE iom,                      ONLY: lk_iomput
+# endif
+      USE par_oce,                  ONLY: jpim1, jpjm1
+      USE sms_medusa,               ONLY: xrfn, xthetapd, xthetapn,      &
+                                          xthetazme, xthetazmi
+      USE trc,                      ONLY: ln_diatrc, med_diag 
 # if defined key_roam
-      USE trcoxy_medusa,     ONLY: oxy_sato
+      USE trcoxy_medusa,            ONLY: oxy_sato
 # endif
 
@@ -55,6 +50,4 @@
 #  include "domzgr_substitute.h90"
 
-      !! time (integer timestep)
-      INTEGER, INTENT( in ) :: kt
       !! level
       INTEGER, INTENT( in ) :: jk
@@ -202,1811 +195,17 @@
       IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
 
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               !! OPEN wet point IF..THEN loop
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !!-------------------------------------------------------
-                  !! Add in XML diagnostics stuff
-                  !!-------------------------------------------------------
-                  !!
-                  !! ** 2D diagnostics
-#   if defined key_debug_medusa
-                  IF (lwp) write (numout,*)                                  &
-                     'trc_bio_medusa: diag in ij-jj-jk loop'
-                  CALL flush(numout)
-#   endif
-                  IF ( med_diag%PRN%dgsave ) THEN
-                      fprn2d(ji,jj) = fprn2d(ji,jj) +                        &
-                                      (fprn(ji,jj)  * zphn(ji,jj) *          &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF ( med_diag%MPN%dgsave ) THEN
-                      fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn(ji,jj) *         &
-                                                       fse3t(ji,jj,jk))
-                  ENDIF
-                  IF ( med_diag%PRD%dgsave ) THEN
-                      fprd2d(ji,jj) = fprd2d(ji,jj) +                        &
-                                      (fprd(ji,jj)  * zphd(ji,jj) *          &
-                                       fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%MPD%dgsave ) THEN
-                      fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd(ji,jj) *         &
-                                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  !  IF( med_diag%DSED%dgsave ) THEN
-                  !      CALL iom_put( "DSED"  , ftot_n )
-                  !  ENDIF
-                  IF( med_diag%OPAL%dgsave ) THEN
-                      fprds2d(ji,jj) = fprds2d(ji,jj) +                      &
-                                       (fprds(ji,jj) * zpds(ji,jj) *         &
-                                        fse3t(ji,jj,jk)) 
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%OPALDISS%dgsave ) THEN
-                      fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss(ji,jj) *   &
-                                                           fse3t(ji,jj,jk))  
-                  ENDIF
-                  IF( med_diag%GMIPn%dgsave ) THEN
-                      fgmipn2d(ji,jj) = fgmipn2d(ji,jj) +                    &
-                                        (fgmipn(ji,jj)  * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%GMID%dgsave ) THEN
-                      fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid(ji,jj) *      &
-                                                         fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%MZMI%dgsave ) THEN
-                      fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi(ji,jj) *      &
-                                                         fse3t(ji,jj,jk)) 
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%GMEPN%dgsave ) THEN
-                      fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn(ji,jj) *   &
-                                                           fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%GMEPD%dgsave ) THEN
-                      fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd(ji,jj) *   &
-                                                           fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%GMEZMI%dgsave ) THEN
-                      fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) +                  &
-                                         (fgmezmi(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%GMED%dgsave ) THEN
-                      fgmed2d(ji,jj) = fgmed2d(ji,jj) +                      &
-                                       (fgmed(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%MZME%dgsave ) THEN
-                      fdzme2d(ji,jj) = fdzme2d(ji,jj) +                      &
-                                       (fdzme(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  !  IF( med_diag%DEXP%dgsave ) THEN
-                  !      CALL iom_put( "DEXP"  , ftot_n )
-                  !  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%DETN%dgsave ) THEN
-                      fslown2d(ji,jj) = fslown2d(ji,jj) +                    &
-                                        (fslown(ji,jj) * fse3t(ji,jj,jk))  
-                  ENDIF
-                  IF( med_diag%MDET%dgsave ) THEN
-                      fdd2d(ji,jj) = fdd2d(ji,jj) +                          &
-                                     (fdd(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%AEOLIAN%dgsave ) THEN
-                      ffetop2d(ji,jj) = ffetop2d(ji,jj) +                    &
-                                        (ffetop(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%BENTHIC%dgsave ) THEN
-                      ffebot2d(ji,jj) = ffebot2d(ji,jj) +                    &
-                                        (ffebot(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%SCAVENGE%dgsave ) THEN
-                      ffescav2d(ji,jj) = ffescav2d(ji,jj) +                  &
-                                         (ffescav(ji,jj) * fse3t(ji,jj,jk))  
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%PN_JLIM%dgsave ) THEN
-                      ! fjln2d(ji,jj) = fjln2d(ji,jj) +                       &
-                      !                 (fjln(ji,jj)  * zphn(ji,jj) *         &
-                      !                  fse3t(ji,jj,jk)) 
-                      fjln2d(ji,jj) = fjln2d(ji,jj) +                        &
-                                      (fjlim_pn(ji,jj) * zphn(ji,jj) *       &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%PN_NLIM%dgsave ) THEN
-                      fnln2d(ji,jj) = fnln2d(ji,jj) +                        &
-                                      (fnln(ji,jj) * zphn(ji,jj) *           &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%PN_FELIM%dgsave ) THEN
-                      ffln2d(ji,jj) = ffln2d(ji,jj) +                        &
-                                      (ffln2(ji,jj) * zphn(ji,jj) *          &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%PD_JLIM%dgsave ) THEN
-                      ! fjld2d(ji,jj) = fjld2d(ji,jj) +                       &
-                      !                 (fjld(ji,jj)  * zphd(ji,jj) *         &
-                      !                  fse3t(ji,jj,jk)) 
-                      fjld2d(ji,jj) = fjld2d(ji,jj) +                        &
-                                      (fjlim_pd(ji,jj) * zphd(ji,jj) *       &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%PD_NLIM%dgsave ) THEN
-                      fnld2d(ji,jj) = fnld2d(ji,jj) +                        &
-                                      (fnld(ji,jj) * zphd(ji,jj) *           &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%PD_FELIM%dgsave ) THEN
-                      ffld2d(ji,jj) = ffld2d(ji,jj) +                        &
-                                      (ffld(ji,jj) * zphd(ji,jj) *           &
-                                       fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%PD_SILIM%dgsave ) THEN
-                      fsld2d2(ji,jj) = fsld2d2(ji,jj) +                      &
-                                       (fsld2(ji,jj) * zphd(ji,jj) *         &
-                                        fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%PDSILIM2%dgsave ) THEN
-                      fsld2d(ji,jj) = fsld2d(ji,jj) +                        &
-                                      (fsld(ji,jj) * zphd(ji,jj) *           &
-                                       fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! 
-                  IF( med_diag%TOTREG_N%dgsave ) THEN
-                      fregen2d(ji,jj) = fregen2d(ji,jj) + fregen(ji,jj)
-                  ENDIF
-                  IF( med_diag%TOTRG_SI%dgsave ) THEN
-                      fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi(ji,jj)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! 
-                  IF( med_diag%FASTN%dgsave ) THEN
-                      ftempn2d(ji,jj) = ftempn2d(ji,jj) +                    &
-                                        (ftempn(ji,jj)  * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%FASTSI%dgsave ) THEN
-                      ftempsi2d(ji,jj) = ftempsi2d(ji,jj) +                  &
-                                         (ftempsi(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%FASTFE%dgsave ) THEN
-                      ftempfe2d(ji,jj) = ftempfe2d(ji,jj) +                  &
-                                         (ftempfe(ji,jj) * fse3t(ji,jj,jk))  
-                  ENDIF
-                  IF( med_diag%FASTC%dgsave ) THEN
-                      ftempc2d(ji,jj) = ftempc2d(ji,jj) +                    &
-                                        (ftempc(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%FASTCA%dgsave ) THEN
-                      ftempca2d(ji,jj) = ftempca2d(ji,jj) +                  &
-                                         (ftempca(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! 
-                  IF( med_diag%REMINN%dgsave ) THEN
-                      freminn2d(ji,jj) = freminn2d(ji,jj) +                  &
-                                         (freminn(ji,jj)  * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%REMINSI%dgsave ) THEN
-                      freminsi2d(ji,jj) = freminsi2d(ji,jj) +                &
-                                          (freminsi(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%REMINFE%dgsave ) THEN
-                      freminfe2d(ji,jj) = freminfe2d(ji,jj) +                &
-                                          (freminfe(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%REMINC%dgsave ) THEN
-                      freminc2d(ji,jj) = freminc2d(ji,jj) +                  &
-                                         (freminc(ji,jj)  * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  IF( med_diag%REMINCA%dgsave ) THEN
-                      freminca2d(ji,jj) = freminca2d(ji,jj) +                &
-                                          (freminca(ji,jj) * fse3t(ji,jj,jk)) 
-                  ENDIF
-                  !!
-               ENDIF
-            ENDDO
-         ENDDO
-
-# if defined key_roam
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !!
-                  !! AXY (09/11/16): CMIP6 diagnostics
-                  IF( med_diag%FD_NIT3%dgsave ) THEN
-                     fd_nit3(ji,jj,jk) = ffastn(ji,jj)
-                  ENDIF
-                  IF( med_diag%FD_SIL3%dgsave ) THEN
-                     fd_sil3(ji,jj,jk) = ffastsi(ji,jj)
-                  ENDIF
-                  IF( med_diag%FD_CAR3%dgsave ) THEN
-                     fd_car3(ji,jj,jk) = ffastc(ji,jj)
-                  ENDIF
-                  IF( med_diag%FD_CAL3%dgsave ) THEN
-                     fd_cal3(ji,jj,jk) = ffastca(ji,jj)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         IF (jk.eq.i0100) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     IF( med_diag%RR_0100%dgsave ) THEN
-                        ffastca2d(ji,jj) =                                   &
-                           ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
-                     ENDIF
-                  ENDIF
-               ENDDO
-            ENDDO
-         ELSE IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     IF( med_diag%RR_0500%dgsave ) THEN
-                        ffastca2d(ji,jj) =                                   &
-                           ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
-                     ENDIF
-                  ENDIF
-               ENDDO
-            ENDDO
-         ELSE IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     IF( med_diag%RR_1000%dgsave ) THEN
-                        ffastca2d(ji,jj) =                                   &
-                           ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
-                     ENDIF
-                  ENDIF
-               ENDDO
-            ENDDO
-         ELSE IF (jk.eq.mbathy(ji,jj)) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     IF( med_diag%IBEN_N%dgsave ) THEN
-                        iben_n2d(ji,jj) = f_sbenin_n(ji,jj) +                &
-                                          f_fbenin_n(ji,jj)
-                     ENDIF
-                     IF( med_diag%IBEN_FE%dgsave ) THEN
-                        iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) +              &
-                                           f_fbenin_fe(ji,jj)
-                     ENDIF
-                     IF( med_diag%IBEN_C%dgsave ) THEN
-                        iben_c2d(ji,jj) = f_sbenin_c(ji,jj) +                &
-                                          f_fbenin_c(ji,jj)
-                     ENDIF
-                     IF( med_diag%IBEN_SI%dgsave ) THEN
-                        iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
-                     ENDIF
-                     IF( med_diag%IBEN_CA%dgsave ) THEN
-                        iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
-                     ENDIF
-                     IF( med_diag%OBEN_N%dgsave ) THEN
-                        oben_n2d(ji,jj) = f_benout_n(ji,jj)
-                     ENDIF
-                     IF( med_diag%OBEN_FE%dgsave ) THEN
-                        oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
-                     ENDIF
-                     IF( med_diag%OBEN_C%dgsave ) THEN
-                        oben_c2d(ji,jj) = f_benout_c(ji,jj)
-                     ENDIF
-                     IF( med_diag%OBEN_SI%dgsave ) THEN
-                        oben_si2d(ji,jj) = f_benout_si(ji,jj)
-                     ENDIF
-                     IF( med_diag%OBEN_CA%dgsave ) THEN
-                        oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
-                     ENDIF
-                     IF( med_diag%SFR_OCAL%dgsave ) THEN
-                        sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
-                     ENDIF
-                     IF( med_diag%SFR_OARG%dgsave ) THEN
-                        sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
-                     ENDIF
-                     IF( med_diag%LYSO_CA%dgsave ) THEN
-                        lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
-                     ENDIF
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-         !! end bathy-1 diags
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !!
-                  IF( med_diag%RIV_N%dgsave ) THEN
-                     rivn2d(ji,jj) = rivn2d(ji,jj) +                         &
-                                     (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%RIV_SI%dgsave ) THEN
-                     rivsi2d(ji,jj) = rivsi2d(ji,jj) +                       &
-                                      (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%RIV_C%dgsave ) THEN
-                     rivc2d(ji,jj) = rivc2d(ji,jj) +                         &
-                                     (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%RIV_ALK%dgsave ) THEN
-                     rivalk2d(ji,jj) = rivalk2d(ji,jj) +                     &
-                                       (f_riv_loc_alk(ji,jj) *               &
-                                        fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%DETC%dgsave ) THEN
-                     fslowc2d(ji,jj) = fslowc2d(ji,jj) +                     &
-                                       (fslowc(ji,jj)  * fse3t(ji,jj,jk))   
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !!
-                  IF( med_diag%PN_LLOSS%dgsave ) THEN
-                     fdpn22d(ji,jj) = fdpn22d(ji,jj) +                       &
-                                      (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%PD_LLOSS%dgsave ) THEN
-                     fdpd22d(ji,jj) = fdpd22d(ji,jj) +                       &
-                                      (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%ZI_LLOSS%dgsave ) THEN
-                     fdzmi22d(ji,jj) = fdzmi22d(ji,jj) +                     &
-                                       (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_LLOSS%dgsave ) THEN
-                     fdzme22d(ji,jj) = fdzme22d(ji,jj) +                     &
-                                       (fdzme2(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%ZI_MES_N%dgsave ) THEN
-                     zimesn2d(ji,jj) = zimesn2d(ji,jj) +                     &
-                                       (xphi * (fgmipn(ji,jj) +              &
-                                                fgmid(ji,jj)) *              &
-                                        fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZI_MES_D%dgsave ) THEN
-                     zimesd2d(ji,jj) = zimesd2d(ji,jj) +                     & 
-                                       ((1. - xbetan) * finmi(ji,jj) *       &
-                                        fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZI_MES_C%dgsave ) THEN
-                     zimesc2d(ji,jj) = zimesc2d(ji,jj) +                     &
-                                       (xphi * ((xthetapn * fgmipn(ji,jj)) + &
-                                                fgmidc(ji,jj)) *             &
-                                                fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZI_MESDC%dgsave ) THEN
-                     zimesdc2d(ji,jj) = zimesdc2d(ji,jj) +                   &
-                                        ((1. - xbetac) * ficmi(ji,jj) *      &
-                                         fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZI_EXCR%dgsave ) THEN
-                     ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +                     &
-                                       (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZI_RESP%dgsave ) THEN
-                     ziresp2d(ji,jj) = ziresp2d(ji,jj) +                     &
-                                       (fmiresp(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZI_GROW%dgsave ) THEN
-                     zigrow2d(ji,jj) = zigrow2d(ji,jj) +                     &
-                                       (fmigrow(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%ZE_MES_N%dgsave ) THEN
-                     zemesn2d(ji,jj) = zemesn2d(ji,jj) +                     &
-                                       (xphi *                               &
-                                        (fgmepn(ji,jj) + fgmepd(ji,jj) +     &
-                                         fgmezmi(ji,jj) + fgmed(ji,jj)) *    &
-                                        fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_MES_D%dgsave ) THEN
-                     zemesd2d(ji,jj) = zemesd2d(ji,jj) +                     &
-                                       ((1. - xbetan) * finme(ji,jj) *       &
-                                        fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_MES_C%dgsave ) THEN
-                     zemesc2d(ji,jj) = zemesc2d(ji,jj) +                     & 
-                                       (xphi *                               &
-                                        ((xthetapn * fgmepn(ji,jj)) +        &
-                                         (xthetapd * fgmepd(ji,jj)) +        &
-                                         (xthetazmi * fgmezmi(ji,jj)) +      &
-                                         fgmedc(ji,jj)) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_MESDC%dgsave ) THEN
-                     zemesdc2d(ji,jj) = zemesdc2d(ji,jj) +                   &
-                                        ((1. - xbetac) * ficme(ji,jj) *      &
-                                         fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_EXCR%dgsave ) THEN
-                     zeexcr2d(ji,jj) = zeexcr2d(ji,jj) +                     &
-                                       (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_RESP%dgsave ) THEN
-                     zeresp2d(ji,jj) = zeresp2d(ji,jj) +                     &
-                                       (fmeresp(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%ZE_GROW%dgsave ) THEN
-                     zegrow2d(ji,jj) = zegrow2d(ji,jj) +                     &
-                                       (fmegrow(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%MDETC%dgsave ) THEN
-                     mdetc2d(ji,jj) = mdetc2d(ji,jj) +                       &
-                                      (fddc(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%GMIDC%dgsave ) THEN
-                     gmidc2d(ji,jj) = gmidc2d(ji,jj) +                       &
-                                      (fgmidc(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-                  IF( med_diag%GMEDC%dgsave ) THEN
-                     gmedc2d(ji,jj) = gmedc2d(ji,jj) +                       &
-                                      (fgmedc(ji,jj) * fse3t(ji,jj,jk))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-# endif                   
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !!
-                  !! ** 3D diagnostics
-                  IF( med_diag%TPP3%dgsave ) THEN
-                     tpp3d(ji,jj,jk) = (fprn(ji,jj) * zphn(ji,jj)) +         &
-                                       (fprd(ji,jj) * zphd(ji,jj))
-                     !CALL iom_put( "TPP3"  , tpp3d )
-                  ENDIF
-                  IF( med_diag%TPPD3%dgsave ) THEN
-                     tppd3(ji,jj,jk) = (fprd(ji,jj) * zphd(ji,jj))
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  
-                  IF( med_diag%REMIN3N%dgsave ) THEN
-                     !! remineralisation
-                     remin3dn(ji,jj,jk) = fregen(ji,jj) +                    &
-                                          (freminn(ji,jj) * fse3t(ji,jj,jk))
-                     !CALL iom_put( "REMIN3N"  , remin3dn )
-                  ENDIF
-                  !! IF( med_diag%PH3%dgsave ) THEN
-                  !!     CALL iom_put( "PH3"  , f3_pH )
-                  !! ENDIF
-                  !! IF( med_diag%OM_CAL3%dgsave ) THEN
-                  !!     CALL iom_put( "OM_CAL3"  , f3_omcal )
-                  !! ENDIF
-		  !! 
-		  !! AXY (09/11/16): CMIP6 diagnostics
-		  IF ( med_diag%DCALC3%dgsave   ) THEN
-                     dcalc3(ji,jj,jk) = freminca(ji,jj)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-		  IF ( med_diag%FEDISS3%dgsave  ) THEN
-                     fediss3(ji,jj,jk) = ffetop(ji,jj)
-                  ENDIF
-		  IF ( med_diag%FESCAV3%dgsave  ) THEN
-                     fescav3(ji,jj,jk) = ffescav(ji,jj)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-		  IF ( med_diag%MIGRAZP3%dgsave ) THEN
-                     migrazp3(ji,jj,jk) = fgmipn(ji,jj) * xthetapn
-                  ENDIF
-		  IF ( med_diag%MIGRAZD3%dgsave ) THEN
-                     migrazd3(ji,jj,jk) = fgmidc(ji,jj)
-                  ENDIF
-		  IF ( med_diag%MEGRAZP3%dgsave ) THEN
-                     megrazp3(ji,jj,jk) = (fgmepn(ji,jj) * xthetapn) +       &
-                                          (fgmepd(ji,jj) * xthetapd)
-                  ENDIF
-		  IF ( med_diag%MEGRAZD3%dgsave ) THEN
-                     megrazd3(ji,jj,jk) = fgmedc(ji,jj)
-                  ENDIF
-		  IF ( med_diag%MEGRAZZ3%dgsave ) THEN
-                     megrazz3(ji,jj,jk) = (fgmezmi(ji,jj) * xthetazmi)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-		  IF ( med_diag%PBSI3%dgsave    ) THEN
-                     pbsi3(ji,jj,jk)    = (fprds(ji,jj) * zpds(ji,jj))
-                  ENDIF
-		  IF ( med_diag%PCAL3%dgsave    ) THEN
-                     pcal3(ji,jj,jk)    = ftempca(ji,jj)
-                  ENDIF
-		  IF ( med_diag%REMOC3%dgsave   ) THEN
-                     remoc3(ji,jj,jk)   = freminc(ji,jj)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-		  IF ( med_diag%PNLIMJ3%dgsave  ) THEN
-                     ! pnlimj3(ji,jj,jk)  = fjln(ji,jj)
-                     pnlimj3(ji,jj,jk)  = fjlim_pn(ji,jj)
-                  ENDIF
-		  IF ( med_diag%PNLIMN3%dgsave  ) THEN
-                     pnlimn3(ji,jj,jk)  = fnln(ji,jj)
-                  ENDIF
-		  IF ( med_diag%PNLIMFE3%dgsave ) THEN
-                     pnlimfe3(ji,jj,jk) = ffln2(ji,jj)
-                  ENDIF
-		  IF ( med_diag%PDLIMJ3%dgsave  ) THEN
-                     ! pdlimj3(ji,jj,jk)  = fjld(ji,jj)
-                     pdlimj3(ji,jj,jk)  = fjlim_pd(ji,jj)
-                  ENDIF
-		  IF ( med_diag%PDLIMN3%dgsave  ) THEN
-                     pdlimn3(ji,jj,jk)  = fnld(ji,jj)
-                  ENDIF
-		  IF ( med_diag%PDLIMFE3%dgsave ) THEN
-                     pdlimfe3(ji,jj,jk) = ffld(ji,jj)
-                  ENDIF
-		  IF ( med_diag%PDLIMSI3%dgsave ) THEN
-                     pdlimsi3(ji,jj,jk) = fsld2(ji,jj)
-                  ENDIF
-               ENDIF
-            ENDDO
-         ENDDO
+         !!---------------------------------------------------------------
+         !! Calculates the diagnostics used with iom_put
+         !!---------------------------------------------------------------
+         CALL bio_med_diag_iomput( jk )
 
       ELSE IF( ln_diatrc ) THEN
 
-         !!
-         !! ** Without using iom_use
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk ln_diatrc'
-         CALL flush(numout)
-#   endif
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) then
-                  !!-------------------------------------------------------
-                  !! Prepare 2D diagnostics
-                  !!-------------------------------------------------------
-                  !!
-                  !! if ((kt / 240*240).eq.kt) then
-                  !!    IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt
-                  !! endif     
-                  !! nitrogen inventory
-                  trc2d(ji,jj,1)  =  ftot_n(ji,jj)
-                  !! silicon  inventory
-                  trc2d(ji,jj,2)  =  ftot_si(ji,jj)
-                  !! iron     inventory
-                  trc2d(ji,jj,3)  =  ftot_fe(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! non-diatom production
-                  trc2d(ji,jj,4)  = trc2d(ji,jj,4)  +                        &
-                                    (fprn(ji,jj)  * zphn(ji,jj) *            &
-                                     fse3t(ji,jj,jk))
-                  !! non-diatom non-grazing losses
-                  trc2d(ji,jj,5)  = trc2d(ji,jj,5)  +                        &
-                                    (fdpn(ji,jj) * fse3t(ji,jj,jk))
-                  !! diatom production
-                  trc2d(ji,jj,6)  = trc2d(ji,jj,6)  +                        &
-                                    (fprd(ji,jj) * zphd(ji,jj) *             &
-                                     fse3t(ji,jj,jk))
-                  !! diatom non-grazing losses
-                  !! diagnostic field  8 is (ostensibly) supplied by trcsed.F
-                  trc2d(ji,jj,7)  = trc2d(ji,jj,7)  +                        &
-                                    (fdpd(ji,jj) * fse3t(ji,jj,jk))
-                  !! diatom silicon production
-                  trc2d(ji,jj,9)  = trc2d(ji,jj,9)  +                        &
-                                    (fprds(ji,jj) * zpds(ji,jj) *            &
-                                     fse3t(ji,jj,jk))
-                  !! diatom silicon dissolution
-                  trc2d(ji,jj,10) = trc2d(ji,jj,10) +                        &
-                                    (fsdiss(ji,jj)  * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! microzoo grazing on non-diatoms
-                  trc2d(ji,jj,11) = trc2d(ji,jj,11) +                        &
-                                    (fgmipn(ji,jj)  * fse3t(ji,jj,jk))
-                  !! microzoo grazing on detrital nitrogen
-                  trc2d(ji,jj,12) = trc2d(ji,jj,12) +                        &
-                                    (fgmid(ji,jj) * fse3t(ji,jj,jk))
-                  !! microzoo non-grazing losses
-                  trc2d(ji,jj,13) = trc2d(ji,jj,13) +                        &
-                                    (fdzmi(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! mesozoo  grazing on non-diatoms
-                  trc2d(ji,jj,14) = trc2d(ji,jj,14) +                        &
-                                    (fgmepn(ji,jj)  * fse3t(ji,jj,jk))
-                  !! mesozoo  grazing on diatoms
-                  trc2d(ji,jj,15) = trc2d(ji,jj,15) +                        &
-                                    (fgmepd(ji,jj)  * fse3t(ji,jj,jk))
-                  !! mesozoo  grazing on microzoo
-                  trc2d(ji,jj,16) = trc2d(ji,jj,16) +                        &
-                                    (fgmezmi(ji,jj) * fse3t(ji,jj,jk))
-                  !! mesozoo  grazing on detrital nitrogen
-                  trc2d(ji,jj,17) = trc2d(ji,jj,17) +                        &
-                                    (fgmed(ji,jj) * fse3t(ji,jj,jk))
-                  !! mesozoo  non-grazing losses
-                  trc2d(ji,jj,18) = trc2d(ji,jj,18) +                        &
-                                    (fdzme(ji,jj)   * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! diagnostic field 19 is (ostensibly) supplied by trcexp.F
-                  !! slow sinking detritus N production
-                  trc2d(ji,jj,20) = trc2d(ji,jj,20) +                        &
-                                    (fslown(ji,jj) * fse3t(ji,jj,jk))
-                  !! detrital remineralisation
-                  trc2d(ji,jj,21) = trc2d(ji,jj,21) +                        &
-                                    (fdd(ji,jj) * fse3t(ji,jj,jk))
-                  !! aeolian  iron addition
-                  trc2d(ji,jj,22) = trc2d(ji,jj,22) +                        &
-                                    (ffetop(ji,jj) * fse3t(ji,jj,jk))
-                  !! seafloor iron addition
-                  trc2d(ji,jj,23) = trc2d(ji,jj,23) +                        &
-                                    (ffebot(ji,jj) * fse3t(ji,jj,jk))
-                  !! "free" iron scavenging
-                  trc2d(ji,jj,24) = trc2d(ji,jj,24) +                        &
-                                    (ffescav(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! non-diatom J  limitation term 
-                  trc2d(ji,jj,25) = trc2d(ji,jj,25) +                        &
-                                    (fjlim_pn(ji,jj) * zphn(ji,jj) *         &
-                                     fse3t(ji,jj,jk))
-                  !! non-diatom N  limitation term 
-                  trc2d(ji,jj,26) = trc2d(ji,jj,26) +                        &
-                                    (fnln(ji,jj) * zphn(ji,jj) *             &
-                                     fse3t(ji,jj,jk))
-                  !! non-diatom Fe limitation term 
-                  trc2d(ji,jj,27) = trc2d(ji,jj,27) +                        &
-                                    (ffln2(ji,jj) * zphn(ji,jj) *            &
-                                     fse3t(ji,jj,jk))
-                  !! diatom     J  limitation term 
-                  trc2d(ji,jj,28) = trc2d(ji,jj,28) +                        &
-                                    (fjlim_pd(ji,jj) * zphd(ji,jj) *         &
-                                     fse3t(ji,jj,jk))
-                  !! diatom     N  limitation term 
-                  trc2d(ji,jj,29) = trc2d(ji,jj,29) +                        &
-                                    (fnld(ji,jj) * zphd(ji,jj) *             &
-                                     fse3t(ji,jj,jk))
-                  !! diatom     Fe limitation term 
-                  trc2d(ji,jj,30) = trc2d(ji,jj,30) +                        &
-                                    (ffld(ji,jj) * zphd(ji,jj) *             &
-                                     fse3t(ji,jj,jk))
-                  !! diatom     Si limitation term 
-                  trc2d(ji,jj,31) = trc2d(ji,jj,31) +                        &
-                                    (fsld2(ji,jj) * zphd(ji,jj) *            &
-                                     fse3t(ji,jj,jk))
-                  !! diatom     Si uptake limitation term
-                  trc2d(ji,jj,32) = trc2d(ji,jj,32) +                        &
-                                    (fsld(ji,jj) * zphd(ji,jj) *             &
-                                     fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         IF (jk.eq.i0100) THEN 
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at  100 m
-                     trc2d(ji,jj,33) = fslownflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at  200 m
-                     trc2d(ji,jj,34) = fslownflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at  500 m
-                     trc2d(ji,jj,35) = fslownflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at 1000 m
-                     trc2d(ji,jj,36) = fslownflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! non-fast N  full column regeneration
-                  trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen(ji,jj)
-                  !! non-fast Si full column regeneration
-                  trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi(ji,jj)
-                  !! non-fast N  regeneration to  100 m
-               ENDIF
-            ENDDO
-         ENDDO
-
-         IF (jk.eq.i0100) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     trc2d(ji,jj,39) = trc2d(ji,jj,37)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! non-fast N  regeneration to  200 m
-                     trc2d(ji,jj,40) = trc2d(ji,jj,37)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! non-fast N  regeneration to  500 m
-                     trc2d(ji,jj,41) = trc2d(ji,jj,37)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! non-fast N  regeneration to 1000 m
-                     trc2d(ji,jj,42) = trc2d(ji,jj,37)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! fast sinking detritus N production
-                  trc2d(ji,jj,43) = trc2d(ji,jj,43) +                        &
-                                    (ftempn(ji,jj) * fse3t(ji,jj,jk))
-                  !! fast sinking detritus Si production
-                  trc2d(ji,jj,44) = trc2d(ji,jj,44) +                        &
-                                    (ftempsi(ji,jj) * fse3t(ji,jj,jk))
-                  !! fast sinking detritus Fe production
-                  trc2d(ji,jj,45) = trc2d(ji,jj,45) +                        &
-                                    (ftempfe(ji,jj) * fse3t(ji,jj,jk))
-                  !! fast sinking detritus C production
-                  trc2d(ji,jj,46) = trc2d(ji,jj,46) +                        &
-                                    (ftempc(ji,jj)  * fse3t(ji,jj,jk))
-                  !! fast sinking detritus CaCO3 production
-                  trc2d(ji,jj,47) = trc2d(ji,jj,47) +                        &
-                                    (ftempca(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         IF (jk.eq.i0100) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus N  flux at  100 m
-                     trc2d(ji,jj,48) = ffastn(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus N  flux at  200 m
-                     trc2d(ji,jj,49) = ffastn(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus N  flux at  500 m
-                     trc2d(ji,jj,50) = ffastn(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus N  flux at 1000 m
-                     trc2d(ji,jj,51) = ffastn(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0100) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! N  regeneration to  100 m
-                     trc2d(ji,jj,52) = fregenfast(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! N  regeneration to  200 m
-                     trc2d(ji,jj,53) = fregenfast(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! N  regeneration to  500 m
-                     trc2d(ji,jj,54) = fregenfast(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! N  regeneration to 1000 m
-                     trc2d(ji,jj,55) = fregenfast(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0100) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus Si flux at  100 m
-                     trc2d(ji,jj,56) = ffastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus Si flux at  200 m
-                     trc2d(ji,jj,57) = ffastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus Si flux at  500 m
-                     trc2d(ji,jj,58) = ffastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! fast detritus Si flux at 1000 m
-                     trc2d(ji,jj,59) = ffastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0100) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! Si regeneration to  100 m
-                     trc2d(ji,jj,60) = fregenfastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! Si regeneration to  200 m
-                     trc2d(ji,jj,61) = fregenfastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! Si regeneration to  500 m
-                     trc2d(ji,jj,62) = fregenfastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! Si regeneration to 1000 m
-                     trc2d(ji,jj,63) = fregenfastsi(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! sum of fast-sinking N  fluxes
-                  trc2d(ji,jj,64) = trc2d(ji,jj,64) +                        &
-                                    (freminn(ji,jj) * fse3t(ji,jj,jk))
-                  !! sum of fast-sinking Si fluxes
-                  trc2d(ji,jj,65) = trc2d(ji,jj,65) +                        &
-                                    (freminsi(ji,jj) * fse3t(ji,jj,jk))
-                  !! sum of fast-sinking Fe fluxes
-                  trc2d(ji,jj,66) = trc2d(ji,jj,66) +                        &
-                                    (freminfe(ji,jj) * fse3t(ji,jj,jk))
-                  !! sum of fast-sinking C  fluxes
-                  trc2d(ji,jj,67) = trc2d(ji,jj,67) +                        &
-                                    (freminc(ji,jj) * fse3t(ji,jj,jk))
-                  !! sum of fast-sinking Ca fluxes
-                  trc2d(ji,jj,68) = trc2d(ji,jj,68) +                        &
-                                    (freminca(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-
-         if (jk.eq.mbathy(ji,jj)) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! N  sedimentation flux
-                     trc2d(ji,jj,69) = fsedn(ji,jj)
-                     !! Si sedimentation flux
-                     trc2d(ji,jj,70) = fsedsi(ji,jj)
-                     !! Fe sedimentation flux
-                     trc2d(ji,jj,71) = fsedfe(ji,jj)
-                     !! C  sedimentation flux
-                     trc2d(ji,jj,72) = fsedc(ji,jj)
-                     !! Ca sedimentation flux
-                     trc2d(ji,jj,73) = fsedca(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk.eq.1) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     trc2d(ji,jj,74) = qsr(ji,jj)
-                     trc2d(ji,jj,75) = xpar(ji,jj,jk)
-                     !! trc2d(ji,jj,75) = real(iters(ji,jj))
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! diagnostic fields 76 to 80 calculated below
-                  !! mixed layer non-diatom production
-                  trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj)
-                  !! mixed layer     diatom production
-                  trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-
-# if defined key_gulf_finland
-         if (jk.eq.1) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! Gulf of Finland check
-                     trc2d(ji,jj,83) = real(ibio_switch)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-# else
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! calcite CCD depth
-                  trc2d(ji,jj,83) = ocal_ccd(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-# endif
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! last model level above calcite CCD depth
-                  trc2d(ji,jj,84) = fccd(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-
-         IF (jk.eq.1) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! surface "free" iron
-                     trc2d(ji,jj,85) = xFree(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-! I'm keeping this the same as before, but it looks like it should 
-! be i0100 and not i0200 - marc 8/5/17
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! "free" iron at  100 m
-                     trc2d(ji,jj,86) = xFree(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-
-         IF (jk.eq.i0200) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! "free" iron at  200 m
-                     trc2d(ji,jj,87) = xFree(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-
-         IF (jk.eq.i0500) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! "free" iron at  500 m
-                     trc2d(ji,jj,88) = xFree(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-
-         IF (jk.eq.i1000) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! "free" iron at 1000 m
-                     trc2d(ji,jj,89) = xFree(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-
-         IF (jk.eq.1) THEN
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! AXY (27/06/12): extract "euphotic depth"
-                     trc2d(ji,jj,90) = xze(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-# if defined key_roam
-         if (jk .eq. 1) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! ROAM provisionally has access to a further 20 2D 
-                     !! diagnostics
-                     !! surface wind
-                     trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + wndm(ji,jj)
-                     !! atmospheric pCO2
-                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm(ji,jj)
-                     !! ocean pH
-                     trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph(ji,jj)
-                     !! ocean pCO2
-                     trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w(ji,jj)
-                     !! ocean H2CO3 conc.
-                     trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3(ji,jj)
-                     !! ocean HCO3 conc.
-                     trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3(ji,jj)
-                     !! ocean CO3 conc.
-                     trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3(ji,jj)
-                     !! air-sea CO2 flux
-                     trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux(ji,jj)
-                  ENDIF
-               ENDDO
-           ENDDO
-
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! ocean omega calcite 
-                     trc2d(ji,jj,99)  = trc2d(ji,jj,99)  + f_omcal(ji,jj)
-                     !! ocean omega aragonite
-                     trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj)
-                     !! ocean TDIC
-                     trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC(ji,jj)
-                     !! ocean TALK
-                     trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK(ji,jj)
-                     !! surface kw660
-                     trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660(ji,jj)
-                     !! surface pressure
-                     trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0(ji,jj)
-                     !! air-sea O2 flux
-                     trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux(ji,jj)
-                     !! ocean O2 saturation
-                     trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat(ji,jj)
-                     !! depth calcite CCD
-                     trc2d(ji,jj,107) = f2_ccd_cal(ji,jj)
-                     !! depth aragonite CCD
-                     trc2d(ji,jj,108) = f2_ccd_arg(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk .eq. mbathy(ji,jj)) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! seafloor omega calcite
-                     trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)
-                     !! seafloor omega aragonite
-                     trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk.eq.i0100) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! diagnostic fields 111 to 117 calculated below
-                     !! rain ratio at  100 m
-                     trc2d(ji,jj,118) =                                      &
-                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk.eq.i0500) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! rain ratio at  500 m
-                     trc2d(ji,jj,119) =                                      &
-                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk.eq.i1000) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! rain ratio at 1000 m
-                     trc2d(ji,jj,120) =                                      &
-                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk.eq.mbathy(ji,jj)) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! AXY (18/01/12): benthic flux diagnostics
-                     trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
-                     trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
-                     trc2d(ji,jj,123) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
-                     trc2d(ji,jj,124) = f_fbenin_si(ji,jj)
-                     trc2d(ji,jj,125) = f_fbenin_ca(ji,jj)
-                     trc2d(ji,jj,126) = f_benout_n(ji,jj)
-                     trc2d(ji,jj,127) = f_benout_fe(ji,jj)
-                     trc2d(ji,jj,128) = f_benout_c(ji,jj)
-                     trc2d(ji,jj,129) = f_benout_si(ji,jj)
-                     trc2d(ji,jj,130) = f_benout_ca(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! diagnostics fields 131 to 135 calculated below
-                  trc2d(ji,jj,136) = f_runoff(ji,jj)
-                  !! AXY (19/07/12): amended to allow for riverine 
-                  !! nutrient addition below surface
-                  trc2d(ji,jj,137) = trc2d(ji,jj,137) +                      &
-                                     (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
-                  trc2d(ji,jj,138) = trc2d(ji,jj,138) +                      &
-                                     (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
-                  trc2d(ji,jj,139) = trc2d(ji,jj,139) +                      &
-                                     (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
-                  trc2d(ji,jj,140) = trc2d(ji,jj,140) +                      &
-                                     (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
-                  !! slow sinking detritus C production
-                  trc2d(ji,jj,141) = trc2d(ji,jj,141) +                      &
-                                     (fslowc(ji,jj)  * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         if (jk.eq.i0100) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at  100 m
-                     trc2d(ji,jj,142) = fslowcflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         if (jk.eq.i0200) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at  200 m
-                     trc2d(ji,jj,143) = fslowcflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-
-         if (jk.eq.i0500) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at  500 m
-                     trc2d(ji,jj,144) = fslowcflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-
-         if (jk.eq.i1000) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! slow detritus flux at 1000 m
-                     trc2d(ji,jj,145) = fslowcflux(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! carbon     inventory
-                  trc2d(ji,jj,146)  = trc2d(ji,jj,146)  + ftot_c(ji,jj)
-                  !! alkalinity inventory
-                  trc2d(ji,jj,147)  = trc2d(ji,jj,147)  + ftot_a(ji,jj)
-                  !! oxygen     inventory
-                  trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-
-         if (jk.eq.mbathy(ji,jj)) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! community respiration
-                  trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-		  !!
-		  !! AXY (14/02/14): a Valentines Day gift to BASIN - a 
-                  !!                 shedload of new diagnostics that 
-                  !!                 they'll most likely never need!
-                  !!                 (actually, as with all such gifts, 
-                  !!                 I'm giving them some things I'd like 
-                  !!                 myself!)
-                  !! 
-                  !! ------------------------------------------------------
-                  !! linear losses
-                  !! non-diatom
-                  trc2d(ji,jj,151) = trc2d(ji,jj,151) +                      &
-                                     (fdpn2(ji,jj) * fse3t(ji,jj,jk))
-                  !! diatom
-                  trc2d(ji,jj,152) = trc2d(ji,jj,152) +                      &
-                                     (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
-                  !! microzooplankton
-                  trc2d(ji,jj,153) = trc2d(ji,jj,153) +                      &
-                                     (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
-                  !! mesozooplankton
-                  trc2d(ji,jj,154) = trc2d(ji,jj,154) +                      &
-                                     (fdzme2(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! ------------------------------------------------------
-                  !! microzooplankton grazing
-                  !! microzooplankton messy -> N
-                  trc2d(ji,jj,155) = trc2d(ji,jj,155) +                      &
-                                     (xphi * (fgmipn(ji,jj) +                &
-                                              fgmid(ji,jj)) * fse3t(ji,jj,jk))
-                  !! microzooplankton messy -> D
-                  trc2d(ji,jj,156) = trc2d(ji,jj,156) +                      &
-                                     ((1. - xbetan) * finmi(ji,jj) *         &
-                                      fse3t(ji,jj,jk))
-                  !! microzooplankton messy -> DIC
-                  trc2d(ji,jj,157) = trc2d(ji,jj,157) +                      &
-                                     (xphi * ((xthetapn * fgmipn(ji,jj)) +   &
-                                              fgmidc(ji,jj)) *               &
-                                      fse3t(ji,jj,jk))
-                  !! microzooplankton messy -> Dc
-                  trc2d(ji,jj,158) = trc2d(ji,jj,158) +                      &
-                                     ((1. - xbetac) * ficmi(ji,jj) *         &
-                                      fse3t(ji,jj,jk))
-                  !! microzooplankton excretion
-                  trc2d(ji,jj,159) = trc2d(ji,jj,159) +                      &
-                                     (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
-                  !! microzooplankton respiration
-                  trc2d(ji,jj,160) = trc2d(ji,jj,160) +                      &
-                                     (fmiresp(ji,jj) * fse3t(ji,jj,jk))
-                  !! microzooplankton growth
-                  trc2d(ji,jj,161) = trc2d(ji,jj,161) +                      &
-                                     (fmigrow(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! ------------------------------------------------------
-                  !! mesozooplankton grazing
-                  !! mesozooplankton messy -> N
-                  trc2d(ji,jj,162) = trc2d(ji,jj,162) +                      &
-                                     (xphi *                                 &
-                                      (fgmepn(ji,jj) + fgmepd(ji,jj) +       &
-                                       fgmezmi(ji,jj) + fgmed(ji,jj)) *      &
-                                      fse3t(ji,jj,jk))
-                  !! mesozooplankton messy -> D
-                  trc2d(ji,jj,163) = trc2d(ji,jj,163) +                      &
-                                     ((1. - xbetan) * finme(ji,jj) *         &
-                                      fse3t(ji,jj,jk))
-                  !! mesozooplankton messy -> DIC
-                  trc2d(ji,jj,164) = trc2d(ji,jj,164) +                      &
-                                     (xphi *                                 &
-                                      ((xthetapn * fgmepn(ji,jj)) +          &
-                                       (xthetapd * fgmepd(ji,jj)) +          &
-                                       (xthetazmi * fgmezmi(ji,jj)) +        &
-                                      fgmedc(ji,jj)) * fse3t(ji,jj,jk))
-                  !! mesozooplankton messy -> Dc
-                  trc2d(ji,jj,165) = trc2d(ji,jj,165) +                      &
-                                     ((1. - xbetac) * ficme(ji,jj) *         &
-                                      fse3t(ji,jj,jk))
-                  !! mesozooplankton excretion
-                  trc2d(ji,jj,166) = trc2d(ji,jj,166) +                      &
-                                     (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
-                  !! mesozooplankton respiration
-                  trc2d(ji,jj,167) = trc2d(ji,jj,167) +                      &
-                                     (fmeresp(ji,jj) * fse3t(ji,jj,jk))
-                  !! mesozooplankton growth
-                  trc2d(ji,jj,168) = trc2d(ji,jj,168) +                      &
-                                     (fmegrow(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! ------------------------------------------------------
-                  !! miscellaneous
-                  !! detrital C remineralisation
-                  trc2d(ji,jj,169) = trc2d(ji,jj,169) +                      &
-                                     (fddc(ji,jj) * fse3t(ji,jj,jk))
-                  !! microzoo grazing on detrital carbon
-                  trc2d(ji,jj,170) = trc2d(ji,jj,170) +                      &
-                                     (fgmidc(ji,jj)  * fse3t(ji,jj,jk))
-                  !! mesozoo  grazing on detrital carbon
-                  trc2d(ji,jj,171) = trc2d(ji,jj,171) +                      &
-                                     (fgmedc(ji,jj)  * fse3t(ji,jj,jk))
-                  !!
-               ENDIF
-            ENDDO
-         ENDDO
-
-         !! ------------------------------------------------------
-	 !!
-         !! AXY (23/10/14): extract primary production related 
-         !!                 surface fields to deal with diel 
-         !!                 cycle issues; hijacking BASIN 150m
-         !!                 diagnostics to do so (see commented 
-         !!                 out diagnostics below this section)
-         !!
-         !! extract relevant BASIN fields at 150m
-         if (jk .eq. i0150) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! Pn PP
-                     trc2d(ji,jj,172) = trc2d(ji,jj,4)
-                     !! Pn linear loss
-                     trc2d(ji,jj,173) = trc2d(ji,jj,151)
-                     !! Pn non-linear loss
-                     trc2d(ji,jj,174) = trc2d(ji,jj,5)
-                     !! Pn grazing to Zmi
-                     trc2d(ji,jj,175) = trc2d(ji,jj,11)
-                     !! Pn grazing to Zme
-                     trc2d(ji,jj,176) = trc2d(ji,jj,14)
-                     !! Pd PP
-                     trc2d(ji,jj,177) = trc2d(ji,jj,6)
-                     !! Pd linear loss
-                     trc2d(ji,jj,178) = trc2d(ji,jj,152)
-                     !! Pd non-linear loss
-                     trc2d(ji,jj,179) = trc2d(ji,jj,7)
-                     !! Pd grazing to Zme
-                     trc2d(ji,jj,180) = trc2d(ji,jj,15)
-                     !! Zmi grazing on D
-                     trc2d(ji,jj,181) = trc2d(ji,jj,12)
-                     !! Zmi grazing on Dc
-                     trc2d(ji,jj,182) = trc2d(ji,jj,170)
-                     !! Zmi messy feeding loss to N
-                     trc2d(ji,jj,183) = trc2d(ji,jj,155)
-                     !! Zmi messy feeding loss to D
-                     trc2d(ji,jj,184) = trc2d(ji,jj,156)
-                     !! Zmi messy feeding loss to DIC
-                     trc2d(ji,jj,185) = trc2d(ji,jj,157)
-                     !! Zmi messy feeding loss to Dc
-                     trc2d(ji,jj,186) = trc2d(ji,jj,158)
-                     !! Zmi excretion
-                     trc2d(ji,jj,187) = trc2d(ji,jj,159)
-                     !! Zmi respiration
-                     trc2d(ji,jj,188) = trc2d(ji,jj,160)
-                     !! Zmi growth
-                     trc2d(ji,jj,189) = trc2d(ji,jj,161)
-                     !! Zmi linear loss
-                     trc2d(ji,jj,190) = trc2d(ji,jj,153)
-                     !! Zmi non-linear loss
-                     trc2d(ji,jj,191) = trc2d(ji,jj,13)
-                     !! Zmi grazing to Zme
-                     trc2d(ji,jj,192) = trc2d(ji,jj,16)
-                     !! Zme grazing on D
-                     trc2d(ji,jj,193) = trc2d(ji,jj,17)
-                     !! Zme grazing on Dc
-                     trc2d(ji,jj,194) = trc2d(ji,jj,171)
-                     !! Zme messy feeding loss to N
-                     trc2d(ji,jj,195) = trc2d(ji,jj,162)
-                     !! Zme messy feeding loss to D
-                     trc2d(ji,jj,196) = trc2d(ji,jj,163)
-                     !! Zme messy feeding loss to DIC
-                     trc2d(ji,jj,197) = trc2d(ji,jj,164)
-                     !! Zme messy feeding loss to Dc
-                     trc2d(ji,jj,198) = trc2d(ji,jj,165)
-                     !! Zme excretion
-                     trc2d(ji,jj,199) = trc2d(ji,jj,166)
-                     !! Zme respiration
-                     trc2d(ji,jj,200) = trc2d(ji,jj,167)
-                     !! Zme growth
-                     trc2d(ji,jj,201) = trc2d(ji,jj,168)
-                     !! Zme linear loss
-                     trc2d(ji,jj,202) = trc2d(ji,jj,154)
-                     !! Zme non-linear loss
-                     trc2d(ji,jj,203) = trc2d(ji,jj,18)
-                     !! Slow detritus production, N
-                     trc2d(ji,jj,204) = trc2d(ji,jj,20)
-                     !! Slow detritus remineralisation, N
-                     trc2d(ji,jj,205) = trc2d(ji,jj,21)
-                     !! Slow detritus production, C
-                     trc2d(ji,jj,206) = trc2d(ji,jj,141)
-                     !! Slow detritus remineralisation, C
-                     trc2d(ji,jj,207) = trc2d(ji,jj,169)
-                     !! Fast detritus production, N
-                     trc2d(ji,jj,208) = trc2d(ji,jj,43)
-                     !! Fast detritus remineralisation, N
-                     trc2d(ji,jj,209) = trc2d(ji,jj,21)
-                     !! Fast detritus production, C
-                     trc2d(ji,jj,210) = trc2d(ji,jj,64)
-                     !! Fast detritus remineralisation, C
-                     trc2d(ji,jj,211) = trc2d(ji,jj,67)
-                     !! Community respiration
-                     trc2d(ji,jj,212) = trc2d(ji,jj,150)
-                     !! Slow detritus N flux at 150 m
-                     trc2d(ji,jj,213) = fslownflux(ji,jj)
-                     !! Slow detritus C flux at 150 m
-                     trc2d(ji,jj,214) = fslowcflux(ji,jj)
-                     !! Fast detritus N flux at 150 m
-                     trc2d(ji,jj,215) = ffastn(ji,jj)
-                     !! Fast detritus C flux at 150 m
-                     trc2d(ji,jj,216) = ffastc(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-
-         !! 
-         !! Jpalm (11-08-2014)
-         !! Add UKESM1 diagnoatics 
-         !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-         if ((jk .eq. 1) .and.( jdms.eq.1)) then
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  IF (tmask(ji,jj,jk) == 1) THEN
-                     !! DMS surface concentration 
-                     trc2d(ji,jj,221) = dms_surf(ji,jj)
-                     !! AXY (13/03/15): add in other DMS estimates
-                     !! DMS surface concentration 
-                     trc2d(ji,jj,222) = dms_andr(ji,jj)
-                     !! DMS surface concentration 
-                     trc2d(ji,jj,223) = dms_simo(ji,jj)
-                     !! DMS surface concentration 
-                     trc2d(ji,jj,224) = dms_aran(ji,jj)
-                     !! DMS surface concentration 
-                     trc2d(ji,jj,225) = dms_hall(ji,jj)
-                  ENDIF
-               ENDDO
-            ENDDO
-         endif
-# endif
-
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! other possible future diagnostics include:
-                  !!   - integrated tracer values (esp. biological)
-                  !!   - mixed layer tracer values
-                  !!   - sub-surface chlorophyll maxima (plus depth)
-                  !!   - different mixed layer depth criteria (T, sigma, 
-                  !!     var. sigma)
-                  !!-------------------------------------------------------
-                  !! Prepare 3D diagnostics
-                  !!-------------------------------------------------------
-                  !!
-                  !! primary production  
-                  trc3d(ji,jj,jk,1)  = ((fprn(ji,jj) + fprd(ji,jj)) *        &
-                                        zphn(ji,jj))
-                  !! detrital flux
-                  trc3d(ji,jj,jk,2)  = fslownflux(ji,jj) + ffastn(ji,jj)
-                  !! remineralisation
-                  trc3d(ji,jj,jk,3)  = fregen(ji,jj) +                       &
-                                       (freminn(ji,jj) * fse3t(ji,jj,jk))
-               ENDIF
-            ENDDO
-         ENDDO
-# if defined key_roam
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! pH
-                  trc3d(ji,jj,jk,4)  = f3_pH(ji,jj,jk)
-                  !! omega calcite
-                  trc3d(ji,jj,jk,5)  = f3_omcal(ji,jj,jk)
-               ENDIF
-            ENDDO
-         ENDDO
-# else
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  !! fast Si flux
-                  trc3d(ji,jj,jk,4)  = ffastsi(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-# endif
-
-      ENDIF   ! end of ln_diatrc option
+         !!---------------------------------------------------------------
+         !! The diagnostics without using iom_use
+         !!---------------------------------------------------------------
+         CALL bio_med_diag_trc( jk )
+
+      ENDIF
 
    END SUBROUTINE bio_medusa_diag
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 8023)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8025)
@@ -633,5 +633,5 @@
          !! Diagnostics
          !!---------------------------------------------------------------
-         CALL bio_medusa_diag( kt, jk )
+         CALL bio_medusa_diag( jk )
 
          !!-------------------------------------------------------
