New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
iron_chem_scav.F90 in branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90 @ 7978

Last change on this file since 7978 was 7978, checked in by marc, 7 years ago

Moved the iron chemistry and scavenging out of trcbio_medusa.F90

File size: 24.1 KB
Line 
1MODULE iron_chem_scav_mod
2   !!======================================================================
3   !!                         ***  MODULE iron_chem_scav_mod  ***
4   !! Calculate the iron chemistry and scavenging.
5   !!======================================================================
6   !! History :
7   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
8   !!----------------------------------------------------------------------
9#if defined key_medusa
10   !!----------------------------------------------------------------------
11   !!                                                   MEDUSA bio-model
12   !!----------------------------------------------------------------------
13
14   IMPLICIT NONE
15   PRIVATE
16     
17   PUBLIC   iron_chem_scav        ! Called in trcbio_medusa.F90
18
19   !!----------------------------------------------------------------------
20   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
21   !! $Id$
22   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
23   !!----------------------------------------------------------------------
24
25CONTAINS
26
27   SUBROUTINE iron_chem_scav( jk )
28      !!-------------------------------------------------------------------
29      !!                     ***  ROUTINE iron_chem_scav  ***
30      !! This called from TRC_BIO_MEDUSA and
31      !!  -
32      !!-------------------------------------------------------------------
33      USE bio_medusa_mod,    ONLY: ffastc, ffastca, ffastsi,              &
34                                   ffetop, ffebot, ffescav, xfree,        & 
35                                   zdet, zfer, zphd, zphn, zzme, zzmi 
36      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,        &
37                                   mbathy, tmask
38      USE par_kind,          ONLY: wp
39      USE par_oce,           ONLY: jpim1, jpjm1
40      USE sms_medusa,        ONLY: i0500, jiron, xfe_sed, xfe_sol,        &
41                                   xk_FeL, xk_sc_Fe, xLgT,                &
42                                   xmassc, xmassca, xmasssi,              &
43                                   xthetad, xthetapd, xthetapn,           &
44                                   xthetazme, xthetazmi,                  &
45                                   zirondep
46
47   !!* Substitution
48#  include "domzgr_substitute.h90"
49
50      !! Level
51      INTEGER, INTENT( in ) :: jk
52
53      !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
54      !! state variables for iron-ligand system
55      REAL(wp) ::    xLgF, xFeT, xFeF, xFeL
56      !! iron-ligand parameters
57      REAL(wp) ::    xb_coef_tmp, xb2M4ac
58      !! max Fe' parameters
59      REAL(wp) ::    xmaxFeF,fdeltaFe
60      !!
61      !! local parameters for Moore et al. (2004) alternative scavenging
62      !! scheme
63      REAL(wp) ::    fbase_scav,fscal_sink,fscal_part,fscal_scav
64      !!
65      !! local parameters for Moore et al. (2008) alternative scavenging
66      !! scheme
67      REAL(wp) ::    fscal_csink,fscal_sisink,fscal_casink
68      !!
69      !! local parameters for Galbraith et al. (2010) alternative
70      !! scavenging scheme.
71      !! organic portion of scavenging
72      REAL(wp) ::    xCscav1, xCscav2, xk_org, xORGscav
73      !! inorganic portion of scavenging
74      REAL(wp) ::    xk_inorg, xINORGscav
75
76      INTEGER :: ji, jj
77
78! I'D LIKE TO PULL THE IF (jiron) statements outside the DO loops - marc
79
80      !!------------------------------------------------------------------
81      !! Iron chemistry and fractionation
82      !! following the Parekh et al. (2004) scheme adopted by the Met.
83      !! Office, Medusa models total iron but considers "free" and
84      !! ligand-bound forms for the purposes of scavenging (only "free"
85      !! iron can be scavenged
86      !!------------------------------------------------------------------
87      DO jj = 2,jpjm1
88         DO ji = 2,jpim1
89            !! OPEN wet point IF..THEN loop
90            if (tmask(ji,jj,1) == 1) then
91               !!
92               !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
93               xFeT        = zfer(ji,jj) * 1.e3
94               !!
95               !! calculate fractionation (based on Diat-HadOCC; in turn
96               !! based on Parekh et al., 2004)
97               xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0
98               xb2M4ac     = max(((xb_coef_tmp * xb_coef_tmp) +              &
99                                  (4.0 * xk_FeL * xLgT)), 0.0)
100               !!
101               !! "free" ligand concentration
102               xLgF        = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL
103               !!
104               !! ligand-bound iron concentration
105               xFeL        = xLgT - xLgF
106               !!
107               !! "free" iron concentration (and convert to mmol Fe / m3)
108               xFeF        = (xFeT - xFeL) * 1.e-3
109               xFree(ji,jj)= xFeF / (zfer(ji,jj) + tiny(zfer(ji,jj)))
110               !!
111               !! scavenging of iron (multiple schemes); I'm only really
112               !! happy with the first one at the moment - the others
113               !! involve assumptions (sometimes guessed at by me) that
114               !! are potentially questionable
115               !!
116               if (jiron.eq.1) then
117                  !!------------------------------------------------------
118                  !! Scheme 1: Dutkiewicz et al. (2005)
119                  !! This scheme includes a single scavenging term based
120                  !! solely on a fixed rate and the availablility of
121                  !! "free" iron
122                  !!------------------------------------------------------
123                  !! = mmol/m3/d
124                  ffescav(ji,jj)     = xk_sc_Fe * xFeF
125                  !!
126                  !!------------------------------------------------------
127                  !!
128                  !! Mick's code contains a further (optional) implicit
129                  !! "scavenging" of iron that sets an upper bound on
130                  !! "free" iron concentration, and essentially caps the
131                  !! concentration of total iron as xFeL + "free" iron;
132                  !! since the former is constrained by a fixed total
133                  !! ligand concentration (= 1.0 umol/m3), and the latter
134                  !! isn't allowed above this upper bound, total iron is
135                  !! constrained to a maximum of ...
136                  !!
137                  !!    xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3
138                  !!                                  = 1.3 umol / m3
139                  !!
140                  !! In Mick's code, the actual value of total iron is
141                  !! reset to this sum (i.e. TFe = FeL + Fe'; but
142                  !! Fe' <= 0.3 umol/m3); this isn't our favoured approach
143                  !! to tracer updating here (not least because of the
144                  !! leapfrog), so here the amount scavenged is augmented
145                  !! by an additional amount that serves to drag total
146                  !! iron back towards that expected from this limitation
147                  !! on iron concentration ...
148                  !!
149                  !! = umol/m3
150                  xmaxFeF     = min((xFeF * 1.e3), 0.3)
151                  !!
152                  !! Here, the difference between current total Fe and
153                  !! (FeL + Fe') is calculated and added to the scavenging
154                  !! flux already calculated above ...
155                  !!
156                  !! = mmol/m3
157                  fdeltaFe    = (xFeT - (xFeL + xmaxFeF)) * 1.e-3
158                  !!
159                  !! This assumes that the "excess" iron is dissipated
160                  !! with a time-scale of 1 day; seems reasonable to me
161                  !! ... (famous last words)
162                  !!
163                  !! = mmol/m3/d
164                  ffescav(ji,jj)     = ffescav(ji,jj) + fdeltaFe
165                  !!
166# if defined key_deep_fe_fix
167                  !! AXY (17/01/13)
168                  !! stop scavenging for iron concentrations below
169                  !! 0.5 umol / m3 at depths greater than 1000 m; this
170                  !! aims to end MEDUSA's continual loss of iron at depth
171                  !! without impacting things at the surface too much; the
172                  !! justification for this is that it appears to be what
173                  !! Mick Follows et al. do in their work (as evidenced by
174                  !! the iron initial condition they supplied me with); to
175                  !! be honest, it looks like Follow et al. do this at
176                  !! shallower depths than 1000 m, but I'll stick with this
177                  !! for now; I suspect that this seemingly arbitrary
178                  !! approach effectively "parameterises" the
179                  !! particle-based scavenging rates that other models use
180                  !! (i.e. at depth there are no sinking particles, so
181                  !! scavenging stops); it might be fun justifying this in
182                  !! a paper though!
183                  !!
184                  if ((fsdepw(ji,jj,jk).gt.1000.) .and. (xFeT.lt.0.5)) then
185                     ffescav(ji,jj) = 0.
186                  endif
187# endif
188                  !!
189               elseif (jiron.eq.2) then
190                  !!------------------------------------------------------
191                  !! Scheme 2: Moore et al. (2004)
192                  !! This scheme includes a single scavenging term that
193                  !! accounts for both suspended and sinking particles in
194                  !! the water column; this term scavenges total iron rather
195                  !! than "free" iron
196                  !!------------------------------------------------------
197                  !!
198                  !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
199                  xFeT        = zfer(ji,jj) * 1.e3
200                  !!
201                  !! this has a base scavenging rate (12% / y) which is
202                  !! modified by local particle concentration and sinking
203                  !! flux (and dust - but I'm ignoring that here for now)
204                  !! and which is accelerated when Fe concentration gets
205                  !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased
206                  !! as concentrations below 0.4 nM (= 0.4 umol/m3 =
207                  !! 0.0004 mmol/m3)
208                  !!
209                  !! base scavenging rate (0.12 / y)
210                  fbase_scav = 0.12 / 365.25
211                  !!
212                  !! calculate sinking particle part of scaling factor
213                  !! this takes local fast sinking carbon (mmol C / m2 / d)
214                  !! and gets it into nmol C / cm3 / s ("rdt" below is the
215                  !! number of seconds in a model timestep)
216                  !!
217                  !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.)
218                  !!
219                  !! ... actually, re-reading Moore et al.'s equations, it
220                  !! looks like he uses his sinking flux directly, without
221                  !! scaling it by time-step or anything, so I'll copy this
222                  !! here ...
223                  !!
224                  fscal_sink = ffastc(ji,jj) * 1.e2
225                  !!
226                  !! calculate particle part of scaling factor
227                  !! this totals up the carbon in suspended particles
228                  !! (Pn, Pd, Zmi, Zme, D),
229                  !! which comes out in mmol C / m3 (= nmol C / cm3), and
230                  !! then multiplies it by a magic factor, 0.002, to get it
231                  !! into nmol C / cm2 / s
232                  !!
233                  fscal_part = ( (xthetapn * zphn(ji,jj)) +                  &
234                                 (xthetapd * zphd(ji,jj)) +                  &
235                                 (xthetazmi * zzmi(ji,jj)) +                 &
236                                 (xthetazme * zzme(ji,jj)) +                 &
237                                 (xthetad * zdet(ji,jj)) ) * 0.002
238                  !!
239                  !! calculate scaling factor for base scavenging rate
240                  !! this uses the (now correctly scaled) sinking flux and
241                  !! standing
242                  !! particle concentration, divides through by some sort
243                  !! of reference value (= 0.0066 nmol C / cm2 / s) and
244                  !! then uses this, or not if its too high, to rescale the
245                  !! base scavenging rate
246                  !!
247                  fscal_scav = fbase_scav *                                  &
248                               min(((fscal_sink + fscal_part) / 0.0066), 4.0)
249                  !!
250                  !! the resulting scavenging rate is then scaled further
251                  !! according to the local iron concentration (i.e.
252                  !! diminished in low iron regions; enhanced in high iron
253                  !! regions; less alone in intermediate iron regions)
254                  !!
255                  if (xFeT.lt.0.4) then
256                     !!
257                     !! low iron region
258                     !!
259                     fscal_scav = fscal_scav * (xFeT / 0.4)
260                     !!
261                  elseif (xFeT.gt.0.6) then
262                     !!
263                     !! high iron region
264                     !!
265                     fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4))
266                     !!
267                  else
268                     !!
269                     !! intermediate iron region: do nothing
270                     !!
271                  endif
272                  !!
273                  !! apply the calculated scavenging rate ...
274                  !!
275                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
276                  !!
277               elseif (jiron.eq.3) then
278                  !!------------------------------------------------------
279                  !! Scheme 3: Moore et al. (2008)
280                  !! This scheme includes a single scavenging term that
281                  !! accounts for sinking particles in the water column,
282                  !! and includes organic C, biogenic opal, calcium
283                  !! carbonate and dust in this (though the latter is
284                  !! ignored here until I work out what units the incoming
285                  !! "dust" flux is in); this term scavenges total iron
286                  !! rather than "free" iron
287                  !!------------------------------------------------------
288                  !!
289                  !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
290                  xFeT        = zfer(ji,jj) * 1.e3
291                  !!
292                  !! this has a base scavenging rate which is modified by
293                  !! local particle sinking flux (including dust - but I'm
294                  !! ignoring that here for now) and which is accelerated
295                  !! when Fe concentration is > 0.6 nM (= 0.6 umol/m3 =
296                  !! 0.0006 mmol/m3), and decreased as concentrations <
297                  !! 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3)
298                  !!
299                  !! base scavenging rate (Fe_b in paper; units may be
300                  !! wrong there)
301                  fbase_scav = 0.00384 ! (ng)^-1 cm
302                  !!
303                  !! calculate sinking particle part of scaling factor;
304                  !! this converts mmol / m2 / d fluxes of organic carbon,
305                  !! silicon and calcium carbonate into ng / cm2 / s
306                  !! fluxes; it is assumed here that the mass conversions
307                  !! simply consider the mass of the main element
308                  !! (C, Si and Ca) and *not* the mass of the molecules
309                  !! that they are part of; Moore et al. (2008) is unclear
310                  !! on the conversion that should be used
311                  !!
312                  !! milli -> nano; mol -> gram; /m2 -> /cm2; /d -> /s
313                  !! ng C  / cm2 / s
314                  fscal_csink  = (ffastc(ji,jj)  * 1.e6 * xmassc  *           &
315                                  1.e-4 / 86400.)
316                  !! ng Si / cm2 / s
317                  fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi *           &
318                                  1.e-4 / 86400.)
319                  !! ng Ca / cm2 / s
320                  fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca *           &
321                                  1.e-4 / 86400.)
322                  !!
323                  !! sum up these sinking fluxes and convert to ng / cm
324                  !! by dividing through by a sinking rate of
325                  !! 100 m / d = 1.157 cm / s
326                  !! ng / cm
327                  fscal_sink   = ((fscal_csink * 6.) + fscal_sisink +        &
328                                  fscal_casink) / (100. * 1.e3 / 86400)
329                  !!
330                  !! now calculate the scavenging rate based upon the base
331                  !! rate and this particle flux scaling; according to the
332                  !! published units, the result actually has *no* units,
333                  !! but as it must be expressed per unit time for it to
334                  !! make any sense, I'm assuming a missing "per second"
335                  !! / s
336                  fscal_scav = fbase_scav * fscal_sink
337                  !!
338                  !! the resulting scavenging rate is then scaled further
339                  !! according to the local iron concentration (i.e.
340                  !! diminished in low iron regions; enhanced in high iron
341                  !! regions; less alone in intermediate iron regions)
342                  !!
343                  if (xFeT.lt.0.5) then
344                     !!
345                     !! low iron region (0.5 instead of the 0.4 in Moore
346                     !! et al., 2004)
347                     !!
348                     fscal_scav = fscal_scav * (xFeT / 0.5)
349                     !!
350                  elseif (xFeT.gt.0.6) then
351                     !!
352                     !! high iron region (functional form different in
353                     !! Moore et al., 2004)
354                     !!
355                     fscal_scav = fscal_scav + ((xFeT - 0.6) * 0.00904)
356                     !!
357                  else
358                     !!
359                     !! intermediate iron region: do nothing
360                     !!
361                  endif
362                  !!
363                  !! apply the calculated scavenging rate ...
364                  !!
365                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
366                  !!
367               elseif (jiron.eq.4) then
368                  !!------------------------------------------------------
369                  !! Scheme 4: Galbraith et al. (2010)
370                  !! This scheme includes two scavenging terms, one for
371                  !! organic, particle-based scavenging, and another for
372                  !! inorganic scavenging; both terms scavenge "free" iron
373                  !! only
374                  !!------------------------------------------------------
375                  !!
376                  !! Galbraith et al. (2010) present a more straightforward
377                  !! outline of the scheme in Parekh et al. (2005) ...
378                  !!
379                  !! sinking particulate carbon available for scavenging
380                  !! this assumes a sinking rate of 100 m / d (Moore &
381                  !! Braucher, 2008),
382                  xCscav1     = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3
383                  !!
384                  !! scale by Honeyman et al. (1981) exponent coefficient
385                  !! multiply by 1.e-3 to express C flux in g C rather than
386                  !! mg C
387                  xCscav2     = (xCscav1 * 1.e-3)**0.58
388                  !!
389                  !! multiply by Galbraith et al. (2010) scavenging rate
390                  xk_org      = 0.5 ! ((g C m/3)^-1) / d
391                  xORGscav    = xk_org * xCscav2 * xFeF
392                  !!
393                  !! Galbraith et al. (2010) also include an inorganic bit ...
394                  !!
395                  !! this occurs at a fixed rate, again based on the
396                  !! availability of "free" iron
397                  !!
398                  !! k_inorg  = 1000 d**-1 nmol Fe**-0.5 kg**-0.5
399                  !!
400                  !! to implement this here, scale xFeF by 1026 to put in
401                  !! units of umol/m3 which approximately equal nmol/kg
402                  !!
403                  xk_inorg    = 1000. ! ((nmol Fe / kg)^1.5)
404                  xINORGscav  = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e-3
405                  !!
406                  !! sum these two terms together
407                  ffescav(ji,jj)     = xORGscav + xINORGscav
408               else
409                  !!------------------------------------------------------
410                  !! No Scheme: you coward!
411                  !! This scheme puts its head in the sand and eskews any
412                  !! decision about how iron is removed from the ocean;
413                  !! prepare to get deluged in iron you fool!
414                  !!------------------------------------------------------
415                  ffescav(ji,jj)     = 0.
416               endif
417
418               !!---------------------------------------------------------
419               !! Other iron cycle processes
420               !!---------------------------------------------------------
421               !!
422               !! aeolian iron deposition
423               if (jk.eq.1) then
424                  !! zirondep   is in mmol-Fe / m2 / day
425                  !! ffetop(ji,jj)     is in mmol-dissolved-Fe / m3 / day
426                  ffetop(ji,jj)  = zirondep(ji,jj) * xfe_sol / fse3t(ji,jj,jk) 
427               else
428                  ffetop(ji,jj)  = 0.0
429               endif
430               !!
431               !! seafloor iron addition
432               !! AXY (10/07/12): amended to only apply sedimentary flux up
433               !! to ~500 m down
434               !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then
435               if ((jk.eq.mbathy(ji,jj)).AND.jk.le.i0500) then
436                  !! Moore et al. (2004) cite a coastal California value of
437                  !! 5 umol/m2/d, but adopt a global value of 2 umol/m2/d
438                  !! for all areas < 1100 m; here we use this latter value
439                  !! but apply it everywhere
440                  !! AXY (21/07/09): actually, let's just apply it below
441                  !! 1100 m (levels 1-37)
442                  ffebot(ji,jj)  = (xfe_sed / fse3t(ji,jj,jk))
443               else
444                  ffebot(ji,jj)  = 0.0
445               endif
446
447               !! AXY (16/12/09): remove iron addition/removal processes
448               !! For the purposes of the quarter degree run, the iron
449               !! cycle is being pegged to the initial condition supplied
450               !! by Mick Follows via restoration with a 30 day period;
451               !! iron addition at the seafloor is still permitted with
452               !! the idea that this extra iron will be removed by the
453               !! restoration away from the source
454               !! ffescav(ji,jj) = 0.0
455               !! ffetop(ji,jj)  = 0.0
456               !! ffebot(ji,jj)  = 0.0
457
458# if defined key_debug_medusa
459               !! report miscellaneous calculations
460               if (idf.eq.1.AND.idfval.eq.1) then
461                  IF (lwp) write (numout,*) '------------------------------'
462                  IF (lwp) write (numout,*) 'xfe_sol  = ', xfe_sol
463                  IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass
464                  IF (lwp) write (numout,*) 'ffetop(',jk,')  = ', ffetop(ji,jj)
465                  IF (lwp) write (numout,*) 'ffebot(',jk,')  = ', ffebot(ji,jj)
466                  IF (lwp) write (numout,*) 'xFree(',jk,')   = ', xFree(ji,jj)
467                  IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav(ji,jj)
468               endif
469# endif
470            ENDIF
471         ENDDO
472      ENDDO
473
474   END SUBROUTINE iron_chem_scav
475
476#else
477   !!======================================================================
478   !!  Dummy module :                                   No MEDUSA bio-model
479   !!======================================================================
480CONTAINS
481   SUBROUTINE iron_chem_scav( )                    ! Empty routine
482      WRITE(*,*) 'iron_chem_scav: You should not have seen this print! error?'
483   END SUBROUTINE iron_chem_scav
484#endif 
485
486   !!======================================================================
487END MODULE iron_chem_scav_mod
Note: See TracBrowser for help on using the repository browser.