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 7886)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7912)
@@ -95,9 +95,11 @@
       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
 
       IMPLICIT NONE
       PRIVATE
       
-      PUBLIC   trc_bio_medusa    ! called in ???
+      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
 
    !!* Substitution
@@ -187,5 +189,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::    fun_Q10
       !! AXY (01/03/10): add in mixed layer PP diagnostics
-      REAL(wp), DIMENSION(jpi,jpj) ::    fprn_ml,fprd_ml
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fprn_ml,fprd_ml
       !!
       !! nutrient limiting factors
@@ -200,5 +202,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::    ffetop,ffebot,ffescav
       REAL(wp) ::    xLgF, xFeT, xFeF, xFeL         !! state variables for iron-ligand system
-      REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system
+!      REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system
       REAL(wp) ::    xb_coef_tmp, xb2M4ac           !! iron-ligand parameters
       REAL(wp) ::    xmaxFeF,fdeltaFe               !! max Fe' parameters
@@ -223,7 +225,7 @@
       !!
       !! mortality/Remineralisation (defunct parameter "fz" removed)
-      REAL(wp, DIMENSION(jpi,jpj)) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
+      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
 # if defined key_roam
-      REAL(wp, DIMENSION(jpi,jpj)) ::    fddc
+      REAL(wp), DIMENSION(jpi,jpj) ::    fddc
 # endif
       REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
@@ -231,9 +233,9 @@
       REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux
       REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi
-      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi
 # if defined key_roam
-! Doesn't look like this is used - marc 10/4/17
-!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenc
-      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfastc
+!! Doesn't look like this is used - marc 10/4/17
+!!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenc
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfastc
 # endif
       !!
@@ -242,8 +244,8 @@
       REAL(WP), DIMENSION(jpi,jpj) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
       REAL(wp), DIMENSION(jpi,jpj) ::    freminn,freminsi,freminfe,freminc,freminca
-      REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca
       REAL(wp), DIMENSION(jpi,jpj) ::    fprotf
-      REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca
-      REAL(wp), DIMENSION(jpi,jpj) ::    fccd
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fccd
       REAL(wp), DIMENSION(jpi,jpj) ::    fccd_dep
       !!
@@ -254,5 +256,5 @@
 ! I don't think ffast2slowfe is used - marc 10/4/17
 !      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowfe,ffast2slowc
-      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slow
+      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowc
       !!
       !! conservation law
@@ -266,21 +268,21 @@
       !!
       !! water column nutrient and flux integrals
-      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_n,ftot_si,ftot_fe
-      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_n,fflx_si,fflx_fe
-      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_n,fifd_si,fifd_fe
-      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_n,fofd_si,fofd_fe
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_n,ftot_si,ftot_fe
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_n,fflx_si,fflx_fe
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_n,fifd_si,fifd_fe
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_n,fofd_si,fofd_fe
 # if defined key_roam
-      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_c,ftot_a,ftot_o2
-      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_c,fflx_a,fflx_o2
-      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_c,fifd_a,fifd_o2
-      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_c,fofd_a,fofd_o2
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_c,ftot_a,ftot_o2
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_c,fflx_a,fflx_o2
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_c,fifd_a,fifd_o2
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_c,fofd_a,fofd_o2
 # endif
       !!
       !! zooplankton grazing integrals
-      REAL(wp), DIMENSION(jpi,jpj) ::    fzmi_i,fzmi_o,fzme_i,fzme_o
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fzmi_i,fzmi_o,fzme_i,fzme_o
       !!
       !! limitation term temporary variables
-      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_pn,ftot_pd
-      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_zmi,ftot_zme,ftot_det,ftot_dtc
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_pn,ftot_pd
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_zmi,ftot_zme,ftot_det,ftot_dtc
       !! use ballast scheme (1) or simple exponential scheme (0; a conservation test)
       INTEGER  ::    iball
@@ -293,5 +295,5 @@
       !! nitrogen and silicon production and consumption
       REAL(wp) ::    fn_prod, fn_cons, fs_prod, fs_cons
-      REAL(wp), DIMENSION(jpi,jpj) ::    fnit_prod, fnit_cons, fsil_prod, fsil_cons
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fnit_prod, fnit_cons, fsil_prod, fsil_cons
 # if defined key_roam
       !!
@@ -304,6 +306,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::    f_TDIC, f_TALK, f_dcf, f_henry
       REAL(wp), DIMENSION(jpi,jpj) ::    f_pp0
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat, f_o2sat3
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg
+      REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat
+      REAL(wp)                     ::    f_o2sat3
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg
       !!
       !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
@@ -323,10 +326,10 @@
       !! carbon, alkalinity production and consumption
       REAL(wp) ::    fc_prod, fc_cons, fa_prod, fa_cons
-      REAL(wp), DIMENSION(jpi,jpj) ::    fcomm_resp
-      REAL(wp), DIMENSION(jpi,jpj) ::    fcar_prod, fcar_cons
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fcomm_resp
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fcar_prod, fcar_cons
       !!
       !! oxygen production and consumption (and non-consumption)
       REAL(wp), DIMENSION(jpi,jpj) ::    fo2_prod, fo2_cons, fo2_ncons, fo2_ccons
-      REAL(wp), DIMENSION(jpi,jpj) ::    foxy_prod, foxy_cons, foxy_anox
+!      REAL(wp), DIMENSION(jpi,jpj) ::    foxy_prod, foxy_cons, foxy_anox
       !! Jpalm (11-08-2014)
       !! add DMS in MEDUSA for UKESM1 model
@@ -338,65 +341,17 @@
       !! 
       !! benthic fluxes
-      INTEGER  ::    ibenthic
-      REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe,              f_sbenin_c
-      REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca
-      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca
+!      INTEGER  ::    ibenthic
+!      REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe,              f_sbenin_c
+!      REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca
+!      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca
       REAL(wp) ::    zfact
       !! 
       !! benthic fluxes of CaCO3 that shouldn't happen because of lysocline
-      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca
+!      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca
       !!
       !! riverine fluxes
-      REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk
+!      REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk
       !! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface
       REAL(wp), DIMENSION(jpi,jpj) :: f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk
-      !!
-      !! Jpalm -- 11-10-2015 -- adapt diag to iom_use
-      !! 2D var for diagnostics.
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fprn2d, fdpn2d, fprd2d, fdpd2d, fprds2d, fsdiss2d, fgmipn2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fgmid2d, fdzmi2d, fgmepn2d, fgmepd2d, fgmezmi2d, fgmed2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fdzme2d, fslown2d, fdd2d, ffetop2d, ffebot2d, ffescav2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fjln2d, fnln2d, ffln2d, fjld2d, fnld2d, ffld2d, fsld2d2
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fsld2d, fregen2d, fregensi2d, ftempn2d, ftempsi2d, ftempfe2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: ftempc2d, ftempca2d, freminn2d, freminsi2d, freminfe2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: freminc2d, freminca2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
-# if defined key_roam
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: ffastca2d, rivn2d, rivsi2d, rivc2d, rivalk2d, fslowc2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fdpn22d, fdpd22d, fdzmi22d, fdzme22d, zimesn2d, zimesd2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zimesc2d, zimesdc2d, ziexcr2d, ziresp2d, zigrow2d, zemesn2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zemesd2d, zemesc2d, zemesdc2d, zeexcr2d, zeresp2d, zegrow2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: mdetc2d, gmidc2d, gmedc2d, f_pco2a2d, f_pco2w2d, f_co2flux2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_TDIC2d, f_TALK2d, f_kw6602d, f_pp02d, f_o2flux2d, f_o2sat2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: dms_andr2d, dms_simo2d, dms_aran2d, dms_hall2d, dms_surf2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: iben_n2d, iben_fe2d, iben_c2d, iben_si2d, iben_ca2d, oben_n2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: oben_fe2d, oben_c2d, oben_si2d, oben_ca2d, sfr_ocal2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: sfr_oarg2d, lyso_ca2d 
-      !! AXY (23/11/16): extra MOCSY diagnostics
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_xco2a_2d, f_fco2w_2d, f_fco2a_2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_ocnrhosw_2d, f_ocnschco2_2d, f_ocnkwco2_2d
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_ocnk0_2d, f_co2starair_2d, f_ocndpco2_2d
-# endif
-      !!
-      !! 3D var for diagnostics.
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: tpp3d, detflux3d, remin3dn
-      !!
-# if defined key_roam
-      !! AXY (04/11/16)
-      !! 2D var for new CMIP6 diagnostics (behind a key_roam ifdef for simplicity)
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fgco2, intdissic, intdissin, intdissisi, inttalk, o2min, zo2min
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: fbddtalk, fbddtdic, fbddtdife, fbddtdin, fbddtdisi
-      !!
-      !! 3D var for new CMIP6 diagnostics
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: tppd3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: bddtalk3, bddtdic3, bddtdife3, bddtdin3, bddtdisi3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: fd_nit3, fd_sil3, fd_car3, fd_cal3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: co33, co3satarag3, co3satcalc3, dcalc3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: expc3, expn3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: fediss3, fescav3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: migrazp3, migrazd3, megrazp3, megrazd3, megrazz3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: o2sat3, pbsi3, pcal3, remoc3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: pnlimj3, pnlimn3, pnlimfe3, pdlimj3, pdlimn3, pdlimfe3, pdlimsi3
-# endif
       !!---------------------------------------------------------------------
 
@@ -468,769 +423,15 @@
 
       !!----------------------------------------------------------------------
-      !! blank fast-sinking detritus 2D fields
+      !! Initialise arrays to zero and set up arrays for diagnostics
       !!----------------------------------------------------------------------
-      !!
-      ffastn(:,:)  = 0.0        !! organic nitrogen
-      ffastsi(:,:) = 0.0        !! biogenic silicon
-      ffastfe(:,:) = 0.0        !! organic iron
-      ffastc(:,:)  = 0.0        !! organic carbon
-      ffastca(:,:) = 0.0        !! biogenic calcium carbonate
-      !!
-      fsedn(:,:)   = 0.0        !! Seafloor flux of N 
-      fsedsi(:,:)  = 0.0        !! Seafloor flux of Si
-      fsedfe(:,:)  = 0.0        !! Seafloor flux of Fe
-      fsedc(:,:)   = 0.0        !! Seafloor flux of C
-      fsedca(:,:)  = 0.0        !! Seafloor flux of CaCO3
-      !!
-      fregenfast(:,:)   = 0.0   !! integrated  N regeneration (fast detritus)
-      fregenfastsi(:,:) = 0.0   !! integrated Si regeneration (fast detritus)
-# if defined key_roam
-      fregenfastc(:,:)  = 0.0   !! integrated  C regeneration (fast detritus)
-# endif
-      !!
-      fccd(:,:)    = 0.0        !! last depth level before CCD
-
-      !!----------------------------------------------------------------------
-      !! blank nutrient/flux inventories
-      !!----------------------------------------------------------------------
-      !!
-      fflx_n(:,:)  = 0.0        !! nitrogen flux total
-      fflx_si(:,:) = 0.0        !! silicon  flux total
-      fflx_fe(:,:) = 0.0        !! iron     flux total
-      fifd_n(:,:)  = 0.0        !! nitrogen fast detritus production
-      fifd_si(:,:) = 0.0        !! silicon  fast detritus production
-      fifd_fe(:,:) = 0.0        !! iron     fast detritus production
-      fofd_n(:,:)  = 0.0        !! nitrogen fast detritus remineralisation
-      fofd_si(:,:) = 0.0        !! silicon  fast detritus remineralisation
-      fofd_fe(:,:) = 0.0        !! iron     fast detritus remineralisation
-# if defined key_roam
-      fflx_c(:,:)  = 0.0        !! carbon     flux total
-      fflx_a(:,:)  = 0.0        !! alkalinity flux total
-      fflx_o2(:,:) = 0.0        !! oxygen     flux total
-      ftot_c(:,:)  = 0.0        !! carbon     inventory
-      ftot_a(:,:)  = 0.0        !! alkalinity inventory
-      ftot_o2(:,:) = 0.0        !! oxygen     inventory
-      fifd_c(:,:)  = 0.0        !! carbon     fast detritus production
-      fifd_a(:,:)  = 0.0        !! alkalinity fast detritus production
-      fifd_o2(:,:) = 0.0        !! oxygen     fast detritus production
-      fofd_c(:,:)  = 0.0        !! carbon     fast detritus remineralisation
-      fofd_a(:,:)  = 0.0        !! alkalinity fast detritus remineralisation
-      fofd_o2(:,:) = 0.0        !! oxygen     fast detritus remineralisation
-      !!
-      fnit_prod(:,:) = 0.0      !! (organic)   nitrogen production
-      fnit_cons(:,:) = 0.0      !! (organic)   nitrogen consumption
-      fsil_prod(:,:) = 0.0      !! (inorganic) silicon production
-      fsil_cons(:,:) = 0.0      !! (inorganic) silicon consumption
-      fcar_prod(:,:) = 0.0      !! (organic)   carbon production
-      fcar_cons(:,:) = 0.0      !! (organic)   carbon consumption
-      !!
-      foxy_prod(:,:) = 0.0      !! oxygen production
-      foxy_cons(:,:) = 0.0      !! oxygen consumption
-      foxy_anox(:,:) = 0.0      !! unrealised oxygen consumption
-      !!
-# endif
-      ftot_n(:,:)   = 0.0       !! N inventory 
-      ftot_si(:,:)  = 0.0       !! Si inventory
-      ftot_fe(:,:)  = 0.0       !! Fe inventory
-      ftot_pn(:,:)  = 0.0       !! integrated non-diatom phytoplankton
-      ftot_pd(:,:)  = 0.0       !! integrated diatom     phytoplankton
-      ftot_zmi(:,:) = 0.0       !! integrated microzooplankton
-      ftot_zme(:,:) = 0.0       !! integrated mesozooplankton
-      ftot_det(:,:) = 0.0       !! integrated slow detritus, nitrogen
-      ftot_dtc(:,:) = 0.0       !! integrated slow detritus, carbon
-      !!
-      fzmi_i(:,:)  = 0.0        !! material grazed by microzooplankton
-      fzmi_o(:,:)  = 0.0        !! ... sum of fate of this material
-      fzme_i(:,:)  = 0.0        !! material grazed by  mesozooplankton
-      fzme_o(:,:)  = 0.0        !! ... sum of fate of this material
-      !!
-      f_sbenin_n(:,:)  = 0.0    !! slow detritus N  -> benthic pool
-      f_sbenin_fe(:,:) = 0.0    !! slow detritus Fe -> benthic pool
-      f_sbenin_c(:,:)  = 0.0    !! slow detritus C  -> benthic pool
-      f_fbenin_n(:,:)  = 0.0    !! fast detritus N  -> benthic pool
-      f_fbenin_fe(:,:) = 0.0    !! fast detritus Fe -> benthic pool
-      f_fbenin_si(:,:) = 0.0    !! fast detritus Si -> benthic pool
-      f_fbenin_c(:,:)  = 0.0    !! fast detritus C  -> benthic pool
-      f_fbenin_ca(:,:) = 0.0    !! fast detritus Ca -> benthic pool
-      !!
-      f_benout_n(:,:)  = 0.0    !! benthic N  pool  -> dissolved
-      f_benout_fe(:,:) = 0.0    !! benthic Fe pool  -> dissolved
-      f_benout_si(:,:) = 0.0    !! benthic Si pool  -> dissolved
-      f_benout_c(:,:)  = 0.0    !! benthic C  pool  -> dissolved
-      f_benout_ca(:,:) = 0.0    !! benthic Ca pool  -> dissolved
-      !!
-      f_benout_lyso_ca(:,:) = 0.0 !! benthic Ca pool  -> dissolved (when it shouldn't!)
-      !!
-      f_runoff(:,:)  = 0.0      !! riverine runoff
-      f_riv_n(:,:)   = 0.0      !! riverine N   input 
-      f_riv_si(:,:)  = 0.0      !! riverine Si  input 
-      f_riv_c(:,:)   = 0.0      !! riverine C   input 
-      f_riv_alk(:,:) = 0.0      !! riverine alk input 
-      !! 
-      !! Jpalm -- 06-03-2017 -- Forgotten var to init
-      f_omarg(:,:) = 0.0        !!
-      f_omcal(:,:) = 0.0 
-      xFree(:,:) = 0.0          !! state variables for iron-ligand system
-      fcomm_resp(:,:) = 0.0 
-      fprn_ml(:,:) = 0.0        !! mixed layer PP diagnostics
-      fprd_ml(:,:) = 0.0        !! mixed layer PP diagnostics
-
-      !!
-      !! allocate and initiate 2D diag
-      !! -----------------------------
-      !! Juju :: add kt condition !!
-      IF ( lk_iomput .AND. .NOT.  ln_diatrc ) THEN 
-         !!
-         if ( kt == nittrc000 )   CALL trc_nam_iom_medusa !! initialise iom_use test
-         !!
-         CALL wrk_alloc( jpi, jpj,      zw2d )
-         zw2d(:,:)      = 0.0   !!
-         IF ( med_diag%PRN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fprn2d    )
-            fprn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF ( med_diag%MPN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fdpn2d    )
-            fdpn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF ( med_diag%PRD%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fprd2d    )
-            fprd2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%MPD%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fdpd2d    )
-            fdpd2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OPAL%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fprds2d    )
-            fprds2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OPALDISS%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fsdiss2d    )
-            fsdiss2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMIPn%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fgmipn2d    )
-            fgmipn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMID%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fgmid2d    )
-            fgmid2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%MZMI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fdzmi2d    )
-            fdzmi2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMEPN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fgmepn2d    )
-            fgmepn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMEPD%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fgmepd2d    )
-            fgmepd2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMEZMI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fgmezmi2d    )
-            fgmezmi2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMED%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fgmed2d    )
-            fgmed2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%MZME%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fdzme2d    )
-            fdzme2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%DETN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fslown2d    )
-            fslown2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%MDET%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fdd2d    )
-            fdd2d(:,:)      = 0.0 !!
-         ENDIF      
-         IF( med_diag%AEOLIAN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ffetop2d    )
-            ffetop2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%BENTHIC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    ffebot2d   )
-            ffebot2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%SCAVENGE%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ffescav2d    )
-            ffescav2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PN_JLIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fjln2d    )
-            fjln2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PN_NLIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fnln2d    )
-            fnln2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PN_FELIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ffln2d    )
-            ffln2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PD_JLIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fjld2d    )
-            fjld2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PD_NLIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fnld2d    )
-            fnld2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PD_FELIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ffld2d    )
-            ffld2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PD_SILIM%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fsld2d2    )
-            fsld2d2(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PDSILIM2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fsld2d    )
-            fsld2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip SDT_XXXX diagnostics here
-!!
-         IF( med_diag%TOTREG_N%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fregen2d    )
-            fregen2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%TOTRG_SI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fregensi2d    )
-            fregensi2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip REG_XXXX diagnostics here
-!!
-         IF( med_diag%FASTN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ftempn2d    )
-            ftempn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%FASTSI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ftempsi2d    )
-            ftempsi2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%FASTFE%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,  ftempfe2d     )
-            ftempfe2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%FASTC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,  ftempc2d     )
-            ftempc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%FASTCA%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   ftempca2d    )
-            ftempca2d(:,:)      = 0.0 !!
-         ENDIF     
-!!
-!! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here
-!!
-         IF( med_diag%REMINN%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    freminn2d   )
-            freminn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%REMINSI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    freminsi2d   )
-            freminsi2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%REMINFE%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    freminfe2d   )
-            freminfe2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%REMINC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   freminc2d    )
-            freminc2d(:,:)      = 0.0 !! 
-         ENDIF
-         IF( med_diag%REMINCA%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   freminca2d    )
-            freminca2d(:,:)      = 0.0 !!
-         ENDIF
-# if defined key_roam
-!!
-!! skip SEAFLRXX, MED_XXXX, INTFLX_XX, INT_XX, ML_XXX, OCAL_XXX, FE_XXXX, MED_XZE, WIND diagnostics here
-!!
-         IF( med_diag%RR_0100%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    ffastca2d   )
-            ffastca2d(:,:)      = 0.0 !!
-         ENDIF
-
-         IF( med_diag%ATM_PCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    f_pco2a2d   )
-            f_pco2a2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip OCN_PH diagnostic here
-!!
-         IF( med_diag%OCN_PCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    f_pco2w2d   )
-            f_pco2w2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here
-!!
-         IF( med_diag%CO2FLUX%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   f_co2flux2d    )
-            f_co2flux2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip OM_XXX diagnostics here
-!!
-         IF( med_diag%TCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   f_TDIC2d    )
-            f_TDIC2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%TALK%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    f_TALK2d   )
-            f_TALK2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%KW660%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    f_kw6602d   )
-            f_kw6602d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ATM_PP0%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    f_pp02d   )
-            f_pp02d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%O2FLUX%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   f_o2flux2d    )
-            f_o2flux2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%O2SAT%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    f_o2sat2d   )
-            f_o2sat2d(:,:)      = 0.0 !!
-         ENDIF 
-!!
-!! skip XXX_CCD diagnostics here
-!! 
-         IF( med_diag%SFR_OCAL%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    sfr_ocal2d  )
-            sfr_ocal2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%SFR_OARG%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    sfr_oarg2d  )
-            sfr_oarg2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip XX_PROD, XX_CONS, O2_ANOX, RR_XXXX diagnostics here
-!! 
-         IF( med_diag%IBEN_N%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    iben_n2d  )
-            iben_n2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%IBEN_FE%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   iben_fe2d   )
-            iben_fe2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%IBEN_C%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   iben_c2d   )
-            iben_c2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%IBEN_SI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   iben_si2d   )
-            iben_si2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%IBEN_CA%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   iben_ca2d   )
-            iben_ca2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OBEN_N%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    oben_n2d  )
-            oben_n2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OBEN_FE%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    oben_fe2d  )
-            oben_fe2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OBEN_C%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    oben_c2d  )
-            oben_c2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OBEN_SI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    oben_si2d  )
-            oben_si2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OBEN_CA%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    oben_ca2d  )
-            oben_ca2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip BEN_XX diagnostics here
-!!
-         IF( med_diag%RIV_N%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    rivn2d   )
-            rivn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%RIV_SI%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    rivsi2d   )
-            rivsi2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%RIV_C%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   rivc2d    )
-            rivc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%RIV_ALK%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    rivalk2d   )
-            rivalk2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%DETC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    fslowc2d   )
-            fslowc2d(:,:)      = 0.0 !!
-         ENDIF 
-!!
-!! skip SDC_XXXX, INVTXXX diagnostics here
-!!
-         IF( med_diag%LYSO_CA%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    lyso_ca2d  )
-            lyso_ca2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip COM_RESP diagnostic here
-!!
-         IF( med_diag%PN_LLOSS%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    fdpn22d   )
-            fdpn22d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%PD_LLOSS%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    fdpd22d   )
-            fdpd22d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_LLOSS%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    fdzmi22d   )
-            fdzmi22d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_LLOSS%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   fdzme22d    )
-            fdzme22d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_MES_N%dgsave ) THEN   
-            CALL wrk_alloc( jpi, jpj,   zimesn2d    )
-            zimesn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_MES_D%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zimesd2d   )
-            zimesd2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_MES_C%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zimesc2d   )
-            zimesc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_MESDC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zimesdc2d   )
-            zimesdc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_EXCR%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,     ziexcr2d  )
-            ziexcr2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_RESP%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    ziresp2d   )
-            ziresp2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZI_GROW%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zigrow2d   )
-            zigrow2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_MES_N%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   zemesn2d    )
-            zemesn2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_MES_D%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zemesd2d   )
-            zemesd2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_MES_C%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zemesc2d   )
-            zemesc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_MESDC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zemesdc2d   )
-            zemesdc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_EXCR%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zeexcr2d   )
-            zeexcr2d(:,:)      = 0.0 !!
-         ENDIF                  
-         IF( med_diag%ZE_RESP%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zeresp2d   )
-            zeresp2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ZE_GROW%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    zegrow2d   )
-            zegrow2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%MDETC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,   mdetc2d    )
-            mdetc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMIDC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    gmidc2d   )
-            gmidc2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%GMEDC%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj,    gmedc2d   )
-            gmedc2d(:,:)      = 0.0 !!
-         ENDIF
-!!
-!! skip INT_XXX diagnostics here
-!!
-         IF (jdms .eq. 1) THEN
-            IF( med_diag%DMS_SURF%dgsave ) THEN
-               CALL wrk_alloc( jpi, jpj,   dms_surf2d    )
-               dms_surf2d(:,:)      = 0.0 !!
-            ENDIF
-            IF( med_diag%DMS_ANDR%dgsave ) THEN
-               CALL wrk_alloc( jpi, jpj,   dms_andr2d    )
-               dms_andr2d(:,:)      = 0.0 !!
-            ENDIF
-            IF( med_diag%DMS_SIMO%dgsave ) THEN
-               CALL wrk_alloc( jpi, jpj,  dms_simo2d     )
-               dms_simo2d(:,:)      = 0.0 !!
-            ENDIF
-            IF( med_diag%DMS_ARAN%dgsave ) THEN
-               CALL wrk_alloc( jpi, jpj,   dms_aran2d    )
-               dms_aran2d(:,:)      = 0.0 !!
-            ENDIF
-            IF( med_diag%DMS_HALL%dgsave ) THEN
-               CALL wrk_alloc( jpi, jpj,   dms_hall2d    )
-               dms_hall2d(:,:)      = 0.0 !!
-            ENDIF
-         ENDIF   
-         !!
-         !! AXY (24/11/16): extra MOCSY diagnostics, 2D
-         IF( med_diag%ATM_XCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_xco2a_2d      )
-            f_xco2a_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OCN_FCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_fco2w_2d      )
-            f_fco2w_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%ATM_FCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_fco2a_2d      )
-            f_fco2a_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OCN_RHOSW%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d   )
-            f_ocnrhosw_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OCN_SCHCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d  )
-            f_ocnschco2_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OCN_KWCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d   )
-            f_ocnkwco2_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OCN_K0%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_ocnk0_2d      )
-            f_ocnk0_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%CO2STARAIR%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_co2starair_2d )
-            f_co2starair_2d(:,:)      = 0.0 !!
-         ENDIF
-         IF( med_diag%OCN_DPCO2%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d   )
-            f_ocndpco2_2d(:,:)      = 0.0 !!
-         ENDIF
-# endif  
-         IF( med_diag%TPP3%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, jpk,       tpp3d )
-            tpp3d(:,:,:)      = 0.0 !! 
-         ENDIF
-         IF( med_diag%DETFLUX3%dgsave ) THEN
-            CALL wrk_alloc( jpi, jpj, jpk,        detflux3d )
-            detflux3d(:,:,:)      = 0.0 !! 
-         ENDIF
-         IF( med_diag%REMIN3N%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk,        remin3dn )
-             remin3dn(:,:,:)      = 0.0 !! 
-          ENDIF
-          !! 
-          !! AXY (10/11/16): CMIP6 diagnostics, 2D
-          !! JPALM -- 17-11-16 -- put fgco2 alloc out of diag request
-          !!                   needed for coupling/passed through restart
-          !! IF( med_diag%FGCO2%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   fgco2    )
-             fgco2(:,:)      = 0.0 !!
-          !! ENDIF
-          IF( med_diag%INTDISSIC%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   intdissic    )
-             intdissic(:,:)  = 0.0 !!
-          ENDIF          
-          IF( med_diag%INTDISSIN%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   intdissin    )
-             intdissin(:,:)  = 0.0 !!
-          ENDIF          
-          IF( med_diag%INTDISSISI%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   intdissisi    )
-             intdissisi(:,:)  = 0.0 !!
-          ENDIF          
-          IF( med_diag%INTTALK%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   inttalk    )
-             inttalk(:,:)  = 0.0 !!
-          ENDIF          
-          IF( med_diag%O2min%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   o2min    )
-             o2min(:,:)  = 1.e3 !! set to high value as we're looking for min(o2)
-          ENDIF          
-          IF( med_diag%ZO2min%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj,   zo2min    )
-             zo2min(:,:)  = 0.0 !!
-          ENDIF          
-          IF( med_diag%FBDDTALK%dgsave  ) THEN
-             CALL wrk_alloc( jpi, jpj, fbddtalk  )
-             fbddtalk(:,:)  = 0.0 !! 
-          ENDIF
-          IF( med_diag%FBDDTDIC%dgsave  ) THEN
-             CALL wrk_alloc( jpi, jpj, fbddtdic  )
-             fbddtdic(:,:)  = 0.0 !! 
-          ENDIF
-          IF( med_diag%FBDDTDIFE%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj, fbddtdife )
-             fbddtdife(:,:) = 0.0 !! 
-          ENDIF
-          IF( med_diag%FBDDTDIN%dgsave  ) THEN
-             CALL wrk_alloc( jpi, jpj, fbddtdin  )
-             fbddtdin(:,:)  = 0.0 !! 
-          ENDIF
-          IF( med_diag%FBDDTDISI%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj, fbddtdisi )
-             fbddtdisi(:,:) = 0.0 !! 
-          ENDIF
-          !! 
-          !! AXY (10/11/16): CMIP6 diagnostics, 3D
-          IF( med_diag%TPPD3%dgsave     ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, tppd3     )
-             tppd3(:,:,:)     = 0.0 !! 
-          ENDIF
-          IF( med_diag%BDDTALK3%dgsave  ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, bddtalk3  )
-             bddtalk3(:,:,:)  = 0.0 !! 
-          ENDIF
-          IF( med_diag%BDDTDIC3%dgsave  ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, bddtdic3  )
-             bddtdic3(:,:,:)  = 0.0 !! 
-          ENDIF
-          IF( med_diag%BDDTDIFE3%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, bddtdife3 )
-             bddtdife3(:,:,:) = 0.0 !! 
-          ENDIF
-          IF( med_diag%BDDTDIN3%dgsave  ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, bddtdin3  )
-             bddtdin3(:,:,:)  = 0.0 !! 
-          ENDIF
-          IF( med_diag%BDDTDISI3%dgsave ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, bddtdisi3 )
-             bddtdisi3(:,:,:) = 0.0 !! 
-          ENDIF
-          IF( med_diag%FD_NIT3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, fd_nit3   )
-             fd_nit3(:,:,:)   = 0.0 !! 
-          ENDIF
-          IF( med_diag%FD_SIL3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, fd_sil3   )
-             fd_sil3(:,:,:)   = 0.0 !! 
-          ENDIF
-          IF( med_diag%FD_CAR3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, fd_car3   )
-             fd_car3(:,:,:)   = 0.0 !! 
-          ENDIF
-          IF( med_diag%FD_CAL3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, fd_cal3   )
-             fd_cal3(:,:,:)   = 0.0 !! 
-          ENDIF
-          IF( med_diag%DCALC3%dgsave    ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, dcalc3    )
-             dcalc3(:,:,: )   = 0.0 !! 
-          ENDIF
-          IF( med_diag%EXPC3%dgsave     ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, expc3   )
-             expc3(:,:,: )    = 0.0 !! 
-          ENDIF
-          IF( med_diag%EXPN3%dgsave     ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, expn3   )
-             expn3(:,:,: )    = 0.0 !! 
-          ENDIF
-          IF( med_diag%FEDISS3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, fediss3   )
-             fediss3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%FESCAV3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, fescav3   )
-             fescav3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%MIGRAZP3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, migrazp3  )
-             migrazp3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%MIGRAZD3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, migrazd3  )
-             migrazd3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%MEGRAZP3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, megrazp3  )
-             megrazp3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%MEGRAZD3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, megrazd3  )
-             megrazd3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%MEGRAZZ3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, megrazz3  )
-             megrazz3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%O2SAT3%dgsave     ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, o2sat3    )
-             o2sat3(:,:,: )    = 0.0 !! 
-          ENDIF
-          IF( med_diag%PBSI3%dgsave      ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pbsi3     )
-             pbsi3(:,:,: )     = 0.0 !! 
-          ENDIF
-          IF( med_diag%PCAL3%dgsave      ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pcal3     )
-             pcal3(:,:,: )     = 0.0 !! 
-          ENDIF
-          IF( med_diag%REMOC3%dgsave     ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, remoc3    )
-             remoc3(:,:,: )    = 0.0 !! 
-          ENDIF
-          IF( med_diag%PNLIMJ3%dgsave    ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pnlimj3   )
-             pnlimj3(:,:,: )   = 0.0 !! 
-          ENDIF
-          IF( med_diag%PNLIMN3%dgsave    ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pnlimn3   )
-             pnlimn3(:,:,: )   = 0.0 !! 
-          ENDIF
-          IF( med_diag%PNLIMFE3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pnlimfe3  )
-             pnlimfe3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%PDLIMJ3%dgsave    ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pdlimj3   )
-             pdlimj3(:,:,: )   = 0.0 !! 
-          ENDIF
-          IF( med_diag%PDLIMN3%dgsave    ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pdlimn3   )
-             pdlimn3(:,:,: )   = 0.0 !! 
-          ENDIF
-          IF( med_diag%PDLIMFE3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pdlimfe3  )
-             pdlimfe3(:,:,: )  = 0.0 !! 
-          ENDIF
-          IF( med_diag%PDLIMSI3%dgsave   ) THEN
-             CALL wrk_alloc( jpi, jpj, jpk, pdlimsi3  )
-             pdlimsi3(:,:,: )  = 0.0 !! 
-          ENDIF
-
-       ENDIF
-       !! lk_iomput                                   
+! tmp - marc
+      write(numout,*) 'bbb13. before call to bio_medusa_init,kt=',kt
+      flush(numout)
+!
+      CALL bio_medusa_init( kt )
+! tmp - marc
+      write(numout,*) 'bbb14. after call to bio_medusa_init,kt=',kt
+      flush(numout)
+!
        !!
 # if defined key_axy_nancheck
@@ -1884,5 +1085,5 @@
                       ENDIF
                       IF( med_diag%CO2FLUX%dgsave ) THEN
-                         f_co2flux(ji,jj)2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk)           !! mmol/m3/d -> mmol/m2/d
+                         f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk)           !! mmol/m3/d -> mmol/m2/d
                       ENDIF
                       IF( med_diag%TCO2%dgsave ) THEN
