Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 7973)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 7975)
@@ -63,5 +63,5 @@
       USE mocsy_wrapper,     ONLY: mocsy_interface
       USE oce,               ONLY: PCO2a_in_cpl
-      USE par_oce,           ONLY: jpim1, jpjm1, jpk
+      USE par_oce,           ONLY: jpim1, jpjm1
       USE sbc_oce,           ONLY: fr_i, lk_oasis, qsr, wndm
       USE sms_medusa,        ONLY: jdms, jdms_input, jdms_model,          &
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 7973)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 7975)
@@ -20,19 +20,59 @@
 
    !! model state variables
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zchn,zchd,zphn,zphd,zpds,zzmi
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zzme,zdet,zdtc,zdin,zsil,zfer
-# if defined key_roam
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zdic, zalk, zoxy
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    ztmp, zsal
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zchn,zchd,zphn,zphd,zpds,zzmi
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zzme,zdet,zdtc,zdin,zsil,zfer
+# if defined key_roam
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdic, zalk, zoxy
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmp, zsal
 # endif
 # if defined key_mocsy
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zpho
-# endif
-
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpho
+# endif
+
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fthetan,faln,fchn1,fchn
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fjln,fprn,frn
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fthetad,fald,fchd1,fchd
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fjld,fprd,frd
+
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fjlim_pn, fjlim_pd
+   !! AXY (03/02/11): add in Liebig terms
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fpnlim, fpdlim
+   !! AXY (16/07/09): add in Eppley curve functionality
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fun_T,xvpnT,xvpdT
+
+   !! AXY (16/05/11): per Katya's prompting, add in new T-dependence
+   !!                 for phytoplankton growth only (i.e. no change
+   !!                 for remineralisation)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fun_Q10
    !! AXY (01/03/10): add in mixed layer PP diagnostics
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fprn_ml,fprd_ml
+   !!
+   !! nutrient limiting factors
+   !! N and Fe (renaming ffln to ffln2 to avoid conflict with
+   !! ffln in module sms_medusa - marc 25/4/17)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fnln,ffln2
+   !! N, Fe and Si
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fnld,ffld,fsld,fsld2
+   !!
+   !! silicon cycle
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fsin,fnsi,fprds,fsdiss
 
    !! Variable for iron-ligand system
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: xFree
+
+   !! Microzooplankton grazing
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fmi1,fmi,fgmipn,fgmid,fgmidc
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: finmi,ficmi,fstarmi,fmith
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fmigrow,fmiexcr,fmiresp
+   !!
+   !! Mesozooplankton grazing
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fme1,fme,fgmepn,fgmepd
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fgmepds,fgmezmi,fgmed,fgmedc
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: finme,ficme,fstarme,fmeth
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fmegrow,fmeexcr,fmeresp
+   !!
+   !! mortality/Remineralisation (defunct parameter "fz" removed)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
 
    !! Mortality/Remineralisation
@@ -44,4 +84,6 @@
 
    !! Particle flux
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdep1,fcaco3
+
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ffastn,ffastsi,ffastfe
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ffastc,ffastca
@@ -67,4 +109,7 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ftot_pn,ftot_pd
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ftot_zmi,ftot_zme,ftot_det,ftot_dtc
+
+   !! diagnose fluxes (should only be used in 1D runs)
+   INTEGER                               :: idf, idfval
 
    !! Nitrogen and silicon production and consumption
@@ -208,6 +253,36 @@
                ztmp(jpi,jpj),zsal(jpi,jpj),                           &
                zpho(jpi,jpj),                                         &
+               fthetan(jpi,jpj),faln(jpi,jpj),fchn1(jpi,jpj),         &
+               fchn(jpi,jpj),fjln(jpi,jpj),fprn(jpi,jpj),             &
+               frn(jpi,jpj),                                          &
+               fthetad(jpi,jpj),fald(jpi,jpj),fchd1(jpi,jpj),         &
+               fchd(jpi,jpj),fjld(jpi,jpj),fprd(jpi,jpj),             &
+               frd(jpi,jpj),                                          &
+               fjlim_pn(jpi,jpj), fjlim_pd(jpi,jpj),                  &
+               fpnlim(jpi,jpj), fpdlim(jpi,jpj),                      &
+               fun_T(jpi,jpj),xvpnT(jpi,jpj),xvpdT(jpi,jpj),          &
+               fun_Q10(jpi,jpj),                                      &
                fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),                     &
+               fnln(jpi,jpj),ffln2(jpi,jpj),                          &
+               fnld(jpi,jpj),ffld(jpi,jpj),fsld(jpi,jpj),             &
+               fsld2(jpi,jpj),                                        &
+               fsin(jpi,jpj),fnsi(jpi,jpj),fprds(jpi,jpj),            &
+               fsdiss(jpi,jpj),                                       &
                xFree(jpi,jpj),                                        &
+               fmi1(jpi,jpj),fmi(jpi,jpj),fgmipn(jpi,jpj),            &
+               fgmid(jpi,jpj),fgmidc(jpi,jpj),                        &
+               finmi(jpi,jpj),ficmi(jpi,jpj),fstarmi(jpi,jpj),        &
+               fmith(jpi,jpj),fmigrow(jpi,jpj),fmiexcr(jpi,jpj),      &
+               fmiresp(jpi,jpj),                                      &
+               fme1(jpi,jpj),fme(jpi,jpj),fgmepn(jpi,jpj),            &
+               fgmepd(jpi,jpj),fgmepds(jpi,jpj),fgmezmi(jpi,jpj),     &
+               fgmed(jpi,jpj),fgmedc(jpi,jpj),                        &
+               finme(jpi,jpj),ficme(jpi,jpj),fstarme(jpi,jpj),        &
+               fmeth(jpi,jpj),fmegrow(jpi,jpj),fmeexcr(jpi,jpj),      &
+               fmeresp(jpi,jpj),                                      &
+               fdpn(jpi,jpj),fdpd(jpi,jpj),fdpds(jpi,jpj),            &
+               fdzmi(jpi,jpj),fdzme(jpi,jpj),fdd(jpi,jpj),            &
+               fdpn2(jpi,jpj),fdpd2(jpi,jpj),fdpds2(jpi,jpj),         &
+               fdzmi2(jpi,jpj),fdzme2(jpi,jpj),                       &
                fslownflux(jpi,jpj),fslowcflux(jpi,jpj),               &
                fregenfast(jpi,jpj),fregenfastsi(jpi,jpj),             &
@@ -215,4 +290,5 @@
                fregenfastc(jpi,jpj),                                  &
 # endif
+               fdep1(jpi,jpj),fcaco3(jpi,jpj),                        &
                ffastn(jpi,jpj),ffastsi(jpi,jpj),ffastfe(jpi,jpj),     &
                ffastc(jpi,jpj),ffastca(jpi,jpj),                      &
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 7973)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 7975)
@@ -40,6 +40,6 @@
                                    f_schmidtco2, f_TALK, f_TDIC, f_xco2a, &
                                    zalk, zdic, zpho, zsal, zsil, ztmp 
-      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_n, gphit,      &
-                                   mbathy, tmask
+      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_0, gdepw_n,    &
+                                   gphit, mbathy, tmask
       USE in_out_manager,    ONLY: lwp, numout
       USE mocsy_wrapper,     ONLY: mocsy_interface
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 7975)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 7975)
@@ -0,0 +1,474 @@
+MODULE phytoplankton_mod
+   !!======================================================================
+   !!                         ***  MODULE phytoplankton_mod  ***
+   !! Calculates the phytoplankton growth
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   phytoplankton        ! Called in plankton.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE phytoplankton( jk )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE phytoplankton  ***
+      !! This called from PLANKTON and calculates the phytoplankton
+      !! growth.
+      !!----------------------------------------------------------------------
+      USE bio_medusa_mod,    ONLY: fald, faln, fchd, fchd1, fchn, fchn1, &
+                                   fdep1, ffld, ffln2, fjld, fjln,       &
+                                   fjlim_pd, fjlim_pn, fjln,             &
+                                   fnld, fnln, fnsi,                     &
+                                   fpdlim, fpnlim, fprd, fprd_ml, fprds, &
+                                   fprn, fprn_ml, frd, frn,              &
+                                   fsin, fsld, fsld2, fthetad, fthetan,  & 
+                                   ftot_det, ftot_dtc, ftot_pd,          &
+                                   ftot_pn, ftot_zme, ftot_zmi,          &
+                                   fun_Q10, fun_T, idf, idfval,          &
+                                   xvpdT, xvpnT,                         &
+                                   zchd, zchn, zdet, zdin, zdtc,         &
+                                   zfer, zpds, zphd, zphn, zsil,         &
+                                   zzme, zzmi
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
+      USE in_out_manager,    ONLY: lwp, numout
+      USE oce,               ONLY: tsn
+      USE par_kind,          ONLY: wp
+      USE par_oce,           ONLY: jp_tem, jpim1, jpjm1
+      USE phycst,            ONLY: rsmall
+      USE sms_medusa,        ONLY: jliebig, jphy, jq10,                  &
+                                   xald, xaln, xfld, xfln,               &
+                                   xnld, xnln, xnsi0, xpar,              &
+                                   xsin0, xsld, xthetam, xthetamd, xuif, &
+                                   xvpd, xvpn, xxi
+      USE zdfmxl,            ONLY: hmld
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+
+      !! Level
+      INTEGER, INTENT( in ) :: jk
+
+      INTEGER :: ji, jj
+
+      REAL(wp)                     ::    fsin1,fnsi1,fnsi2
+      REAL(wp) ::    fq0
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
+               !!----------------------------------------------------------
+               !! Chlorophyll calculations
+               !!----------------------------------------------------------
+               !!
+               !! non-diatoms
+	       if (zphn(ji,jj).GT.rsmall) then
+                  fthetan(ji,jj) = max(tiny(zchn(ji,jj)),                    &
+                                       (zchn(ji,jj) * xxi) /                 &
+                                       (zphn(ji,jj) + tiny(zphn(ji,jj))))
+                  faln(ji,jj)    = xaln * fthetan(ji,jj)
+               else
+                  fthetan(ji,jj) = 0.
+                  faln(ji,jj)    = 0.
+               endif
+               !!
+               !! diatoms
+	       if (zphd(ji,jj).GT.rsmall) then
+                  fthetad(ji,jj) = max(tiny(zchd(ji,jj)),                   &
+                                       (zchd(ji,jj) * xxi) /                &
+                                       (zphd(ji,jj) + tiny(zphd(ji,jj))))
+                  fald(ji,jj)    = xald * fthetad(ji,jj)
+               else
+                  fthetad(ji,jj) = 0.
+                  fald(ji,jj)    = 0.
+               endif
+
+# if defined key_debug_medusa
+               !! report biological calculations
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'faln(',jk,') = ', faln(ji,jj)
+                  IF (lwp) write (numout,*) 'fald(',jk,') = ', fald(ji,jj)
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Phytoplankton light limitation
+               !!----------------------------------------------------------
+               !!
+               !! It is assumed xpar is the depth-averaged (vertical layer) PAR 
+               !! Light limitation (check self-shading) in W/m2
+               !!
+               !! Note that there is no temperature dependence in phytoplankton
+               !! growth rate or any other function. 
+               !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to 
+               !! avoid NaNs in case of Phy==0.  
+               !!
+               !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and 
+               !! non-diat: 
+               !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012
+               !!
+               !! AXY (16/07/09)
+               !! temperature for new Eppley style phytoplankton growth
+               fun_T(ji,jj)   = 1.066**(1.0 * tsn(ji,jj,jk,jp_tem))
+               !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for
+               !!                 phytoplankton growth; remin. unaffected
+               fun_Q10(ji,jj) = jq10**((tsn(ji,jj,jk,jp_tem) - 0.0) / 10.0)
+               if (jphy.eq.1) then
+                  xvpnT(ji,jj) = xvpn * fun_T(ji,jj)
+                  xvpdT(ji,jj) = xvpd * fun_T(ji,jj)
+               elseif (jphy.eq.2) then
+                  xvpnT(ji,jj) = xvpn * fun_Q10(ji,jj)
+                  xvpdT(ji,jj) = xvpd * fun_Q10(ji,jj)
+               else
+                  xvpnT(ji,jj) = xvpn
+                  xvpdT(ji,jj) = xvpd
+               endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then           
+               !!
+               !! non-diatoms
+               fchn1(ji,jj)   = (xvpnT(ji,jj) * xvpnT(ji,jj)) +              &
+                                 (faln(ji,jj) * faln(ji,jj) *                &
+                                  xpar(ji,jj,jk) * xpar(ji,jj,jk))
+               if (fchn1(ji,jj).GT.rsmall) then
+                  fchn(ji,jj)    = xvpnT(ji,jj) / (sqrt(fchn1(ji,jj)) +      &
+                                                   tiny(fchn1(ji,jj)))
+               else
+                  fchn(ji,jj)    = 0.
+               endif
+               !! non-diatom J term
+               fjln(ji,jj)     = fchn(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk)
+               fjlim_pn(ji,jj) = fjln(ji,jj) / xvpnT(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! diatoms
+               fchd1(ji,jj)   = (xvpdT(ji,jj) * xvpdT(ji,jj)) +              &
+                                (fald(ji,jj) * fald(ji,jj) *                 &
+                                 xpar(ji,jj,jk) * xpar(ji,jj,jk))
+               if (fchd1(ji,jj).GT.rsmall) then
+                  fchd(ji,jj)    = xvpdT(ji,jj) / (sqrt(fchd1(ji,jj)) +      &
+                                                   tiny(fchd1(ji,jj)))
+               else
+                  fchd(ji,jj)    = 0.
+               endif
+               !! diatom J term
+               fjld(ji,jj)    = fchd(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk)
+               fjlim_pd(ji,jj) = fjld(ji,jj) / xvpdT(ji,jj)
+      
+# if defined key_debug_medusa
+               !! report phytoplankton light limitation
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn(ji,jj)
+                  IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd(ji,jj)
+                  IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln(ji,jj)
+                  IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld(ji,jj)
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Phytoplankton nutrient limitation
+               !!----------------------------------------------------------
+               !!
+               !! non-diatoms (N, Fe).
+               !! non-diatom Qn term
+               fnln(ji,jj)  = zdin(ji,jj) / (zdin(ji,jj) + xnln)
+               !! non-diatom Qf term
+               ffln2(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfln)
+               !!
+               !! diatoms (N, Si, Fe).
+               !! diatom Qn term
+               fnld(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnld)
+               !! diatom Qs term
+               fsld(ji,jj) = zsil(ji,jj) / (zsil(ji,jj) + xsld)
+               !! diatom Qf term
+               ffld(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfld)
+
+# if defined key_debug_medusa
+               !! report phytoplankton nutrient limitation
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln(ji,jj)
+                  IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld(ji,jj)
+                  IF (lwp) write (numout,*) 'ffln2(',jk,') = ', ffln2(ji,jj)
+                  IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld(ji,jj)
+                  IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld(ji,jj)
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Primary production (non-diatoms)
+               !! (note: still needs multiplying by phytoplankton 
+               !! concentration)
+               !!----------------------------------------------------------
+               !!
+               if (jliebig .eq. 0) then
+                  !! multiplicative nutrient limitation
+                  fpnlim(ji,jj) = fnln(ji,jj) * ffln2(ji,jj)
+               elseif (jliebig .eq. 1) then
+                  !! Liebig Law (= most limiting) nutrient limitation
+                  fpnlim(ji,jj) = min(fnln(ji,jj), ffln2(ji,jj))
+               endif
+               fprn(ji,jj) = fjln(ji,jj) * fpnlim(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Primary production (diatoms)
+               !! (note: still needs multiplying by phytoplankton 
+               !! concentration)
+               !!
+               !! Production here is split between nitrogen production and 
+               !! that of silicon; depending upon the "intracellular" ratio 
+               !! of Si:N, model diatoms will uptake nitrogen/silicon 
+               !! differentially; this borrows from the diatom model of 
+               !! Mongin et al. (2006)
+               !!----------------------------------------------------------
+               !!
+               if (jliebig .eq. 0) then
+                  !! multiplicative nutrient limitation
+                  fpdlim(ji,jj) = fnld(ji,jj) * ffld(ji,jj)
+               elseif (jliebig .eq. 1) then
+                  !! Liebig Law (= most limiting) nutrient limitation
+                  fpdlim(ji,jj) = min(fnld(ji,jj), ffld(ji,jj))
+               endif
+               !!
+	       if (zphd(ji,jj).GT.rsmall .AND. zpds(ji,jj).GT.rsmall) then
+                  !! "intracellular" elemental ratios
+                  ! fsin(ji,jj)  = zpds(ji,jj) / (zphd(ji,jj) +              &
+                  !                               tiny(zphd(ji,jj)))
+                  ! fnsi(ji,jj)  = zphd(ji,jj) / (zpds(ji,jj) +              &
+                  !                               tiny(zpds(ji,jj)))
+                  fsin(ji,jj) = 0.0
+                  IF( zphd(ji,jj) .GT. rsmall) fsin(ji,jj)  = zpds(ji,jj) /  &
+                                                              zphd(ji,jj)
+                  fnsi(ji,jj) = 0.0
+                  IF( zpds(ji,jj) .GT. rsmall) fnsi(ji,jj)  = zphd(ji,jj) /  &
+                                                              zpds(ji,jj)
+                  !! AXY (23/02/10): these next variables derive from 
+                  !! Mongin et al. (2003)
+                  fsin1 = 3.0 * xsin0 !! = 0.6
+                  fnsi1 = 1.0 / fsin1 !! = 1.667
+                  fnsi2 = 1.0 / xsin0 !! = 5.0
+                  !!
+                  !! conditionalities based on ratios
+                  !! nitrogen (and iron and carbon)
+                  if (fsin(ji,jj).le.xsin0) then
+                     fprd(ji,jj)  = 0.0
+                     fsld2(ji,jj) = 0.0
+                  elseif (fsin(ji,jj).lt.fsin1) then
+                     fprd(ji,jj)  = xuif * ((fsin(ji,jj) - xsin0) /          &
+                                            (fsin(ji,jj) +                   &
+                                             tiny(fsin(ji,jj)))) *           &
+                                    (fjld(ji,jj) * fpdlim(ji,jj))
+                     fsld2(ji,jj) = xuif * ((fsin(ji,jj) - xsin0) /          &
+                                            (fsin(ji,jj) +                   &
+                                             tiny(fsin(ji,jj))))
+                  elseif (fsin(ji,jj).ge.fsin1) then
+                     fprd(ji,jj)  = (fjld(ji,jj) * fpdlim(ji,jj))
+                     fsld2(ji,jj) = 1.0
+                  endif
+                  !!
+                  !! silicon
+                  if (fsin(ji,jj).lt.fnsi1) then
+                     fprds(ji,jj) = (fjld(ji,jj) * fsld(ji,jj))
+                  elseif (fsin(ji,jj).lt.fnsi2) then
+                     fprds(ji,jj) = xuif * ((fnsi(ji,jj) - xnsi0) /          &
+                                            (fnsi(ji,jj) +                   &
+                                             tiny(fnsi(ji,jj)))) *           &
+                                    (fjld(ji,jj) * fsld(ji,jj))
+                  else
+                     fprds(ji,jj) = 0.0
+                  endif     
+               else
+                  fsin(ji,jj)  = 0.0
+                  fnsi(ji,jj)  = 0.0
+                  fprd(ji,jj)  = 0.0
+                  fsld2(ji,jj) = 0.0
+                  fprds(ji,jj) = 0.0
+               endif
+
+# if defined key_debug_medusa
+               !! report phytoplankton growth (including diatom silicon 
+               !! submodel)
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin(ji,jj)
+                  IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi(ji,jj)
+                  IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2(ji,jj)
+                  IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn(ji,jj)
+                  IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd(ji,jj)
+                  IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds(ji,jj)
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Mixed layer primary production
+               !! this block calculates the amount of primary production 
+               !! that occurs within the upper mixed layer; this allows the 
+               !! separate diagnosis of "sub-surface" primary production; it 
+               !! does assume that short-term variability in mixed layer 
+               !! depth doesn't mess with things though
+               !!----------------------------------------------------------
+               !!
+               if (fdep1(ji,jj).le.hmld(ji,jj)) then
+                  !! this level is entirely in the mixed layer
+                  fq0 = 1.0
+               elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
+                  !! this level is entirely below the mixed layer
+                  fq0 = 0.0
+               else
+                  !! this level straddles the mixed layer
+                  fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
+               endif
+               !!
+               fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn(ji,jj) * zphn(ji,jj) * &
+                                                  fse3t(ji,jj,jk) * fq0)
+               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * &
+                                                  fse3t(ji,jj,jk) * fq0)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Vertical Integral --
+               !!----------------------------------------------------------
+               !! vertical integral non-diatom phytoplankton
+               ftot_pn(ji,jj)  = ftot_pn(ji,jj)  + (zphn(ji,jj) *             &
+                                                    fse3t(ji,jj,jk))
+               !! vertical integral diatom phytoplankton
+               ftot_pd(ji,jj)  = ftot_pd(ji,jj)  + (zphd(ji,jj) *             &
+                                                    fse3t(ji,jj,jk))
+               !! vertical integral microzooplankton
+               ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi(ji,jj) *             &
+                                                    fse3t(ji,jj,jk))
+               !! vertical integral mesozooplankton
+               ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme(ji,jj) *             &
+                                                    fse3t(ji,jj,jk))
+               !! vertical integral slow detritus, nitrogen
+               ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet(ji,jj) *             &
+                                                    fse3t(ji,jj,jk))
+               !! vertical integral slow detritus, carbon
+               ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc(ji,jj) *             &
+                                                    fse3t(ji,jj,jk))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! More chlorophyll calculations
+               !!----------------------------------------------------------
+               !!
+               !! frn(ji,jj) = (xthetam / fthetan(ji,jj)) *                   &
+               !!              (fprn(ji,jj) / (fthetan(ji,jj) * xpar(ji,jj,jk)))
+               !! frd(ji,jj) = (xthetam / fthetad(ji,jj)) *                   &
+               !!              (fprd(ji,jj) / (fthetad(ji,jj) * xpar(ji,jj,jk)))
+               frn(ji,jj) = (xthetam * fchn(ji,jj) * fnln(ji,jj) *            &
+                             ffln2(ji,jj)) / (fthetan(ji,jj) +                &
+                                             tiny(fthetan(ji,jj)))
+               !! AXY (12/05/09): there's potentially a problem here; fsld, 
+               !!   silicic acid limitation, is used in the following line 
+               !!   to regulate chlorophyll growth in a manner that is 
+               !!   inconsistent with its use in the regulation of biomass 
+               !!   growth; the Mongin term term used in growth is more
+               !!   complex than the simple multiplicative function used
+               !!   below
+               !! frd(ji,jj) = (xthetam * fchd(ji,jj) * fnld(ji,jj) *        &
+               !!               ffld(ji,jj) * fsld(ji,jj)) /                 &
+               !!               (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
+               !! AXY (12/05/09): this replacement line uses the new 
+               !!   variable, fsld2, to regulate chlorophyll growth
+               frd(ji,jj) = (xthetamd * fchd(ji,jj) * fnld(ji,jj) *          &
+                             ffld(ji,jj) * fsld2(ji,jj)) /                   &
+                             (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
+
+# if defined key_debug_medusa
+               !! report chlorophyll calculations
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan(ji,jj)
+                  IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad(ji,jj)
+                  IF (lwp) write (numout,*) 'frn(',jk,')     = ', frn(ji,jj)
+                  IF (lwp) write (numout,*) 'frd(',jk,')     = ', frd(ji,jj)
+               endif
+# endif
+
+            ENDIF
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE phytoplankton
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE phytoplankton( )                    ! Empty routine
+      WRITE(*,*) 'phytoplankton: You should not have seen this print! error?'
+   END SUBROUTINE phytoplankton
+#endif 
+
+   !!======================================================================
+END MODULE phytoplankton_mod
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 7975)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 7975)
@@ -0,0 +1,193 @@
+MODULE plankton_mod
+   !!======================================================================
+   !!                         ***  MODULE plankton_mod  ***
+   !! Calculate the carbon chemistry for the whole ocean
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   plankton        ! 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 plankton( jk )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE plankton  ***
+      !! This called from TRC_BIO_MEDUSA and 
+      !!  - Calculates phytoplankton growth
+      !!  - Zooplankton grazing
+      !!  - Plankton losses
+      !!-------------------------------------------------------------------
+      USE bio_medusa_mod,    ONLY: fdpd, fdpd2, fdpds, fdpds2,             &
+                                   fdpn, fdpn2, fdzme, fdzme2,             &
+                                   fdzmi, fdzmi2, fsin,                    &
+                                   zphd, zphn, zpds, zzme, zzmi
+      USE dom_oce,           ONLY: tmask
+      USE par_oce,           ONLY: jpim1, jpjm1
+      USE phytoplankton_mod, ONLY: phytoplankton
+      USE sms_medusa,        ONLY: jmpd, jmpn, jmzme, jmzmi,               &
+                                   xkphd, xkphn, xkzme, xkzmi,             &
+                                   xmetapd, xmetapn, xmetazme, xmetazmi,   &
+                                   xmpd, xmpn, xmzme, xmzmi
+      USE zooplankton_mod,   ONLY: zooplankton
+
+      !! Level
+      INTEGER, INTENT( in ) :: jk
+
+      INTEGER :: ji, jj
+
+      !!-------------------------------------------------------------------
+      !! Calculate phytoplankton growth
+      !!-------------------------------------------------------------------
+      CALL phytoplankton( jk )
+
+      !!-------------------------------------------------------------------
+      !! Calculate zooplankton grazing
+      !!-------------------------------------------------------------------
+      CALL zooplankton( jk )
+
+      !!-------------------------------------------------------------------
+      !! Miscellaneous plankton losses
+      !!-------------------------------------------------------------------
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
+               !!----------------------------------------------------------
+               !! Plankton metabolic losses
+               !! Linear loss processes assumed to be metabolic in origin
+               !!----------------------------------------------------------
+               !!
+               fdpn2(ji,jj)  = xmetapn  * zphn(ji,jj)
+               fdpd2(ji,jj)  = xmetapd  * zphd(ji,jj)
+               fdpds2(ji,jj) = xmetapd  * zpds(ji,jj)
+               fdzmi2(ji,jj) = xmetazmi * zzmi(ji,jj)
+               fdzme2(ji,jj) = xmetazme * zzme(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Plankton mortality losses
+               !! EKP (26/02/09): phytoplankton hyperbolic mortality term 
+               !! introduced 
+               !! to improve performance in gyres
+               !!----------------------------------------------------------
+               !!
+               !! non-diatom phytoplankton
+               !! linear
+               if (jmpn.eq.1) fdpn(ji,jj) = xmpn * zphn(ji,jj)
+               !! quadratic
+               if (jmpn.eq.2) fdpn(ji,jj) = xmpn * zphn(ji,jj) * zphn(ji,jj)
+               !! hyperbolic
+               if (jmpn.eq.3) fdpn(ji,jj) = xmpn * zphn(ji,jj) *             &
+                                            (zphn(ji,jj) /                   &
+                                             (xkphn + zphn(ji,jj)))
+               !! sigmoid
+               if (jmpn.eq.4) fdpn(ji,jj) = xmpn * zphn(ji,jj) *             &
+                                            ((zphn(ji,jj) * zphn(ji,jj)) /   &
+                                             (xkphn + (zphn(ji,jj) *         &
+                                                       zphn(ji,jj))))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! diatom phytoplankton
+               !! linear
+               if (jmpd.eq.1) fdpd(ji,jj) = xmpd * zphd(ji,jj)
+               !! quadratic
+               if (jmpd.eq.2) fdpd(ji,jj) = xmpd * zphd(ji,jj) * zphd(ji,jj)
+               !! hyperbolic
+               if (jmpd.eq.3) fdpd(ji,jj) = xmpd * zphd(ji,jj) *             &
+                                            (zphd(ji,jj) / (xkphd +          &
+                                                            zphd(ji,jj)))
+               !! sigmoid
+               if (jmpd.eq.4) fdpd(ji,jj) = xmpd * zphd(ji,jj) *             &
+                                            ((zphd(ji,jj) * zphd(ji,jj)) /   &
+                                             (xkphd + (zphd(ji,jj) *         &
+                                                       zphd(ji,jj))))
+               fdpds(ji,jj) = fdpd(ji,jj) * fsin(ji,jj)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! microzooplankton
+               !! linear
+               if (jmzmi.eq.1) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj)
+               !! quadratic
+               if (jmzmi.eq.2) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) *          &
+                                              zzmi(ji,jj)
+               !! hyperbolic
+               if (jmzmi.eq.3) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) *          &
+                                              (zzmi(ji,jj) / (xkzmi +        &
+                                                              zzmi(ji,jj)))
+               !! sigmoid
+               if (jmzmi.eq.4) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * &
+                  ((zzmi(ji,jj) * zzmi(ji,jj)) / (xkzmi + (zzmi(ji,jj) *     &
+                                                           zzmi(ji,jj))))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! mesozooplankton
+               !! linear
+               if (jmzme.eq.1) fdzme(ji,jj) = xmzme * zzme(ji,jj)
+               !! quadratic
+               if (jmzme.eq.2) fdzme(ji,jj) = xmzme * zzme(ji,jj) *          &
+                                              zzme(ji,jj)
+               !! hyperbolic
+               if (jmzme.eq.3) fdzme(ji,jj) = xmzme * zzme(ji,jj) *          &
+                                              (zzme(ji,jj) / (xkzme +        &
+                                                              zzme(ji,jj)))
+               !! sigmoid
+               if (jmzme.eq.4) fdzme(ji,jj) = xmzme * zzme(ji,jj) *          &
+                                              ((zzme(ji,jj) * zzme(ji,jj)) / &
+                                               (xkzme + (zzme(ji,jj) *       &
+                                                         zzme(ji,jj))))
+            ENDIF
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE plankton
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE plankton( )                    ! Empty routine
+      WRITE(*,*) 'plankton: You should not have seen this print! error?'
+   END SUBROUTINE plankton
+#endif 
+
+   !!======================================================================
+END MODULE plankton_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 7973)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7975)
@@ -100,4 +100,5 @@
       USE carb_chem_mod,              ONLY: carb_chem
       USE air_sea_mod,                ONLY: air_sea
+      USE plankton_mod,               ONLY: plankton
       USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
       USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
@@ -181,26 +182,26 @@
       !!
       !! primary production and chl related quantities      
-      REAL(wp), DIMENSION(jpi,jpj) ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
-      REAL(wp), DIMENSION(jpi,jpj) ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
       !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range)
-      REAL(wp), DIMENSION(jpi,jpj) ::    fjlim_pn, fjlim_pd
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fjlim_pn, fjlim_pd
       !! AXY (03/02/11): add in Liebig terms
-      REAL(wp), DIMENSION(jpi,jpj) ::    fpnlim, fpdlim
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fpnlim, fpdlim
       !! AXY (16/07/09): add in Eppley curve functionality
-      REAL(wp), DIMENSION(jpi,jpj) ::    fun_T,xvpnT,xvpdT
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fun_T,xvpnT,xvpdT
       INTEGER  ::    ieppley
       !! AXY (16/05/11): per Katya's prompting, add in new T-dependence
       !!                 for phytoplankton growth only (i.e. no change
       !!                 for remineralisation)
-      REAL(wp), DIMENSION(jpi,jpj) ::    fun_Q10
+!      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
       !!
       !! nutrient limiting factors
-      REAL(wp), DIMENSION(jpi,jpj) ::    fnln,ffln            !! N and Fe
-      REAL(wp), DIMENSION(jpi,jpj) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fnln,ffln2            !! N and Fe
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si
       !!
       !! silicon cycle
-      REAL(wp), DIMENSION(jpi,jpj) ::    fsin,fnsi,fprds,fsdiss
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fsin,fnsi,fprds,fsdiss
       REAL(wp)                     ::    fsin1,fnsi1,fnsi2
       !!
@@ -223,17 +224,17 @@
       !!
       !! microzooplankton grazing
-      REAL(wp), DIMENSION(jpi,jpj) ::    fmi1,fmi,fgmipn,fgmid,fgmidc
-      REAL(wp), DIMENSION(jpi,jpj) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fmi1,fmi,fgmipn,fgmid,fgmidc
+!      REAL(wp), DIMENSION(jpi,jpj) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
       !!
       !! mesozooplankton grazing
-      REAL(wp), DIMENSION(jpi,jpj) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
-      REAL(wp), DIMENSION(jpi,jpj) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
+!      REAL(wp), DIMENSION(jpi,jpj) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
       !!
       !! 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
 # endif
-      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
       REAL(wp), DIMENSION(jpi,jpj) ::    fslown, fslowc
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux
@@ -247,5 +248,5 @@
       !!
       !! particle flux
-      REAL(WP), DIMENSION(jpi,jpj) ::    fdep1,fcaco3
+!      REAL(WP), DIMENSION(jpi,jpj) ::    fdep1,fcaco3
       REAL(WP), DIMENSION(jpi,jpj) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
       REAL(wp), DIMENSION(jpi,jpj) ::    freminn,freminsi,freminfe,freminc,freminca
@@ -297,5 +298,5 @@
       !!
       !! diagnose fluxes (should only be used in 1D runs)
-      INTEGER  ::    idf, idfval
+!      INTEGER  ::    idf, idfval
       !!
       !! nitrogen and silicon production and consumption
@@ -844,484 +845,16 @@
 # endif
 
+         !!------------------------------------------------------------------
+         !! Phytoplankton growth, zooplankton grazing and miscellaneous
+         !! plankton losses. 
+         !!------------------------------------------------------------------
+         CALL plankton( jk )
+
+! Detritus process - marc
+
          DO jj = 2,jpjm1
          DO ji = 2,jpim1
             !! OPEN wet point IF..THEN loop
             if (tmask(ji,jj,jk) == 1) then
-
-               !!----------------------------------------------------------------------
-               !! Chlorophyll calculations
-               !!----------------------------------------------------------------------
-               !!
-               !! non-diatoms
-	       if (zphn(ji,jj).GT.rsmall) then
-                  fthetan(ji,jj) = max(tiny(zchn(ji,jj)), (zchn(ji,jj) * xxi) / (zphn(ji,jj) + tiny(zphn(ji,jj))))
-                  faln(ji,jj)    = xaln * fthetan(ji,jj)
-               else
-                  fthetan(ji,jj) = 0.
-                  faln(ji,jj)    = 0.
-               endif
-               !!
-               !! diatoms
-	       if (zphd(ji,jj).GT.rsmall) then
-                  fthetad(ji,jj) = max(tiny(zchd(ji,jj)), (zchd(ji,jj) * xxi) / (zphd(ji,jj) + tiny(zphd(ji,jj))))
-                  fald(ji,jj)    = xald * fthetad(ji,jj)
-               else
-                  fthetad(ji,jj) = 0.
-                  fald(ji,jj)    = 0.
-               endif
-
-# if defined key_debug_medusa
-               !! report biological calculations
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'faln(',jk,') = ', faln(ji,jj)
-                  IF (lwp) write (numout,*) 'fald(',jk,') = ', fald(ji,jj)
-               endif
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Phytoplankton light limitation
-               !!----------------------------------------------------------------------
-               !!
-               !! It is assumed xpar is the depth-averaged (vertical layer) PAR 
-               !! Light limitation (check self-shading) in W/m2
-               !!
-               !! Note that there is no temperature dependence in phytoplankton
-               !! growth rate or any other function. 
-               !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid
-               !! NaNs in case of Phy==0.  
-               !!
-               !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and non-diat: 
-               !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012
-               !!
-               !! AXY (16/07/09)
-               !! temperature for new Eppley style phytoplankton growth
-               fun_T(ji,jj)   = 1.066**(1.0 * tsn(ji,jj,jk,jp_tem))
-               !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for
-               !phytoplankton
-               !!                 growth; remin. unaffected
-               fun_Q10(ji,jj) = jq10**((tsn(ji,jj,jk,jp_tem) - 0.0) / 10.0)
-               if (jphy.eq.1) then
-                  xvpnT(ji,jj) = xvpn * fun_T(ji,jj)
-                  xvpdT(ji,jj) = xvpd * fun_T(ji,jj)
-               elseif (jphy.eq.2) then
-                  xvpnT(ji,jj) = xvpn * fun_Q10(ji,jj)
-                  xvpdT(ji,jj) = xvpd * fun_Q10(ji,jj)
-               else
-                  xvpnT(ji,jj) = xvpn
-                  xvpdT(ji,jj) = xvpd
-               endif
-               !!
-               !! non-diatoms
-               fchn1(ji,jj)   = (xvpnT(ji,jj) * xvpnT(ji,jj)) + (faln(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk) * xpar(ji,jj,jk))
-               if (fchn1(ji,jj).GT.rsmall) then
-                  fchn(ji,jj)    = xvpnT(ji,jj) / (sqrt(fchn1(ji,jj)) + tiny(fchn1(ji,jj)))
-               else
-                  fchn(ji,jj)    = 0.
-               endif
-               fjln(ji,jj)    = fchn(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk) !! non-diatom J term
-               fjlim_pn(ji,jj) = fjln(ji,jj) / xvpnT(ji,jj)
-               !!
-               !! diatoms
-               fchd1(ji,jj)   = (xvpdT(ji,jj) * xvpdT(ji,jj)) + (fald(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk) * xpar(ji,jj,jk))
-               if (fchd1(ji,jj).GT.rsmall) then
-                  fchd(ji,jj)    = xvpdT(ji,jj) / (sqrt(fchd1(ji,jj)) + tiny(fchd1(ji,jj)))
-               else
-                  fchd(ji,jj)    = 0.
-               endif
-               fjld(ji,jj)    = fchd(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk) !! diatom J term
-               fjlim_pd(ji,jj) = fjld(ji,jj) / xvpdT(ji,jj)
-      
-# if defined key_debug_medusa
-               !! report phytoplankton light limitation
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn(ji,jj)
-                  IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd(ji,jj)
-                  IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln(ji,jj)
-                  IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld(ji,jj)
-               endif
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Phytoplankton nutrient limitation
-               !!----------------------------------------------------------------------
-               !!
-               !! non-diatoms (N, Fe)
-               fnln(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnln) !! non-diatom Qn term
-               ffln(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfln) !! non-diatom Qf term
-               !!
-               !! diatoms (N, Si, Fe)
-               fnld(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnld) !! diatom Qn term
-               fsld(ji,jj) = zsil(ji,jj) / (zsil(ji,jj) + xsld) !! diatom Qs term
-               ffld(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfld) !! diatom Qf term
-
-# if defined key_debug_medusa
-               !! report phytoplankton nutrient limitation
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln(ji,jj)
-                  IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld(ji,jj)
-                  IF (lwp) write (numout,*) 'ffln(',jk,') = ', ffln(ji,jj)
-                  IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld(ji,jj)
-                  IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld(ji,jj)
-               endif
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Primary production (non-diatoms)
-               !! (note: still needs multiplying by phytoplankton concentration)
-               !!----------------------------------------------------------------------
-               !!
-               if (jliebig .eq. 0) then
-                  !! multiplicative nutrient limitation
-                  fpnlim(ji,jj) = fnln(ji,jj) * ffln(ji,jj)
-               elseif (jliebig .eq. 1) then
-                  !! Liebig Law (= most limiting) nutrient limitation
-                  fpnlim(ji,jj) = min(fnln(ji,jj), ffln(ji,jj))
-               endif
-               fprn(ji,jj) = fjln(ji,jj) * fpnlim(ji,jj)
-
-               !!----------------------------------------------------------------------
-               !! Primary production (diatoms)
-               !! (note: still needs multiplying by phytoplankton concentration)
-               !!
-               !! production here is split between nitrogen production and that of
-               !! silicon; depending upon the "intracellular" ratio of Si:N, model
-               !! diatoms will uptake nitrogen/silicon differentially; this borrows
-               !! from the diatom model of Mongin et al. (2006)
-               !!----------------------------------------------------------------------
-               !!
-               if (jliebig .eq. 0) then
-                  !! multiplicative nutrient limitation
-                  fpdlim(ji,jj) = fnld(ji,jj) * ffld(ji,jj)
-               elseif (jliebig .eq. 1) then
-                  !! Liebig Law (= most limiting) nutrient limitation
-                  fpdlim(ji,jj) = min(fnld(ji,jj), ffld(ji,jj))
-               endif
-               !!
-	       if (zphd(ji,jj).GT.rsmall .AND. zpds(ji,jj).GT.rsmall) then
-                  !! "intracellular" elemental ratios
-                  ! fsin(ji,jj)  = zpds(ji,jj) / (zphd(ji,jj) + tiny(zphd(ji,jj)))
-                  ! fnsi(ji,jj)  = zphd(ji,jj) / (zpds(ji,jj) + tiny(zpds(ji,jj)))
-                  fsin(ji,jj) = 0.0
-                  IF( zphd(ji,jj) .GT. rsmall) fsin(ji,jj)  = zpds(ji,jj) / zphd(ji,jj)
-                  fnsi(ji,jj) = 0.0
-                  IF( zpds(ji,jj) .GT. rsmall) fnsi(ji,jj)  = zphd(ji,jj) / zpds(ji,jj)
-                  !! AXY (23/02/10): these next variables derive from Mongin et al. (2003)
-                  fsin1 = 3.0 * xsin0 !! = 0.6
-                  fnsi1 = 1.0 / fsin1 !! = 1.667
-                  fnsi2 = 1.0 / xsin0 !! = 5.0
-                  !!
-                  !! conditionalities based on ratios
-                  !! nitrogen (and iron and carbon)
-                  if (fsin(ji,jj).le.xsin0) then
-                     fprd(ji,jj)  = 0.0
-                     fsld2(ji,jj) = 0.0
-                  elseif (fsin(ji,jj).lt.fsin1) then
-                     fprd(ji,jj)  = xuif * ((fsin(ji,jj) - xsin0) / (fsin(ji,jj) + tiny(fsin(ji,jj)))) * (fjld(ji,jj) * fpdlim(ji,jj))
-                     fsld2(ji,jj) = xuif * ((fsin(ji,jj) - xsin0) / (fsin(ji,jj) + tiny(fsin(ji,jj))))
-                  elseif (fsin(ji,jj).ge.fsin1) then
-                     fprd(ji,jj)  = (fjld(ji,jj) * fpdlim(ji,jj))
-                     fsld2(ji,jj) = 1.0
-                  endif
-                  !!
-                  !! silicon
-                  if (fsin(ji,jj).lt.fnsi1) then
-                     fprds(ji,jj) = (fjld(ji,jj) * fsld(ji,jj))
-                  elseif (fsin(ji,jj).lt.fnsi2) then
-                     fprds(ji,jj) = xuif * ((fnsi(ji,jj) - xnsi0) / (fnsi(ji,jj) + tiny(fnsi(ji,jj)))) * (fjld(ji,jj) * fsld(ji,jj))
-                  else
-                     fprds(ji,jj) = 0.0
-                  endif     
-               else
-                  fsin(ji,jj)  = 0.0
-                  fnsi(ji,jj)  = 0.0
-                  fprd(ji,jj)  = 0.0
-                  fsld2(ji,jj) = 0.0
-                  fprds(ji,jj) = 0.0
-               endif
-
-# if defined key_debug_medusa
-               !! report phytoplankton growth (including diatom silicon submodel)
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin(ji,jj)
-                  IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi(ji,jj)
-                  IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2(ji,jj)
-                  IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn(ji,jj)
-                  IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd(ji,jj)
-                  IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds(ji,jj)
-               endif
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Mixed layer primary production
-               !! this block calculates the amount of primary production that occurs
-               !! within the upper mixed layer; this allows the separate diagnosis
-               !! of "sub-surface" primary production; it does assume that short-
-               !! term variability in mixed layer depth doesn't mess with things
-               !! though
-               !!----------------------------------------------------------------------
-               !!
-               if (fdep1(ji,jj).le.hmld(ji,jj)) then
-                  !! this level is entirely in the mixed layer
-                  fq0 = 1.0
-               elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
-                  !! this level is entirely below the mixed layer
-                  fq0 = 0.0
-               else
-                  !! this level straddles the mixed layer
-                  fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
-               endif
-               !!
-               fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk) * fq0)
-               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk) * fq0)
-               
-               !!----------------------------------------------------------------------
-               !! Vertical Integral --
-               !!----------------------------------------------------------------------
-               ftot_pn(ji,jj)  = ftot_pn(ji,jj)  + (zphn(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral non-diatom phytoplankton
-               ftot_pd(ji,jj)  = ftot_pd(ji,jj)  + (zphd(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral diatom phytoplankton
-               ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral microzooplankton
-               ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral mesozooplankton
-               ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral slow detritus, nitrogen
-               ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral slow detritus, carbon
-               
-               !!----------------------------------------------------------------------
-               !! More chlorophyll calculations
-               !!----------------------------------------------------------------------
-               !!
-               !! frn(ji,jj) = (xthetam / fthetan(ji,jj)) * (fprn(ji,jj) / (fthetan(ji,jj) * xpar(ji,jj,jk)))
-               !! frd(ji,jj) = (xthetam / fthetad(ji,jj)) * (fprd(ji,jj) / (fthetad(ji,jj) * xpar(ji,jj,jk)))
-               frn(ji,jj) = (xthetam * fchn(ji,jj) * fnln(ji,jj) * ffln(ji,jj)       ) / (fthetan(ji,jj) + tiny(fthetan(ji,jj)))
-               !! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid 
-               !!   limitation, is used in the following line to regulate chlorophyll 
-               !!   growth in a manner that is inconsistent with its use in the regulation 
-               !!   of biomass growth; the Mongin term term used in growth is more complex
-               !!   than the simple multiplicative function used below
-               !! frd(ji,jj) = (xthetam * fchd(ji,jj) * fnld(ji,jj) * ffld(ji,jj) * fsld(ji,jj)) / (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
-               !! AXY (12/05/09): this replacement line uses the new variable, fsld2, to
-               !!   regulate chlorophyll growth
-               frd(ji,jj) = (xthetamd * fchd(ji,jj) * fnld(ji,jj) * ffld(ji,jj) * fsld2(ji,jj)) / (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
-
-# if defined key_debug_medusa
-               !! report chlorophyll calculations
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan(ji,jj)
-                  IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad(ji,jj)
-                  IF (lwp) write (numout,*) 'frn(',jk,')     = ', frn(ji,jj)
-                  IF (lwp) write (numout,*) 'frd(',jk,')     = ', frd(ji,jj)
-               endif
-# endif
-
-! MAYBE BUT A BREAK IN HERE, ZOOPLANKTON GRAZING - marc 20/4/17 
-! (plankton growth is 281 lines)
-
-               !!----------------------------------------------------------------------
-               !! Zooplankton Grazing 
-               !! this code supplements the base grazing model with one that
-               !! considers the C:N ratio of grazed food and balances this against
-               !! the requirements of zooplankton growth; this model is derived 
-               !! from that of Anderson & Pondaven (2003)
-               !!
-               !! the current version of the code assumes a fixed C:N ratio for
-               !! detritus (in contrast to Anderson & Pondaven, 2003), though the
-               !! full equations are retained for future extension
-               !!----------------------------------------------------------------------
-               !!
-               !!----------------------------------------------------------------------
-               !! Microzooplankton first
-               !!----------------------------------------------------------------------
-               !!
-               fmi1(ji,jj)    = (xkmi * xkmi) + (xpmipn * zphn(ji,jj) * zphn(ji,jj)) + (xpmid * zdet(ji,jj) * zdet(ji,jj))
-               fmi(ji,jj)     = xgmi * zzmi(ji,jj) / fmi1(ji,jj)
-               fgmipn(ji,jj)  = fmi(ji,jj) * xpmipn * zphn(ji,jj) * zphn(ji,jj)   !! grazing on non-diatoms
-               fgmid(ji,jj)   = fmi(ji,jj) * xpmid  * zdet(ji,jj) * zdet(ji,jj)   !! grazing on detrital nitrogen
-# if defined key_roam
-               fgmidc(ji,jj)  = rsmall !acc
-               IF ( zdet(ji,jj) .GT. rsmall ) fgmidc(ji,jj)  = (zdtc(ji,jj) / (zdet(ji,jj) + tiny(zdet(ji,jj)))) * fgmid(ji,jj)  !! grazing on detrital carbon
-# else
-               !! AXY (26/11/08): implicit detrital carbon change
-               fgmidc(ji,jj)  = xthetad * fgmid(ji,jj)              !! grazing on detrital carbon
-# endif
-               !!
-               !! which translates to these incoming N and C fluxes
-               finmi(ji,jj)   = (1.0 - xphi) * (fgmipn(ji,jj) + fgmid(ji,jj))
-               ficmi(ji,jj)   = (1.0 - xphi) * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))
-               !!
-               !! the ideal food C:N ratio for microzooplankton
-               !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
-               fstarmi(ji,jj) = (xbetan * xthetazmi) / (xbetac * xkc)
-               !!
-               !! process these to determine proportioning of grazed N and C
-               !! (since there is no explicit consideration of respiration,
-               !! only growth and excretion are calculated here)
-               fmith(ji,jj)   = (ficmi(ji,jj) / (finmi(ji,jj) + tiny(finmi(ji,jj))))
-               if (fmith(ji,jj).ge.fstarmi(ji,jj)) then
-                  fmigrow(ji,jj) = xbetan * finmi(ji,jj)
-                  fmiexcr(ji,jj) = 0.0
-               else
-                  fmigrow(ji,jj) = (xbetac * xkc * ficmi(ji,jj)) / xthetazmi
-                  fmiexcr(ji,jj) = ficmi(ji,jj) * ((xbetan / (fmith(ji,jj) + tiny(fmith(ji,jj)))) - ((xbetac * xkc) / xthetazmi))
-               endif
-# if defined key_roam
-               fmiresp(ji,jj) = (xbetac * ficmi(ji,jj)) - (xthetazmi * fmigrow(ji,jj))
-# endif
-
-# if defined key_debug_medusa
-               !! report microzooplankton grazing
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fmi1(',jk,')    = ', fmi1(ji,jj)
-                  IF (lwp) write (numout,*) 'fmi(',jk,')     = ', fmi(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmipn(',jk,')  = ', fgmipn(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmid(',jk,')   = ', fgmid(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmidc(',jk,')  = ', fgmidc(ji,jj)
-                  IF (lwp) write (numout,*) 'finmi(',jk,')   = ', finmi(ji,jj)
-                  IF (lwp) write (numout,*) 'ficmi(',jk,')   = ', ficmi(ji,jj)
-                  IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi(ji,jj)
-                  IF (lwp) write (numout,*) 'fmith(',jk,')   = ', fmith(ji,jj)
-                  IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow(ji,jj)
-                  IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr(ji,jj)
-#  if defined key_roam
-                  IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp(ji,jj)
-#  endif
-               endif
-# endif
-
-               !!----------------------------------------------------------------------
-               !! Mesozooplankton second
-               !!----------------------------------------------------------------------
-               !!
-               fme1(ji,jj)    = (xkme * xkme) + (xpmepn * zphn(ji,jj) * zphn(ji,jj)) + (xpmepd * zphd(ji,jj) * zphd(ji,jj)) + & 
-                         (xpmezmi * zzmi(ji,jj) * zzmi(ji,jj)) + (xpmed * zdet(ji,jj) * zdet(ji,jj))
-               fme(ji,jj)     = xgme * zzme(ji,jj) / fme1(ji,jj)
-               fgmepn(ji,jj)  = fme(ji,jj) * xpmepn  * zphn(ji,jj) * zphn(ji,jj)  !! grazing on non-diatoms
-               fgmepd(ji,jj)  = fme(ji,jj) * xpmepd  * zphd(ji,jj) * zphd(ji,jj)  !! grazing on diatoms
-               fgmepds(ji,jj) = fsin(ji,jj) * fgmepd(ji,jj)                !! grazing on diatom silicon
-               fgmezmi(ji,jj) = fme(ji,jj) * xpmezmi * zzmi(ji,jj) * zzmi(ji,jj)  !! grazing on microzooplankton
-               fgmed(ji,jj)   = fme(ji,jj) * xpmed   * zdet(ji,jj) * zdet(ji,jj)  !! grazing on detrital nitrogen
-# if defined key_roam
-               fgmedc(ji,jj)  = rsmall !acc
-               IF ( zdet(ji,jj) .GT. rsmall ) fgmedc(ji,jj)  = (zdtc(ji,jj) / (zdet(ji,jj) + tiny(zdet(ji,jj)))) * fgmed(ji,jj)  !! grazing on detrital carbon
-# else
-               !! AXY (26/11/08): implicit detrital carbon change
-               fgmedc(ji,jj)  = xthetad * fgmed(ji,jj)              !! grazing on detrital carbon
-# endif
-               !!
-               !! which translates to these incoming N and C fluxes
-               finme(ji,jj)   = (1.0 - xphi) * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))
-               ficme(ji,jj)   = (1.0 - xphi) * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) + &
-                        (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))
-               !!
-               !! the ideal food C:N ratio for mesozooplankton
-               !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
-               fstarme(ji,jj) = (xbetan * xthetazme) / (xbetac * xkc)
-               !!
-               !! process these to determine proportioning of grazed N and C
-               !! (since there is no explicit consideration of respiration,
-               !! only growth and excretion are calculated here)
-               fmeth(ji,jj)   = (ficme(ji,jj) / (finme(ji,jj) + tiny(finme(ji,jj))))
-               if (fmeth(ji,jj).ge.fstarme(ji,jj)) then
-                  fmegrow(ji,jj) = xbetan * finme(ji,jj)
-                  fmeexcr(ji,jj) = 0.0
-               else
-                  fmegrow(ji,jj) = (xbetac * xkc * ficme(ji,jj)) / xthetazme
-                  fmeexcr(ji,jj) = ficme(ji,jj) * ((xbetan / (fmeth(ji,jj) + tiny(fmeth(ji,jj)))) - ((xbetac * xkc) / xthetazme))
-               endif
-# if defined key_roam
-               fmeresp(ji,jj) = (xbetac * ficme(ji,jj)) - (xthetazme * fmegrow(ji,jj))
-# endif
-
-# if defined key_debug_medusa
-               !! report mesozooplankton grazing
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fme1(',jk,')    = ', fme1(ji,jj)
-                  IF (lwp) write (numout,*) 'fme(',jk,')     = ', fme(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmepn(',jk,')  = ', fgmepn(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmepd(',jk,')  = ', fgmepd(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmed(',jk,')   = ', fgmed(ji,jj)
-                  IF (lwp) write (numout,*) 'fgmedc(',jk,')  = ', fgmedc(ji,jj)
-                  IF (lwp) write (numout,*) 'finme(',jk,')   = ', finme(ji,jj)
-                  IF (lwp) write (numout,*) 'ficme(',jk,')   = ', ficme(ji,jj)
-                  IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme(ji,jj)
-                  IF (lwp) write (numout,*) 'fmeth(',jk,')   = ', fmeth(ji,jj)
-                  IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow(ji,jj)
-                  IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr(ji,jj)
-#  if defined key_roam
-                  IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp(ji,jj)
-#  endif
-               endif
-# endif
-
-               fzmi_i(ji,jj)  = fzmi_i(ji,jj)  + fse3t(ji,jj,jk) * (  &
-                  fgmipn(ji,jj) + fgmid(ji,jj) )
-               fzmi_o(ji,jj)  = fzmi_o(ji,jj)  + fse3t(ji,jj,jk) * (  &
-                  fmigrow(ji,jj) + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + fmiexcr(ji,jj) + ((1.0 - xbetan) * finmi(ji,jj)) )
-               fzme_i(ji,jj)  = fzme_i(ji,jj)  + fse3t(ji,jj,jk) * (  &
-                  fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj) )
-               fzme_o(ji,jj)  = fzme_o(ji,jj)  + fse3t(ji,jj,jk) * (  &
-                  fmegrow(ji,jj) + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) + fmeexcr(ji,jj) + ((1.0 - xbetan) * finme(ji,jj)) )
-
-! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PLANKTON LOSSES - marc 20/4/17 
-! (zoo plankton grazing is 152 lines)
-
-               !!----------------------------------------------------------------------
-               !! Plankton metabolic losses
-               !! Linear loss processes assumed to be metabolic in origin
-               !!----------------------------------------------------------------------
-               !!
-               fdpn2(ji,jj)  = xmetapn  * zphn(ji,jj)
-               fdpd2(ji,jj)  = xmetapd  * zphd(ji,jj)
-               fdpds2(ji,jj) = xmetapd  * zpds(ji,jj)
-               fdzmi2(ji,jj) = xmetazmi * zzmi(ji,jj)
-               fdzme2(ji,jj) = xmetazme * zzme(ji,jj)
-
-               !!----------------------------------------------------------------------
-               !! Plankton mortality losses
-               !! EKP (26/02/09): phytoplankton hyperbolic mortality term introduced 
-               !! to improve performance in gyres
-               !!----------------------------------------------------------------------
-               !!
-               !! non-diatom phytoplankton
-               if (jmpn.eq.1) fdpn(ji,jj) = xmpn * zphn(ji,jj)               !! linear
-               if (jmpn.eq.2) fdpn(ji,jj) = xmpn * zphn(ji,jj) * zphn(ji,jj)        !! quadratic
-               if (jmpn.eq.3) fdpn(ji,jj) = xmpn * zphn(ji,jj) * &           !! hyperbolic
-                  (zphn(ji,jj) / (xkphn + zphn(ji,jj)))
-               if (jmpn.eq.4) fdpn(ji,jj) = xmpn * zphn(ji,jj) * &           !! sigmoid
-                  ((zphn(ji,jj) * zphn(ji,jj)) / (xkphn + (zphn(ji,jj) * zphn(ji,jj))))
-               !!
-               !! diatom phytoplankton
-               if (jmpd.eq.1) fdpd(ji,jj) = xmpd * zphd(ji,jj)               !! linear
-               if (jmpd.eq.2) fdpd(ji,jj) = xmpd * zphd(ji,jj) * zphd(ji,jj)        !! quadratic
-               if (jmpd.eq.3) fdpd(ji,jj) = xmpd * zphd(ji,jj) * &           !! hyperbolic
-                  (zphd(ji,jj) / (xkphd + zphd(ji,jj)))
-               if (jmpd.eq.4) fdpd(ji,jj) = xmpd * zphd(ji,jj) * &           !! sigmoid
-                  ((zphd(ji,jj) * zphd(ji,jj)) / (xkphd + (zphd(ji,jj) * zphd(ji,jj))))
-               fdpds(ji,jj) = fdpd(ji,jj) * fsin(ji,jj)
-               !!
-               !! microzooplankton
-               if (jmzmi.eq.1) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj)            !! linear
-               if (jmzmi.eq.2) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * zzmi(ji,jj)     !! quadratic
-               if (jmzmi.eq.3) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * &        !! hyperbolic
-                  (zzmi(ji,jj) / (xkzmi + zzmi(ji,jj)))
-               if (jmzmi.eq.4) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * &        !! sigmoid
-                  ((zzmi(ji,jj) * zzmi(ji,jj)) / (xkzmi + (zzmi(ji,jj) * zzmi(ji,jj))))
-               !!
-               !! mesozooplankton
-               if (jmzme.eq.1) fdzme(ji,jj) = xmzme * zzme(ji,jj)            !! linear
-               if (jmzme.eq.2) fdzme(ji,jj) = xmzme * zzme(ji,jj) * zzme(ji,jj)     !! quadratic
-               if (jmzme.eq.3) fdzme(ji,jj) = xmzme * zzme(ji,jj) * &        !! hyperbolic
-                  (zzme(ji,jj) / (xkzme + zzme(ji,jj)))
-               if (jmzme.eq.4) fdzme(ji,jj) = xmzme * zzme(ji,jj) * &        !! sigmoid
-                  ((zzme(ji,jj) * zzme(ji,jj)) / (xkzme + (zzme(ji,jj) * zzme(ji,jj))))
-
-! MAYBE BUT A BREAK IN HERE, DETRITUS PROCESS- marc 20/4/17 
-! (misc plankton loses is 53 lines - maybe we can link with another section)
 
                !!----------------------------------------------------------------------
@@ -1400,4 +933,12 @@
 ! MAYBE BUT A BREAK IN HERE, IRON CHEMISTRY AND SCAVENGING - marc 20/4/17 
 ! (detritus processes is 74 lines)
+            ENDIF
+         ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
                !!----------------------------------------------------------------------
@@ -1738,4 +1279,12 @@
 ! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PROCESSES - marc 20/4/17 
 ! (iron chemistry and scavenging is 340 lines)
+            ENDIF
+         ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
                !!----------------------------------------------------------------------
@@ -1817,4 +1366,12 @@
 ! MAYBE BUT A BREAK IN HERE, FAST-SINKINIG DETRITUS - marc 20/4/17 
 ! (miscellaneous processes is 79 lines)
+            ENDIF
+         ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
                !!----------------------------------------------------------------------
@@ -2394,4 +1951,12 @@
 ! MAYBE BUT A BREAK IN HERE, BUSINESS AND UPDATING - marc 20/4/17 
 ! (fast sinking detritus is 576 lines)
+            ENDIF
+         ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
                !!======================================================================
@@ -2840,4 +2405,12 @@
 ! (this would make the previous section about 470 lines and the one below
 ! about 700 lines)
+            ENDIF
+         ENDDO
+         ENDDO
+
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
 # if defined key_trc_diabio
@@ -2939,5 +2512,5 @@
                   ENDIF
                   IF( med_diag%PN_FELIM%dgsave ) THEN
-                      ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
+                      ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln2(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PD_JLIM%dgsave ) THEN
@@ -3226,5 +2799,5 @@
                   ENDIF
 		  IF ( med_diag%PNLIMFE3%dgsave ) THEN
-                     pnlimfe3(ji,jj,jk) = ffln(ji,jj)
+                     pnlimfe3(ji,jj,jk) = ffln2(ji,jj)
                   ENDIF
 		  IF ( med_diag%PDLIMJ3%dgsave  ) THEN
@@ -3281,5 +2854,5 @@
                   trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) !! non-diatom J  limitation term 
                   trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom N  limitation term 
-                  trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom Fe limitation term 
+                  trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln2(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom Fe limitation term 
                   trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) !! diatom     J  limitation term 
                   trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     N  limitation term 
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 7975)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 7975)
@@ -0,0 +1,310 @@
+MODULE zooplankton_mod
+   !!======================================================================
+   !!                         ***  MODULE zooplankton_mod  ***
+   !! Calculates the zooplankton grazing
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   zooplankton        ! Called in plankton.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE zooplankton( jk )
+      !!------------------------------------------------------------------
+      !!                     ***  ROUTINE zooplankton  ***
+      !! This called from PLANKTON and calculates the zooplankton 
+      !! grazing
+      !!------------------------------------------------------------------
+      USE bio_medusa_mod,    ONLY: fgmed, fgmedc, fgmepd, fgmepds,       &
+                                   fgmepn, fgmezmi,                      & 
+                                   fgmid, fgmidc, fgmipn,                &
+                                   ficme, ficmi, finme, finmi,           &
+                                   fme, fme1, fmeexcr, fmegrow,          &
+                                   fmeresp, fmeth,                       &
+                                   fmi, fmi1, fmiexcr, fmigrow,          &
+                                   fmiresp, fmith,                       & 
+                                   fsin, fstarme, fstarmi,               &
+                                   fzme_i, fzme_o, fzmi_i, fzmi_o,       &
+                                   idf, idfval,                          &
+                                   zdet, zdtc, zphd, zphn, zzme, zzmi
+      USE dom_oce,           ONLY: e3t_0, e3t_n, tmask
+      USE par_oce,           ONLY: jpim1, jpjm1
+      USE phycst,            ONLY: rsmall
+      USE sms_medusa,        ONLY: xbetac, xbetan, xgme, xgmi,           &
+                                   xkc, xkme, xkmi, xphi,                &
+                                   xpmed, xpmepd, xpmepn, xpmezmi,       &
+                                   xpmid, xpmipn,                        &
+                                   xthetapd, xthetapn,                   &
+                                   xthetazme, xthetazmi
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+
+      !! Level
+      INTEGER, INTENT( in ) :: jk
+
+      INTEGER :: ji, jj
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
+
+               !!----------------------------------------------------------
+               !! Zooplankton Grazing 
+               !! this code supplements the base grazing model with one that
+               !! considers the C:N ratio of grazed food and balances this
+               !! against the requirements of zooplankton growth; this model 
+               !! is derived from that of Anderson & Pondaven (2003)
+               !!
+               !! The current version of the code assumes a fixed C:N ratio 
+               !! for detritus (in contrast to Anderson & Pondaven, 2003), 
+               !! though the full equations are retained for future extension
+               !!----------------------------------------------------------
+               !!
+               !!----------------------------------------------------------
+               !! Microzooplankton first
+               !!----------------------------------------------------------
+               !!
+               fmi1(ji,jj)    = (xkmi * xkmi) + (xpmipn * zphn(ji,jj) *      &
+                                                 zphn(ji,jj)) +              &
+                                (xpmid * zdet(ji,jj) * zdet(ji,jj))
+               fmi(ji,jj)     = xgmi * zzmi(ji,jj) / fmi1(ji,jj)
+               !! grazing on non-diatoms
+               fgmipn(ji,jj)  = fmi(ji,jj) * xpmipn * zphn(ji,jj) *          &
+                                zphn(ji,jj)
+               !! grazing on detrital nitrogen
+               fgmid(ji,jj)   = fmi(ji,jj) * xpmid  * zdet(ji,jj) *          &
+                                zdet(ji,jj)
+# if defined key_roam   
+               ! acc            
+               fgmidc(ji,jj)  = rsmall
+               !! grazing on detrital carbon
+               IF ( zdet(ji,jj) .GT. rsmall ) fgmidc(ji,jj)  =               &
+                  (zdtc(ji,jj) / (zdet(ji,jj) + tiny(zdet(ji,jj)))) *        &
+                  fgmid(ji,jj)
+# else
+               !! AXY (26/11/08): implicit detrital carbon change
+               !! grazing on detrital carbon
+               fgmidc(ji,jj)  = xthetad * fgmid(ji,jj)
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! which translates to these incoming N and C fluxes
+               finmi(ji,jj)   = (1.0 - xphi) * (fgmipn(ji,jj) + fgmid(ji,jj))
+               ficmi(ji,jj)   = (1.0 - xphi) * ((xthetapn * fgmipn(ji,jj)) + &
+                                                fgmidc(ji,jj))
+               !!
+               !! the ideal food C:N ratio for microzooplankton
+               !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
+               fstarmi(ji,jj) = (xbetan * xthetazmi) / (xbetac * xkc)
+               !!
+               !! process these to determine proportioning of grazed N and C
+               !! (since there is no explicit consideration of respiration,
+               !! only growth and excretion are calculated here)
+               fmith(ji,jj)   = (ficmi(ji,jj) / (finmi(ji,jj) +              &
+                                 tiny(finmi(ji,jj))))
+               if (fmith(ji,jj).ge.fstarmi(ji,jj)) then
+                  fmigrow(ji,jj) = xbetan * finmi(ji,jj)
+                  fmiexcr(ji,jj) = 0.0
+               else
+                  fmigrow(ji,jj) = (xbetac * xkc * ficmi(ji,jj)) / xthetazmi
+                  fmiexcr(ji,jj) = ficmi(ji,jj) *                            &
+                                   ((xbetan / (fmith(ji,jj) +                &
+                                               tiny(fmith(ji,jj)))) -        &
+                                    ((xbetac * xkc) / xthetazmi))
+               endif
+# if defined key_roam
+               fmiresp(ji,jj) = (xbetac * ficmi(ji,jj)) -                    &
+                                (xthetazmi * fmigrow(ji,jj))
+# endif
+
+# if defined key_debug_medusa
+               !! report microzooplankton grazing
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'fmi1(',jk,')    = ', fmi1(ji,jj)
+                  IF (lwp) write (numout,*) 'fmi(',jk,')     = ', fmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmipn(',jk,')  = ', fgmipn(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmid(',jk,')   = ', fgmid(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmidc(',jk,')  = ', fgmidc(ji,jj)
+                  IF (lwp) write (numout,*) 'finmi(',jk,')   = ', finmi(ji,jj)
+                  IF (lwp) write (numout,*) 'ficmi(',jk,')   = ', ficmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fmith(',jk,')   = ', fmith(ji,jj)
+                  IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow(ji,jj)
+                  IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr(ji,jj)
+#  if defined key_roam
+                  IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp(ji,jj)
+#  endif
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!----------------------------------------------------------
+               !! Mesozooplankton second
+               !!----------------------------------------------------------
+               !!
+               fme1(ji,jj)    = (xkme * xkme) + (xpmepn * zphn(ji,jj) *       &
+                                                 zphn(ji,jj)) +               &
+                                (xpmepd * zphd(ji,jj) * zphd(ji,jj)) +        & 
+                                (xpmezmi * zzmi(ji,jj) * zzmi(ji,jj)) +       &
+                                (xpmed * zdet(ji,jj) * zdet(ji,jj))
+               fme(ji,jj)     = xgme * zzme(ji,jj) / fme1(ji,jj)
+               !! grazing on non-diatoms
+               fgmepn(ji,jj)  = fme(ji,jj) * xpmepn  * zphn(ji,jj) *          &
+                                zphn(ji,jj)
+               !! grazing on diatoms
+               fgmepd(ji,jj)  = fme(ji,jj) * xpmepd  * zphd(ji,jj) *          &
+                                zphd(ji,jj)
+               !! grazing on diatom silicon
+               fgmepds(ji,jj) = fsin(ji,jj) * fgmepd(ji,jj)
+               !! grazing on microzooplankton
+               fgmezmi(ji,jj) = fme(ji,jj) * xpmezmi * zzmi(ji,jj) *          &
+                                zzmi(ji,jj)
+               !! grazing on detrital nitrogen
+               fgmed(ji,jj)   = fme(ji,jj) * xpmed   * zdet(ji,jj) *          &
+                                zdet(ji,jj)
+# if defined key_roam
+               !! acc
+               fgmedc(ji,jj)  = rsmall
+               !! grazing on detrital carbon
+               IF ( zdet(ji,jj) .GT. rsmall ) fgmedc(ji,jj)  = (zdtc(ji,jj) / &
+                  (zdet(ji,jj) + tiny(zdet(ji,jj)))) * fgmed(ji,jj)
+# else
+               !! AXY (26/11/08): implicit detrital carbon change
+               !! grazing on detrital carbon
+               fgmedc(ji,jj)  = xthetad * fgmed(ji,jj)
+# endif
+               !!
+               !! which translates to these incoming N and C fluxes
+               finme(ji,jj)   = (1.0 - xphi) *                               &
+                                (fgmepn(ji,jj) + fgmepd(ji,jj) +             &
+                                 fgmezmi(ji,jj) + fgmed(ji,jj))
+               ficme(ji,jj)   = (1.0 - xphi) *                               &
+                                ((xthetapn * fgmepn(ji,jj)) +                &
+                                (xthetapd * fgmepd(ji,jj)) +                 &
+                                (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! the ideal food C:N ratio for mesozooplankton
+               !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
+               fstarme(ji,jj) = (xbetan * xthetazme) / (xbetac * xkc)
+               !!
+               !! process these to determine proportioning of grazed N and C
+               !! (since there is no explicit consideration of respiration,
+               !! only growth and excretion are calculated here)
+               fmeth(ji,jj)   = (ficme(ji,jj) / (finme(ji,jj) +              &
+                                                 tiny(finme(ji,jj))))
+               if (fmeth(ji,jj).ge.fstarme(ji,jj)) then
+                  fmegrow(ji,jj) = xbetan * finme(ji,jj)
+                  fmeexcr(ji,jj) = 0.0
+               else
+                  fmegrow(ji,jj) = (xbetac * xkc * ficme(ji,jj)) / xthetazme
+                  fmeexcr(ji,jj) = ficme(ji,jj) *                            &
+                                   ((xbetan / (fmeth(ji,jj) +                &
+                                               tiny(fmeth(ji,jj)))) -        &
+                                    ((xbetac * xkc) / xthetazme))
+               endif
+# if defined key_roam
+               fmeresp(ji,jj) = (xbetac * ficme(ji,jj)) - (xthetazme *       &
+                                                           fmegrow(ji,jj))
+# endif
+
+# if defined key_debug_medusa
+               !! report mesozooplankton grazing
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'fme1(',jk,')    = ', fme1(ji,jj)
+                  IF (lwp) write (numout,*) 'fme(',jk,')     = ', fme(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmepn(',jk,')  = ', fgmepn(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmepd(',jk,')  = ', fgmepd(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmed(',jk,')   = ', fgmed(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmedc(',jk,')  = ', fgmedc(ji,jj)
+                  IF (lwp) write (numout,*) 'finme(',jk,')   = ', finme(ji,jj)
+                  IF (lwp) write (numout,*) 'ficme(',jk,')   = ', ficme(ji,jj)
+                  IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme(ji,jj)
+                  IF (lwp) write (numout,*) 'fmeth(',jk,')   = ', fmeth(ji,jj)
+                  IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow(ji,jj)
+                  IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr(ji,jj)
+#  if defined key_roam
+                  IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp(ji,jj)
+#  endif
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               fzmi_i(ji,jj)  = fzmi_i(ji,jj)  + fse3t(ji,jj,jk) *          &
+                                ( fgmipn(ji,jj) + fgmid(ji,jj) )
+               fzmi_o(ji,jj)  = fzmi_o(ji,jj)  + fse3t(ji,jj,jk) *          &
+                                ( fmigrow(ji,jj) +                          &
+                                  (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + &
+                                  fmiexcr(ji,jj) + ((1.0 - xbetan) *        &
+                                                    finmi(ji,jj)) )
+               fzme_i(ji,jj)  = fzme_i(ji,jj)  + fse3t(ji,jj,jk) *          &
+                                ( fgmepn(ji,jj) + fgmepd(ji,jj) +           &
+                                  fgmezmi(ji,jj) + fgmed(ji,jj) )
+               fzme_o(ji,jj)  = fzme_o(ji,jj)  + fse3t(ji,jj,jk) *          &
+                                ( fmegrow(ji,jj) +                          &
+                                  (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) +  &
+                                  fgmezmi(ji,jj) + fgmed(ji,jj))) +         &
+                                  fmeexcr(ji,jj) + ((1.0 - xbetan) *        &
+                                                    finme(ji,jj)) )
+
+            ENDIF
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE zooplankton
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE zooplankton( )                    ! Empty routine
+      WRITE(*,*) 'zooplankton: You should not have seen this print! error?'
+   END SUBROUTINE zooplankton
+#endif 
+
+   !!======================================================================
+END MODULE zooplankton_mod
