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 7977)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 7978)
@@ -58,4 +58,6 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fsin,fnsi,fprds,fsdiss
 
+   !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ffetop,ffebot,ffescav
    !! Variable for iron-ligand system
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: xFree
@@ -74,4 +76,7 @@
    !! mortality/Remineralisation (defunct parameter "fz" removed)
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
+# if defined key_roam
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fddc
+# endif
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
 
@@ -85,9 +90,14 @@
    !! Particle flux
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdep1,fcaco3
-
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ftempn,ftempsi,ftempfe
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ftempc,ftempca
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: freminn,freminsi,freminfe
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: freminc,freminca
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ffastn,ffastsi,ffastfe
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ffastc,ffastca
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fprotf
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fsedn,fsedsi,fsedfe,fsedc,fsedca
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fccd
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fccd_dep
 
    !! water column nutrient and flux integrals
@@ -269,4 +279,5 @@
                fsin(jpi,jpj),fnsi(jpi,jpj),fprds(jpi,jpj),            &
                fsdiss(jpi,jpj),                                       &
+               ffetop(jpi,jpj),ffebot(jpi,jpj),ffescav(jpi,jpj),      &
                xFree(jpi,jpj),                                        &
                fmi1(jpi,jpj),fmi(jpi,jpj),fgmipn(jpi,jpj),            &
@@ -283,4 +294,5 @@
                fdpn(jpi,jpj),fdpd(jpi,jpj),fdpds(jpi,jpj),            &
                fdzmi(jpi,jpj),fdzme(jpi,jpj),fdd(jpi,jpj),            &
+               fddc(jpi,jpj),                                         &
                fdpn2(jpi,jpj),fdpd2(jpi,jpj),fdpds2(jpi,jpj),         &
                fdzmi2(jpi,jpj),fdzme2(jpi,jpj),                       &
@@ -291,9 +303,15 @@
 # endif
                fdep1(jpi,jpj),fcaco3(jpi,jpj),                        &
+               ftempn(jpi,jpj),ftempsi(jpi,jpj),ftempfe(jpi,jpj),     &
+               ftempc(jpi,jpj),ftempca(jpi,jpj),                      &
+               freminn(jpi,jpj),freminsi(jpi,jpj),freminfe(jpi,jpj),  &
+               freminc(jpi,jpj),freminca(jpi,jpj),                    &
                ffastn(jpi,jpj),ffastsi(jpi,jpj),ffastfe(jpi,jpj),     &
                ffastc(jpi,jpj),ffastca(jpi,jpj),                      &
+               fprotf(jpi,jpj),                                       &
                fsedn(jpi,jpj),fsedsi(jpi,jpj),fsedfe(jpi,jpj),        &
                fsedc(jpi,jpj),fsedca(jpi,jpj),                        &
                fccd(jpi,jpj),                                         &
+               fccd_dep(jpi,jpj),                                     &
                ftot_n(jpi,jpj),ftot_si(jpi,jpj),ftot_fe(jpi,jpj),     &
                fflx_n(jpi,jpj),fflx_si(jpi,jpj),fflx_fe(jpi,jpj),     &
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90	(revision 7978)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90	(revision 7978)
@@ -0,0 +1,487 @@
+MODULE iron_chem_scav_mod
+   !!======================================================================
+   !!                         ***  MODULE iron_chem_scav_mod  ***
+   !! Calculate the iron chemistry and scavenging.
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   iron_chem_scav        ! 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 iron_chem_scav( jk )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE iron_chem_scav  ***
+      !! This called from TRC_BIO_MEDUSA and 
+      !!  - 
+      !!-------------------------------------------------------------------
+      USE bio_medusa_mod,    ONLY: ffastc, ffastca, ffastsi,              &
+                                   ffetop, ffebot, ffescav, xfree,        & 
+                                   zdet, zfer, zphd, zphn, zzme, zzmi 
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,        &
+                                   mbathy, tmask
+      USE par_kind,          ONLY: wp
+      USE par_oce,           ONLY: jpim1, jpjm1
+      USE sms_medusa,        ONLY: i0500, jiron, xfe_sed, xfe_sol,        &
+                                   xk_FeL, xk_sc_Fe, xLgT,                &
+                                   xmassc, xmassca, xmasssi,              &
+                                   xthetad, xthetapd, xthetapn,           &
+                                   xthetazme, xthetazmi,                  &
+                                   zirondep
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+
+      !! Level
+      INTEGER, INTENT( in ) :: jk
+
+      !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
+      !! state variables for iron-ligand system
+      REAL(wp) ::    xLgF, xFeT, xFeF, xFeL
+      !! iron-ligand parameters
+      REAL(wp) ::    xb_coef_tmp, xb2M4ac
+      !! max Fe' parameters
+      REAL(wp) ::    xmaxFeF,fdeltaFe
+      !!
+      !! local parameters for Moore et al. (2004) alternative scavenging 
+      !! scheme
+      REAL(wp) ::    fbase_scav,fscal_sink,fscal_part,fscal_scav
+      !!
+      !! local parameters for Moore et al. (2008) alternative scavenging 
+      !! scheme
+      REAL(wp) ::    fscal_csink,fscal_sisink,fscal_casink
+      !!
+      !! local parameters for Galbraith et al. (2010) alternative 
+      !! scavenging scheme.
+      !! organic portion of scavenging
+      REAL(wp) ::    xCscav1, xCscav2, xk_org, xORGscav
+      !! inorganic portion of scavenging
+      REAL(wp) ::    xk_inorg, xINORGscav
+
+      INTEGER :: ji, jj
+
+! I'D LIKE TO PULL THE IF (jiron) statements outside the DO loops - marc
+
+      !!------------------------------------------------------------------
+      !! Iron chemistry and fractionation
+      !! following the Parekh et al. (2004) scheme adopted by the Met.
+      !! Office, Medusa models total iron but considers "free" and
+      !! ligand-bound forms for the purposes of scavenging (only "free"
+      !! iron can be scavenged
+      !!------------------------------------------------------------------
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
+               xFeT        = zfer(ji,jj) * 1.e3
+               !!
+               !! calculate fractionation (based on Diat-HadOCC; in turn 
+               !! based on Parekh et al., 2004)
+               xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0
+               xb2M4ac     = max(((xb_coef_tmp * xb_coef_tmp) +              &
+                                  (4.0 * xk_FeL * xLgT)), 0.0)
+               !!
+               !! "free" ligand concentration
+               xLgF        = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL
+               !!
+               !! ligand-bound iron concentration
+               xFeL        = xLgT - xLgF
+               !!
+               !! "free" iron concentration (and convert to mmol Fe / m3)
+               xFeF        = (xFeT - xFeL) * 1.e-3
+               xFree(ji,jj)= xFeF / (zfer(ji,jj) + tiny(zfer(ji,jj)))
+               !!
+               !! scavenging of iron (multiple schemes); I'm only really 
+               !! happy with the first one at the moment - the others 
+               !! involve assumptions (sometimes guessed at by me) that 
+               !! are potentially questionable
+               !!
+               if (jiron.eq.1) then
+                  !!------------------------------------------------------
+                  !! Scheme 1: Dutkiewicz et al. (2005)
+                  !! This scheme includes a single scavenging term based 
+                  !! solely on a fixed rate and the availablility of 
+                  !! "free" iron
+                  !!------------------------------------------------------
+                  !! = mmol/m3/d
+                  ffescav(ji,jj)     = xk_sc_Fe * xFeF
+                  !!
+                  !!------------------------------------------------------
+                  !!
+                  !! Mick's code contains a further (optional) implicit 
+                  !! "scavenging" of iron that sets an upper bound on 
+                  !! "free" iron concentration, and essentially caps the 
+                  !! concentration of total iron as xFeL + "free" iron; 
+                  !! since the former is constrained by a fixed total 
+                  !! ligand concentration (= 1.0 umol/m3), and the latter 
+                  !! isn't allowed above this upper bound, total iron is 
+                  !! constrained to a maximum of ...
+                  !!
+                  !!    xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 
+                  !!                                  = 1.3 umol / m3
+                  !! 
+                  !! In Mick's code, the actual value of total iron is 
+                  !! reset to this sum (i.e. TFe = FeL + Fe'; but 
+                  !! Fe' <= 0.3 umol/m3); this isn't our favoured approach 
+                  !! to tracer updating here (not least because of the 
+                  !! leapfrog), so here the amount scavenged is augmented 
+                  !! by an additional amount that serves to drag total 
+                  !! iron back towards that expected from this limitation 
+                  !! on iron concentration ...
+                  !!
+                  !! = umol/m3
+                  xmaxFeF     = min((xFeF * 1.e3), 0.3)
+                  !!
+                  !! Here, the difference between current total Fe and 
+                  !! (FeL + Fe') is calculated and added to the scavenging 
+                  !! flux already calculated above ...
+                  !!
+                  !! = mmol/m3
+                  fdeltaFe    = (xFeT - (xFeL + xmaxFeF)) * 1.e-3
+                  !!
+                  !! This assumes that the "excess" iron is dissipated 
+                  !! with a time-scale of 1 day; seems reasonable to me 
+                  !! ... (famous last words)
+                  !!
+                  !! = mmol/m3/d
+                  ffescav(ji,jj)     = ffescav(ji,jj) + fdeltaFe
+                  !!
+# if defined key_deep_fe_fix
+                  !! AXY (17/01/13)
+                  !! stop scavenging for iron concentrations below 
+                  !! 0.5 umol / m3 at depths greater than 1000 m; this 
+                  !! aims to end MEDUSA's continual loss of iron at depth 
+                  !! without impacting things at the surface too much; the 
+                  !! justification for this is that it appears to be what 
+                  !! Mick Follows et al. do in their work (as evidenced by 
+                  !! the iron initial condition they supplied me with); to 
+                  !! be honest, it looks like Follow et al. do this at
+                  !! shallower depths than 1000 m, but I'll stick with this
+                  !! for now; I suspect that this seemingly arbitrary
+                  !! approach effectively "parameterises" the 
+                  !! particle-based scavenging rates that other models use 
+                  !! (i.e. at depth there are no sinking particles, so 
+                  !! scavenging stops); it might be fun justifying this in 
+                  !! a paper though!
+                  !!
+                  if ((fsdepw(ji,jj,jk).gt.1000.) .and. (xFeT.lt.0.5)) then
+                     ffescav(ji,jj) = 0.
+                  endif
+# endif
+                  !!
+               elseif (jiron.eq.2) then
+                  !!------------------------------------------------------
+                  !! Scheme 2: Moore et al. (2004)
+                  !! This scheme includes a single scavenging term that 
+                  !! accounts for both suspended and sinking particles in 
+                  !! the water column; this term scavenges total iron rather 
+                  !! than "free" iron
+                  !!------------------------------------------------------
+                  !!
+                  !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
+                  xFeT        = zfer(ji,jj) * 1.e3
+                  !!
+                  !! this has a base scavenging rate (12% / y) which is 
+                  !! modified by local particle concentration and sinking 
+                  !! flux (and dust - but I'm ignoring that here for now) 
+                  !! and which is accelerated when Fe concentration gets
+                  !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased 
+                  !! as concentrations below 0.4 nM (= 0.4 umol/m3 = 
+                  !! 0.0004 mmol/m3)
+                  !!
+                  !! base scavenging rate (0.12 / y)
+                  fbase_scav = 0.12 / 365.25
+                  !!
+                  !! calculate sinking particle part of scaling factor
+                  !! this takes local fast sinking carbon (mmol C / m2 / d)
+                  !! and gets it into nmol C / cm3 / s ("rdt" below is the 
+                  !! number of seconds in a model timestep)
+                  !!
+                  !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.)
+                  !!
+                  !! ... actually, re-reading Moore et al.'s equations, it 
+                  !! looks like he uses his sinking flux directly, without 
+                  !! scaling it by time-step or anything, so I'll copy this 
+                  !! here ...
+                  !!
+                  fscal_sink = ffastc(ji,jj) * 1.e2
+                  !!
+                  !! calculate particle part of scaling factor
+                  !! this totals up the carbon in suspended particles 
+                  !! (Pn, Pd, Zmi, Zme, D),
+                  !! which comes out in mmol C / m3 (= nmol C / cm3), and 
+                  !! then multiplies it by a magic factor, 0.002, to get it 
+                  !! into nmol C / cm2 / s
+                  !!
+                  fscal_part = ( (xthetapn * zphn(ji,jj)) +                  &
+                                 (xthetapd * zphd(ji,jj)) +                  &
+                                 (xthetazmi * zzmi(ji,jj)) +                 &
+                                 (xthetazme * zzme(ji,jj)) +                 &
+                                 (xthetad * zdet(ji,jj)) ) * 0.002
+                  !!
+                  !! calculate scaling factor for base scavenging rate
+                  !! this uses the (now correctly scaled) sinking flux and 
+                  !! standing
+                  !! particle concentration, divides through by some sort 
+                  !! of reference value (= 0.0066 nmol C / cm2 / s) and 
+                  !! then uses this, or not if its too high, to rescale the 
+                  !! base scavenging rate
+                  !!
+                  fscal_scav = fbase_scav *                                  &
+                               min(((fscal_sink + fscal_part) / 0.0066), 4.0)
+                  !!
+                  !! the resulting scavenging rate is then scaled further 
+                  !! according to the local iron concentration (i.e. 
+                  !! diminished in low iron regions; enhanced in high iron 
+                  !! regions; less alone in intermediate iron regions)
+                  !!
+                  if (xFeT.lt.0.4) then
+                     !!
+                     !! low iron region
+                     !!
+                     fscal_scav = fscal_scav * (xFeT / 0.4)
+                     !!
+                  elseif (xFeT.gt.0.6) then
+                     !!
+                     !! high iron region
+                     !!
+                     fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4))
+                     !!
+                  else
+                     !!
+                     !! intermediate iron region: do nothing
+                     !!
+                  endif
+                  !! 
+                  !! apply the calculated scavenging rate ...
+                  !!
+                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
+                  !!
+               elseif (jiron.eq.3) then
+                  !!------------------------------------------------------
+                  !! Scheme 3: Moore et al. (2008)
+                  !! This scheme includes a single scavenging term that 
+                  !! accounts for sinking particles in the water column, 
+                  !! and includes organic C, biogenic opal, calcium 
+                  !! carbonate and dust in this (though the latter is 
+                  !! ignored here until I work out what units the incoming
+                  !! "dust" flux is in); this term scavenges total iron 
+                  !! rather than "free" iron
+                  !!------------------------------------------------------
+                  !!
+                  !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
+                  xFeT        = zfer(ji,jj) * 1.e3
+                  !!
+                  !! this has a base scavenging rate which is modified by 
+                  !! local particle sinking flux (including dust - but I'm 
+                  !! ignoring that here for now) and which is accelerated 
+                  !! when Fe concentration is > 0.6 nM (= 0.6 umol/m3 = 
+                  !! 0.0006 mmol/m3), and decreased as concentrations < 
+                  !! 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3)
+                  !!
+                  !! base scavenging rate (Fe_b in paper; units may be 
+                  !! wrong there)
+                  fbase_scav = 0.00384 ! (ng)^-1 cm
+                  !!
+                  !! calculate sinking particle part of scaling factor; 
+                  !! this converts mmol / m2 / d fluxes of organic carbon, 
+                  !! silicon and calcium carbonate into ng / cm2 / s 
+                  !! fluxes; it is assumed here that the mass conversions 
+                  !! simply consider the mass of the main element
+                  !! (C, Si and Ca) and *not* the mass of the molecules 
+                  !! that they are part of; Moore et al. (2008) is unclear 
+                  !! on the conversion that should be used
+                  !!
+                  !! milli -> nano; mol -> gram; /m2 -> /cm2; /d -> /s
+                  !! ng C  / cm2 / s
+                  fscal_csink  = (ffastc(ji,jj)  * 1.e6 * xmassc  *           &
+                                  1.e-4 / 86400.)
+                  !! ng Si / cm2 / s
+                  fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi *           &
+                                  1.e-4 / 86400.)
+                  !! ng Ca / cm2 / s
+                  fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca *           &
+                                  1.e-4 / 86400.)
+                  !! 
+                  !! sum up these sinking fluxes and convert to ng / cm 
+                  !! by dividing through by a sinking rate of 
+                  !! 100 m / d = 1.157 cm / s
+                  !! ng / cm
+                  fscal_sink   = ((fscal_csink * 6.) + fscal_sisink +        &
+                                  fscal_casink) / (100. * 1.e3 / 86400)
+                  !!
+                  !! now calculate the scavenging rate based upon the base 
+                  !! rate and this particle flux scaling; according to the 
+                  !! published units, the result actually has *no* units, 
+                  !! but as it must be expressed per unit time for it to 
+                  !! make any sense, I'm assuming a missing "per second"
+                  !! / s
+                  fscal_scav = fbase_scav * fscal_sink
+                  !!
+                  !! the resulting scavenging rate is then scaled further 
+                  !! according to the local iron concentration (i.e. 
+                  !! diminished in low iron regions; enhanced in high iron 
+                  !! regions; less alone in intermediate iron regions)
+                  !!
+                  if (xFeT.lt.0.5) then
+                     !!
+                     !! low iron region (0.5 instead of the 0.4 in Moore 
+                     !! et al., 2004)
+                     !!
+                     fscal_scav = fscal_scav * (xFeT / 0.5)
+                     !!
+                  elseif (xFeT.gt.0.6) then
+                     !!
+                     !! high iron region (functional form different in 
+                     !! Moore et al., 2004)
+                     !!
+                     fscal_scav = fscal_scav + ((xFeT - 0.6) * 0.00904)
+                     !!
+                  else
+                     !!
+                     !! intermediate iron region: do nothing
+                     !!
+                  endif
+                  !! 
+                  !! apply the calculated scavenging rate ...
+                  !!
+                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
+                  !!
+               elseif (jiron.eq.4) then
+                  !!------------------------------------------------------
+                  !! Scheme 4: Galbraith et al. (2010)
+                  !! This scheme includes two scavenging terms, one for 
+                  !! organic, particle-based scavenging, and another for 
+                  !! inorganic scavenging; both terms scavenge "free" iron 
+                  !! only
+                  !!------------------------------------------------------
+                  !!
+                  !! Galbraith et al. (2010) present a more straightforward 
+                  !! outline of the scheme in Parekh et al. (2005) ...
+                  !! 
+                  !! sinking particulate carbon available for scavenging
+                  !! this assumes a sinking rate of 100 m / d (Moore & 
+                  !! Braucher, 2008),
+                  xCscav1     = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3
+                  !! 
+                  !! scale by Honeyman et al. (1981) exponent coefficient
+                  !! multiply by 1.e-3 to express C flux in g C rather than 
+                  !! mg C
+                  xCscav2     = (xCscav1 * 1.e-3)**0.58
+                  !!
+                  !! multiply by Galbraith et al. (2010) scavenging rate
+                  xk_org      = 0.5 ! ((g C m/3)^-1) / d
+                  xORGscav    = xk_org * xCscav2 * xFeF
+                  !!
+                  !! Galbraith et al. (2010) also include an inorganic bit ...
+                  !! 
+                  !! this occurs at a fixed rate, again based on the 
+                  !! availability of "free" iron
+                  !!
+                  !! k_inorg  = 1000 d**-1 nmol Fe**-0.5 kg**-0.5
+                  !!
+                  !! to implement this here, scale xFeF by 1026 to put in 
+                  !! units of umol/m3 which approximately equal nmol/kg
+                  !!
+                  xk_inorg    = 1000. ! ((nmol Fe / kg)^1.5)
+                  xINORGscav  = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e-3
+                  !!
+                  !! sum these two terms together
+                  ffescav(ji,jj)     = xORGscav + xINORGscav
+               else
+                  !!------------------------------------------------------
+                  !! No Scheme: you coward!
+                  !! This scheme puts its head in the sand and eskews any 
+                  !! decision about how iron is removed from the ocean; 
+                  !! prepare to get deluged in iron you fool!
+                  !!------------------------------------------------------
+                  ffescav(ji,jj)     = 0.
+               endif
+
+               !!---------------------------------------------------------
+               !! Other iron cycle processes
+               !!---------------------------------------------------------
+               !!
+               !! aeolian iron deposition
+               if (jk.eq.1) then
+                  !! zirondep   is in mmol-Fe / m2 / day
+                  !! ffetop(ji,jj)     is in mmol-dissolved-Fe / m3 / day
+                  ffetop(ji,jj)  = zirondep(ji,jj) * xfe_sol / fse3t(ji,jj,jk) 
+               else
+                  ffetop(ji,jj)  = 0.0
+               endif
+               !!
+               !! seafloor iron addition
+               !! AXY (10/07/12): amended to only apply sedimentary flux up 
+               !! to ~500 m down
+               !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then
+               if ((jk.eq.mbathy(ji,jj)).AND.jk.le.i0500) then
+                  !! Moore et al. (2004) cite a coastal California value of 
+                  !! 5 umol/m2/d, but adopt a global value of 2 umol/m2/d 
+                  !! for all areas < 1100 m; here we use this latter value
+                  !! but apply it everywhere
+                  !! AXY (21/07/09): actually, let's just apply it below 
+                  !! 1100 m (levels 1-37)
+                  ffebot(ji,jj)  = (xfe_sed / fse3t(ji,jj,jk))
+               else
+                  ffebot(ji,jj)  = 0.0
+               endif
+
+               !! AXY (16/12/09): remove iron addition/removal processes
+               !! For the purposes of the quarter degree run, the iron 
+               !! cycle is being pegged to the initial condition supplied 
+               !! by Mick Follows via restoration with a 30 day period;
+               !! iron addition at the seafloor is still permitted with 
+               !! the idea that this extra iron will be removed by the 
+               !! restoration away from the source
+               !! ffescav(ji,jj) = 0.0
+               !! ffetop(ji,jj)  = 0.0
+               !! ffebot(ji,jj)  = 0.0
+
+# if defined key_debug_medusa
+               !! report miscellaneous calculations
+               if (idf.eq.1.AND.idfval.eq.1) then
+                  IF (lwp) write (numout,*) '------------------------------'
+                  IF (lwp) write (numout,*) 'xfe_sol  = ', xfe_sol
+                  IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass
+                  IF (lwp) write (numout,*) 'ffetop(',jk,')  = ', ffetop(ji,jj)
+                  IF (lwp) write (numout,*) 'ffebot(',jk,')  = ', ffebot(ji,jj)
+                  IF (lwp) write (numout,*) 'xFree(',jk,')   = ', xFree(ji,jj)
+                  IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav(ji,jj)
+               endif
+# endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE iron_chem_scav
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE iron_chem_scav( )                    ! Empty routine
+      WRITE(*,*) 'iron_chem_scav: You should not have seen this print! error?'
+   END SUBROUTINE iron_chem_scav
+#endif 
+
+   !!======================================================================
+END MODULE iron_chem_scav_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 7977)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7978)
@@ -101,4 +101,5 @@
       USE air_sea_mod,                ONLY: air_sea
       USE plankton_mod,               ONLY: plankton
+      USE iron_chem_scav_mod,         ONLY: iron_chem_scav
       USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
       USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
@@ -207,5 +208,5 @@
       !!
       !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
-      REAL(wp), DIMENSION(jpi,jpj) ::    ffetop,ffebot,ffescav
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ffetop,ffebot,ffescav
       REAL(wp) ::    xLgF, xFeT, xFeF, xFeL         !! state variables for iron-ligand system
 !      REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system
@@ -234,5 +235,5 @@
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
 # if defined key_roam
-      REAL(wp), DIMENSION(jpi,jpj) ::    fddc
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fddc
 # endif
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
@@ -249,11 +250,11 @@
       !! particle flux
 !      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
+!      REAL(WP), DIMENSION(jpi,jpj) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
+!      REAL(wp), DIMENSION(jpi,jpj) ::    freminn,freminsi,freminfe,freminc,freminca
 !      REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca
-      REAL(wp), DIMENSION(jpi,jpj) ::    fprotf
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fprotf
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca
 !      REAL(wp), DIMENSION(jpi,jpj) ::    fccd
-      REAL(wp), DIMENSION(jpi,jpj) ::    fccd_dep
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fccd_dep
       !!
       !! AXY (06/07/11): alternative fast detritus schemes
@@ -937,349 +938,10 @@
          ENDDO
 
-         DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,jk) == 1) then
-
-               !!----------------------------------------------------------------------
-               !! Iron chemistry and fractionation
-               !! following the Parekh et al. (2004) scheme adopted by the Met.
-               !! Office, Medusa models total iron but considers "free" and
-               !! ligand-bound forms for the purposes of scavenging (only "free"
-               !! iron can be scavenged
-               !!----------------------------------------------------------------------
-               !!
-               !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
-               xFeT        = zfer(ji,jj) * 1.e3
-               !!
-               !! calculate fractionation (based on Diat-HadOCC; in turn based on Parekh et al., 2004)
-               xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0
-               xb2M4ac     = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0)
-               !!
-               !! "free" ligand concentration
-               xLgF        = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL
-               !!
-               !! ligand-bound iron concentration
-               xFeL        = xLgT - xLgF
-               !!
-               !! "free" iron concentration (and convert to mmol Fe / m3)
-               xFeF        = (xFeT - xFeL) * 1.e-3
-               xFree(ji,jj)= xFeF / (zfer(ji,jj) + tiny(zfer(ji,jj)))
-               !!
-               !! scavenging of iron (multiple schemes); I'm only really happy with the 
-               !! first one at the moment - the others involve assumptions (sometimes
-               !! guessed at by me) that are potentially questionable
-               !!
-               if (jiron.eq.1) then
-                  !!----------------------------------------------------------------------
-                  !! Scheme 1: Dutkiewicz et al. (2005)
-                  !! This scheme includes a single scavenging term based solely on a
-                  !! fixed rate and the availablility of "free" iron
-                  !!----------------------------------------------------------------------
-                  !!
-                  ffescav(ji,jj)     = xk_sc_Fe * xFeF                     ! = mmol/m3/d
-                  !!
-                  !!----------------------------------------------------------------------
-                  !!
-                  !! Mick's code contains a further (optional) implicit "scavenging" of 
-                  !! iron that sets an upper bound on "free" iron concentration, and 
-                  !! essentially caps the concentration of total iron as xFeL + "free" 
-                  !! iron; since the former is constrained by a fixed total ligand 
-                  !! concentration (= 1.0 umol/m3), and the latter isn't allowed above 
-                  !! this upper bound, total iron is constrained to a maximum of ...
-                  !!
-                  !!    xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3
-                  !! 
-                  !! In Mick's code, the actual value of total iron is reset to this
-                  !! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't
-                  !! our favoured approach to tracer updating here (not least because
-                  !! of the leapfrog), so here the amount scavenged is augmented by an
-                  !! additional amount that serves to drag total iron back towards that
-                  !! expected from this limitation on iron concentration ...
-                  !!
-                  xmaxFeF     = min((xFeF * 1.e3), 0.3)             ! = umol/m3
-                  !!
-                  !! Here, the difference between current total Fe and (FeL + Fe') is
-                  !! calculated and added to the scavenging flux already calculated
-                  !! above ...
-                  !!
-                  fdeltaFe    = (xFeT - (xFeL + xmaxFeF)) * 1.e-3   ! = mmol/m3
-                  !!
-                  !! This assumes that the "excess" iron is dissipated with a time-
-                  !! scale of 1 day; seems reasonable to me ... (famous last words)
-                  !!
-                  ffescav(ji,jj)     = ffescav(ji,jj) + fdeltaFe                  ! = mmol/m3/d
-                  !!
-# if defined key_deep_fe_fix
-                  !! AXY (17/01/13)
-                  !! stop scavenging for iron concentrations below 0.5 umol / m3
-                  !! at depths greater than 1000 m; this aims to end MEDUSA's
-                  !! continual loss of iron at depth without impacting things
-                  !! at the surface too much; the justification for this is that
-                  !! it appears to be what Mick Follows et al. do in their work
-                  !! (as evidenced by the iron initial condition they supplied
-                  !! me with); to be honest, it looks like Follow et al. do this
-                  !! at shallower depths than 1000 m, but I'll stick with this
-                  !! for now; I suspect that this seemingly arbitrary approach
-                  !! effectively "parameterises" the particle-based scavenging
-                  !! rates that other models use (i.e. at depth there are no
-                  !! sinking particles, so scavenging stops); it might be fun
-                  !! justifying this in a paper though!
-                  !!
-                  if ((fsdepw(ji,jj,jk).gt.1000.) .and. (xFeT.lt.0.5)) then
-                     ffescav(ji,jj) = 0.
-                  endif
-# endif
-                  !!
-               elseif (jiron.eq.2) then
-                  !!----------------------------------------------------------------------
-                  !! Scheme 2: Moore et al. (2004)
-                  !! This scheme includes a single scavenging term that accounts for
-                  !! both suspended and sinking particles in the water column; this
-                  !! term scavenges total iron rather than "free" iron
-                  !!----------------------------------------------------------------------
-                  !!
-                  !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
-                  xFeT        = zfer(ji,jj) * 1.e3
-                  !!
-                  !! this has a base scavenging rate (12% / y) which is modified by local
-                  !! particle concentration and sinking flux (and dust - but I'm ignoring
-                  !! that here for now) and which is accelerated when Fe concentration gets
-                  !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as concentrations
-                  !! below 0.4 nM (= 0.4 umol/m3 = 0.0004 mmol/m3)
-                  !!
-                  !! base scavenging rate (0.12 / y)
-                  fbase_scav = 0.12 / 365.25
-                  !!
-                  !! calculate sinking particle part of scaling factor
-                  !! this takes local fast sinking carbon (mmol C / m2 / d) and
-                  !! gets it into nmol C / cm3 / s ("rdt" below is the number of seconds in
-                  !! a model timestep)
-                  !!
-                  !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.)
-                  !!
-                  !! ... actually, re-reading Moore et al.'s equations, it looks like he uses
-                  !! his sinking flux directly, without scaling it by time-step or anything,
-                  !! so I'll copy this here ...
-                  !!
-                  fscal_sink = ffastc(ji,jj) * 1.e2
-                  !!
-                  !! calculate particle part of scaling factor
-                  !! this totals up the carbon in suspended particles (Pn, Pd, Zmi, Zme, D),
-                  !! which comes out in mmol C / m3 (= nmol C / cm3), and then multiplies it
-                  !! by a magic factor, 0.002, to get it into nmol C / cm2 / s
-                  !!
-                  fscal_part = ((xthetapn * zphn(ji,jj)) + (xthetapd * zphd(ji,jj)) + (xthetazmi * zzmi(ji,jj)) + &
-                  (xthetazme * zzme(ji,jj)) + (xthetad * zdet(ji,jj))) * 0.002
-                  !!
-                  !! calculate scaling factor for base scavenging rate
-                  !! this uses the (now correctly scaled) sinking flux and standing
-                  !! particle concentration, divides through by some sort of reference
-                  !! value (= 0.0066 nmol C / cm2 / s) and then uses this, or not if its
-                  !! too high, to rescale the base scavenging rate
-                  !!
-                  fscal_scav = fbase_scav * min(((fscal_sink + fscal_part) / 0.0066), 4.0)
-                  !!
-                  !! the resulting scavenging rate is then scaled further according to the
-                  !! local iron concentration (i.e. diminished in low iron regions; enhanced
-                  !! in high iron regions; less alone in intermediate iron regions)
-                  !!
-                  if (xFeT.lt.0.4) then
-                     !!
-                     !! low iron region
-                     !!
-                     fscal_scav = fscal_scav * (xFeT / 0.4)
-                     !!
-                  elseif (xFeT.gt.0.6) then
-                     !!
-                     !! high iron region
-                     !!
-                     fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4))
-                     !!
-                  else
-                     !!
-                     !! intermediate iron region: do nothing
-                     !!
-                  endif
-                  !! 
-                  !! apply the calculated scavenging rate ...
-                  !!
-                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
-                  !!
-               elseif (jiron.eq.3) then
-                  !!----------------------------------------------------------------------
-                  !! Scheme 3: Moore et al. (2008)
-                  !! This scheme includes a single scavenging term that accounts for
-                  !! sinking particles in the water column, and includes organic C,
-                  !! biogenic opal, calcium carbonate and dust in this (though the
-                  !! latter is ignored here until I work out what units the incoming
-                  !! "dust" flux is in); this term scavenges total iron rather than 
-                  !! "free" iron
-                  !!----------------------------------------------------------------------
-                  !!
-                  !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
-                  xFeT        = zfer(ji,jj) * 1.e3
-                  !!
-                  !! this has a base scavenging rate which is modified by local
-                  !! particle sinking flux (including dust - but I'm ignoring that 
-                  !! here for now) and which is accelerated when Fe concentration
-                  !! is > 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as 
-                  !! concentrations < 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3)
-                  !!
-                  !! base scavenging rate (Fe_b in paper; units may be wrong there)
-                  fbase_scav = 0.00384 ! (ng)^-1 cm
-                  !!
-                  !! calculate sinking particle part of scaling factor; this converts
-                  !! mmol / m2 / d fluxes of organic carbon, silicon and calcium
-                  !! carbonate into ng / cm2 / s fluxes; it is assumed here that the
-                  !! mass conversions simply consider the mass of the main element
-                  !! (C, Si and Ca) and *not* the mass of the molecules that they are
-                  !! part of; Moore et al. (2008) is unclear on the conversion that
-                  !! should be used
-                  !!
-                  !! milli -> nano; mol -> gram; /m2 -> /cm2; /d -> /s
-                  fscal_csink  = (ffastc(ji,jj)  * 1.e6 * xmassc  * 1.e-4 / 86400.)      ! ng C  / cm2 / s
-                  fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi * 1.e-4 / 86400.)      ! ng Si / cm2 / s
-                  fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca * 1.e-4 / 86400.)      ! ng Ca / cm2 / s
-                  !! 
-                  !! sum up these sinking fluxes and convert to ng / cm by dividing
-                  !! through by a sinking rate of 100 m / d = 1.157 cm / s
-                  fscal_sink   = ((fscal_csink * 6.) + fscal_sisink + fscal_casink) / &
-                  (100. * 1.e3 / 86400)                                                  ! ng / cm
-                  !!
-                  !! now calculate the scavenging rate based upon the base rate and
-                  !! this particle flux scaling; according to the published units,
-                  !! the result actually has *no* units, but as it must be expressed
-                  !! per unit time for it to make any sense, I'm assuming a missing
-                  !! "per second"
-                  fscal_scav = fbase_scav * fscal_sink                                   ! / s
-                  !!
-                  !! the resulting scavenging rate is then scaled further according to the
-                  !! local iron concentration (i.e. diminished in low iron regions; enhanced
-                  !! in high iron regions; less alone in intermediate iron regions)
-                  !!
-                  if (xFeT.lt.0.5) then
-                     !!
-                     !! low iron region (0.5 instead of the 0.4 in Moore et al., 2004)
-                     !!
-                     fscal_scav = fscal_scav * (xFeT / 0.5)
-                     !!
-                  elseif (xFeT.gt.0.6) then
-                     !!
-                     !! high iron region (functional form different in Moore et al., 2004)
-                     !!
-                     fscal_scav = fscal_scav + ((xFeT - 0.6) * 0.00904)
-                     !!
-                  else
-                     !!
-                     !! intermediate iron region: do nothing
-                     !!
-                  endif
-                  !! 
-                  !! apply the calculated scavenging rate ...
-                  !!
-                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
-                  !!
-               elseif (jiron.eq.4) then
-                  !!----------------------------------------------------------------------
-                  !! Scheme 4: Galbraith et al. (2010)
-                  !! This scheme includes two scavenging terms, one for organic,
-                  !! particle-based scavenging, and another for inorganic scavenging;
-                  !! both terms scavenge "free" iron only
-                  !!----------------------------------------------------------------------
-                  !!
-                  !! Galbraith et al. (2010) present a more straightforward outline of 
-                  !! the scheme in Parekh et al. (2005) ...
-                  !! 
-                  !! sinking particulate carbon available for scavenging
-                  !! this assumes a sinking rate of 100 m / d (Moore & Braucher, 2008),
-                  xCscav1     = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3
-                  !! 
-                  !! scale by Honeyman et al. (1981) exponent coefficient
-                  !! multiply by 1.e-3 to express C flux in g C rather than mg C
-                  xCscav2     = (xCscav1 * 1.e-3)**0.58
-                  !!
-                  !! multiply by Galbraith et al. (2010) scavenging rate
-                  xk_org      = 0.5 ! ((g C m/3)^-1) / d
-                  xORGscav    = xk_org * xCscav2 * xFeF
-                  !!
-                  !! Galbraith et al. (2010) also include an inorganic bit ...
-                  !! 
-                  !! this occurs at a fixed rate, again based on the availability of
-                  !! "free" iron
-                  !!
-                  !! k_inorg  = 1000 d**-1 nmol Fe**-0.5 kg**-0.5
-                  !!
-                  !! to implement this here, scale xFeF by 1026 to put in units of
-                  !! umol/m3 which approximately equal nmol/kg
-                  !!
-                  xk_inorg    = 1000. ! ((nmol Fe / kg)^1.5)
-                  xINORGscav  = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e-3
-                  !!
-                  !! sum these two terms together
-                  ffescav(ji,jj)     = xORGscav + xINORGscav
-               else
-                  !!----------------------------------------------------------------------
-                  !! No Scheme: you coward!
-                  !! This scheme puts its head in the sand and eskews any decision about
-                  !! how iron is removed from the ocean; prepare to get deluged in iron
-                  !! you fool!
-                  !!----------------------------------------------------------------------
-                  ffescav(ji,jj)     = 0.
-               endif
-
-               !!----------------------------------------------------------------------
-               !! Other iron cycle processes
-               !!----------------------------------------------------------------------
-               !!
-               !! aeolian iron deposition
-               if (jk.eq.1) then
-                  !! zirondep   is in mmol-Fe / m2 / day
-                  !! ffetop(ji,jj)     is in mmol-dissolved-Fe / m3 / day
-                  ffetop(ji,jj)  = zirondep(ji,jj) * xfe_sol / fse3t(ji,jj,jk) 
-               else
-                  ffetop(ji,jj)  = 0.0
-               endif
-               !!
-               !! seafloor iron addition
-               !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down
-               !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then
-               if ((jk.eq.mbathy(ji,jj)).AND.jk.le.i0500) then
-                  !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a
-                  !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value
-                  !! but apply it everywhere
-                  !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 1-37)
-                  ffebot(ji,jj)  = (xfe_sed / fse3t(ji,jj,jk))
-               else
-                  ffebot(ji,jj)  = 0.0
-               endif
-
-               !! AXY (16/12/09): remove iron addition/removal processes
-               !! For the purposes of the quarter degree run, the iron cycle is being pegged to the
-               !! initial condition supplied by Mick Follows via restoration with a 30 day period;
-               !! iron addition at the seafloor is still permitted with the idea that this extra
-               !! iron will be removed by the restoration away from the source
-               !! ffescav(ji,jj) = 0.0
-               !! ffetop(ji,jj)  = 0.0
-               !! ffebot(ji,jj)  = 0.0
-
-# if defined key_debug_medusa
-               !! report miscellaneous calculations
-               if (idf.eq.1.AND.idfval.eq.1) then
-                  IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'xfe_sol  = ', xfe_sol
-                  IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass
-                  IF (lwp) write (numout,*) 'ffetop(',jk,')  = ', ffetop(ji,jj)
-                  IF (lwp) write (numout,*) 'ffebot(',jk,')  = ', ffebot(ji,jj)
-                  IF (lwp) write (numout,*) 'xFree(',jk,')   = ', xFree(ji,jj)
-                  IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav(ji,jj)
-               endif
-# endif
-
-! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PROCESSES - marc 20/4/17 
-! (iron chemistry and scavenging is 340 lines)
-            ENDIF
-         ENDDO
-         ENDDO
+         !!----------------------------------------------------------------
+         !! Iron chemistry and scavenging
+         !!----------------------------------------------------------------
+         CALL iron_chem_scav( jk )
+
+! Miscellaneous processes - marc
 
          DO jj = 2,jpjm1
