Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/AGE/trcsms_age.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/AGE/trcsms_age.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/AGE/trcsms_age.F90	(revision 8442)
@@ -57,7 +57,9 @@
       IF( nn_timing == 1 )  CALL timing_start('trc_sms_age')
       !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' trc_sms_age:  AGE model'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      IF( kt == nittrc000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' trc_sms_age:  AGE model'
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      ENDIF
 
       IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrage )
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/C14b/trcnam_c14b.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/C14b/trcnam_c14b.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/C14b/trcnam_c14b.F90	(revision 8442)
@@ -49,9 +49,6 @@
       ! definition of additional diagnostic as a structure
       INTEGER :: jl, jn
-      TYPE(DIAG), DIMENSION(jp_c14b_2d) :: c14dia2d
-      TYPE(DIAG), DIMENSION(jp_c14b_3d) :: c14dia3d
       !!
       NAMELIST/namc14date/ ndate_beg_b, nyear_res_b
-      NAMELIST/namc14dia/  c14dia2d, c14dia3d     ! additional diagnostics
       !!-------------------------------------------------------------------
       !                             ! Open namelist file
@@ -77,50 +74,4 @@
       IF(lwp) WRITE(numout,*) '    initial year (aa)                  nyear_beg_b = ', nyear_beg_b
       !
-      IF( .NOT.lk_iomput .AND. ln_diatrc ) THEN
-         !
-         ! Namelist namc14dia
-         ! -------------------
-         REWIND( numnatb_ref )              ! Namelist namc14dia in reference namelist : c14b diagnostics
-         READ  ( numnatb_ref, namc14dia, IOSTAT = ios, ERR = 903)
-903      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc14dia in reference namelist', lwp )
-
-         REWIND( numnatb_cfg )              ! Namelist namc14dia in configuration namelist : c14b diagnostics
-         READ  ( numnatb_cfg, namc14dia, IOSTAT = ios, ERR = 904 )
-904      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc14dia in configuration namelist', lwp )
-         IF(lwm) WRITE ( numonb, namc14dia )
-
-         DO jl = 1, jp_c14b_2d
-            jn = jp_c14b0_2d + jl - 1
-            ctrc2d(jn) = c14dia2d(jl)%sname
-            ctrc2l(jn) = c14dia2d(jl)%lname
-            ctrc2u(jn) = c14dia2d(jl)%units
-         END DO
-
-         DO jl = 1, jp_c14b_3d
-            jn = jp_c14b0_3d + jl - 1
-            ctrc3d(jn) = c14dia3d(jl)%sname
-            ctrc3l(jn) = c14dia3d(jl)%lname
-            ctrc3u(jn) = c14dia3d(jl)%units
-         END DO
-
-         IF(lwp) THEN                   ! control print
-            WRITE(numout,*)
-            WRITE(numout,*) ' Namelist : natadd'
-            DO jl = 1, jp_c14b_3d
-               jn = jp_c14b0_3d + jl - 1
-               WRITE(numout,*) '  3d diag nb : ', jn, '    short name : ', ctrc3d(jn), &
-                 &             '  long name  : ', ctrc3l(jn), '   unit : ', ctrc3u(jn)
-            END DO
-            WRITE(numout,*) ' '
-
-            DO jl = 1, jp_c14b_2d
-               jn = jp_c14b0_2d + jl - 1
-               WRITE(numout,*) '  2d diag nb : ', jn, '    short name : ', ctrc2d(jn), &
-                 &             '  long name  : ', ctrc2l(jn), '   unit : ', ctrc2u(jn)
-            END DO
-            WRITE(numout,*) ' '
-         ENDIF
-         !
-      ENDIF
 
    IF(lwm) CALL FLUSH ( numonb )     ! flush output namelist C14b
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcnam_cfc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcnam_cfc.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcnam_cfc.F90	(revision 8442)
@@ -47,8 +47,6 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       INTEGER :: jl, jn
-      TYPE(DIAG), DIMENSION(jp_cfc_2d) :: cfcdia2d
       !!
       NAMELIST/namcfcdate/ ndate_beg, nyear_res, simu_type 
-      NAMELIST/namcfcdia/  cfcdia2d     ! additional diagnostics
       !!----------------------------------------------------------------------
       !                             ! Open namelist files
@@ -82,37 +80,4 @@
       !
 
-      IF( .NOT.lk_iomput .AND. ln_diatrc ) THEN
-         !
-         ! Namelist namcfcdia
-         ! -------------------
-         REWIND( numnatc_ref )              ! Namelist namcfcdia in reference namelist : CFC diagnostics
-         READ  ( numnatc_ref, namcfcdia, IOSTAT = ios, ERR = 903)
-903      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfcdia in reference namelist', lwp )
-
-         REWIND( numnatc_cfg )              ! Namelist namcfcdia in configuration namelist : CFC diagnostics
-         READ  ( numnatc_cfg, namcfcdia, IOSTAT = ios, ERR = 904 )
-904      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfcdia in configuration namelist', lwp )
-         IF(lwm) WRITE ( numonc, namcfcdia )
-
-         DO jl = 1, jp_cfc_2d
-            jn = jp_cfc0_2d + jl - 1
-            ctrc2d(jn) = TRIM( cfcdia2d(jl)%sname )
-            ctrc2l(jn) = TRIM( cfcdia2d(jl)%lname )
-            ctrc2u(jn) = TRIM( cfcdia2d(jl)%units )
-         END DO
-
-         IF(lwp) THEN                   ! control print
-            WRITE(numout,*)
-            WRITE(numout,*) ' Namelist : natadd'
-            DO jl = 1, jp_cfc_2d
-               jn = jp_cfc0_2d + jl - 1
-               WRITE(numout,*) '  2d diag nb : ', jn, '    short name : ', ctrc2d(jn), &
-                 &             '  long name  : ', ctrc2l(jn), '   unit : ', ctrc2u(jn)
-            END DO
-            WRITE(numout,*) ' '
-         ENDIF
-         !
-      ENDIF
-
    IF(lwm) CALL FLUSH ( numonc )     ! flush output namelist CFC
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 8442)
@@ -257,21 +257,10 @@
       !ENDIF                                            
       !
-      IF( lk_iomput ) THEN
-         IF  (iom_use("qtrCFC11"))  CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
-         IF  (iom_use("qintCFC11")) CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
-         IF  (iom_use("qtrCFC12"))  CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) )
-         IF  (iom_use("qintCFC12")) CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) )
-         IF  (iom_use("qtrSF6"))    CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) )
-         IF  (iom_use("qintSF6"))   CALL iom_put( "qintSF6"   , qint_cfc(:,:,3) )
-      ELSE
-         IF( ln_diatrc ) THEN
-            trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1)
-            trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1)
-            trc2d(:,:,jp_cfc0_2d + 2) = qtr_cfc (:,:,2)
-            trc2d(:,:,jp_cfc0_2d + 3) = qint_cfc(:,:,2)
-            trc2d(:,:,jp_cfc0_2d + 4) = qtr_cfc (:,:,3)
-            trc2d(:,:,jp_cfc0_2d + 5) = qint_cfc(:,:,3)
-         END IF
-      END IF
+      IF  (iom_use("qtrCFC11"))  CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
+      IF  (iom_use("qintCFC11")) CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
+      IF  (iom_use("qtrCFC12"))  CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) )
+      IF  (iom_use("qintCFC12")) CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) )
+      IF  (iom_use("qtrSF6"))    CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) )
+      IF  (iom_use("qintSF6"))   CALL iom_put( "qintSF6"   , qint_cfc(:,:,3) )
       !
       IF( l_trdtrc ) THEN
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/IDTRA/trcsms_idtra.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/IDTRA/trcsms_idtra.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/IDTRA/trcsms_idtra.F90	(revision 8442)
@@ -165,15 +165,7 @@
       !ENDIF
       !
-      IF( lk_iomput ) THEN
          CALL iom_put( "qtrIDTRA"  , qtr_idtra (:,:,1) )
          CALL iom_put( "qintIDTRA" , qint_idtra(:,:,1) )
          CALL iom_put( "invIDTRA" , inv_idtra(:,:,1) )
-      ELSE
-         IF( ln_diatrc ) THEN
-            trc2d(:,:,jp_idtra0_2d    ) = qtr_idtra (:,:,1)
-            trc2d(:,:,jp_idtra0_2d + 1) = qint_idtra(:,:,1)
-            trc2d(:,:,jp_idtra0_2d + 2) = inv_idtra(:,:,1)
-         END IF
-      END IF
       !
 # if defined key_debug_medusa
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8442)
@@ -124,4 +124,10 @@
       !! Air-sea gas exchange
       !!-----------------------------------------------------------
+
+#   if defined key_debug_medusa
+               IF (lwp) write (numout,*)                     & 
+               'air-sea: gas_transfer kt = ', kt
+               CALL flush(numout)
+#   endif
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
@@ -145,18 +151,15 @@
                !!                 Wanninkhof (2014), option 7
                !!
+               CALL gas_transfer( wndm(ji,jj), 1, 7,         &  ! inputs
+                                  f_kw660(ji,jj) )              ! outputs
+            ENDIF
+         ENDDO
+      ENDDO
+
 #   if defined key_debug_medusa
-               IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer'
+               IF (lwp) write (numout,*)                     &
+               'air-sea: carb-chem kt = ', kt
                CALL flush(numout)
 #   endif
-               CALL gas_transfer( wndm(ji,jj), 1, 7,         &  ! inputs
-                                  f_kw660(ji,jj) )              ! outputs
-#   if defined key_debug_medusa
-               IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer'
-               CALL flush(numout)
-#   endif
-            ENDIF
-         ENDDO
-      ENDDO
-
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
@@ -228,5 +231,5 @@
                !!                 failure position can be determined
                if (iters .eq. 25) then
-                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', &
+                  IF(lwp) WRITE(numout,*) 'air-sea: ITERS WARNING, ',       &
                      iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
                endif
@@ -335,7 +338,6 @@
                                          hmld(ji,jj),qsr(ji,jj),              &
                                          zdin(ji,jj), dms_nlim(ji,jj),        &
-                                         dms_andr(ji,jj),dms_simo(ji,jj),     &
-                                         dms_aran(ji,jj),dms_hall(ji,jj),     &
-                                         dms_andm(ji,jj))
+                                         dms_andr,dms_simo,dms_aran,dms_hall, & 
+                                         dms_andm)
                   else
                      !! use diel-average inputs
@@ -346,53 +348,49 @@
                                          zn_dms_mld(ji,jj),zn_dms_qsr(ji,jj), &
                                          zn_dms_din(ji,jj),dms_nlim(ji,jj),   &
-                                         dms_andr(ji,jj),dms_simo(ji,jj),     &
-                                         dms_aran(ji,jj),dms_hall(ji,jj),     &
-                                         dms_andm(ji,jj))
+                                         dms_andr,dms_simo,dms_aran,dms_hall, & 
+                                         dms_andm)
                   endif
                   !!
                   !! assign correct output to variable passed to atmosphere
                   if (jdms_model .eq. 1) then
-                     dms_surf(ji,jj) = dms_andr(ji,jj)
+                     dms_surf = dms_andr
                   elseif (jdms_model .eq. 2) then
-                     dms_surf(ji,jj) = dms_simo(ji,jj)
+                     dms_surf = dms_simo
                   elseif (jdms_model .eq. 3) then
-                     dms_surf(ji,jj) = dms_aran(ji,jj)
+                     dms_surf = dms_aran
                   elseif (jdms_model .eq. 4) then
-                     dms_surf(ji,jj) = dms_hall(ji,jj)
+                     dms_surf = dms_hall
                   elseif (jdms_model .eq. 5) then
-                     dms_surf(ji,jj) = dms_andm(ji,jj)
+                     dms_surf = dms_andm
                   endif
                   !!
                   !! 2D diag through iom_use
-                  IF( lk_iomput ) THEN
-                     IF( med_diag%DMS_SURF%dgsave ) THEN
-                        dms_surf2d(ji,jj) = dms_surf(ji,jj)
-                     ENDIF
-                     IF( med_diag%DMS_ANDR%dgsave ) THEN
-                        dms_andr2d(ji,jj) = dms_andr(ji,jj)
-                     ENDIF
-                     IF( med_diag%DMS_SIMO%dgsave ) THEN
-                        dms_simo2d(ji,jj) = dms_simo(ji,jj)
-                     ENDIF
-                     IF( med_diag%DMS_ARAN%dgsave ) THEN
-                        dms_aran2d(ji,jj) = dms_aran(ji,jj)
-                     ENDIF
-                     IF( med_diag%DMS_HALL%dgsave ) THEN
-                        dms_hall2d(ji,jj) = dms_hall(ji,jj)
-                     ENDIF 
-                     IF( med_diag%DMS_ANDM%dgsave ) THEN 
-                        dms_andm2d(ji,jj) = dms_andm(ji,jj) 
-                     ENDIF
-#   if defined key_debug_medusa
-                     IF (lwp) write (numout,*)                                &
-                        'trc_bio_medusa: finish calculating dms'
-                     CALL flush(numout)
-#   endif 
-                  ENDIF      !! End iom
+                  IF( med_diag%DMS_SURF%dgsave ) THEN
+                     dms_surf2d(ji,jj) = dms_surf
+                  ENDIF
+                  IF( med_diag%DMS_ANDR%dgsave ) THEN
+                     dms_andr2d(ji,jj) = dms_andr
+                  ENDIF
+                  IF( med_diag%DMS_SIMO%dgsave ) THEN
+                     dms_simo2d(ji,jj) = dms_simo
+                  ENDIF
+                  IF( med_diag%DMS_ARAN%dgsave ) THEN
+                     dms_aran2d(ji,jj) = dms_aran
+                  ENDIF
+                  IF( med_diag%DMS_HALL%dgsave ) THEN
+                     dms_hall2d(ji,jj) = dms_hall
+                  ENDIF 
+                  IF( med_diag%DMS_ANDM%dgsave ) THEN 
+                     dms_andm2d(ji,jj) = dms_andm
+                  ENDIF
                ENDIF
             ENDDO
          ENDDO
+#   if defined key_debug_medusa
+         IF (lwp) write (numout,*)                                &
+            'air-sea: finish calculating dms kt = ',kt
+            CALL flush(numout)
+#   endif 
       ENDIF                  !! End IF (jdms == 1)
-
 
       !!
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 8442)
@@ -61,5 +61,5 @@
 #   if defined key_debug_medusa
                IF (lwp) write (numout,*)                                     &
-                  'trc_bio_medusa: diag in ij-jj-jk loop'
+                  'bio_med_diag_iomput: in ij-jj loop jk = ', jk
                CALL flush(numout)
 #   endif
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90	(revision 8441)
+++ 	(revision )
@@ -1,1203 +1,0 @@
-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
-
-
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            IF (tmask(ji,jj,jk) == 1) THEN
-               if (jk.eq.mbathy(ji,jj)) 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
-            ENDIF
-         ENDDO
-      ENDDO
-
-      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 for reproducibility, 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
-
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            IF (tmask(ji,jj,jk) == 1) THEN
-               if (jk .eq. mbathy(ji,jj)) 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
-            ENDIF
-         ENDDO
-      ENDDO
-
-      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
-
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            IF (tmask(ji,jj,jk) == 1) THEN
-               if (jk.eq.mbathy(ji,jj)) 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
-            ENDIF
-         ENDDO
-      ENDDO
-
-      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)
-
-               if (jk.eq.mbathy(ji,jj)) then
-                  trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
-               endif
-
-               !! 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_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8442)
@@ -31,5 +31,4 @@
       !!-------------------------------------------------------------------
       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,                  &
@@ -42,5 +41,5 @@
       USE sms_medusa,               ONLY: xrfn, xthetapd, xthetapn,      &
                                           xthetazme, xthetazmi
-      USE trc,                      ONLY: ln_diatrc, med_diag 
+      USE trc,                      ONLY: med_diag 
 # if defined key_roam
       USE trcoxy_medusa,            ONLY: oxy_sato
@@ -193,19 +192,8 @@
 # endif
 
-      IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
-
-         !!---------------------------------------------------------------
-         !! Calculates the diagnostics used with iom_put
-         !!---------------------------------------------------------------
-         CALL bio_med_diag_iomput( jk )
-
-      ELSE IF( ln_diatrc ) THEN
-
-         !!---------------------------------------------------------------
-         !! The diagnostics without using iom_use
-         !!---------------------------------------------------------------
-         CALL bio_med_diag_trc( jk )
-
-      ENDIF
+      !!---------------------------------------------------------------
+      !! Calculates the diagnostics used with iom_put
+      !!---------------------------------------------------------------
+      CALL bio_med_diag_iomput( jk )
 
    END SUBROUTINE bio_medusa_diag
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 8442)
@@ -60,9 +60,10 @@
       !!
       !!-----------------------------------------
-      IF (jk.eq.1) THEN
 #   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 1'
+         IF (lwp) write (numout,*) 'bio_medusa_diag_slice: start jk = ', jk
          CALL flush(numout)
 #   endif
+      !!
+      IF (jk.eq.1) THEN
          !! JPALM -- 02-06-2017 --
          !! add Chl surf coupling
@@ -245,8 +246,4 @@
 # endif                     
       ELSE IF (jk.eq.i0100) THEN 
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 100'
-         CALL flush(numout)
-#   endif
          IF( med_diag%SDT__100%dgsave ) THEN
             zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
@@ -293,15 +290,6 @@
             CALL iom_put( "epSI100"   , ffastsi )
          ENDIF		     
-      ELSE IF (jk.eq.i0150) THEN
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 150'
-         CALL flush(numout)
-#   endif
 # endif                     
       ELSE IF (jk.eq.i0200) THEN
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 200'
-         CALL flush(numout)
-#   endif
          IF( med_diag%SDT__200%dgsave ) THEN
             zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
@@ -333,8 +321,4 @@
 # endif                     
       ELSE IF (jk.eq.i0500) THEN
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 500'
-         CALL flush(numout)
-#   endif
          IF( med_diag%SDT__500%dgsave ) THEN
             zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
@@ -369,8 +353,4 @@
 # endif                      
       ELSE IF (jk.eq.i1000) THEN
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: diag jk = 1000'
-         CALL flush(numout)
-#   endif
          IF( med_diag%SDT_1000%dgsave ) THEN
             zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
@@ -420,4 +400,8 @@
       ENDIF		     
 # endif		     
+#   if defined key_debug_medusa
+         IF (lwp) write (numout,*) 'bio_medusa_diag_slice: end jk = ', jk
+         CALL flush(numout)
+#   endif
 
    END SUBROUTINE bio_medusa_diag_slice
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 8442)
@@ -51,7 +51,5 @@
                                    zn_sed_c, zn_sed_ca, zn_sed_fe,      &
                                    zn_sed_n, zn_sed_si
-      USE trc,               ONLY: ctrc2d, ctrc3d, ln_diatrc,           &
-                                   med_diag, nittrc000,                 &
-                                   trbio, trc2d, trc3d
+      USE trc,               ONLY: med_diag, nittrc000 
       USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
  
@@ -100,23 +98,4 @@
          zn_sed_ca(:,:) = za_sed_ca(:,:)
       endif
-      IF( ln_diatrc ) THEN
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               trc2d(ji,jj,131) = za_sed_n(ji,jj)
-               trc2d(ji,jj,132) = za_sed_fe(ji,jj)
-               trc2d(ji,jj,133) = za_sed_c(ji,jj)
-               trc2d(ji,jj,134) = za_sed_si(ji,jj)
-               trc2d(ji,jj,135) = za_sed_ca(ji,jj)
-            ENDDO
-         ENDDO
-         !! AXY (07/07/15): temporary hijacking
-# if defined key_roam
-  !!       trc2d(:,:,126) = zn_dms_chn(:,:)
-  !!       trc2d(:,:,127) = zn_dms_chd(:,:)
-  !!       trc2d(:,:,128) = zn_dms_mld(:,:)
-  !!       trc2d(:,:,129) = zn_dms_qsr(:,:)
-  !!       trc2d(:,:,130) = zn_dms_din(:,:)
-# endif
-      ENDIF 
       !!
       if (ibenthic.eq.2) then
@@ -238,5 +217,5 @@
          ENDIF
       endif
-     
+
 #  if defined key_debug_medusa
          !! AXY (12/07/17)
@@ -279,5 +258,5 @@
                ENDIF
             ENDDO
-         ENDDO
+         ENDDO   
          !! silicon
          DO jj = 2,jpjm1
@@ -292,5 +271,5 @@
                ENDIF
             ENDDO
-         ENDDO
+         ENDDO   
          !! carbon
          DO jj = 2,jpjm1
@@ -306,5 +285,5 @@
                ENDIF
             ENDDO
-         ENDDO
+         ENDDO   
          !! alkalinity
          DO jj = 2,jpjm1
@@ -319,176 +298,7 @@
                ENDIF
             ENDDO
-         ENDDO
+         ENDDO   
 #  endif
- 
-      IF( ln_diatrc ) THEN
-         !!-----------------------------------------------------------------
-         !! Output several accumulated diagnostics
-         !!   - biomass-average phytoplankton limitation terms
-         !!   - integrated tendency terms
-         !!-----------------------------------------------------------------
-         !! 
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               !! non-diatom phytoplankton limitations
-               trc2d(ji,jj,25)  = trc2d(ji,jj,25) / MAX(ftot_pn(ji,jj), rsmall)
-               trc2d(ji,jj,26)  = trc2d(ji,jj,26) / MAX(ftot_pn(ji,jj), rsmall)
-               trc2d(ji,jj,27)  = trc2d(ji,jj,27) / MAX(ftot_pn(ji,jj), rsmall)
-               !! diatom phytoplankton limitations
-               trc2d(ji,jj,28)  = trc2d(ji,jj,28) / MAX(ftot_pd(ji,jj), rsmall)
-               trc2d(ji,jj,29)  = trc2d(ji,jj,29) / MAX(ftot_pd(ji,jj), rsmall)
-               trc2d(ji,jj,30)  = trc2d(ji,jj,30) / MAX(ftot_pd(ji,jj), rsmall)
-               trc2d(ji,jj,31)  = trc2d(ji,jj,31) / MAX(ftot_pd(ji,jj), rsmall)
-               trc2d(ji,jj,32)  = trc2d(ji,jj,32) / MAX(ftot_pd(ji,jj), rsmall)
-               !! tendency terms
-               trc2d(ji,jj,76)  = fflx_n(ji,jj)
-               trc2d(ji,jj,77)  = fflx_si(ji,jj)
-               trc2d(ji,jj,78)  = fflx_fe(ji,jj)
-               !! Integrated biomass.
-               !! integrated non-diatom phytoplankton
-               trc2d(ji,jj,79)  = ftot_pn(ji,jj)
-               !! integrated diatom phytoplankton
-               trc2d(ji,jj,80)  = ftot_pd(ji,jj)
-               !! Integrated microzooplankton
-               trc2d(ji,jj,217) = ftot_zmi(ji,jj)
-               !! Integrated mesozooplankton
-               trc2d(ji,jj,218) = ftot_zme(ji,jj)
-               !! Integrated slow detritus, nitrogen
-               trc2d(ji,jj,219) = ftot_det(ji,jj)
-               !! Integrated slow detritus, carbon
-               trc2d(ji,jj,220) = ftot_dtc(ji,jj)
-# if defined key_roam
-               !! The balance of nitrogen production/consumption.
-               !! integrated nitrogen production
-               trc2d(ji,jj,111) = fnit_prod(ji,jj)
-               !! integrated nitrogen consumption
-               trc2d(ji,jj,112) = fnit_cons(ji,jj)
-               !! The balance of carbon production/consumption.
-               !! integrated carbon production
-               trc2d(ji,jj,113) = fcar_prod(ji,jj)
-               !! integrated carbon consumption
-               trc2d(ji,jj,114) = fcar_cons(ji,jj)
-               !! The balance of oxygen production/consumption.
-               !! integrated oxygen production
-               trc2d(ji,jj,115) = foxy_prod(ji,jj)
-               !! integrated oxygen consumption
-               trc2d(ji,jj,116) = foxy_cons(ji,jj)
-               !! integrated unrealised oxygen consumption
-               trc2d(ji,jj,117) = foxy_anox(ji,jj)
-# endif
-            ENDDO
-         ENDDO
-         
-# if defined key_roam
-#  if defined key_axy_nancheck
-         !!---------------------------------------------------------------
-         !! Check for NaNs in diagnostic outputs
-         !!---------------------------------------------------------------
-         !! 
-         !! 2D diagnostics
-         DO jn = 1,150
-            fq0 = SUM(trc2d(:,:,jn))
-            !! AXY (30/01/14): "isnan" problem on HECTOR
-            !! if (fq0 /= fq0 ) then
-            if ( ieee_is_nan( fq0 ) ) then
-               !! there's a NaN here
-               if (lwp) write(numout,*)                                     &
-                  'NAN detected in 2D diagnostic field', jn, 'at time',     &
-                   kt, 'at position:'
-               DO jj = 1,jpj
-                  DO ji = 1,jpi
-                     if ( ieee_is_nan( trc2d(ji,jj,jn) ) ) then
-                        if (lwp) write (numout,'(a,3i6)') 'NAN-CHECK',      &
-                                 ji, jj, jn
-                     endif
-                  ENDDO
-               ENDDO
-	       CALL ctl_stop( 'trcbio_medusa, NAN in 2D diagnostic field' )
-            endif
-         ENDDO
-         !!
-         !! 3D diagnostics
-         DO jn = 1,5
-            fq0 = SUM(trc3d(:,:,:,jn))
-            !! AXY (30/01/14): "isnan" problem on HECTOR
-            !! if (fq0 /= fq0 ) then
-            if ( ieee_is_nan( fq0 ) ) then
-               !! there's a NaN here
-               if (lwp) write(numout,*)                                     &
-                  'NAN detected in 3D diagnostic field', jn, 'at time',     &
-                  kt, 'at position:'
-               DO jk = 1,jpk
-                  DO jj = 1,jpj
-                     DO ji = 1,jpi
-                        if ( ieee_is_nan( trc3d(ji,jj,jk,jn) ) ) then
-                           if (lwp) write (numout,'(a,4i6)') 'NAN-CHECK',   &
-                                    ji, jj, jk, jn
-                        endif
-                     ENDDO
-                  ENDDO
-               ENDDO
-	       CALL ctl_stop( 'trcbio_medusa, NAN in 3D diagnostic field' )
-            endif
-         ENDDO
-	 CALL flush(numout)
-#  endif
-# endif
-
-         !!----------------------------------------------------------------
-         !! Don't know what this does; belongs to someone else ...
-         !!----------------------------------------------------------------
-         !! 
-         !! Lateral boundary conditions on trc2d
-         DO jn=1,jp_medusa_2d
-             CALL lbc_lnk(trc2d(:,:,jn),'T',1. )
-         ENDDO 
-
-         !! Lateral boundary conditions on trc3d
-         DO jn=1,jp_medusa_3d
-             CALL lbc_lnk(trc3d(:,:,1,jn),'T',1. )
-         ENDDO 
-
-
-# if defined key_axy_nodiag
-         !!----------------------------------------------------------------
-         !! Blank diagnostics as a NaN-trap
-         !!----------------------------------------------------------------
-         !! 
-         !! blank 2D diagnostic array
-         trc2d(:,:,:) = 0.e0
-         !!
-         !! blank 3D diagnostic array
-         trc3d(:,:,:,:) = 0.e0
-# endif
-
-
-         !!----------------------------------------------------------------
-         !! Add in XML diagnostics stuff
-         !!----------------------------------------------------------------
-         !!
-         !! ** 2D diagnostics
-         DO jn=1,jp_medusa_2d
-            CALL iom_put(TRIM(ctrc2d(jn)), trc2d(:,:,jn))
-         END DO
-!! AXY (17/02/14): don't think I need this if I modify the above for 
-!! all diagnostics
-!! #  if defined key_roam
-!!          DO jn=91,jp_medusa_2d
-!!             CALL iom_put(TRIM(ctrc2d(jn)), trc2d(:,:,jn))
-!!          END DO      
-!! #  endif
-         !!
-         !! ** 3D diagnostics
-         DO jn=1,jp_medusa_3d
-            CALL iom_put(TRIM(ctrc3d(jn)), trc3d(:,:,:,jn))
-         END DO
-!! AXY (17/02/14): don't think I need this if I modify the above for 
-!! all diagnostics
-!! #  if defined key_roam
-!!          CALL iom_put(TRIM(ctrc3d(5)), trc3d(:,:,:,5))
-!! #  endif
-
-
-      ELSE IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN
+
          !!!---------------------------------------------------------------
          !! Add very last diag calculations 
@@ -529,5 +339,5 @@
          !! ** 2D diagnostics
 #   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: export all diag.'
+         IF (lwp) write (numout,*) 'bio_medusa_fin: export all diag kt = ', kt
          CALL flush(numout)
 #   endif
@@ -1181,6 +991,4 @@
           DEALLOCATE( zw2d )
 
-       ENDIF                    ! end of ln_diatrc option
-
    END SUBROUTINE bio_medusa_fin
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 8442)
@@ -35,6 +35,6 @@
       USE par_oce,           ONLY: jpi, jpj, jpk
       USE sms_medusa,        ONLY: jdms
-      USE trc,               ONLY: ln_diatrc, med_diag, nittrc000,     &
-                                   trc2d, trc3d
+      USE trc,               ONLY: ln_diatrc, med_diag, nittrc000 
+      USE in_out_manager,    ONLY: lwp
 
 # if defined key_iomput
@@ -47,9 +47,6 @@
 
       IF( ln_diatrc ) THEN
-         !! blank 2D diagnostic array
-         trc2d(:,:,:) = 0.e0
-         !!
-         !! blank 3D diagnostic array
-         trc3d(:,:,:,:) = 0.e0
+         IF (lwp) write (numout,*) 'Diagnostics are now ALL through XIOS (key_xios)'
+         IF (lwp) write (numout,*) 'No more key_diatrc anymore.'
       ENDIF
 
@@ -171,5 +168,5 @@
       !! -----------------------------
       !! Juju :: add kt condition !!
-      IF ( lk_iomput .AND. .NOT.  ln_diatrc ) THEN 
+      IF ( lk_iomput ) THEN 
 
          !! initialise iom_use test
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 8442)
@@ -161,8 +161,8 @@
 
    !! Add DMS in MEDUSA for UKESM1 model
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_surf
+   REAL(wp)                              :: dms_surf,dms_andm
    !! AXY (13/03/15): add in other DMS calculations
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_andr,dms_simo,dms_aran,dms_hall
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_andm, dms_nlim, dms_wtkn
+   REAL(wp)                              :: dms_andr,dms_simo,dms_aran,dms_hall
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_nlim, dms_wtkn
 # endif
 
@@ -360,7 +360,4 @@
                foxy_prod(jpi,jpj), foxy_cons(jpi,jpj),                &
                foxy_anox(jpi,jpj),                                    &
-               dms_surf(jpi,jpj),                                     &
-               dms_andr(jpi,jpj),dms_simo(jpi,jpj),                   &
-               dms_aran(jpi,jpj),dms_hall(jpi,jpj),dms_andm(jpi,jpj), &
                dms_nlim(jpi,jpj),dms_wtkn(jpi,jpj),                   &
 # endif
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 8442)
@@ -781,6 +781,6 @@
 #   if defined key_debug_medusa
       IF (lwp) write (numout,*) '------'
-      IF (lwp) write (numout,*) 'trc_bio_medusa: end all calculations'
-      IF (lwp) write (numout,*) 'trc_bio_medusa: now outputs'
+      IF (lwp) write (numout,*) 'bio_medusa_update: end all calculations'
+      IF (lwp) write (numout,*) 'bio_medusa_update: now outputs kt = ', kt
       CALL flush(numout)
 #   endif
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8442)
@@ -99,7 +99,5 @@
       USE sbc_oce,                    ONLY: lk_oasis
       USE sms_medusa,                 ONLY: hist_pco2
-      USE trc,                        ONLY: ln_diatrc, ln_rsttr,            &
-                                            nittrc000, trn
- 
+      USE trc,                        ONLY: ln_rsttr, nittrc000, trn
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
       USE carb_chem_mod,              ONLY: carb_chem
@@ -640,5 +638,5 @@
          !! 2d specific k level diags
          !!-------------------------------------------------------
-         IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
+         IF( lk_iomput ) THEN
             CALL bio_medusa_diag_slice( jk )
          ENDIF
@@ -651,11 +649,4 @@
       !!------------------------------------------------------------------
       CALL bio_medusa_fin( kt )
-
-# if defined key_trc_diabio
-       !! Lateral boundary conditions on trcbio
-       DO jn=1,jp_medusa_trd
-          CALL lbc_lnk(trbio(:,:,1,jn),'T',1. )
-       ENDDO 
-# endif
 
 # if defined key_debug_medusa
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 8442)
@@ -93,11 +93,6 @@
       INTEGER :: jl, jn
       INTEGER :: ios                 ! Local integer output status for namelist read
-      TYPE(DIAG), DIMENSION(jp_medusa_2d)  :: meddia2d
-      TYPE(DIAG), DIMENSION(jp_medusa_3d)  :: meddia3d
-      TYPE(DIAG), DIMENSION(jp_medusa_trd) :: meddiabio
       CHARACTER(LEN=32)   ::   clname
       !!
-      NAMELIST/nammeddia/ meddia3d, meddia2d     ! additional diagnostics
-
       !!----------------------------------------------------------------------
 
@@ -126,94 +121,4 @@
 # if defined key_debug_medusa
       CALL flush(numout)
-# endif
-      !
-# if defined key_debug_medusa
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'Just before reading namelist_medusa :: nammeddia'
-      IF (lwp) write (numout,*) ' '
-      CALL flush(numout)
-# endif
-
-     IF( ( .NOT.lk_iomput .AND. ln_diatrc ) .OR. ( ln_diatrc .AND. lk_medusa ) ) THEN
-         !
-         ! Namelist nammeddia
-         ! -------------------
-         REWIND( numnatp_ref )              ! Namelist nammeddia in reference namelist : MEDUSA diagnostics
-         READ  ( numnatp_ref, nammeddia, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in reference namelist', lwp )
-
-         REWIND( numnatp_cfg )              ! Namelist nammeddia in configuration namelist : MEDUSA diagnostics
-         READ  ( numnatp_cfg, nammeddia, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in configuration namelist', lwp )
-         IF(lwm) WRITE ( numonp, nammeddia )
-
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) 'reading namelist_medusa :: nammeddia OK'
-         IF (lwp) write (numout,*) 'Check number of variable in nammeddia:'
-         IF (lwp) write (numout,*) 'jp_medusa_2d: ',jp_medusa_2d ,'jp_medusa_3d: ',jp_medusa_3d
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
-         DO jl = 1, jp_medusa_2d
-            jn = jp_msa0_2d + jl - 1
-# if defined key_debug_medusa
-            IF (lwp) write (numout,*) 'Check what is readden in nammeddia: -- 2D'
-            IF (lwp) write (numout,*) jl,'meddia2d-sname: ',meddia2d(jl)%sname 
-            IF (lwp) write (numout,*) jl,'meddia2d-lname: ',meddia2d(jl)%lname 
-            IF (lwp) write (numout,*) jl,'meddia2d-units: ',meddia2d(jl)%units 
-            CALL flush(numout)
-# endif
-            ctrc2d(jn) = meddia2d(jl)%sname
-            ctrc2l(jn) = meddia2d(jl)%lname
-            ctrc2u(jn) = meddia2d(jl)%units
-         END DO
-
-         DO jl = 1, jp_medusa_3d
-            jn = jp_msa0_3d + jl - 1
-# if defined key_debug_medusa
-            IF (lwp) write (numout,*) 'Check what is readden in nammeddia: -- 3D'
-            IF (lwp) write (numout,*) jl,'meddia3d-sname: ',meddia3d(jl)%sname 
-            IF (lwp) write (numout,*) jl,'meddia3d-lname: ',meddia3d(jl)%lname
-            IF (lwp) write (numout,*) jl,'meddia3d-units: ',meddia3d(jl)%units
-            CALL flush(numout)
-# endif
-            ctrc3d(jn) = meddia3d(jl)%sname
-            ctrc3l(jn) = meddia3d(jl)%lname
-            ctrc3u(jn) = meddia3d(jl)%units
-         END DO
-
-         IF(lwp) THEN                   ! control print
-# if defined key_debug_medusa
-            IF (lwp) write (numout,*) '------------------------------'
-            IF (lwp) write (numout,*) 'Jpalm - debug'
-            IF (lwp) write (numout,*) 'Var name assignation OK'
-            IF (lwp) write (numout,*) 'next check var names'
-            IF (lwp) write (numout,*) ' '
-            CALL flush(numout)
-# endif
-            WRITE(numout,*)
-            WRITE(numout,*) ' Namelist : natadd'
-            DO jl = 1, jp_medusa_3d
-               jn = jp_msa0_3d + jl - 1
-               WRITE(numout,*) '  3d diag nb : ', jn, '    short name : ', ctrc3d(jn), &
-                 &             '  long name  : ', ctrc3l(jn), '   unit : ', ctrc3u(jn)
-            END DO
-            WRITE(numout,*) ' '
-
-            DO jl = 1, jp_medusa_2d
-               jn = jp_msa0_2d + jl - 1
-               WRITE(numout,*) '  2d diag nb : ', jn, '    short name : ', ctrc2d(jn), &
-                 &             '  long name  : ', ctrc2l(jn), '   unit : ', ctrc2u(jn)
-            END DO
-            WRITE(numout,*) ' '
-         ENDIF
-         !
-      ENDIF   
-         !
-# if defined key_debug_medusa
-            CALL flush(numout)
 # endif
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcrst_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcrst_medusa.F90	(revision 8441)
+++ 	(revision )
@@ -1,306 +1,0 @@
-MODULE trcrst_medusa
-   !!======================================================================
-   !!                       ***  MODULE trcrst_medusa  ***
-   !! TOP :   create, write, read the restart files of MEDUSA tracer
-   !!======================================================================
-   !! History :   1.0  !  2010-01 (C. Ethe) Original
-   !!             1.1  !  2014-07 (A. Yool; J. Palmieri) converted for MEDUSA
-   !!             1.2  !  2015-07 (A. Yool) add averaged fields for DMS
-   !!----------------------------------------------------------------------
-#if defined key_medusa
-   !!----------------------------------------------------------------------
-   !!   'key_medusa'                                               medusa tracers
-   !!----------------------------------------------------------------------
-   !!   trc_rst_read_medusa   : read  restart file
-   !!   trc_rst_wri_medusa    : write restart file
-   !!----------------------------------------------------------------------
-   USE oce_trc         ! Ocean variables
-   USE par_trc         ! TOP parameters
-   USE trc             ! TOP variables
-   USE trcsms_medusa   ! MEDUSA sms trends
-   USE sms_medusa      ! MEDUSA sms trends
-   USE iom
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC  trc_rst_read_medusa  ! called by trcini.F90 module (actually trcrst.F90)
-   PUBLIC  trc_rst_wri_medusa   ! called by trcini.F90 module (actually trcrst.F90)
-
-CONTAINS
-   
-   SUBROUTINE trc_rst_read_medusa( knum ) 
-      !!----------------------------------------------------------------------
-      !!                     ***  trc_rst_read_medusa  ***  
-      !!
-      !! ** Purpose : Read in restart file specific variables from medusa model
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)  :: knum  ! unit of the restart file
-      !! AXY (07/07/14): temporary variables
-      REAL(wp) ::    fq0,fq1,fq2
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' trc_rst_read_medusa : Read specific variables from medusa model '
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
-
-      !! AXY (07/07/14): remove LOBSTER calls
-      !! CALL iom_get( knum, jpdom_autoglo, 'SEDB'//ctrcnm(jp_lob_det), sedpocb(:,:) ) 
-      !! CALL iom_get( knum, jpdom_autoglo, 'SEDN'//ctrcnm(jp_lob_det), sedpocn(:,:) ) 
-
-      !! AXY (07/07/14): as well as passive tracers, the restart files
-      !!                 contain 2D fields of sediments; these need to
-      !!                 be read in specially; hence this routine
-      !!
-      !! AXY (13/01/12): check if the restart contains sediment fields;
-      !!                 this is only relevant for simulations that include
-      !!                 biogeochemistry and are restarted from earlier runs
-      !!                 in which there was no sediment component
-      !!
-      IF( iom_varid( knum, 'B_SED_N', ldstop = .FALSE. ) > 0 ) THEN
-         !! YES; in which case read them
-         !!
-         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields present - reading in ...'
-         CALL iom_get( knum, jpdom_autoglo, 'B_SED_N',  zb_sed_n(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_SED_N',  zn_sed_n(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'B_SED_FE', zb_sed_fe(:,:) )
-         CALL iom_get( knum, jpdom_autoglo, 'N_SED_FE', zn_sed_fe(:,:) )
-         CALL iom_get( knum, jpdom_autoglo, 'B_SED_SI', zb_sed_si(:,:) )
-         CALL iom_get( knum, jpdom_autoglo, 'N_SED_SI', zn_sed_si(:,:) )
-         CALL iom_get( knum, jpdom_autoglo, 'B_SED_C',  zb_sed_c(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_SED_C',  zn_sed_c(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'B_SED_CA', zb_sed_ca(:,:) )
-         CALL iom_get( knum, jpdom_autoglo, 'N_SED_CA', zn_sed_ca(:,:) )
-      ELSE
-         !! NO; in which case set them to zero
-         !!
-         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields absent - setting to zero ...'
-         zb_sed_n(:,:)  = 0.0   !! organic N
-         zn_sed_n(:,:)  = 0.0
-         zb_sed_fe(:,:) = 0.0   !! organic Fe
-         zn_sed_fe(:,:) = 0.0
-         zb_sed_si(:,:) = 0.0   !! inorganic Si
-         zn_sed_si(:,:) = 0.0
-         zb_sed_c(:,:)  = 0.0   !! organic C
-         zn_sed_c(:,:)  = 0.0
-         zb_sed_ca(:,:) = 0.0   !! inorganic C
-         zn_sed_ca(:,:) = 0.0
-      ENDIF
-      !!
-      !! calculate stats on these fields
-      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
-      fq0 = MINVAL(zn_sed_n(:,:))
-      fq1 = MAXVAL(zn_sed_n(:,:))
-      fq2 = SUM(zn_sed_n(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  N ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_fe(:,:))
-      fq1 = MAXVAL(zn_sed_fe(:,:))
-      fq2 = SUM(zn_sed_fe(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Fe ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_si(:,:))
-      fq1 = MAXVAL(zn_sed_si(:,:))
-      fq2 = SUM(zn_sed_si(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Si ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_c(:,:))
-      fq1 = MAXVAL(zn_sed_c(:,:))
-      fq2 = SUM(zn_sed_c(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  C ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_ca(:,:))
-      fq1 = MAXVAL(zn_sed_ca(:,:))
-      fq2 = SUM(zn_sed_ca(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Ca ', &
-         &        fq0, fq1, fq2
-
-      !! AXY (07/07/15): read in temporally averaged fields for DMS
-      !!                 calculations
-      !!
-      IF( iom_varid( knum, 'B_DMS_CHN', ldstop = .FALSE. ) > 0 ) THEN
-         !! YES; in which case read them
-         !!
-         IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS present - reading in ...'
-         CALL iom_get( knum, jpdom_autoglo, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'B_DMS_DIN',  zb_dms_din(:,:)  )
-         CALL iom_get( knum, jpdom_autoglo, 'N_DMS_DIN',  zn_dms_din(:,:)  )
-      ELSE
-         !! NO; in which case set them to zero
-         !!
-         IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS absent - setting to zero ...'
-         zb_dms_chn(:,:)  = 0.0   !! CHN
-         zn_dms_chn(:,:)  = 0.0
-         zb_dms_chd(:,:)  = 0.0   !! CHD
-         zn_dms_chd(:,:)  = 0.0
-         zb_dms_mld(:,:)  = 0.0   !! MLD
-         zn_dms_mld(:,:)  = 0.0
-         zb_dms_qsr(:,:)  = 0.0   !! QSR
-         zn_dms_qsr(:,:)  = 0.0
-         zb_dms_din(:,:)  = 0.0   !! DIN
-         zn_dms_din(:,:)  = 0.0
-      ENDIF
-      !!
-      !! calculate stats on these fields
-      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
-      fq0 = MINVAL(zn_dms_chn(:,:))
-      fq1 = MAXVAL(zn_dms_chn(:,:))
-      fq2 = SUM(zn_dms_chn(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, CHN ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_chd(:,:))
-      fq1 = MAXVAL(zn_dms_chd(:,:))
-      fq2 = SUM(zn_dms_chd(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, CHD ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_mld(:,:))
-      fq1 = MAXVAL(zn_dms_mld(:,:))
-      fq2 = SUM(zn_dms_mld(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, MLD ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_qsr(:,:))
-      fq1 = MAXVAL(zn_dms_qsr(:,:))
-      fq2 = SUM(zn_dms_qsr(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, QSR ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_din(:,:))
-      fq1 = MAXVAL(zn_dms_din(:,:))
-      fq2 = SUM(zn_dms_din(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, DIN ', fq0, fq1, fq2
-
-   END SUBROUTINE trc_rst_read_medusa
-
-   SUBROUTINE trc_rst_wri_medusa( kt, kitrst, knum )
-      !!----------------------------------------------------------------------
-      !!                     ***  trc_rst_read_medusa  ***
-      !!
-      !! ** Purpose : Read in restart file specific variables from medusa model
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)  :: kt      ! time step
-      INTEGER, INTENT(in)  :: kitrst  ! time step of restart write
-      INTEGER, INTENT(in)  :: knum    ! unit of the restart file
-      !! AXY (07/07/14): temporary variables
-      REAL(wp) ::    fq0,fq1,fq2
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' trc_rst_wri_medusa : Write specific variables from medusa model '
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
-
-      !! AXY (07/07/14): remove LOBSTER calls
-      !! CALL iom_rstput( kt, kitrst, knum, 'SEDB'//ctrcnm(jp_lob_det), sedpocb(:,:) )
-      !! CALL iom_rstput( kt, kitrst, knum, 'SEDN'//ctrcnm(jp_lob_det), sedpocn(:,:) )
-
-      !! AXY (07/07/14): as well as passive tracers, the restart files
-      !!                 contain 2D fields of sediments; these need to
-      !!                 be written out specially; hence this routine
-      !!
-      !! AXY (13/01/12): write out "before" and "now" state of seafloor
-      !!                 sediment pools into restart; this happens
-      !!                 whether or not the pools are to be used by
-      !!                 MEDUSA (which is controlled by a switch in the
-      !!                 namelist_top file)
-      !!
-      IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields - writing out ...'
-      CALL iom_rstput( kt, kitrst, knum, 'B_SED_N',  zb_sed_n(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_SED_N',  zn_sed_n(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'B_SED_FE', zb_sed_fe(:,:) )
-      CALL iom_rstput( kt, kitrst, knum, 'N_SED_FE', zn_sed_fe(:,:) )
-      CALL iom_rstput( kt, kitrst, knum, 'B_SED_SI', zb_sed_si(:,:) )
-      CALL iom_rstput( kt, kitrst, knum, 'N_SED_SI', zn_sed_si(:,:) )
-      CALL iom_rstput( kt, kitrst, knum, 'B_SED_C',  zb_sed_c(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_SED_C',  zn_sed_c(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'B_SED_CA', zb_sed_ca(:,:) )
-      CALL iom_rstput( kt, kitrst, knum, 'N_SED_CA', zn_sed_ca(:,:) )
-      !!
-      !! calculate stats on these fields
-      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
-      fq0 = MINVAL(zn_sed_n(:,:))
-      fq1 = MAXVAL(zn_sed_n(:,:))
-      fq2 = SUM(zn_sed_n(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  N ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_fe(:,:))
-      fq1 = MAXVAL(zn_sed_fe(:,:))
-      fq2 = SUM(zn_sed_fe(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Fe ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_si(:,:))
-      fq1 = MAXVAL(zn_sed_si(:,:))
-      fq2 = SUM(zn_sed_si(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Si ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_c(:,:))
-      fq1 = MAXVAL(zn_sed_c(:,:))
-      fq2 = SUM(zn_sed_c(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  C ', &
-         &        fq0, fq1, fq2
-      fq0 = MINVAL(zn_sed_ca(:,:))
-      fq1 = MAXVAL(zn_sed_ca(:,:))
-      fq2 = SUM(zn_sed_ca(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Ca ', &
-         &        fq0, fq1, fq2
-
-      !! AXY (07/07/15): write out temporally averaged fields for DMS
-      !!                 calculations
-      !!
-      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS - writing out ...'
-      CALL iom_rstput( kt, kitrst, knum, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'B_DMS_DIN',  zb_dms_din(:,:)  )
-      CALL iom_rstput( kt, kitrst, knum, 'N_DMS_DIN',  zn_dms_din(:,:)  )
-      !!
-      !! calculate stats on these fields
-      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
-      fq0 = MINVAL(zn_dms_chn(:,:))
-      fq1 = MAXVAL(zn_dms_chn(:,:))
-      fq2 = SUM(zn_dms_chn(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, CHN ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_chd(:,:))
-      fq1 = MAXVAL(zn_dms_chd(:,:))
-      fq2 = SUM(zn_dms_chd(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, CHD ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_mld(:,:))
-      fq1 = MAXVAL(zn_dms_mld(:,:))
-      fq2 = SUM(zn_dms_mld(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, MLD ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_qsr(:,:))
-      fq1 = MAXVAL(zn_dms_qsr(:,:))
-      fq2 = SUM(zn_dms_qsr(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, QSR ', fq0, fq1, fq2
-      fq0 = MINVAL(zn_dms_din(:,:))
-      fq1 = MAXVAL(zn_dms_din(:,:))
-      fq2 = SUM(zn_dms_din(:,:))
-      if (lwp) write (numout,'(a,3f15.5)') 'DMS, DIN ', fq0, fq1, fq2
-
-   END SUBROUTINE trc_rst_wri_medusa
-
-#else
-   !!----------------------------------------------------------------------
-   !!  Dummy module :                                     No passive tracer
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_rst_read_medusa( knum )
-     INTEGER, INTENT(in)  :: knum
-     WRITE(*,*) 'trc_rst_wri_medusa: You should not have seen this print! error?',knum
-   END SUBROUTINE trc_rst_read_medusa
-
-   SUBROUTINE trc_rst_wri_medusa( kt, kitrst, knum )
-     INTEGER, INTENT(in)  :: kt, kitrst, knum
-     WRITE(*,*) 'trc_rst_wri_medusa: You should not have seen this print! error?', kt, kitrst, knum
-   END SUBROUTINE trc_rst_wri_medusa
-#endif
-
-   !!======================================================================
-END MODULE trcrst_medusa
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90	(revision 8442)
@@ -45,7 +45,7 @@
 
    !! * Module variables
-   INTEGER ::                   &
-     ryyss,                     &  !: number of seconds per year
-     rmtss                         !: number of seconds per month
+   !! INTEGER ::                   &
+     !! ryyss,                     &  !: number of seconds per year
+     !! rmtss                         !: number of seconds per month
 
    !! AXY (10/02/09)
@@ -123,6 +123,6 @@
 
       ! Number of seconds per year and per month
-      ryyss = nyear_len(1) * rday
-      rmtss = ryyss / raamo
+      !! ryyss = nyear_len(1) * rday
+      !! rmtss = ryyss / raamo
 
       !! AXY (20/11/14): alter this to report on first MEDUSA call
@@ -173,14 +173,7 @@
                ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
                tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + ztra
-# if defined key_trc_diabio
-               trbio(ji,jj,jk,8) = ztra
-# endif
-               IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
-                     IF( med_diag%DSED%dgsave ) THEN
-                         zw2d(ji,jj) = zw2d(ji,jj) + ztra * fse3t(ji,jj,jk) * 86400.
-                      ENDIF   
-               ELSE IF( ln_diatrc )  THEN
-                    trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
-               ENDIF    
+               IF( med_diag%DSED%dgsave ) THEN
+                   zw2d(ji,jj) = zw2d(ji,jj) + ztra * fse3t(ji,jj,jk) * 86400.
+               ENDIF   
                 
             END DO
@@ -188,16 +181,7 @@
       END DO
       !
-# if defined key_trc_diabio
-      CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
-# endif
-      IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
-      !!
-      IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
-           IF( med_diag%DSED%dgsave ) THEN
-                CALL iom_put( "DSED"  ,  zw2d)
-                CALL wrk_dealloc( jpi, jpj,    zw2d  )
-            ENDIF
-      ELSE IF (lk_iomput .AND. ln_diatrc)  THEN    
-          CALL iom_put( "DSED",trc2d(:,:,8) )
+      IF( med_diag%DSED%dgsave ) THEN
+           CALL iom_put( "DSED"  ,  zw2d)
+           CALL wrk_dealloc( jpi, jpj,    zw2d  )
       ENDIF
       !!
@@ -229,20 +213,8 @@
                ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
                tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + ztra
-!! #  if defined key_trc_diabio
-!!                trbio(ji,jj,jk,8) = ztra
-!! #  endif
-!!             IF( ln_diatrc ) &
-!!                &  trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
             END DO
          END DO
       END DO
       !
-!! #  if defined key_trc_diabio
-!!       CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
-!! #  endif
-!!       IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
-!! #  if defined key_iomput
-!!       CALL iom_put( "DSED",trc2d(:,:,8) )
-!! #  endif
 
 # endif
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 8442)
@@ -29,5 +29,4 @@
    USE trdtra
    USE prtctl_trc      ! Print control
-   !! USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
@@ -109,11 +108,4 @@
       zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
       zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
-      ! 
-      !! Jpalm -- 14-01-2016 -- restart and proc pb - try this... 
-      !! DO jn = 1, jptra
-      !!   CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
-      !!   CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
-      !! END DO
-      !
 
       IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary)
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 8442)
@@ -77,10 +77,4 @@
          IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends
                                 CALL trc_adv( kstp )            ! horizontal & vertical advection 
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA trc_trp after trc_adv at kt =', kstp
-         CALL trc_rst_tra_stat
-         CALL flush(numout)
-# endif
-
          IF( ln_zps ) THEN
            IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kstp, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
@@ -95,9 +89,4 @@
 #endif
                                 CALL trc_zdf( kstp )            ! vertical mixing and after tracer fields
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA trc_trp after trc_zdf at kt =', kstp
-         CALL trc_rst_tra_stat
-         CALL flush(numout)
-# endif
                                 CALL trc_nxt( kstp )            ! tracer fields at next time step     
 # if defined key_debug_medusa
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/par_trc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/par_trc.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/par_trc.F90	(revision 8442)
@@ -15,8 +15,8 @@
    USE par_c14b      ! C14 bomb tracer
    USE par_cfc       ! CFC 11 and 12 tracers
+   USE par_age       ! AGE  tracer
    USE par_my_trc    ! user defined passive tracers
+   USE par_idtra     ! Idealize tracer
    USE par_medusa    ! MEDUSA model
-   USE par_idtra     ! Idealize tracer
-   USE par_age       ! AGE  tracer
 
    IMPLICIT NONE
@@ -28,9 +28,9 @@
    ! Passive tracers : Total size
    ! ---------------               ! total number of passive tracers, of 2d and 3d output and trend arrays
-   INTEGER, PUBLIC,  PARAMETER ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc    + jp_medusa    + jp_idtra     + jp_age
-   INTEGER, PUBLIC,  PARAMETER ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d + jp_medusa_2d + jp_idtra_2d  + jp_age_2d
-   INTEGER, PUBLIC,  PARAMETER ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d + jp_medusa_3d + jp_idtra_3d  + jp_age_3d
+   INTEGER, PUBLIC,  PARAMETER ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_age    + jp_my_trc    + jp_idtra     + jp_medusa   
+   INTEGER, PUBLIC,  PARAMETER ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_age_2d + jp_my_trc_2d + jp_idtra_2d  + jp_medusa_2d
+   INTEGER, PUBLIC,  PARAMETER ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_age_3d + jp_my_trc_3d + jp_idtra_3d  + jp_medusa_3d
    !                     ! total number of sms diagnostic arrays
-   INTEGER, PUBLIC,  PARAMETER ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd + jp_medusa_trd + jp_idtra_trd + jp_age_trd
+   INTEGER, PUBLIC,  PARAMETER ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_age_trd + jp_my_trc_trd + jp_idtra_trd + jp_medusa_trd 
    
    !  1D configuration ("key_c1d")
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 8442)
@@ -24,8 +24,8 @@
    USE trcini_pisces   ! PISCES   initialisation
    USE trcini_c14b     ! C14 bomb initialisation
+   USE trcini_age      ! AGE      initialisation
    USE trcini_my_trc   ! MY_TRC   initialisation
+   USE trcini_idtra    ! idealize tracer initialisation
    USE trcini_medusa   ! MEDUSA   initialisation
-   USE trcini_idtra    ! idealize tracer initialisation
-   USE trcini_age      ! AGE      initialisation
    USE trcdta          ! initialisation from files
    USE daymod          ! calendar manager
@@ -79,7 +79,5 @@
          &   CALL ctl_warn(' Coupling with passive tracers and used of diurnal cycle. &
          & Computation of a daily mean shortwave for some biogeochemical models) ')
-          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-          !!!!! CHECK For MEDUSA
-          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
       IF( nn_cla == 1 )   &
          &  CALL ctl_stop( ' Cross Land Advection not yet implemented with passive tracer ; nn_cla must be 0' )
@@ -102,22 +100,12 @@
 
       IF( lk_pisces  )       CALL trc_ini_pisces       ! PISCES  bio-model
-      IF( lk_medusa  )       CALL trc_ini_medusa       ! MEDUSA  tracers
-      IF( lk_idtra   )       CALL trc_ini_idtra        ! Idealize tracers
       IF( lk_cfc     )       CALL trc_ini_cfc          ! CFC     tracers
       IF( lk_c14b    )       CALL trc_ini_c14b         ! C14 bomb  tracer
       IF( lk_age     )       CALL trc_ini_age          ! AGE       tracer
       IF( lk_my_trc  )       CALL trc_ini_my_trc       ! MY_TRC  tracers
+      IF( lk_idtra   )       CALL trc_ini_idtra        ! Idealize tracers
+      IF( lk_medusa  )       CALL trc_ini_medusa       ! MEDUSA  tracers
 
       CALL trc_ice_ini                                 ! Tracers in sea ice
-
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' in trc_init'
-         IF (lwp) write (numout,*) ' sms init OK'
-         IF (lwp) write (numout,*) ' next: open tracer.stat'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
 
       IF( ln_ctl ) THEN
@@ -133,19 +121,5 @@
       ENDIF
 
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' in trc_init'
-         IF (lwp) write (numout,*) 'open tracer.stat -- OK'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
-
-
       IF( ln_trcdta ) THEN
-#if defined key_medusa
-         IF(lwp) WRITE(numout,*) 'AXY: calling trc_dta_init'
-         IF(lwp) CALL flush(numout)
-#endif
          CALL trc_dta_init(jptra)
       ENDIF
@@ -153,19 +127,7 @@
       IF( ln_rsttr ) THEN
         !
-#if defined key_medusa
-        IF(lwp) WRITE(numout,*) 'AXY: calling trc_rst_read'
-        IF(lwp) CALL flush(numout)
-#endif
         CALL trc_rst_read              ! restart from a file
         !
       ELSE
-        !
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' Init from file -- will call trc_dta'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
         !
         IF( ln_trcdta .AND. nb_trcdta > 0 ) THEN  ! Initialisation of tracer from a file that may also be used for damping
@@ -188,48 +150,13 @@
         ENDIF
         !
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' in trc_init'
-         IF (lwp) write (numout,*) ' before trb = trn'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
-        !
         trb(:,:,:,:) = trn(:,:,:,:)
         ! 
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' in trc_init'
-         IF (lwp) write (numout,*) ' trb = trn -- OK'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
-        ! 
       ENDIF
  
       tra(:,:,:,:) = 0._wp
       !
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' in trc_init'
-         IF (lwp) write (numout,*) ' partial step -- OK'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
-      !
       IF( nn_dttrc /= 1 )        CALL trc_sub_ini      ! Initialize variables for substepping passive tracers
       !
-# if defined key_debug_medusa
-         IF (lwp) write (numout,*) '------------------------------'
-         IF (lwp) write (numout,*) 'Jpalm - debug'
-         IF (lwp) write (numout,*) ' in trc_init'
-         IF (lwp) write (numout,*) ' before initiate tracer contents'
-         IF (lwp) write (numout,*) ' '
-         CALL flush(numout)
-# endif
-      !
+
       trai(:) = 0._wp                                                   ! initial content of all tracers
       DO jn = 1, jptra
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 8442)
@@ -25,8 +25,8 @@
    USE trcnam_cfc        ! CFC SMS namelist
    USE trcnam_c14b       ! C14 SMS namelist
+   USE trcnam_age        ! AGE SMS namelist
    USE trcnam_my_trc     ! MY_TRC SMS namelist
+   USE trcnam_idtra      ! Idealise tracer namelist
    USE trcnam_medusa     ! MEDUSA namelist
-   USE trcnam_idtra      ! Idealise tracer namelist
-   USE trcnam_age        ! AGE SMS namelist
    USE trd_oce       
    USE trdtrc_oce
@@ -65,49 +65,11 @@
       
       !                                        !  passive tracer informations
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_trc'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
       CALL trc_nam_trc
       
       !                                        !   Parameters of additional diagnostics
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_trc -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_dia'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-
       CALL trc_nam_dia
 
       !                                        !   namelist of transport
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_dia -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_trp'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
       CALL trc_nam_trp
-      !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_trp -- OK'
-      IF (lwp) write (numout,*) 'continue trc_nam '
-      IF (lwp) write (numout,*) ' '
-      CALL flush(numout)
-# endif
-      !
 
 
@@ -131,7 +93,4 @@
          END DO
          WRITE(numout,*) ' '
-# if defined key_debug_medusa
-      CALL flush(numout)
-# endif
       ENDIF
 
@@ -152,7 +111,4 @@
             WRITE(numout,*)
          ENDIF
-# if defined key_debug_medusa
-      CALL flush(numout)
-# endif
       ENDIF
 
@@ -170,7 +126,4 @@
         WRITE(numout,*) '    Passive Tracer  time step    rdttrc  = ', rdttrc(1)
         WRITE(numout,*) 
-# if defined key_debug_medusa
-      CALL flush(numout)
-# endif
       ENDIF
 
@@ -200,30 +153,11 @@
                IF( ln_trdtrc(jn) ) WRITE(numout,*) '    compute ML trends for tracer number :', jn
             END DO
-         WRITE(numout,*) ' '
-         CALL flush(numout)
          ENDIF
 #endif
 
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'just before ice module for tracers call : '
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
 
       ! Call the ice module for tracers
       ! -------------------------------
       CALL trc_nam_ice
-
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'Will now read SMS namelists : '
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
 
       ! namelist of SMS
@@ -232,84 +166,29 @@
       ELSE                    ;   IF(lwp) WRITE(numout,*) '          PISCES not used'
       ENDIF
-      !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_pisces  -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_medusa'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
+
+      IF( lk_cfc     ) THEN   ;   CALL trc_nam_cfc         ! CFC     tracers
+      ELSE                    ;   IF(lwp) WRITE(numout,*) '          CFC not used'
+      ENDIF
+
+      IF( lk_c14b     ) THEN   ;   CALL trc_nam_c14b         ! C14 bomb     tracers
+      ELSE                    ;   IF(lwp) WRITE(numout,*) '          C14 not used'
+      ENDIF
+
+      IF( lk_age     ) THEN  ;   CALL trc_nam_age         ! AGE     tracer
+      ELSE                   ;   IF(lwp) WRITE(numout,*) '          AGE not used'
+      ENDIF
+
+      IF( lk_my_trc  ) THEN   ;   CALL trc_nam_my_trc      ! MY_TRC  tracers
+      ELSE                    ;   IF(lwp) WRITE(numout,*) '          MY_TRC not used'
+      ENDIF
+
+      IF( lk_idtra   ) THEN   ;   CALL trc_nam_idtra       ! Idealize tracers
+      ELSE                    ;   IF(lwp) WRITE(numout,*) '          Idealize tracers not used'
+      ENDIF
+
       IF( lk_medusa  ) THEN   ;   CALL trc_nam_medusa      ! MEDUSA  tracers
       ELSE                    ;   IF(lwp) WRITE(numout,*) '          MEDUSA not used'
       ENDIF
       !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_medusa -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_idtra'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-      IF( lk_idtra   ) THEN   ;   CALL trc_nam_idtra       ! Idealize tracers
-      ELSE                    ;   IF(lwp) WRITE(numout,*) '          Idealize tracers not used'
-      ENDIF
-      !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_idtra -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_cfc'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-      IF( lk_cfc     ) THEN   ;   CALL trc_nam_cfc         ! CFC     tracers
-      ELSE                    ;   IF(lwp) WRITE(numout,*) '          CFC not used'
-      ENDIF
-      !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_cfc -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_c14'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-      IF( lk_c14b     ) THEN   ;   CALL trc_nam_c14b         ! C14 bomb     tracers
-      ELSE                    ;   IF(lwp) WRITE(numout,*) '          C14 not used'
-      ENDIF
-      !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_c14 -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_age'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-      IF( lk_age     ) THEN  ;   CALL trc_nam_age         ! AGE     tracer
-      ELSE                   ;   IF(lwp) WRITE(numout,*) '          AGE not used'
-      ENDIF
-      !
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'CALL trc_nam_age -- OK'
-      IF (lwp) write (numout,*) 'in trc_nam - CALL trc_nam -- OK'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-      IF( lk_my_trc  ) THEN   ;   CALL trc_nam_my_trc      ! MY_TRC  tracers
-      ELSE                    ;   IF(lwp) WRITE(numout,*) '          MY_TRC not used'
-      ENDIF
-       
-      IF(lwp)   CALL flush(numout)
    END SUBROUTINE trc_nam
 
@@ -450,6 +329,5 @@
          ln_trc_wri(jn) =       sn_tracer(jn)%llsave
       END DO
-      IF(lwp)  CALL flush(numout)      
-
+      
     END SUBROUTINE trc_nam_trc
 
@@ -504,13 +382,5 @@
          CALL flush(numout)
       ENDIF
-!!
-!! JPALM -- 17-07-2015 --
-!! MEDUSA is not yet up-to-date with the iom server.
-!! we use it for the main tracer, but not fully with diagnostics.
-!! will have to adapt it properly when visiting Christian Ethee
-!! for now, we change 
-!! IF( ln_diatrc .AND. .NOT. lk_iomput ) THEN
-!! to :
-!!
+
       IF( ( ln_diatrc .AND. .NOT. lk_iomput ) .OR. ( ln_diatrc .AND. lk_medusa ) ) THEN 
          ALLOCATE( trc2d(jpi,jpj,jpdia2d), trc3d(jpi,jpj,jpk,jpdia3d),  &
@@ -522,6 +392,4 @@
          trc3d(:,:,:,:) = 0._wp  ;   ctrc3d(:) = ' '   ;   ctrc3l(:) = ' '    ;    ctrc3u(:) = ' ' 
          !
-      !! ELSE IF  ( lk_iomput .AND. lk_medusa .AND. .NOT. ln_diatrc) THEN
-      !!    CALL trc_nam_iom_medusa
       ENDIF
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcsms.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcsms.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcsms.F90	(revision 8442)
@@ -16,10 +16,10 @@
    USE trc                !
    USE trcsms_pisces      ! PISCES biogeo-model
-   USE trcsms_medusa      ! MEDUSA tracers
-   USE trcsms_idtra       ! Idealize Tracer
    USE trcsms_cfc         ! CFC 11 & 12
    USE trcsms_c14b        ! C14b tracer 
    USE trcsms_age         ! AGE tracer 
    USE trcsms_my_trc      ! MY_TRC  tracers
+   USE trcsms_idtra       ! Idealize Tracer
+   USE trcsms_medusa      ! MEDUSA tracers
    USE prtctl_trc         ! Print control for debbuging
 
@@ -46,5 +46,4 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
       !!
-      INTEGER            ::  jn
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -53,30 +52,10 @@
       !
       IF( lk_pisces  )   CALL trc_sms_pisces ( kt )    ! main program of PISCES 
+      IF( lk_cfc     )   CALL trc_sms_cfc    ( kt )    ! surface fluxes of CFC
+      IF( lk_c14b    )   CALL trc_sms_c14b   ( kt )    ! surface fluxes of C14
+      IF( lk_age     )   CALL trc_sms_age    ( kt )    ! AGE tracer
+      IF( lk_my_trc  )   CALL trc_sms_my_trc ( kt )    ! MY_TRC  tracers
+      IF( lk_idtra   )   CALL trc_sms_idtra  ( kt )    ! radioactive decay of Id. tracer
       IF( lk_medusa  )   CALL trc_sms_medusa ( kt )    ! MEDUSA  tracers
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) '--trcsms : MEDUSA OK --  next IDTRA -- '
-      CALL flush(numout)
-# endif
-      IF( lk_idtra   )   CALL trc_sms_idtra  ( kt )    ! radioactive decay of Id. tracer
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) '--trcsms : IDTRA OK --  next CFC -- '
-      CALL flush(numout)
-# endif
-      IF( lk_cfc     )   CALL trc_sms_cfc    ( kt )    ! surface fluxes of CFC
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) '--trcsms : CFC OK --  next C14 -- '
-      CALL flush(numout)
-# endif
-      IF( lk_c14b    )   CALL trc_sms_c14b   ( kt )    ! surface fluxes of C14
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) '--trcsms : C14 OK --  next C14 -- '
-      CALL flush(numout)
-# endif
-      IF( lk_age     )   CALL trc_sms_age    ( kt )    ! AGE tracer
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) '--trcsms : Age OK --  Continue  -- '
-      CALL flush(numout)
-# endif
-      IF( lk_my_trc  )   CALL trc_sms_my_trc ( kt )    ! MY_TRC  tracers
 
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 8442)
@@ -89,13 +89,5 @@
          tra(:,:,:,:) = 0.e0
          !
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA trc_stp begins at kt =', kt
-         CALL flush(numout)
-# endif
                                    CALL trc_rst_opn  ( kt )       ! Open tracer restart file 
-# if defined key_debug_medusa
-                                   CALL trc_rst_stat 
-                                   CALL trc_rst_tra_stat
-# endif
          IF( lrst_trc )            CALL trc_rst_cal  ( kt, 'WRITE' )   ! calendar
          IF( lk_iomput ) THEN  ;   CALL trc_wri      ( kt )       ! output of passive tracers with iom I/O manager
@@ -124,8 +116,4 @@
          !
          IF( nn_dttrc /= 1   )     CALL trc_sub_reset( kt )       ! resetting physical variables when sub-stepping
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA trc_stp ends at kt =', kt
-         CALL flush(numout)
-# endif
          !
       ENDIF
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcwri.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 8441)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 8442)
@@ -21,8 +21,8 @@
    USE trcwri_cfc
    USE trcwri_c14b
+   USE trcwri_age
    USE trcwri_my_trc
+   USE trcwri_idtra
    USE trcwri_medusa
-   USE trcwri_idtra
-   USE trcwri_age
 
    IMPLICIT NONE
@@ -61,10 +61,10 @@
       ! ---------------------------------------
       IF( lk_pisces  )   CALL trc_wri_pisces     ! PISCES 
-      IF( lk_medusa  )   CALL trc_wri_medusa     ! MESDUSA
-      IF( lk_idtra   )   CALL trc_wri_idtra       ! Idealize tracers
       IF( lk_cfc     )   CALL trc_wri_cfc        ! surface fluxes of CFC
       IF( lk_c14b    )   CALL trc_wri_c14b       ! surface fluxes of C14
       IF( lk_age     )   CALL trc_wri_age        ! AGE tracer
       IF( lk_my_trc  )   CALL trc_wri_my_trc     ! MY_TRC  tracers
+      IF( lk_idtra   )   CALL trc_wri_idtra       ! Idealize tracers
+      IF( lk_medusa  )   CALL trc_wri_medusa     ! MESDUSA
       !
       IF( nn_timing == 1 )  CALL timing_stop('trc_wri')
