Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 7920)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 7920)
@@ -0,0 +1,1111 @@
+MODULE bio_medusa_fin_mod
+   !!======================================================================
+   !!                         ***  MODULE bio_medusa_fin_mod  ***
+   !! Finalisation for TRC_BIO_MEDUSA
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   bio_medusa_fin     ! Called in trcbio_medusa.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bio_medusa_fin( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE bio_medusa_fin  ***
+      !! This called from TRC_BIO_MEDUSA and 
+      !!  - ...
+      !!----------------------------------------------------------------------
+      USE bio_medusa_mod
+      USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt
+      USE in_out_manager,    ONLY: lwp, numout
+      USE iom,               ONLY: iom_put, lk_iomput
+      USE lbclnk,            ONLY: lbc_lnk
+      USE par_medusa,        ONLY: jp_medusa_2d, jp_medusa_3d,          &
+                                   jp_medusa_trd
+      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk
+      USE phycst,            ONLY: rsmall
+      USE sms_medusa,        ONLY: jinorgben, jorgben,                  &
+                                   f3_co3, f3_h2co3, f3_hco3,           &
+                                   f3_omarg, f3_omcal, f3_pH,           &
+                                   za_sed_c, za_sed_ca, za_sed_fe,      &
+                                   za_sed_n, za_sed_si,                 &
+                                   zb_sed_c, zb_sed_ca, zb_sed_fe,      &
+                                   zb_sed_n, zb_sed_si,                 &
+                                   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 trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
+      USE wrk_nemo,          ONLY: wrk_dealloc
+ 
+      !! time (integer timestep)
+      INTEGER, INTENT( in ) ::    kt
+
+      INTEGER :: ji, jj
+      INTEGER :: jn
+      !!----------------------------------------------------------------------
+      !! Process benthic in/out fluxes
+      !! These can be handled outside of the 3D calculations since the
+      !! benthic pools (and fluxes) are 2D in nature; this code is
+      !! (shamelessly) borrowed from corresponding code in the LOBSTER
+      !! model
+      !!----------------------------------------------------------------------
+      !!
+      !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt
+      if (jorgben.eq.1) then
+         za_sed_n(:,:)  = zn_sed_n(:,:)  +                                 & 
+                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -          &
+                            f_benout_n(:,:)  ) * (rdt / 86400.)
+         zn_sed_n(:,:)  = za_sed_n(:,:)
+         !!
+         za_sed_fe(:,:) = zn_sed_fe(:,:) +                                 &
+                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -          &
+                            f_benout_fe(:,:) ) * (rdt / 86400.)
+         zn_sed_fe(:,:) = za_sed_fe(:,:)
+         !!
+         za_sed_c(:,:)  = zn_sed_c(:,:)  +                                 &
+                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -          &
+                            f_benout_c(:,:)  ) * (rdt / 86400.)
+         zn_sed_c(:,:)  = za_sed_c(:,:)
+      endif
+      if (jinorgben.eq.1) then
+         za_sed_si(:,:) = zn_sed_si(:,:) +                                 & 
+                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        &
+                          (rdt / 86400.)
+         zn_sed_si(:,:) = za_sed_si(:,:)
+         !!
+         za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 &
+                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        &
+                          (rdt / 86400.)
+         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
+         !! The code below (in this if ... then ... endif loop) is
+         !! effectively commented out because it does not work as 
+         !! anticipated; it can be deleted at a later date
+         if (jorgben.eq.1) then
+            za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -       &
+                               f_benout_n(:,:)  ) * rdt
+            za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -       &
+                               f_benout_fe(:,:) ) * rdt
+            za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -       &
+                               f_benout_c(:,:)  ) * rdt
+         endif
+         if (jinorgben.eq.1) then
+            za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt
+            za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt
+         endif
+         !!
+         !! Leap-frog scheme - only in explicit case, otherwise the 
+         !! time stepping is already being done in trczdf
+         !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN
+         !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc )
+         !!    IF( neuler == 0 .AND. kt == nittrc000 )  zfact = rdttra(jk) * 
+         !!                                             FLOAT(ndttrc)
+         !!    if (jorgben.eq.1) then
+         !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  )
+         !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) )
+         !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  )
+         !!    endif
+         !!    if (jinorgben.eq.1) then
+         !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) )
+         !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) )
+         !!    endif
+         !! ENDIF
+         !! 
+         !! Time filter and swap of arrays
+         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme
+            IF( neuler == 0 .AND. kt == nittrc000 ) THEN
+               if (jorgben.eq.1) then
+                  zb_sed_n(:,:)  = zn_sed_n(:,:)
+                  zn_sed_n(:,:)  = za_sed_n(:,:)
+                  za_sed_n(:,:)  = 0.0
+                  !!
+                  zb_sed_fe(:,:) = zn_sed_fe(:,:)
+                  zn_sed_fe(:,:) = za_sed_fe(:,:)
+                  za_sed_fe(:,:) = 0.0
+                  !!
+                  zb_sed_c(:,:)  = zn_sed_c(:,:)
+                  zn_sed_c(:,:)  = za_sed_c(:,:)
+                  za_sed_c(:,:)  = 0.0
+               endif
+               if (jinorgben.eq.1) then
+                  zb_sed_si(:,:) = zn_sed_si(:,:)
+                  zn_sed_si(:,:) = za_sed_si(:,:)
+                  za_sed_si(:,:) = 0.0
+                  !!
+                  zb_sed_ca(:,:) = zn_sed_ca(:,:)
+                  zn_sed_ca(:,:) = za_sed_ca(:,:)
+                  za_sed_ca(:,:) = 0.0
+               endif
+            ELSE
+               if (jorgben.eq.1) then
+                  zb_sed_n(:,:)  = (atfp  *                                 &
+                                    ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) +  &
+                                      (atfp1 * zn_sed_n(:,:) )
+                  zn_sed_n(:,:)  = za_sed_n(:,:)
+                  za_sed_n(:,:)  = 0.0
+                  !!
+                  zb_sed_fe(:,:) = (atfp  *                                 &
+                                    ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) +  &
+                                      (atfp1 * zn_sed_fe(:,:))
+                  zn_sed_fe(:,:) = za_sed_fe(:,:)
+                  za_sed_fe(:,:) = 0.0
+                  !!
+                  zb_sed_c(:,:)  = (atfp  *                                 &
+                                    ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) +  &
+                                      (atfp1 * zn_sed_c(:,:) )
+                  zn_sed_c(:,:)  = za_sed_c(:,:)
+                  za_sed_c(:,:)  = 0.0
+               endif
+               if (jinorgben.eq.1) then
+                  zb_sed_si(:,:) = (atfp  *                                 &
+                                    ( zb_sed_si(:,:) + za_sed_si(:,:) )) +  &
+                                      (atfp1 * zn_sed_si(:,:))
+                  zn_sed_si(:,:) = za_sed_si(:,:)
+                  za_sed_si(:,:) = 0.0
+                  !!
+                  zb_sed_ca(:,:) = (atfp  *                                 &
+                                    ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) +  &
+                                      (atfp1 * zn_sed_ca(:,:))
+                  zn_sed_ca(:,:) = za_sed_ca(:,:)
+                  za_sed_ca(:,:) = 0.0
+               endif
+            ENDIF
+         ELSE                   !  case of smolar scheme or muscl
+            if (jorgben.eq.1) then
+               zb_sed_n(:,:)  = za_sed_n(:,:)
+               zn_sed_n(:,:)  = za_sed_n(:,:)
+               za_sed_n(:,:)  = 0.0
+               !!
+               zb_sed_fe(:,:) = za_sed_fe(:,:)
+               zn_sed_fe(:,:) = za_sed_fe(:,:)
+               za_sed_fe(:,:) = 0.0
+               !!
+               zb_sed_c(:,:)  = za_sed_c(:,:)
+               zn_sed_c(:,:)  = za_sed_c(:,:)
+               za_sed_c(:,:)  = 0.0
+            endif
+            if (jinorgben.eq.1) then
+               zb_sed_si(:,:) = za_sed_si(:,:)
+               zn_sed_si(:,:) = za_sed_si(:,:)
+               za_sed_si(:,:) = 0.0
+               !!
+               zb_sed_ca(:,:) = za_sed_ca(:,:)
+               zn_sed_ca(:,:) = za_sed_ca(:,:)
+               za_sed_ca(:,:) = 0.0
+            endif
+         ENDIF
+      endif
+      
+      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 
+         !!!---------------------------------------------------------------
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               !!         
+               IF( med_diag%PN_JLIM%dgsave ) THEN
+                  fjln2d(ji,jj) = fjln2d(ji,jj)   / MAX(ftot_pn(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PN_NLIM%dgsave ) THEN
+                  fnln2d(ji,jj) = fnln2d(ji,jj)   / MAX(ftot_pn(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PN_FELIM%dgsave ) THEN
+                  ffln2d(ji,jj) = ffln2d(ji,jj)   / MAX(ftot_pn(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PD_JLIM%dgsave ) THEN
+                  fjld2d(ji,jj) = fjld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PD_NLIM%dgsave ) THEN
+                  fnld2d(ji,jj) = fnld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PD_FELIM%dgsave ) THEN
+                  ffld2d(ji,jj) = ffld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PD_SILIM%dgsave ) THEN
+                  fsld2d2(ji,jj) = fsld2d2(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
+               ENDIF
+               IF( med_diag%PDSILIM2%dgsave ) THEN
+                  fsld2d(ji,jj) = fsld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
+               ENDIF
+            ENDDO
+         ENDDO
+         !!----------------------------------------------------------------
+         !! Add in XML diagnostics stuff
+         !!----------------------------------------------------------------
+         !!
+         !! ** 2D diagnostics
+#   if defined key_debug_medusa
+         IF (lwp) write (numout,*) 'trc_bio_medusa: export all diag.'
+         CALL flush(numout)
+#   endif
+         IF ( med_diag%INVTN%dgsave ) THEN
+            CALL iom_put( "INVTN"  , ftot_n )
+         ENDIF
+         IF ( med_diag%INVTSI%dgsave ) THEN
+            CALL iom_put( "INVTSI"  , ftot_si )
+         ENDIF
+         IF ( med_diag%INVTFE%dgsave ) THEN
+            CALL iom_put( "INVTFE"  , ftot_fe )
+         ENDIF                           
+         IF ( med_diag%ML_PRN%dgsave ) THEN
+            CALL iom_put( "ML_PRN"  , fprn_ml )
+         ENDIF
+         IF ( med_diag%ML_PRD%dgsave ) THEN
+            CALL iom_put( "ML_PRD"  , fprd_ml )
+         ENDIF
+         IF ( med_diag%OCAL_LVL%dgsave ) THEN
+            CALL iom_put( "OCAL_LVL"  , fccd )
+         ENDIF
+         IF ( med_diag%PN_JLIM%dgsave ) THEN
+            CALL iom_put( "PN_JLIM"  , fjln2d )
+            CALL wrk_dealloc( jpi, jpj,   fjln2d   )
+         ENDIF
+         IF ( med_diag%PN_NLIM%dgsave ) THEN
+            CALL iom_put( "PN_NLIM"  , fnln2d )
+            CALL wrk_dealloc( jpi, jpj,   fnln2d   )
+         ENDIF
+         IF ( med_diag%PN_FELIM%dgsave ) THEN
+            CALL iom_put( "PN_FELIM"  , ffln2d )
+            CALL wrk_dealloc( jpi, jpj,   ffln2d   )
+         ENDIF
+         IF ( med_diag%PD_JLIM%dgsave ) THEN
+            CALL iom_put( "PD_JLIM"  , fjld2d )
+            CALL wrk_dealloc( jpi, jpj,  fjld2d    )
+         ENDIF
+         IF ( med_diag%PD_NLIM%dgsave ) THEN
+            CALL iom_put( "PD_NLIM"  , fnld2d )
+            CALL wrk_dealloc( jpi, jpj,   fnld2d  )
+         ENDIF
+         IF ( med_diag%PD_FELIM%dgsave ) THEN
+            CALL iom_put( "PD_FELIM"  , ffld2d )
+            CALL wrk_dealloc( jpi, jpj,  ffld2d    )
+         ENDIF
+         IF ( med_diag%PD_SILIM%dgsave ) THEN
+            CALL iom_put( "PD_SILIM"  , fsld2d2 )
+            CALL wrk_dealloc( jpi, jpj,   fsld2d2   )
+         ENDIF
+         IF ( med_diag%PDSILIM2%dgsave ) THEN
+            CALL iom_put( "PDSILIM2"  , fsld2d )
+            CALL wrk_dealloc( jpi, jpj,   fsld2d   )
+         ENDIF
+         IF ( med_diag%INTFLX_N%dgsave ) THEN
+            CALL iom_put( "INTFLX_N"  , fflx_n )
+         ENDIF
+         IF ( med_diag%INTFLX_SI%dgsave ) THEN
+            CALL iom_put( "INTFLX_SI"  , fflx_si )
+         ENDIF
+         IF ( med_diag%INTFLX_FE%dgsave ) THEN
+            CALL iom_put( "INTFLX_FE"  , fflx_fe )
+         ENDIF        
+         IF ( med_diag%INT_PN%dgsave ) THEN
+            CALL iom_put( "INT_PN"  , ftot_pn )
+         ENDIF
+         IF ( med_diag%INT_PD%dgsave ) THEN
+            CALL iom_put( "INT_PD"  , ftot_pd )
+         ENDIF         
+         IF ( med_diag%INT_ZMI%dgsave ) THEN
+            CALL iom_put( "INT_ZMI"  , ftot_zmi )
+         ENDIF
+         IF ( med_diag%INT_ZME%dgsave ) THEN
+            CALL iom_put( "INT_ZME"  , ftot_zme )
+         ENDIF
+         IF ( med_diag%INT_DET%dgsave ) THEN
+            CALL iom_put( "INT_DET"  , ftot_det )
+         ENDIF
+         IF ( med_diag%INT_DTC%dgsave ) THEN
+            CALL iom_put( "INT_DTC"  , ftot_dtc )
+         ENDIF
+         IF ( med_diag%BEN_N%dgsave ) THEN
+            CALL iom_put( "BEN_N"  , za_sed_n )
+         ENDIF
+         IF ( med_diag%BEN_FE%dgsave ) THEN
+            CALL iom_put( "BEN_FE"  , za_sed_fe )
+         ENDIF
+         IF ( med_diag%BEN_C%dgsave ) THEN
+            CALL iom_put( "BEN_C"  , za_sed_c )
+         ENDIF
+         IF ( med_diag%BEN_SI%dgsave ) THEN
+            CALL iom_put( "BEN_SI"  , za_sed_si )
+         ENDIF
+         IF ( med_diag%BEN_CA%dgsave ) THEN
+            CALL iom_put( "BEN_CA"  , za_sed_ca )
+         ENDIF
+         IF ( med_diag%RUNOFF%dgsave ) THEN
+            CALL iom_put( "RUNOFF"  , f_runoff )
+         ENDIF 
+# if defined key_roam        
+         IF ( med_diag%N_PROD%dgsave ) THEN
+            CALL iom_put( "N_PROD"  , fnit_prod )
+         ENDIF
+         IF ( med_diag%N_CONS%dgsave ) THEN
+            CALL iom_put( "N_CONS"  , fnit_cons )
+         ENDIF
+         IF ( med_diag%C_PROD%dgsave ) THEN
+            CALL iom_put( "C_PROD"  , fcar_prod )
+         ENDIF
+         IF ( med_diag%C_CONS%dgsave ) THEN
+            CALL iom_put( "C_CONS"  , fcar_cons )
+         ENDIF
+         IF ( med_diag%O2_PROD%dgsave ) THEN
+            CALL iom_put( "O2_PROD"  , foxy_prod )
+         ENDIF
+         IF ( med_diag%O2_CONS%dgsave ) THEN
+            CALL iom_put( "O2_CONS"  , foxy_cons )
+         ENDIF
+         IF ( med_diag%O2_ANOX%dgsave ) THEN
+            CALL iom_put( "O2_ANOX"  , foxy_anox )
+         ENDIF
+         IF ( med_diag%INVTC%dgsave ) THEN
+            CALL iom_put( "INVTC"  , ftot_c )
+         ENDIF
+         IF ( med_diag%INVTALK%dgsave ) THEN
+            CALL iom_put( "INVTALK"  , ftot_a )
+         ENDIF
+         IF ( med_diag%INVTO2%dgsave ) THEN
+            CALL iom_put( "INVTO2"  , ftot_o2 )
+         ENDIF
+         IF ( med_diag%COM_RESP%dgsave ) THEN
+            CALL iom_put( "COM_RESP"  , fcomm_resp )
+         ENDIF         
+# endif      
+         !!
+         !! diagnostic filled in the i-j-k main loop
+         !!--------------------------------------------
+         IF ( med_diag%PRN%dgsave ) THEN
+            CALL iom_put( "PRN"  , fprn2d )
+            CALL wrk_dealloc( jpi, jpj,   fprn2d   )
+         ENDIF
+         IF ( med_diag%MPN%dgsave ) THEN
+            CALL iom_put( "MPN"  ,fdpn2d  )
+            CALL wrk_dealloc( jpi, jpj,    fdpn2d  )
+         ENDIF
+         IF ( med_diag%PRD%dgsave ) THEN
+            CALL iom_put( "PRD"  ,fprd2d  )
+            CALL wrk_dealloc( jpi, jpj,   fprd2d  )
+         ENDIF
+         IF( med_diag%MPD%dgsave ) THEN
+            CALL iom_put( "MPD"  , fdpd2d )
+            CALL wrk_dealloc( jpi, jpj,    fdpd2d )
+         ENDIF
+         !  IF( med_diag%DSED%dgsave ) THEN
+         !      CALL iom_put( "DSED"  , ftot_n )
+         !  ENDIF
+         IF( med_diag%OPAL%dgsave ) THEN
+            CALL iom_put( "OPAL"  , fprds2d )
+            CALL wrk_dealloc( jpi, jpj,   fprds2d  )
+         ENDIF
+         IF( med_diag%OPALDISS%dgsave ) THEN
+            CALL iom_put( "OPALDISS"  , fsdiss2d )
+            CALL wrk_dealloc( jpi, jpj,   fsdiss2d  )
+         ENDIF
+         IF( med_diag%GMIPn%dgsave ) THEN
+            CALL iom_put( "GMIPn"  , fgmipn2d )
+            CALL wrk_dealloc( jpi, jpj,   fgmipn2d  )
+         ENDIF
+         IF( med_diag%GMID%dgsave ) THEN
+            CALL iom_put( "GMID"  , fgmid2d )
+            CALL wrk_dealloc( jpi, jpj,  fgmid2d  )
+         ENDIF
+         IF( med_diag%MZMI%dgsave ) THEN
+            CALL iom_put( "MZMI"  , fdzmi2d )
+            CALL wrk_dealloc( jpi, jpj,   fdzmi2d   )
+         ENDIF
+         IF( med_diag%GMEPN%dgsave ) THEN
+            CALL iom_put( "GMEPN"  , fgmepn2d )
+            CALL wrk_dealloc( jpi, jpj,   fgmepn2d  )
+         ENDIF
+         IF( med_diag%GMEPD%dgsave ) THEN
+            CALL iom_put( "GMEPD"  , fgmepd2d )
+            CALL wrk_dealloc( jpi, jpj,   fgmepd2d   )
+         ENDIF
+         IF( med_diag%GMEZMI%dgsave ) THEN
+            CALL iom_put( "GMEZMI"  , fgmezmi2d )
+            CALL wrk_dealloc( jpi, jpj,   fgmezmi2d   )
+         ENDIF
+         IF( med_diag%GMED%dgsave ) THEN
+            CALL iom_put( "GMED"  , fgmed2d )
+            CALL wrk_dealloc( jpi, jpj,    fgmed2d  )
+         ENDIF
+         IF( med_diag%MZME%dgsave ) THEN
+            CALL iom_put( "MZME"  , fdzme2d )
+            CALL wrk_dealloc( jpi, jpj,  fdzme2d    )
+         ENDIF
+         !  IF( med_diag%DEXP%dgsave ) THEN
+         !      CALL iom_put( "DEXP"  , ftot_n )
+         !  ENDIF
+         IF( med_diag%DETN%dgsave ) THEN
+            CALL iom_put( "DETN"  , fslown2d )
+            CALL wrk_dealloc( jpi, jpj,  fslown2d    )
+         ENDIF
+         IF( med_diag%MDET%dgsave ) THEN
+            CALL iom_put( "MDET"  , fdd2d )
+            CALL wrk_dealloc( jpi, jpj,   fdd2d   )
+         ENDIF
+         IF( med_diag%AEOLIAN%dgsave ) THEN
+            CALL iom_put( "AEOLIAN"  , ffetop2d )
+            CALL wrk_dealloc( jpi, jpj,   ffetop2d   )
+         ENDIF
+         IF( med_diag%BENTHIC%dgsave ) THEN
+            CALL iom_put( "BENTHIC"  , ffebot2d )
+            CALL wrk_dealloc( jpi, jpj,   ffebot2d   )
+         ENDIF
+         IF( med_diag%SCAVENGE%dgsave ) THEN
+            CALL iom_put( "SCAVENGE"  , ffescav2d )
+            CALL wrk_dealloc( jpi, jpj,   ffescav2d  )
+         ENDIF
+         !! 
+         IF( med_diag%TOTREG_N%dgsave ) THEN
+            CALL iom_put( "TOTREG_N"  , fregen2d )
+            CALL wrk_dealloc( jpi, jpj,   fregen2d   )
+         ENDIF
+         IF( med_diag%TOTRG_SI%dgsave ) THEN
+            CALL iom_put( "TOTRG_SI"  , fregensi2d )
+            CALL wrk_dealloc( jpi, jpj,    fregensi2d  )
+         ENDIF
+         !! 
+         IF( med_diag%FASTN%dgsave ) THEN
+            CALL iom_put( "FASTN"  , ftempn2d )
+            CALL wrk_dealloc( jpi, jpj,   ftempn2d   )
+         ENDIF
+         IF( med_diag%FASTSI%dgsave ) THEN
+            CALL iom_put( "FASTSI"  , ftempsi2d )
+            CALL wrk_dealloc( jpi, jpj,   ftempsi2d   )
+         ENDIF
+         IF( med_diag%FASTFE%dgsave ) THEN
+            CALL iom_put( "FASTFE"  , ftempfe2d )
+            CALL wrk_dealloc( jpi, jpj,    ftempfe2d  )
+         ENDIF
+         IF( med_diag%FASTC%dgsave ) THEN
+            CALL iom_put( "FASTC"  , ftempc2d )
+            CALL wrk_dealloc( jpi, jpj,  ftempc2d    )
+         ENDIF
+         IF( med_diag%FASTCA%dgsave ) THEN
+            CALL iom_put( "FASTCA"  , ftempca2d )
+            CALL wrk_dealloc( jpi, jpj,  ftempca2d   )
+         ENDIF
+         !! 
+         IF( med_diag%REMINN%dgsave ) THEN
+            CALL iom_put( "REMINN"  , freminn2d )
+            CALL wrk_dealloc( jpi, jpj,   freminn2d   )
+         ENDIF
+         IF( med_diag%REMINSI%dgsave ) THEN
+            CALL iom_put( "REMINSI"  , freminsi2d )
+            CALL wrk_dealloc( jpi, jpj,   freminsi2d   )
+         ENDIF
+         IF( med_diag%REMINFE%dgsave ) THEN
+            CALL iom_put( "REMINFE"  , freminfe2d )
+            CALL wrk_dealloc( jpi, jpj,  freminfe2d    )
+         ENDIF
+         IF( med_diag%REMINC%dgsave ) THEN
+            CALL iom_put( "REMINC"  , freminc2d )
+            CALL wrk_dealloc( jpi, jpj,    freminc2d  )
+         ENDIF
+         IF( med_diag%REMINCA%dgsave ) THEN
+            CALL iom_put( "REMINCA"  , freminca2d )
+            CALL wrk_dealloc( jpi, jpj,   freminca2d  )
+         ENDIF
+         IF( med_diag%SEAFLRN%dgsave ) THEN
+            CALL iom_put( "SEAFLRN"  , fsedn )
+         ENDIF
+         IF( med_diag%SEAFLRSI%dgsave ) THEN
+            CALL iom_put( "SEAFLRSI"  , fsedsi )
+         ENDIF
+         IF( med_diag%SEAFLRFE%dgsave ) THEN
+            CALL iom_put( "SEAFLRFE"  , fsedfe )
+         ENDIF
+         IF( med_diag%SEAFLRC%dgsave ) THEN
+            CALL iom_put( "SEAFLRC"  , fsedc )
+         ENDIF
+         IF( med_diag%SEAFLRCA%dgsave ) THEN
+            CALL iom_put( "SEAFLRCA"  , fsedca )
+         ENDIF
+         !!
+# if defined key_roam            
+         !!
+         IF( med_diag%RIV_N%dgsave ) THEN
+            CALL iom_put( "RIV_N"  , rivn2d )
+            CALL wrk_dealloc( jpi, jpj,    rivn2d  )
+         ENDIF
+         IF( med_diag%RIV_SI%dgsave ) THEN
+            CALL iom_put( "RIV_SI"  , rivsi2d )
+            CALL wrk_dealloc( jpi, jpj,   rivsi2d   )
+         ENDIF
+         IF( med_diag%RIV_C%dgsave ) THEN
+            CALL iom_put( "RIV_C"  , rivc2d )
+            CALL wrk_dealloc( jpi, jpj,    rivc2d  )
+         ENDIF
+         IF( med_diag%RIV_ALK%dgsave ) THEN
+            CALL iom_put( "RIV_ALK"  , rivalk2d )
+            CALL wrk_dealloc( jpi, jpj,  rivalk2d    )
+         ENDIF
+         IF( med_diag%DETC%dgsave ) THEN
+            CALL iom_put( "DETC"  , fslowc2d )
+            CALL wrk_dealloc( jpi, jpj,   fslowc2d   )
+         ENDIF
+         !!
+         IF( med_diag%PN_LLOSS%dgsave ) THEN
+            CALL iom_put( "PN_LLOSS"  , fdpn22d )
+            CALL wrk_dealloc( jpi, jpj,   fdpn22d   )
+         ENDIF
+         IF( med_diag%PD_LLOSS%dgsave ) THEN
+            CALL iom_put( "PD_LLOSS"  , fdpd22d )
+            CALL wrk_dealloc( jpi, jpj,   fdpd22d   )
+         ENDIF
+         IF( med_diag%ZI_LLOSS%dgsave ) THEN
+            CALL iom_put( "ZI_LLOSS"  , fdzmi22d )
+             CALL wrk_dealloc( jpi, jpj,    fdzmi22d  )
+          ENDIF
+          IF( med_diag%ZE_LLOSS%dgsave ) THEN
+             CALL iom_put( "ZE_LLOSS"  , fdzme22d )
+             CALL wrk_dealloc( jpi, jpj,   fdzme22d   )
+          ENDIF
+          IF( med_diag%ZI_MES_N%dgsave ) THEN
+             CALL iom_put( "ZI_MES_N"  , zimesn2d )
+             CALL wrk_dealloc( jpi, jpj,    zimesn2d  )
+          ENDIF
+          IF( med_diag%ZI_MES_D%dgsave ) THEN
+             CALL iom_put( "ZI_MES_D"  , zimesd2d )
+             CALL wrk_dealloc( jpi, jpj,    zimesd2d  )
+          ENDIF
+          IF( med_diag%ZI_MES_C%dgsave ) THEN
+             CALL iom_put( "ZI_MES_C"  , zimesc2d )
+             CALL wrk_dealloc( jpi, jpj,    zimesc2d  )
+          ENDIF
+          IF( med_diag%ZI_MESDC%dgsave ) THEN
+             CALL iom_put( "ZI_MESDC"  ,zimesdc2d  )
+             CALL wrk_dealloc( jpi, jpj,    zimesdc2d  )
+          ENDIF
+          IF( med_diag%ZI_EXCR%dgsave ) THEN
+             CALL iom_put( "ZI_EXCR"  , ziexcr2d )
+             CALL wrk_dealloc( jpi, jpj,    ziexcr2d )
+          ENDIF
+          IF( med_diag%ZI_RESP%dgsave ) THEN
+             CALL iom_put( "ZI_RESP"  , ziresp2d )
+             CALL wrk_dealloc( jpi, jpj,   ziresp2d   )
+          ENDIF
+          IF( med_diag%ZI_GROW%dgsave ) THEN
+             CALL iom_put( "ZI_GROW"  , zigrow2d )
+             CALL wrk_dealloc( jpi, jpj,   zigrow2d   )
+          ENDIF
+          IF( med_diag%ZE_MES_N%dgsave ) THEN
+             CALL iom_put( "ZE_MES_N"  , zemesn2d )
+             CALL wrk_dealloc( jpi, jpj,    zemesn2d  )
+          ENDIF
+          IF( med_diag%ZE_MES_D%dgsave ) THEN
+             CALL iom_put( "ZE_MES_D"  , zemesd2d )
+             CALL wrk_dealloc( jpi, jpj,    zemesd2d  )
+          ENDIF
+          IF( med_diag%ZE_MES_C%dgsave ) THEN
+             CALL iom_put( "ZE_MES_C"  , zemesc2d )
+             CALL wrk_dealloc( jpi, jpj,   zemesc2d   )
+          ENDIF
+          IF( med_diag%ZE_MESDC%dgsave ) THEN
+             CALL iom_put( "ZE_MESDC"  , zemesdc2d )
+             CALL wrk_dealloc( jpi, jpj,   zemesdc2d   )
+          ENDIF
+          IF( med_diag%ZE_EXCR%dgsave ) THEN
+             CALL iom_put( "ZE_EXCR"  , zeexcr2d )
+             CALL wrk_dealloc( jpi, jpj,   zeexcr2d   )
+          ENDIF
+          IF( med_diag%ZE_RESP%dgsave ) THEN
+             CALL iom_put( "ZE_RESP"  , zeresp2d )
+             CALL wrk_dealloc( jpi, jpj,    zeresp2d  )
+          ENDIF
+          IF( med_diag%ZE_GROW%dgsave ) THEN
+             CALL iom_put( "ZE_GROW"  , zegrow2d )
+             CALL wrk_dealloc( jpi, jpj,   zegrow2d   )
+          ENDIF
+          IF( med_diag%MDETC%dgsave ) THEN
+             CALL iom_put( "MDETC"  , mdetc2d )
+             CALL wrk_dealloc( jpi, jpj,   mdetc2d   )
+          ENDIF
+          IF( med_diag%GMIDC%dgsave ) THEN
+             CALL iom_put( "GMIDC"  , gmidc2d )
+             CALL wrk_dealloc( jpi, jpj,    gmidc2d  )
+          ENDIF
+          IF( med_diag%GMEDC%dgsave ) THEN
+             CALL iom_put( "GMEDC"  , gmedc2d )
+             CALL wrk_dealloc( jpi, jpj,    gmedc2d  )
+          ENDIF
+          IF( med_diag%IBEN_N%dgsave ) THEN
+             CALL iom_put( "IBEN_N"  , iben_n2d )
+             CALL wrk_dealloc( jpi, jpj,    iben_n2d  )
+          ENDIF
+          IF( med_diag%IBEN_FE%dgsave ) THEN
+             CALL iom_put( "IBEN_FE"  , iben_fe2d )
+             CALL wrk_dealloc( jpi, jpj,   iben_fe2d   )
+          ENDIF
+          IF( med_diag%IBEN_C%dgsave ) THEN
+             CALL iom_put( "IBEN_C"  , iben_c2d )
+             CALL wrk_dealloc( jpi, jpj,   iben_c2d   )
+          ENDIF
+          IF( med_diag%IBEN_SI%dgsave ) THEN
+             CALL iom_put( "IBEN_SI"  , iben_si2d )
+             CALL wrk_dealloc( jpi, jpj,   iben_si2d   )
+          ENDIF
+          IF( med_diag%IBEN_CA%dgsave ) THEN
+             CALL iom_put( "IBEN_CA"  , iben_ca2d )
+             CALL wrk_dealloc( jpi, jpj,   iben_ca2d   )
+          ENDIF
+          IF( med_diag%OBEN_N%dgsave ) THEN
+             CALL iom_put( "OBEN_N"  , oben_n2d )
+             CALL wrk_dealloc( jpi, jpj,    oben_n2d  )
+          ENDIF
+          IF( med_diag%OBEN_FE%dgsave ) THEN
+             CALL iom_put( "OBEN_FE"  , oben_fe2d )
+             CALL wrk_dealloc( jpi, jpj,    oben_fe2d  )
+          ENDIF
+          IF( med_diag%OBEN_C%dgsave ) THEN
+             CALL iom_put( "OBEN_C"  , oben_c2d )
+             CALL wrk_dealloc( jpi, jpj,    oben_c2d  )
+          ENDIF
+          IF( med_diag%OBEN_SI%dgsave ) THEN
+             CALL iom_put( "OBEN_SI"  , oben_si2d )
+             CALL wrk_dealloc( jpi, jpj,    oben_si2d  )
+          ENDIF
+          IF( med_diag%OBEN_CA%dgsave ) THEN
+             CALL iom_put( "OBEN_CA"  , oben_ca2d )
+             CALL wrk_dealloc( jpi, jpj, oben_ca2d     )
+          ENDIF
+          IF( med_diag%SFR_OCAL%dgsave ) THEN
+             CALL iom_put( "SFR_OCAL"  , sfr_ocal2d )
+             CALL wrk_dealloc( jpi, jpj,    sfr_ocal2d  )
+          ENDIF
+          IF( med_diag%SFR_OARG%dgsave ) THEN
+             CALL iom_put( "SFR_OARG"  , sfr_oarg2d )
+             CALL wrk_dealloc( jpi, jpj,    sfr_oarg2d  )
+          ENDIF
+          IF( med_diag%LYSO_CA%dgsave ) THEN
+             CALL iom_put( "LYSO_CA"  , lyso_ca2d )
+             CALL wrk_dealloc( jpi, jpj,    lyso_ca2d  )
+          ENDIF
+# endif                   
+          !!
+          !! ** 3D diagnostics
+          IF( med_diag%TPP3%dgsave ) THEN
+             CALL iom_put( "TPP3"  , tpp3d )
+             CALL wrk_dealloc( jpi, jpj, jpk,   tpp3d  )
+          ENDIF
+          IF( med_diag%DETFLUX3%dgsave ) THEN
+             CALL iom_put( "DETFLUX3"  , detflux3d )
+             CALL wrk_dealloc( jpi, jpj, jpk,    detflux3d )
+          ENDIF
+          IF( med_diag%REMIN3N%dgsave ) THEN
+             CALL iom_put( "REMIN3N"  , remin3dn )
+             CALL wrk_dealloc( jpi, jpj, jpk,   remin3dn  )
+          ENDIF
+# if defined key_roam          
+          IF( med_diag%PH3%dgsave ) THEN
+             CALL iom_put( "PH3"  , f3_pH )
+          ENDIF
+          IF( med_diag%OM_CAL3%dgsave ) THEN
+             CALL iom_put( "OM_CAL3"  , f3_omcal )
+          ENDIF
+          !!
+          !! AXY (09/11/16): 2D CMIP6 diagnostics
+          IF( med_diag%INTDISSIC%dgsave ) THEN
+             CALL iom_put( "INTDISSIC"  , intdissic )
+             CALL wrk_dealloc( jpi, jpj, intdissic   )
+          ENDIF          
+          IF( med_diag%INTDISSIN%dgsave ) THEN
+             CALL iom_put( "INTDISSIN"  , intdissin )
+             CALL wrk_dealloc( jpi, jpj, intdissin   )
+          ENDIF          
+          IF( med_diag%INTDISSISI%dgsave ) THEN
+             CALL iom_put( "INTDISSISI"  , intdissisi )
+             CALL wrk_dealloc( jpi, jpj, intdissisi  )
+          ENDIF          
+          IF( med_diag%INTTALK%dgsave ) THEN
+             CALL iom_put( "INTTALK"  , inttalk )
+             CALL wrk_dealloc( jpi, jpj, inttalk     )
+          ENDIF          
+          IF( med_diag%O2min%dgsave ) THEN
+             CALL iom_put( "O2min"  , o2min )
+             CALL wrk_dealloc( jpi, jpj, o2min       )
+          ENDIF          
+          IF( med_diag%ZO2min%dgsave ) THEN
+             CALL iom_put( "ZO2min"  , zo2min )
+             CALL wrk_dealloc( jpi, jpj, zo2min      )
+          ENDIF          
+          IF( med_diag%FBDDTALK%dgsave ) THEN
+             CALL iom_put( "FBDDTALK"  , fbddtalk   )
+             CALL wrk_dealloc( jpi, jpj, fbddtalk   )
+          ENDIF          
+          IF( med_diag%FBDDTDIC%dgsave ) THEN
+             CALL iom_put( "FBDDTDIC"  , fbddtdic   )
+             CALL wrk_dealloc( jpi, jpj, fbddtdic   )
+          ENDIF          
+          IF( med_diag%FBDDTDIFE%dgsave ) THEN
+             CALL iom_put( "FBDDTDIFE" , fbddtdife  )
+             CALL wrk_dealloc( jpi, jpj, fbddtdife  )
+          ENDIF          
+          IF( med_diag%FBDDTDIN%dgsave ) THEN
+             CALL iom_put( "FBDDTDIN"  , fbddtdin   )
+             CALL wrk_dealloc( jpi, jpj, fbddtdin   )
+          ENDIF          
+          IF( med_diag%FBDDTDISI%dgsave ) THEN
+             CALL iom_put( "FBDDTDISI" , fbddtdisi  )
+             CALL wrk_dealloc( jpi, jpj, fbddtdisi  )
+          ENDIF    
+          !!
+          !! AXY (09/11/16): 3D CMIP6 diagnostics
+          IF( med_diag%TPPD3%dgsave ) THEN
+             CALL iom_put( "TPPD3"     , tppd3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, tppd3      )
+          ENDIF          
+          IF( med_diag%BDDTALK3%dgsave ) THEN
+             CALL iom_put( "BDDTALK3"  , bddtalk3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, bddtalk3   )
+          ENDIF          
+          IF( med_diag%BDDTDIC3%dgsave ) THEN
+             CALL iom_put( "BDDTDIC3"  , bddtdic3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, bddtdic3   )
+          ENDIF          
+          IF( med_diag%BDDTDIFE3%dgsave ) THEN
+             CALL iom_put( "BDDTDIFE3" , bddtdife3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, bddtdife3  )
+          ENDIF          
+          IF( med_diag%BDDTDIN3%dgsave ) THEN
+             CALL iom_put( "BDDTDIN3"  , bddtdin3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, bddtdin3   )
+          ENDIF          
+          IF( med_diag%BDDTDISI3%dgsave ) THEN
+             CALL iom_put( "BDDTDISI3" , bddtdisi3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, bddtdisi3  )
+          ENDIF    
+          IF( med_diag%FD_NIT3%dgsave ) THEN
+             CALL iom_put( "FD_NIT3"  , fd_nit3 )
+             CALL wrk_dealloc( jpi, jpj, jpk,   fd_nit3  )
+          ENDIF
+          IF( med_diag%FD_SIL3%dgsave ) THEN
+             CALL iom_put( "FD_SIL3"  , fd_sil3 )
+             CALL wrk_dealloc( jpi, jpj, jpk,   fd_sil3  )
+          ENDIF
+          IF( med_diag%FD_CAL3%dgsave ) THEN
+             CALL iom_put( "FD_CAL3"  , fd_cal3 )
+             CALL wrk_dealloc( jpi, jpj, jpk,   fd_cal3  )
+          ENDIF
+          IF( med_diag%FD_CAR3%dgsave ) THEN
+             CALL iom_put( "FD_CAR3"  , fd_car3 )
+             CALL wrk_dealloc( jpi, jpj, jpk,   fd_car3  )
+          ENDIF
+          IF( med_diag%CO33%dgsave ) THEN
+             CALL iom_put( "CO33"  , f3_co3 )
+          ENDIF          	        
+          IF( med_diag%CO3SATARAG3%dgsave ) THEN
+             CALL iom_put( "CO3SATARAG3"  , f3_omarg )
+          ENDIF          	        
+          IF( med_diag%CO3SATCALC3%dgsave ) THEN
+             CALL iom_put( "CO3SATCALC3"  , f3_omcal )
+          ENDIF          	        
+          IF( med_diag%EXPC3%dgsave ) THEN
+             CALL iom_put( "EXPC3"  , expc3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, expc3  )
+          ENDIF          	        
+          IF( med_diag%EXPN3%dgsave ) THEN
+             CALL iom_put( "EXPN3"  , expn3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, expn3  )
+          ENDIF          	        
+          IF( med_diag%DCALC3%dgsave ) THEN
+             CALL iom_put( "DCALC3"  , dcalc3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, dcalc3  )
+          ENDIF          	        
+          IF( med_diag%FEDISS3%dgsave ) THEN
+             CALL iom_put( "FEDISS3"  , fediss3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, fediss3  )
+          ENDIF          	        
+          IF( med_diag%FESCAV3%dgsave ) THEN
+             CALL iom_put( "FESCAV3"  , fescav3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, fescav3  )
+          ENDIF          	        
+          IF( med_diag%MIGRAZP3%dgsave ) THEN
+             CALL iom_put( "MIGRAZP3"  , migrazp3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, migrazp3  )
+          ENDIF          	        
+          IF( med_diag%MIGRAZD3%dgsave ) THEN
+             CALL iom_put( "MIGRAZD3"  , migrazd3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, migrazd3  )
+          ENDIF          	        
+          IF( med_diag%MEGRAZP3%dgsave ) THEN
+             CALL iom_put( "MEGRAZP3"  , megrazp3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, megrazp3  )
+          ENDIF          	        
+          IF( med_diag%MEGRAZD3%dgsave ) THEN
+             CALL iom_put( "MEGRAZD3"  , megrazd3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, megrazd3  )
+          ENDIF          	        
+          IF( med_diag%MEGRAZZ3%dgsave ) THEN
+             CALL iom_put( "MEGRAZZ3"  , megrazz3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, megrazz3  )
+          ENDIF          	        
+          IF( med_diag%O2SAT3%dgsave ) THEN
+             CALL iom_put( "O2SAT3"  , o2sat3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, o2sat3 )
+          ENDIF          	        
+          IF( med_diag%PBSI3%dgsave ) THEN
+             CALL iom_put( "PBSI3"  , pbsi3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pbsi3  )
+          ENDIF          	        
+          IF( med_diag%PCAL3%dgsave ) THEN
+             CALL iom_put( "PCAL3"  , pcal3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pcal3  )
+          ENDIF          	        
+          IF( med_diag%REMOC3%dgsave ) THEN
+             CALL iom_put( "REMOC3"  , remoc3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, remoc3 )
+          ENDIF          	        
+          IF( med_diag%PNLIMJ3%dgsave ) THEN
+             CALL iom_put( "PNLIMJ3" , pnlimj3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pnlimj3  )
+          ENDIF          	        
+          IF( med_diag%PNLIMN3%dgsave ) THEN
+             CALL iom_put( "PNLIMN3" , pnlimn3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pnlimn3  )
+          ENDIF          	        
+          IF( med_diag%PNLIMFE3%dgsave ) THEN
+             CALL iom_put( "PNLIMFE3" , pnlimfe3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pnlimfe3 )
+          ENDIF          	        
+          IF( med_diag%PDLIMJ3%dgsave ) THEN
+             CALL iom_put( "PDLIMJ3" , pdlimj3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pdlimj3  )
+          ENDIF          	        
+          IF( med_diag%PDLIMN3%dgsave ) THEN
+             CALL iom_put( "PDLIMN3" , pdlimn3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pdlimn3  )
+          ENDIF          	        
+          IF( med_diag%PDLIMFE3%dgsave ) THEN
+             CALL iom_put( "PDLIMFE3" , pdlimfe3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pdlimfe3 )
+          ENDIF          	        
+          IF( med_diag%PDLIMSI3%dgsave ) THEN
+             CALL iom_put( "PDLIMSI3" , pdlimsi3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, pdlimsi3 )
+          ENDIF          	        
+          
+# endif         
+
+          CALL wrk_dealloc( jpi, jpj,   zw2d   )
+
+       ENDIF                    ! end of ln_diatrc option
+
+   END SUBROUTINE bio_medusa_fin
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE bio_medusa_fin( )                    ! Empty routine
+      WRITE(*,*) 'bio_medusa_fin: You should not have seen this print! error?'
+   END SUBROUTINE bio_medusa_fin
+#endif 
+
+   !!======================================================================
+END MODULE bio_medusa_fin_mod
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7919)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7920)
@@ -93,8 +93,10 @@
 
       !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
-      USE sbc_oce, ONLY: lk_oasis
-      USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, PCO2a_in_cpl
+      USE sbc_oce,                ONLY: lk_oasis
+      USE oce,                    ONLY: CO2Flux_out_cpl, DMS_out_cpl,       &
+                                        PCO2a_in_cpl
       USE bio_medusa_mod
-      USE bio_medusa_init_mod, ONLY: bio_medusa_init
+      USE bio_medusa_init_mod,    ONLY: bio_medusa_init
+      USE bio_medusa_fin_mod,     ONLY: bio_medusa_fin
 
       IMPLICIT NONE
@@ -4269,1000 +4271,7 @@
 
       !!----------------------------------------------------------------------
-      !! Process benthic in/out fluxes
-      !! These can be handled outside of the 3D calculations since the
-      !! benthic pools (and fluxes) are 2D in nature; this code is
-      !! (shamelessly) borrowed from corresponding code in the LOBSTER
-      !! model
+      !! Final calculations for diagnostics
       !!----------------------------------------------------------------------
-      !!
-      !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt
-      if (jorgben.eq.1) then
-         za_sed_n(:,:)  = zn_sed_n(:,:)  + & 
-         &                ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ) * (rdt / 86400.)
-         zn_sed_n(:,:)  = za_sed_n(:,:)
-         !!
-         za_sed_fe(:,:) = zn_sed_fe(:,:) + &
-         &                ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * (rdt / 86400.)
-         zn_sed_fe(:,:) = za_sed_fe(:,:)
-         !!
-         za_sed_c(:,:)  = zn_sed_c(:,:)  + &
-         &                ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ) * (rdt / 86400.)
-         zn_sed_c(:,:)  = za_sed_c(:,:)
-      endif
-      if (jinorgben.eq.1) then
-         za_sed_si(:,:) = zn_sed_si(:,:) + & 
-         &                ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * (rdt / 86400.)
-         zn_sed_si(:,:) = za_sed_si(:,:)
-         !!
-         za_sed_ca(:,:) = zn_sed_ca(:,:) + &
-         &                ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * (rdt / 86400.)
-         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
-         !! The code below (in this if ... then ... endif loop) is
-         !! effectively commented out because it does not work as 
-         !! anticipated; it can be deleted at a later date
-         if (jorgben.eq.1) then
-            za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ) * rdt
-            za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt
-            za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ) * rdt
-         endif
-         if (jinorgben.eq.1) then
-            za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt
-            za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt
-         endif
-         !!
-         !! Leap-frog scheme - only in explicit case, otherwise the time stepping
-         !! is already being done in trczdf
-         !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN
-         !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc )
-         !!    IF( neuler == 0 .AND. kt == nittrc000 )   zfact = rdttra(jk) * FLOAT(ndttrc)
-         !!    if (jorgben.eq.1) then
-         !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  )
-         !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) )
-         !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  )
-         !!    endif
-         !!    if (jinorgben.eq.1) then
-         !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) )
-         !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) )
-         !!    endif
-         !! ENDIF
-         !! 
-         !! Time filter and swap of arrays
-         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme
-            IF( neuler == 0 .AND. kt == nittrc000 ) THEN
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = zn_sed_n(:,:)
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = zn_sed_fe(:,:)
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = zn_sed_c(:,:)
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = zn_sed_si(:,:)
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = zn_sed_ca(:,:)
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ELSE
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = (atfp  * ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) + (atfp1 * zn_sed_n(:,:) )
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = (atfp  * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:))
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = (atfp  * ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) + (atfp1 * zn_sed_c(:,:) )
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = (atfp  * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:))
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = (atfp  * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:))
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ENDIF
-         ELSE                   !  case of smolar scheme or muscl
-            if (jorgben.eq.1) then
-               zb_sed_n(:,:)  = za_sed_n(:,:)
-               zn_sed_n(:,:)  = za_sed_n(:,:)
-               za_sed_n(:,:)  = 0.0
-               !!
-               zb_sed_fe(:,:) = za_sed_fe(:,:)
-               zn_sed_fe(:,:) = za_sed_fe(:,:)
-               za_sed_fe(:,:) = 0.0
-               !!
-               zb_sed_c(:,:)  = za_sed_c(:,:)
-               zn_sed_c(:,:)  = za_sed_c(:,:)
-               za_sed_c(:,:)  = 0.0
-            endif
-            if (jinorgben.eq.1) then
-               zb_sed_si(:,:) = za_sed_si(:,:)
-               zn_sed_si(:,:) = za_sed_si(:,:)
-               za_sed_si(:,:) = 0.0
-               !!
-               zb_sed_ca(:,:) = za_sed_ca(:,:)
-               zn_sed_ca(:,:) = za_sed_ca(:,:)
-               za_sed_ca(:,:) = 0.0
-            endif
-         ENDIF
-      endif
-      
-      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
-               trc2d(ji,jj,79)  = ftot_pn(ji,jj)       !! integrated non-diatom phytoplankton
-               trc2d(ji,jj,80)  = ftot_pd(ji,jj)       !! integrated diatom phytoplankton
-               trc2d(ji,jj,217) = ftot_zmi(ji,jj)      !! Integrated microzooplankton
-               trc2d(ji,jj,218) = ftot_zme(ji,jj)      !! Integrated mesozooplankton
-               trc2d(ji,jj,219) = ftot_det(ji,jj)      !! Integrated slow detritus, nitrogen
-               trc2d(ji,jj,220) = ftot_dtc(ji,jj)      !! Integrated slow detritus, carbon
-# if defined key_roam
-               !! the balance of nitrogen production/consumption
-               trc2d(ji,jj,111) = fnit_prod(ji,jj)  !! integrated nitrogen production
-               trc2d(ji,jj,112) = fnit_cons(ji,jj)  !! integrated nitrogen consumption
-               !! the balance of carbon production/consumption
-               trc2d(ji,jj,113) = fcar_prod(ji,jj)  !! integrated carbon production
-               trc2d(ji,jj,114) = fcar_cons(ji,jj)  !! integrated carbon consumption
-               !! the balance of oxygen production/consumption
-               trc2d(ji,jj,115) = foxy_prod(ji,jj)  !! integrated oxygen production
-               trc2d(ji,jj,116) = foxy_cons(ji,jj)  !! integrated oxygen consumption
-               trc2d(ji,jj,117) = foxy_anox(ji,jj)  !! integrated unrealised oxygen consumption
-# 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 
-         !!!----------------------------------------------------------------------
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               !!         
-               IF( med_diag%PN_JLIM%dgsave ) THEN
-                  fjln2d(ji,jj) = fjln2d(ji,jj)   / MAX(ftot_pn(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PN_NLIM%dgsave ) THEN
-                  fnln2d(ji,jj) = fnln2d(ji,jj)   / MAX(ftot_pn(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PN_FELIM%dgsave ) THEN
-                  ffln2d(ji,jj) = ffln2d(ji,jj)   / MAX(ftot_pn(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PD_JLIM%dgsave ) THEN
-                  fjld2d(ji,jj) = fjld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PD_NLIM%dgsave ) THEN
-                  fnld2d(ji,jj) = fnld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PD_FELIM%dgsave ) THEN
-                  ffld2d(ji,jj) = ffld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PD_SILIM%dgsave ) THEN
-                  fsld2d2(ji,jj) = fsld2d2(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
-               ENDIF
-               IF( med_diag%PDSILIM2%dgsave ) THEN
-                  fsld2d(ji,jj) = fsld2d(ji,jj)   / MAX(ftot_pd(ji,jj), rsmall)
-               ENDIF
-            ENDDO
-         ENDDO
-         !!----------------------------------------------------------------------
-         !! Add in XML diagnostics stuff
-         !!----------------------------------------------------------------------
-         !!
-         !! ** 2D diagnostics
-#   if defined key_debug_medusa
-         IF (lwp) write (numout,*) 'trc_bio_medusa: export all diag.'
-         CALL flush(numout)
-#   endif
-         IF ( med_diag%INVTN%dgsave ) THEN
-            CALL iom_put( "INVTN"  , ftot_n )
-         ENDIF
-         IF ( med_diag%INVTSI%dgsave ) THEN
-            CALL iom_put( "INVTSI"  , ftot_si )
-         ENDIF
-         IF ( med_diag%INVTFE%dgsave ) THEN
-            CALL iom_put( "INVTFE"  , ftot_fe )
-         ENDIF                           
-         IF ( med_diag%ML_PRN%dgsave ) THEN
-            CALL iom_put( "ML_PRN"  , fprn_ml )
-         ENDIF
-         IF ( med_diag%ML_PRD%dgsave ) THEN
-            CALL iom_put( "ML_PRD"  , fprd_ml )
-         ENDIF
-         IF ( med_diag%OCAL_LVL%dgsave ) THEN
-            CALL iom_put( "OCAL_LVL"  , fccd )
-         ENDIF
-         IF ( med_diag%PN_JLIM%dgsave ) THEN
-            CALL iom_put( "PN_JLIM"  , fjln2d )
-            CALL wrk_dealloc( jpi, jpj,   fjln2d   )
-         ENDIF
-         IF ( med_diag%PN_NLIM%dgsave ) THEN
-            CALL iom_put( "PN_NLIM"  , fnln2d )
-            CALL wrk_dealloc( jpi, jpj,   fnln2d   )
-         ENDIF
-         IF ( med_diag%PN_FELIM%dgsave ) THEN
-            CALL iom_put( "PN_FELIM"  , ffln2d )
-            CALL wrk_dealloc( jpi, jpj,   ffln2d   )
-         ENDIF
-         IF ( med_diag%PD_JLIM%dgsave ) THEN
-            CALL iom_put( "PD_JLIM"  , fjld2d )
-            CALL wrk_dealloc( jpi, jpj,  fjld2d    )
-         ENDIF
-         IF ( med_diag%PD_NLIM%dgsave ) THEN
-            CALL iom_put( "PD_NLIM"  , fnld2d )
-            CALL wrk_dealloc( jpi, jpj,   fnld2d  )
-         ENDIF
-         IF ( med_diag%PD_FELIM%dgsave ) THEN
-            CALL iom_put( "PD_FELIM"  , ffld2d )
-            CALL wrk_dealloc( jpi, jpj,  ffld2d    )
-         ENDIF
-         IF ( med_diag%PD_SILIM%dgsave ) THEN
-            CALL iom_put( "PD_SILIM"  , fsld2d2 )
-            CALL wrk_dealloc( jpi, jpj,   fsld2d2   )
-         ENDIF
-         IF ( med_diag%PDSILIM2%dgsave ) THEN
-            CALL iom_put( "PDSILIM2"  , fsld2d )
-            CALL wrk_dealloc( jpi, jpj,   fsld2d   )
-         ENDIF
-         IF ( med_diag%INTFLX_N%dgsave ) THEN
-            CALL iom_put( "INTFLX_N"  , fflx_n )
-         ENDIF
-         IF ( med_diag%INTFLX_SI%dgsave ) THEN
-            CALL iom_put( "INTFLX_SI"  , fflx_si )
-         ENDIF
-         IF ( med_diag%INTFLX_FE%dgsave ) THEN
-            CALL iom_put( "INTFLX_FE"  , fflx_fe )
-         ENDIF        
-         IF ( med_diag%INT_PN%dgsave ) THEN
-            CALL iom_put( "INT_PN"  , ftot_pn )
-         ENDIF
-         IF ( med_diag%INT_PD%dgsave ) THEN
-            CALL iom_put( "INT_PD"  , ftot_pd )
-         ENDIF         
-         IF ( med_diag%INT_ZMI%dgsave ) THEN
-            CALL iom_put( "INT_ZMI"  , ftot_zmi )
-         ENDIF
-         IF ( med_diag%INT_ZME%dgsave ) THEN
-            CALL iom_put( "INT_ZME"  , ftot_zme )
-         ENDIF
-         IF ( med_diag%INT_DET%dgsave ) THEN
-            CALL iom_put( "INT_DET"  , ftot_det )
-         ENDIF
-         IF ( med_diag%INT_DTC%dgsave ) THEN
-            CALL iom_put( "INT_DTC"  , ftot_dtc )
-         ENDIF
-         IF ( med_diag%BEN_N%dgsave ) THEN
-            CALL iom_put( "BEN_N"  , za_sed_n )
-         ENDIF
-         IF ( med_diag%BEN_FE%dgsave ) THEN
-            CALL iom_put( "BEN_FE"  , za_sed_fe )
-         ENDIF
-         IF ( med_diag%BEN_C%dgsave ) THEN
-            CALL iom_put( "BEN_C"  , za_sed_c )
-         ENDIF
-         IF ( med_diag%BEN_SI%dgsave ) THEN
-            CALL iom_put( "BEN_SI"  , za_sed_si )
-         ENDIF
-         IF ( med_diag%BEN_CA%dgsave ) THEN
-            CALL iom_put( "BEN_CA"  , za_sed_ca )
-         ENDIF
-         IF ( med_diag%RUNOFF%dgsave ) THEN
-            CALL iom_put( "RUNOFF"  , f_runoff )
-         ENDIF 
-# if defined key_roam        
-         IF ( med_diag%N_PROD%dgsave ) THEN
-            CALL iom_put( "N_PROD"  , fnit_prod )
-         ENDIF
-         IF ( med_diag%N_CONS%dgsave ) THEN
-            CALL iom_put( "N_CONS"  , fnit_cons )
-         ENDIF
-         IF ( med_diag%C_PROD%dgsave ) THEN
-            CALL iom_put( "C_PROD"  , fcar_prod )
-         ENDIF
-         IF ( med_diag%C_CONS%dgsave ) THEN
-            CALL iom_put( "C_CONS"  , fcar_cons )
-         ENDIF
-         IF ( med_diag%O2_PROD%dgsave ) THEN
-            CALL iom_put( "O2_PROD"  , foxy_prod )
-         ENDIF
-         IF ( med_diag%O2_CONS%dgsave ) THEN
-            CALL iom_put( "O2_CONS"  , foxy_cons )
-         ENDIF
-         IF ( med_diag%O2_ANOX%dgsave ) THEN
-            CALL iom_put( "O2_ANOX"  , foxy_anox )
-         ENDIF
-         IF ( med_diag%INVTC%dgsave ) THEN
-            CALL iom_put( "INVTC"  , ftot_c )
-         ENDIF
-         IF ( med_diag%INVTALK%dgsave ) THEN
-            CALL iom_put( "INVTALK"  , ftot_a )
-         ENDIF
-         IF ( med_diag%INVTO2%dgsave ) THEN
-            CALL iom_put( "INVTO2"  , ftot_o2 )
-         ENDIF
-         IF ( med_diag%COM_RESP%dgsave ) THEN
-            CALL iom_put( "COM_RESP"  , fcomm_resp )
-         ENDIF         
-# endif      
-         !!
-         !! diagnostic filled in the i-j-k main loop
-         !!--------------------------------------------
-         IF ( med_diag%PRN%dgsave ) THEN
-            CALL iom_put( "PRN"  , fprn2d )
-            CALL wrk_dealloc( jpi, jpj,   fprn2d   )
-         ENDIF
-         IF ( med_diag%MPN%dgsave ) THEN
-            CALL iom_put( "MPN"  ,fdpn2d  )
-            CALL wrk_dealloc( jpi, jpj,    fdpn2d  )
-         ENDIF
-         IF ( med_diag%PRD%dgsave ) THEN
-            CALL iom_put( "PRD"  ,fprd2d  )
-            CALL wrk_dealloc( jpi, jpj,   fprd2d  )
-         ENDIF
-         IF( med_diag%MPD%dgsave ) THEN
-            CALL iom_put( "MPD"  , fdpd2d )
-            CALL wrk_dealloc( jpi, jpj,    fdpd2d )
-         ENDIF
-         !  IF( med_diag%DSED%dgsave ) THEN
-         !      CALL iom_put( "DSED"  , ftot_n )
-         !  ENDIF
-         IF( med_diag%OPAL%dgsave ) THEN
-            CALL iom_put( "OPAL"  , fprds2d )
-            CALL wrk_dealloc( jpi, jpj,   fprds2d  )
-         ENDIF
-         IF( med_diag%OPALDISS%dgsave ) THEN
-            CALL iom_put( "OPALDISS"  , fsdiss2d )
-            CALL wrk_dealloc( jpi, jpj,   fsdiss2d  )
-         ENDIF
-         IF( med_diag%GMIPn%dgsave ) THEN
-            CALL iom_put( "GMIPn"  , fgmipn2d )
-            CALL wrk_dealloc( jpi, jpj,   fgmipn2d  )
-         ENDIF
-         IF( med_diag%GMID%dgsave ) THEN
-            CALL iom_put( "GMID"  , fgmid2d )
-            CALL wrk_dealloc( jpi, jpj,  fgmid2d  )
-         ENDIF
-         IF( med_diag%MZMI%dgsave ) THEN
-            CALL iom_put( "MZMI"  , fdzmi2d )
-            CALL wrk_dealloc( jpi, jpj,   fdzmi2d   )
-         ENDIF
-         IF( med_diag%GMEPN%dgsave ) THEN
-            CALL iom_put( "GMEPN"  , fgmepn2d )
-            CALL wrk_dealloc( jpi, jpj,   fgmepn2d  )
-         ENDIF
-         IF( med_diag%GMEPD%dgsave ) THEN
-            CALL iom_put( "GMEPD"  , fgmepd2d )
-            CALL wrk_dealloc( jpi, jpj,   fgmepd2d   )
-         ENDIF
-         IF( med_diag%GMEZMI%dgsave ) THEN
-            CALL iom_put( "GMEZMI"  , fgmezmi2d )
-            CALL wrk_dealloc( jpi, jpj,   fgmezmi2d   )
-         ENDIF
-         IF( med_diag%GMED%dgsave ) THEN
-            CALL iom_put( "GMED"  , fgmed2d )
-            CALL wrk_dealloc( jpi, jpj,    fgmed2d  )
-         ENDIF
-         IF( med_diag%MZME%dgsave ) THEN
-            CALL iom_put( "MZME"  , fdzme2d )
-            CALL wrk_dealloc( jpi, jpj,  fdzme2d    )
-         ENDIF
-         !  IF( med_diag%DEXP%dgsave ) THEN
-         !      CALL iom_put( "DEXP"  , ftot_n )
-         !  ENDIF
-         IF( med_diag%DETN%dgsave ) THEN
-            CALL iom_put( "DETN"  , fslown2d )
-            CALL wrk_dealloc( jpi, jpj,  fslown2d    )
-         ENDIF
-         IF( med_diag%MDET%dgsave ) THEN
-            CALL iom_put( "MDET"  , fdd2d )
-            CALL wrk_dealloc( jpi, jpj,   fdd2d   )
-         ENDIF
-         IF( med_diag%AEOLIAN%dgsave ) THEN
-            CALL iom_put( "AEOLIAN"  , ffetop2d )
-            CALL wrk_dealloc( jpi, jpj,   ffetop2d   )
-         ENDIF
-         IF( med_diag%BENTHIC%dgsave ) THEN
-            CALL iom_put( "BENTHIC"  , ffebot2d )
-            CALL wrk_dealloc( jpi, jpj,   ffebot2d   )
-         ENDIF
-         IF( med_diag%SCAVENGE%dgsave ) THEN
-            CALL iom_put( "SCAVENGE"  , ffescav2d )
-            CALL wrk_dealloc( jpi, jpj,   ffescav2d  )
-         ENDIF
-         !! 
-         IF( med_diag%TOTREG_N%dgsave ) THEN
-            CALL iom_put( "TOTREG_N"  , fregen2d )
-            CALL wrk_dealloc( jpi, jpj,   fregen2d   )
-         ENDIF
-         IF( med_diag%TOTRG_SI%dgsave ) THEN
-            CALL iom_put( "TOTRG_SI"  , fregensi2d )
-            CALL wrk_dealloc( jpi, jpj,    fregensi2d  )
-         ENDIF
-         !! 
-         IF( med_diag%FASTN%dgsave ) THEN
-            CALL iom_put( "FASTN"  , ftempn2d )
-            CALL wrk_dealloc( jpi, jpj,   ftempn2d   )
-         ENDIF
-         IF( med_diag%FASTSI%dgsave ) THEN
-            CALL iom_put( "FASTSI"  , ftempsi2d )
-            CALL wrk_dealloc( jpi, jpj,   ftempsi2d   )
-         ENDIF
-         IF( med_diag%FASTFE%dgsave ) THEN
-            CALL iom_put( "FASTFE"  , ftempfe2d )
-            CALL wrk_dealloc( jpi, jpj,    ftempfe2d  )
-         ENDIF
-         IF( med_diag%FASTC%dgsave ) THEN
-            CALL iom_put( "FASTC"  , ftempc2d )
-            CALL wrk_dealloc( jpi, jpj,  ftempc2d    )
-         ENDIF
-         IF( med_diag%FASTCA%dgsave ) THEN
-            CALL iom_put( "FASTCA"  , ftempca2d )
-            CALL wrk_dealloc( jpi, jpj,  ftempca2d   )
-         ENDIF
-         !! 
-         IF( med_diag%REMINN%dgsave ) THEN
-            CALL iom_put( "REMINN"  , freminn2d )
-            CALL wrk_dealloc( jpi, jpj,   freminn2d   )
-         ENDIF
-         IF( med_diag%REMINSI%dgsave ) THEN
-            CALL iom_put( "REMINSI"  , freminsi2d )
-            CALL wrk_dealloc( jpi, jpj,   freminsi2d   )
-         ENDIF
-         IF( med_diag%REMINFE%dgsave ) THEN
-            CALL iom_put( "REMINFE"  , freminfe2d )
-            CALL wrk_dealloc( jpi, jpj,  freminfe2d    )
-         ENDIF
-         IF( med_diag%REMINC%dgsave ) THEN
-            CALL iom_put( "REMINC"  , freminc2d )
-            CALL wrk_dealloc( jpi, jpj,    freminc2d  )
-         ENDIF
-         IF( med_diag%REMINCA%dgsave ) THEN
-            CALL iom_put( "REMINCA"  , freminca2d )
-            CALL wrk_dealloc( jpi, jpj,   freminca2d  )
-         ENDIF
-         IF( med_diag%SEAFLRN%dgsave ) THEN
-            CALL iom_put( "SEAFLRN"  , fsedn )
-         ENDIF
-         IF( med_diag%SEAFLRSI%dgsave ) THEN
-            CALL iom_put( "SEAFLRSI"  , fsedsi )
-         ENDIF
-         IF( med_diag%SEAFLRFE%dgsave ) THEN
-            CALL iom_put( "SEAFLRFE"  , fsedfe )
-         ENDIF
-         IF( med_diag%SEAFLRC%dgsave ) THEN
-            CALL iom_put( "SEAFLRC"  , fsedc )
-         ENDIF
-         IF( med_diag%SEAFLRCA%dgsave ) THEN
-            CALL iom_put( "SEAFLRCA"  , fsedca )
-         ENDIF
-         !!
-# if defined key_roam            
-         !!
-         IF( med_diag%RIV_N%dgsave ) THEN
-            CALL iom_put( "RIV_N"  , rivn2d )
-            CALL wrk_dealloc( jpi, jpj,    rivn2d  )
-         ENDIF
-         IF( med_diag%RIV_SI%dgsave ) THEN
-            CALL iom_put( "RIV_SI"  , rivsi2d )
-            CALL wrk_dealloc( jpi, jpj,   rivsi2d   )
-         ENDIF
-         IF( med_diag%RIV_C%dgsave ) THEN
-            CALL iom_put( "RIV_C"  , rivc2d )
-            CALL wrk_dealloc( jpi, jpj,    rivc2d  )
-         ENDIF
-         IF( med_diag%RIV_ALK%dgsave ) THEN
-            CALL iom_put( "RIV_ALK"  , rivalk2d )
-            CALL wrk_dealloc( jpi, jpj,  rivalk2d    )
-         ENDIF
-         IF( med_diag%DETC%dgsave ) THEN
-            CALL iom_put( "DETC"  , fslowc2d )
-            CALL wrk_dealloc( jpi, jpj,   fslowc2d   )
-         ENDIF
-         !!
-         IF( med_diag%PN_LLOSS%dgsave ) THEN
-            CALL iom_put( "PN_LLOSS"  , fdpn22d )
-            CALL wrk_dealloc( jpi, jpj,   fdpn22d   )
-         ENDIF
-         IF( med_diag%PD_LLOSS%dgsave ) THEN
-            CALL iom_put( "PD_LLOSS"  , fdpd22d )
-            CALL wrk_dealloc( jpi, jpj,   fdpd22d   )
-         ENDIF
-         IF( med_diag%ZI_LLOSS%dgsave ) THEN
-            CALL iom_put( "ZI_LLOSS"  , fdzmi22d )
-             CALL wrk_dealloc( jpi, jpj,    fdzmi22d  )
-          ENDIF
-          IF( med_diag%ZE_LLOSS%dgsave ) THEN
-             CALL iom_put( "ZE_LLOSS"  , fdzme22d )
-             CALL wrk_dealloc( jpi, jpj,   fdzme22d   )
-          ENDIF
-          IF( med_diag%ZI_MES_N%dgsave ) THEN
-             CALL iom_put( "ZI_MES_N"  , zimesn2d )
-             CALL wrk_dealloc( jpi, jpj,    zimesn2d  )
-          ENDIF
-          IF( med_diag%ZI_MES_D%dgsave ) THEN
-             CALL iom_put( "ZI_MES_D"  , zimesd2d )
-             CALL wrk_dealloc( jpi, jpj,    zimesd2d  )
-          ENDIF
-          IF( med_diag%ZI_MES_C%dgsave ) THEN
-             CALL iom_put( "ZI_MES_C"  , zimesc2d )
-             CALL wrk_dealloc( jpi, jpj,    zimesc2d  )
-          ENDIF
-          IF( med_diag%ZI_MESDC%dgsave ) THEN
-             CALL iom_put( "ZI_MESDC"  ,zimesdc2d  )
-             CALL wrk_dealloc( jpi, jpj,    zimesdc2d  )
-          ENDIF
-          IF( med_diag%ZI_EXCR%dgsave ) THEN
-             CALL iom_put( "ZI_EXCR"  , ziexcr2d )
-             CALL wrk_dealloc( jpi, jpj,    ziexcr2d )
-          ENDIF
-          IF( med_diag%ZI_RESP%dgsave ) THEN
-             CALL iom_put( "ZI_RESP"  , ziresp2d )
-             CALL wrk_dealloc( jpi, jpj,   ziresp2d   )
-          ENDIF
-          IF( med_diag%ZI_GROW%dgsave ) THEN
-             CALL iom_put( "ZI_GROW"  , zigrow2d )
-             CALL wrk_dealloc( jpi, jpj,   zigrow2d   )
-          ENDIF
-          IF( med_diag%ZE_MES_N%dgsave ) THEN
-             CALL iom_put( "ZE_MES_N"  , zemesn2d )
-             CALL wrk_dealloc( jpi, jpj,    zemesn2d  )
-          ENDIF
-          IF( med_diag%ZE_MES_D%dgsave ) THEN
-             CALL iom_put( "ZE_MES_D"  , zemesd2d )
-             CALL wrk_dealloc( jpi, jpj,    zemesd2d  )
-          ENDIF
-          IF( med_diag%ZE_MES_C%dgsave ) THEN
-             CALL iom_put( "ZE_MES_C"  , zemesc2d )
-             CALL wrk_dealloc( jpi, jpj,   zemesc2d   )
-          ENDIF
-          IF( med_diag%ZE_MESDC%dgsave ) THEN
-             CALL iom_put( "ZE_MESDC"  , zemesdc2d )
-             CALL wrk_dealloc( jpi, jpj,   zemesdc2d   )
-          ENDIF
-          IF( med_diag%ZE_EXCR%dgsave ) THEN
-             CALL iom_put( "ZE_EXCR"  , zeexcr2d )
-             CALL wrk_dealloc( jpi, jpj,   zeexcr2d   )
-          ENDIF
-          IF( med_diag%ZE_RESP%dgsave ) THEN
-             CALL iom_put( "ZE_RESP"  , zeresp2d )
-             CALL wrk_dealloc( jpi, jpj,    zeresp2d  )
-          ENDIF
-          IF( med_diag%ZE_GROW%dgsave ) THEN
-             CALL iom_put( "ZE_GROW"  , zegrow2d )
-             CALL wrk_dealloc( jpi, jpj,   zegrow2d   )
-          ENDIF
-          IF( med_diag%MDETC%dgsave ) THEN
-             CALL iom_put( "MDETC"  , mdetc2d )
-             CALL wrk_dealloc( jpi, jpj,   mdetc2d   )
-          ENDIF
-          IF( med_diag%GMIDC%dgsave ) THEN
-             CALL iom_put( "GMIDC"  , gmidc2d )
-             CALL wrk_dealloc( jpi, jpj,    gmidc2d  )
-          ENDIF
-          IF( med_diag%GMEDC%dgsave ) THEN
-             CALL iom_put( "GMEDC"  , gmedc2d )
-             CALL wrk_dealloc( jpi, jpj,    gmedc2d  )
-          ENDIF
-          IF( med_diag%IBEN_N%dgsave ) THEN
-             CALL iom_put( "IBEN_N"  , iben_n2d )
-             CALL wrk_dealloc( jpi, jpj,    iben_n2d  )
-          ENDIF
-          IF( med_diag%IBEN_FE%dgsave ) THEN
-             CALL iom_put( "IBEN_FE"  , iben_fe2d )
-             CALL wrk_dealloc( jpi, jpj,   iben_fe2d   )
-          ENDIF
-          IF( med_diag%IBEN_C%dgsave ) THEN
-             CALL iom_put( "IBEN_C"  , iben_c2d )
-             CALL wrk_dealloc( jpi, jpj,   iben_c2d   )
-          ENDIF
-          IF( med_diag%IBEN_SI%dgsave ) THEN
-             CALL iom_put( "IBEN_SI"  , iben_si2d )
-             CALL wrk_dealloc( jpi, jpj,   iben_si2d   )
-          ENDIF
-          IF( med_diag%IBEN_CA%dgsave ) THEN
-             CALL iom_put( "IBEN_CA"  , iben_ca2d )
-             CALL wrk_dealloc( jpi, jpj,   iben_ca2d   )
-          ENDIF
-          IF( med_diag%OBEN_N%dgsave ) THEN
-             CALL iom_put( "OBEN_N"  , oben_n2d )
-             CALL wrk_dealloc( jpi, jpj,    oben_n2d  )
-          ENDIF
-          IF( med_diag%OBEN_FE%dgsave ) THEN
-             CALL iom_put( "OBEN_FE"  , oben_fe2d )
-             CALL wrk_dealloc( jpi, jpj,    oben_fe2d  )
-          ENDIF
-          IF( med_diag%OBEN_C%dgsave ) THEN
-             CALL iom_put( "OBEN_C"  , oben_c2d )
-             CALL wrk_dealloc( jpi, jpj,    oben_c2d  )
-          ENDIF
-          IF( med_diag%OBEN_SI%dgsave ) THEN
-             CALL iom_put( "OBEN_SI"  , oben_si2d )
-             CALL wrk_dealloc( jpi, jpj,    oben_si2d  )
-          ENDIF
-          IF( med_diag%OBEN_CA%dgsave ) THEN
-             CALL iom_put( "OBEN_CA"  , oben_ca2d )
-             CALL wrk_dealloc( jpi, jpj, oben_ca2d     )
-          ENDIF
-          IF( med_diag%SFR_OCAL%dgsave ) THEN
-             CALL iom_put( "SFR_OCAL"  , sfr_ocal2d )
-             CALL wrk_dealloc( jpi, jpj,    sfr_ocal2d  )
-          ENDIF
-          IF( med_diag%SFR_OARG%dgsave ) THEN
-             CALL iom_put( "SFR_OARG"  , sfr_oarg2d )
-             CALL wrk_dealloc( jpi, jpj,    sfr_oarg2d  )
-          ENDIF
-          IF( med_diag%LYSO_CA%dgsave ) THEN
-             CALL iom_put( "LYSO_CA"  , lyso_ca2d )
-             CALL wrk_dealloc( jpi, jpj,    lyso_ca2d  )
-          ENDIF
-# endif                   
-          !!
-          !! ** 3D diagnostics
-          IF( med_diag%TPP3%dgsave ) THEN
-             CALL iom_put( "TPP3"  , tpp3d )
-             CALL wrk_dealloc( jpi, jpj, jpk,   tpp3d  )
-          ENDIF
-          IF( med_diag%DETFLUX3%dgsave ) THEN
-             CALL iom_put( "DETFLUX3"  , detflux3d )
-             CALL wrk_dealloc( jpi, jpj, jpk,    detflux3d )
-          ENDIF
-          IF( med_diag%REMIN3N%dgsave ) THEN
-             CALL iom_put( "REMIN3N"  , remin3dn )
-             CALL wrk_dealloc( jpi, jpj, jpk,   remin3dn  )
-          ENDIF
-# if defined key_roam          
-          IF( med_diag%PH3%dgsave ) THEN
-             CALL iom_put( "PH3"  , f3_pH )
-          ENDIF
-          IF( med_diag%OM_CAL3%dgsave ) THEN
-             CALL iom_put( "OM_CAL3"  , f3_omcal )
-          ENDIF
-          !!
-          !! AXY (09/11/16): 2D CMIP6 diagnostics
-          IF( med_diag%INTDISSIC%dgsave ) THEN
-             CALL iom_put( "INTDISSIC"  , intdissic )
-             CALL wrk_dealloc( jpi, jpj, intdissic   )
-          ENDIF          
-          IF( med_diag%INTDISSIN%dgsave ) THEN
-             CALL iom_put( "INTDISSIN"  , intdissin )
-             CALL wrk_dealloc( jpi, jpj, intdissin   )
-          ENDIF          
-          IF( med_diag%INTDISSISI%dgsave ) THEN
-             CALL iom_put( "INTDISSISI"  , intdissisi )
-             CALL wrk_dealloc( jpi, jpj, intdissisi  )
-          ENDIF          
-          IF( med_diag%INTTALK%dgsave ) THEN
-             CALL iom_put( "INTTALK"  , inttalk )
-             CALL wrk_dealloc( jpi, jpj, inttalk     )
-          ENDIF          
-          IF( med_diag%O2min%dgsave ) THEN
-             CALL iom_put( "O2min"  , o2min )
-             CALL wrk_dealloc( jpi, jpj, o2min       )
-          ENDIF          
-          IF( med_diag%ZO2min%dgsave ) THEN
-             CALL iom_put( "ZO2min"  , zo2min )
-             CALL wrk_dealloc( jpi, jpj, zo2min      )
-          ENDIF          
-          IF( med_diag%FBDDTALK%dgsave ) THEN
-             CALL iom_put( "FBDDTALK"  , fbddtalk   )
-             CALL wrk_dealloc( jpi, jpj, fbddtalk   )
-          ENDIF          
-          IF( med_diag%FBDDTDIC%dgsave ) THEN
-             CALL iom_put( "FBDDTDIC"  , fbddtdic   )
-             CALL wrk_dealloc( jpi, jpj, fbddtdic   )
-          ENDIF          
-          IF( med_diag%FBDDTDIFE%dgsave ) THEN
-             CALL iom_put( "FBDDTDIFE" , fbddtdife  )
-             CALL wrk_dealloc( jpi, jpj, fbddtdife  )
-          ENDIF          
-          IF( med_diag%FBDDTDIN%dgsave ) THEN
-             CALL iom_put( "FBDDTDIN"  , fbddtdin   )
-             CALL wrk_dealloc( jpi, jpj, fbddtdin   )
-          ENDIF          
-          IF( med_diag%FBDDTDISI%dgsave ) THEN
-             CALL iom_put( "FBDDTDISI" , fbddtdisi  )
-             CALL wrk_dealloc( jpi, jpj, fbddtdisi  )
-          ENDIF    
-          !!
-          !! AXY (09/11/16): 3D CMIP6 diagnostics
-          IF( med_diag%TPPD3%dgsave ) THEN
-             CALL iom_put( "TPPD3"     , tppd3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, tppd3      )
-          ENDIF          
-          IF( med_diag%BDDTALK3%dgsave ) THEN
-             CALL iom_put( "BDDTALK3"  , bddtalk3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, bddtalk3   )
-          ENDIF          
-          IF( med_diag%BDDTDIC3%dgsave ) THEN
-             CALL iom_put( "BDDTDIC3"  , bddtdic3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, bddtdic3   )
-          ENDIF          
-          IF( med_diag%BDDTDIFE3%dgsave ) THEN
-             CALL iom_put( "BDDTDIFE3" , bddtdife3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, bddtdife3  )
-          ENDIF          
-          IF( med_diag%BDDTDIN3%dgsave ) THEN
-             CALL iom_put( "BDDTDIN3"  , bddtdin3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, bddtdin3   )
-          ENDIF          
-          IF( med_diag%BDDTDISI3%dgsave ) THEN
-             CALL iom_put( "BDDTDISI3" , bddtdisi3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, bddtdisi3  )
-          ENDIF    
-          IF( med_diag%FD_NIT3%dgsave ) THEN
-             CALL iom_put( "FD_NIT3"  , fd_nit3 )
-             CALL wrk_dealloc( jpi, jpj, jpk,   fd_nit3  )
-          ENDIF
-          IF( med_diag%FD_SIL3%dgsave ) THEN
-             CALL iom_put( "FD_SIL3"  , fd_sil3 )
-             CALL wrk_dealloc( jpi, jpj, jpk,   fd_sil3  )
-          ENDIF
-          IF( med_diag%FD_CAL3%dgsave ) THEN
-             CALL iom_put( "FD_CAL3"  , fd_cal3 )
-             CALL wrk_dealloc( jpi, jpj, jpk,   fd_cal3  )
-          ENDIF
-          IF( med_diag%FD_CAR3%dgsave ) THEN
-             CALL iom_put( "FD_CAR3"  , fd_car3 )
-             CALL wrk_dealloc( jpi, jpj, jpk,   fd_car3  )
-          ENDIF
-          IF( med_diag%CO33%dgsave ) THEN
-             CALL iom_put( "CO33"  , f3_co3 )
-          ENDIF          	        
-          IF( med_diag%CO3SATARAG3%dgsave ) THEN
-             CALL iom_put( "CO3SATARAG3"  , f3_omarg )
-          ENDIF          	        
-          IF( med_diag%CO3SATCALC3%dgsave ) THEN
-             CALL iom_put( "CO3SATCALC3"  , f3_omcal )
-          ENDIF          	        
-          IF( med_diag%EXPC3%dgsave ) THEN
-             CALL iom_put( "EXPC3"  , expc3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, expc3  )
-          ENDIF          	        
-          IF( med_diag%EXPN3%dgsave ) THEN
-             CALL iom_put( "EXPN3"  , expn3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, expn3  )
-          ENDIF          	        
-          IF( med_diag%DCALC3%dgsave ) THEN
-             CALL iom_put( "DCALC3"  , dcalc3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, dcalc3  )
-          ENDIF          	        
-          IF( med_diag%FEDISS3%dgsave ) THEN
-             CALL iom_put( "FEDISS3"  , fediss3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, fediss3  )
-          ENDIF          	        
-          IF( med_diag%FESCAV3%dgsave ) THEN
-             CALL iom_put( "FESCAV3"  , fescav3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, fescav3  )
-          ENDIF          	        
-          IF( med_diag%MIGRAZP3%dgsave ) THEN
-             CALL iom_put( "MIGRAZP3"  , migrazp3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, migrazp3  )
-          ENDIF          	        
-          IF( med_diag%MIGRAZD3%dgsave ) THEN
-             CALL iom_put( "MIGRAZD3"  , migrazd3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, migrazd3  )
-          ENDIF          	        
-          IF( med_diag%MEGRAZP3%dgsave ) THEN
-             CALL iom_put( "MEGRAZP3"  , megrazp3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, megrazp3  )
-          ENDIF          	        
-          IF( med_diag%MEGRAZD3%dgsave ) THEN
-             CALL iom_put( "MEGRAZD3"  , megrazd3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, megrazd3  )
-          ENDIF          	        
-          IF( med_diag%MEGRAZZ3%dgsave ) THEN
-             CALL iom_put( "MEGRAZZ3"  , megrazz3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, megrazz3  )
-          ENDIF          	        
-          IF( med_diag%O2SAT3%dgsave ) THEN
-             CALL iom_put( "O2SAT3"  , o2sat3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, o2sat3 )
-          ENDIF          	        
-          IF( med_diag%PBSI3%dgsave ) THEN
-             CALL iom_put( "PBSI3"  , pbsi3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pbsi3  )
-          ENDIF          	        
-          IF( med_diag%PCAL3%dgsave ) THEN
-             CALL iom_put( "PCAL3"  , pcal3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pcal3  )
-          ENDIF          	        
-          IF( med_diag%REMOC3%dgsave ) THEN
-             CALL iom_put( "REMOC3"  , remoc3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, remoc3 )
-          ENDIF          	        
-          IF( med_diag%PNLIMJ3%dgsave ) THEN
-             CALL iom_put( "PNLIMJ3" , pnlimj3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pnlimj3  )
-          ENDIF          	        
-          IF( med_diag%PNLIMN3%dgsave ) THEN
-             CALL iom_put( "PNLIMN3" , pnlimn3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pnlimn3  )
-          ENDIF          	        
-          IF( med_diag%PNLIMFE3%dgsave ) THEN
-             CALL iom_put( "PNLIMFE3" , pnlimfe3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pnlimfe3 )
-          ENDIF          	        
-          IF( med_diag%PDLIMJ3%dgsave ) THEN
-             CALL iom_put( "PDLIMJ3" , pdlimj3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pdlimj3  )
-          ENDIF          	        
-          IF( med_diag%PDLIMN3%dgsave ) THEN
-             CALL iom_put( "PDLIMN3" , pdlimn3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pdlimn3  )
-          ENDIF          	        
-          IF( med_diag%PDLIMFE3%dgsave ) THEN
-             CALL iom_put( "PDLIMFE3" , pdlimfe3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pdlimfe3 )
-          ENDIF          	        
-          IF( med_diag%PDLIMSI3%dgsave ) THEN
-             CALL iom_put( "PDLIMSI3" , pdlimsi3 )
-             CALL wrk_dealloc( jpi, jpj, jpk, pdlimsi3 )
-          ENDIF          	        
-          
-# endif         
-
-          CALL wrk_dealloc( jpi, jpj,   zw2d   )
-
-       ENDIF                    ! end of ln_diatrc option
+      CALL bio_medusa_fin( kt )
 
 # if defined key_trc_diabio
