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.
detritus_fast_sink.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90 @ 9309

Last change on this file since 9309 was 8441, checked in by frrh, 7 years ago

Commit changes relating to Met Office GMED ticket 339 for the modularisation of
of trcbio_medusa.F90.

Branch http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
from revisions 8394 to 8423 inclusive refer.

File size: 50.9 KB
RevLine 
[8395]1MODULE detritus_fast_sink_mod
2   !!======================================================================
3   !!                         ***  MODULE detritus_fast_sink_mod  ***
4   !! Calculates fast-sinking detritus
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   detritus_fast_sink        ! Called in detritus.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 detritus_fast_sink( jk, iball )
28      !!-------------------------------------------------------------------
29      !!                     ***  ROUTINE detritus_fast_sink  ***
30      !! This called from DETRITUS and calculates the fast-sinking detritus
31      !!-------------------------------------------------------------------
32      USE bio_medusa_mod,    ONLY: b0,                                     &
33                                   f_benout_c, f_benout_ca, f_benout_fe,   &
34                                   f_benout_lyso_ca, f_benout_n,           &
35                                   f_benout_si,                            &
36                                   f_fbenin_c, f_fbenin_ca, f_fbenin_fe,   &
37                                   f_fbenin_n, f_fbenin_si, f_omcal,       &
38                                   fccd, fdep1, fdd,                       &
39                                   fdpd, fdpd2, fdpds, fdpds2,             &
40                                   fdpn, fdpn2,                            &
41                                   fdzme, fdzme2, fdzmi, fdzmi2,           &
42                                   ffast2slowc, ffast2slown,               &
43                                   ffastc, ffastca, ffastfe, ffastn,       &
44                                   ffastsi,                                &
45                                   fgmed, fgmepd, fgmepds, fgmepn,         &
46                                   fgmezmi,                                &
47                                   fgmid, fgmipn,                          &
48                                   ficme, ficmi,                           &
49                                   fifd_fe, fifd_n, fifd_si,               &
50                                   finme, finmi,                           &
51                                   fmeexcr, fmiexcr,                       &
52                                   fofd_fe, fofd_n, fofd_si,               &
53                                   fregen, fregenfast, fregenfastsi,       &
54                                   fregensi,                               &
55                                   freminc, freminca, freminfe,            &
56                                   freminn, freminsi,                      &
57                                   fsdiss,                                 &
58                                   fsedc, fsedca, fsedn, fsedfe, fsedsi,   &
59                                   fslowc, fslowcflux,                     &
60                                   fslown, fslownflux,                     &
61                                   ftempc, ftempca, ftempfe, ftempn,       &
62                                   ftempsi,                                &
63# if defined key_roam
64                                   fifd_c, fofd_c, fregenfastc,            &
65# endif
66                                   idf, idfval,                            &
67                                   zdet, zdtc
68      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,         &
69                                   gphit, mbathy, tmask
70      USE in_out_manager,    ONLY: lwp, numout
71      USE oce,               ONLY: tsn
72      USE par_kind,          ONLY: wp
73      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
74      USE sms_medusa,        ONLY: f2_ccd_cal, f3_omcal,                   &
75                                   jexport, jfdfate, jinorgben, jocalccd,  &
76                                   jorgben, jp_tem, jrratio,               &
77                                   ocal_ccd, vsed,                         &
78                                   xbetac, xbetan,                         &
79                                   xcaco3a, xcaco3b,                       &
80                                   xfastc, xfastca, xfastsi,               &
81                                   xfdfrac1, xfdfrac2, xfdfrac3,           &
82                                   xmassc, xmassca, xmasssi,               &
83                                   xphi, xprotca, xprotsi,                 &
84                                   xrfn, xridg_r0,                         &
85                                   xsedc, xsedca, xsedfe,xsedn, xsedsi,    &
86                                   xthetapd, xthetapn,                     &
87                                   xthetazme, xthetazmi,                   &
88                                   zn_sed_c, zn_sed_ca, zn_sed_fe,         &
89                                   zn_sed_n, zn_sed_si
90
91   !!* Substitution
92#  include "domzgr_substitute.h90"
93
94      !! Level
95      INTEGER, INTENT( in ) :: jk
96      !! Fast detritus ballast scheme (0 = no; 1 = yes)
97      INTEGER, INTENT( in ) :: iball
98
99      !! Loop variables
100      INTEGER :: ji, jj
101
102      REAL(wp) :: fb_val, fl_sst
103      !! Particle flux
104      REAL(wp) :: fcaco3
105      REAL(wp) :: fprotf
106      REAL(wp), DIMENSION(jpi,jpj) :: fccd_dep
107      !! temporary variables
108      REAL(wp) :: fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8
109
110      !! The two sections below, slow detritus creation and Nutrient
111      !! regeneration are moved from just above the CALL to DETRITUS
112      !! in trcbio_medusa.F90.
113      !!---------------------------------------------------------
114      !! Slow detritus creation
115      !!---------------------------------------------------------
116      DO jj = 2,jpjm1
117         DO ji = 2,jpim1
118            IF (tmask(ji,jj,jk) == 1) THEN
119               !!
120               !! this variable integrates the creation of slow sinking
121               !! detritus to allow the split between fast and slow
122               !! detritus to be diagnosed
123               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) +                 &
124                                ((1.0 - xfdfrac1) * fdpd(ji,jj)) +           &
125                                ((1.0 - xfdfrac2) * fdzme(ji,jj)) +          &
126                                ((1.0 - xbetan) *                            &
127                                 (finmi(ji,jj) + finme(ji,jj)))
128               !!
129               !! this variable records the slow detrital sinking flux at
130               !! this particular depth; it is used in the output of this
131               !! flux at standard depths in the diagnostic outputs;
132               !! needs to be adjusted from per second to per day because
133               !! of parameter vsed
134               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400.
135# if defined key_roam
136               !!
137               !! and the same for detrital carbon
138               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) +                   &
139                                (xthetazmi * fdzmi(ji,jj)) +                 &
140                                (xthetapd * (1.0 - xfdfrac1) *               &
141                                 fdpd(ji,jj)) +                              &
142                                (xthetazme * (1.0 - xfdfrac2) *              &
143                                 fdzme(ji,jj)) +                             &
144                                ((1.0 - xbetac) * (ficmi(ji,jj) +            &
145                                                   ficme(ji,jj)))
146               !!
147               !! this variable records the slow detrital sinking flux
148               !! at this particular depth; it is used in the output of
149               !! this flux at standard depths in the diagnostic
150               !! outputs; needs to be adjusted from per second to per
151               !! day because of parameter vsed
152               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400.
153# endif
154            ENDIF
155         ENDDO
156      ENDDO
157
158      !!---------------------------------------------------------
159      !! Nutrient regeneration
160      !! this variable integrates total nitrogen regeneration down the
161      !! watercolumn; its value is stored and output as a 2D diagnostic;
162      !! the corresponding dissolution flux of silicon (from sources
163      !! other than fast detritus) is also integrated; note that,
164      !! confusingly, the linear loss terms from plankton compartments
165      !! are labelled as fdX2 when one might have expected fdX or fdX1
166      !!---------------------------------------------------------
167      DO jj = 2,jpjm1
168         DO ji = 2,jpim1
169            IF (tmask(ji,jj,jk) == 1) THEN
170               !!
171               !! nitrogen
172               fregen(ji,jj) =                                             &
173                                     ! messy feeding
174                               (((xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + &
175                                 (xphi *                                   &
176                                  (fgmepn(ji,jj) + fgmepd(ji,jj) +         &
177                                   fgmezmi(ji,jj) + fgmed(ji,jj))) +       &
178                                     ! excretion + D remin.
179                                 fmiexcr(ji,jj) + fmeexcr(ji,jj) +         &
180                                 fdd(ji,jj) +                              &
181                                     ! linear mortality
182                                 fdpn2(ji,jj) + fdpd2(ji,jj) +             &
183                                 fdzmi2(ji,jj) + fdzme2(ji,jj)) *          &
184                                fse3t(ji,jj,jk))
185               !!
186               !! silicon
187               fregensi(ji,jj) =                                           &
188                                     ! dissolution + non-lin. mortality
189                                 ((fsdiss(ji,jj) +                         &
190                                   ((1.0 - xfdfrac1) * fdpds(ji,jj)) +     &
191                                     ! egestion by zooplankton
192                                   ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +   &
193                                     ! linear mortality
194                                   fdpds2(ji,jj)) * fse3t(ji,jj,jk))
195            ENDIF
196         ENDDO
197      ENDDO
198
199      !!-------------------------------------------------------------------
200      !! Fast-sinking detritus terms
201      !! "local" variables declared so that conservation can be checked;
202      !! the calculated terms are added to the fast-sinking flux later on
203      !! only after the flux entering this level has experienced some
204      !! remineralisation
205      !! note: these fluxes need to be scaled by the level thickness
206      !!-------------------------------------------------------------------
207      DO jj = 2,jpjm1
208         DO ji = 2,jpim1
209            !! OPEN wet point IF..THEN loop
210            if (tmask(ji,jj,jk) == 1) then
211
212               !! nitrogen:   diatom and mesozooplankton mortality
213               ftempn(ji,jj)  = b0 * ((xfdfrac1 * fdpd(ji,jj))  +            &
214                                      (xfdfrac2 * fdzme(ji,jj)))
215               !!
216               !! silicon:    diatom mortality and grazed diatoms
217               ftempsi(ji,jj) = b0 * ((xfdfrac1 * fdpds(ji,jj)) +            &
218                                      (xfdfrac3 * fgmepds(ji,jj)))
219               !!
220               !! iron:       diatom and mesozooplankton mortality
221               ftempfe(ji,jj) = b0 * (((xfdfrac1 * fdpd(ji,jj)) +            &
222                                       (xfdfrac2 * fdzme(ji,jj))) * xrfn)
223               !!
224               !! carbon:     diatom and mesozooplankton mortality
225               ftempc(ji,jj)  = b0 * ((xfdfrac1 * xthetapd * fdpd(ji,jj)) +  & 
226                                      (xfdfrac2 * xthetazme * fdzme(ji,jj)))
227               !!
228            ENDIF
229         ENDDO
230      ENDDO
231
232# if defined key_roam
233      DO jj = 2,jpjm1
234         DO ji = 2,jpim1
235            if (tmask(ji,jj,jk) == 1) then
236               if (jrratio.eq.0) then
237                  !! CaCO3:      latitudinally-based fraction of total
238                  !!             primary production
239                  !!               0.10 at equator; 0.02 at pole
240                  fcaco3 = xcaco3a + ((xcaco3b - xcaco3a) *                  &
241                                      ((90.0 - abs(gphit(ji,jj))) / 90.0))
242               elseif (jrratio.eq.1) then
243                  !! CaCO3:      Ridgwell et al. (2007) submodel, version 1
244                  !!             this uses SURFACE omega calcite to regulate
245                  !!             rain ratio
246                  if (f_omcal(ji,jj).ge.1.0) then
247                     fq1 = (f_omcal(ji,jj) - 1.0)**0.81
248                  else
249                     fq1 = 0.
250                  endif
251                  fcaco3 = xridg_r0 * fq1
252               elseif (jrratio.eq.2) then
253                  !! CaCO3:      Ridgwell et al. (2007) submodel, version 2
254                  !!             this uses FULL 3D omega calcite to regulate
255                  !!             rain ratio
256                  if (f3_omcal(ji,jj,jk).ge.1.0) then
257                     fq1 = (f3_omcal(ji,jj,jk) - 1.0)**0.81
258                  else
259                     fq1 = 0.
260                  endif
261                  fcaco3 = xridg_r0 * fq1
262               endif
263# else
264               !! CaCO3:      latitudinally-based fraction of total primary
265               !!              production
266               !!               0.10 at equator; 0.02 at pole
267               fcaco3 = xcaco3a + ((xcaco3b - xcaco3a) *                      &
268                                   ((90.0 - abs(gphit(ji,jj))) / 90.0))
269# endif
270               !! AXY (09/03/09): convert CaCO3 production from function of
271               !! primary production into a function of fast-sinking material;
272               !! technically, this is what Dunne et al. (2007) do anyway; they
273               !! convert total primary production estimated from surface
274               !! chlorophyll to an export flux for which they apply conversion
275               !! factors to estimate the various elemental fractions (Si, Ca)
276               ftempca(ji,jj) = ftempc(ji,jj) * fcaco3
277
278# if defined key_debug_medusa
279               !! integrate total fast detritus production
280               if (idf.eq.1) then
281                  fifd_n(ji,jj)  = fifd_n(ji,jj)  + (ftempn(ji,jj)  *         &
282                                                     fse3t(ji,jj,jk))
283                  fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi(ji,jj) *         &
284                                                     fse3t(ji,jj,jk))
285                  fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe(ji,jj) *         &
286                                                     fse3t(ji,jj,jk))
287#  if defined key_roam
288                  fifd_c(ji,jj)  = fifd_c(ji,jj)  + (ftempc(ji,jj)  *         &
289                                                     fse3t(ji,jj,jk))
290#  endif
291               endif
292
293               !! report quantities of fast-sinking detritus for each component
294               if (idf.eq.1.AND.idfval.eq.1) then
295                  IF (lwp) write (numout,*) '------------------------------'
296! These variables are not in this routine - marc 28/4/17
297!                  IF (lwp) write (numout,*) 'fdpd(',jk,')    = ', fdpd(ji,jj)
298!                  IF (lwp) write (numout,*) 'fdzme(',jk,')   = ', fdzme(ji,jj)
299                  IF (lwp) write (numout,*) 'ftempn(',jk,')  = ', ftempn(ji,jj)
300                  IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi(ji,jj)
301                  IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe(ji,jj)
302                  IF (lwp) write (numout,*) 'ftempc(',jk,')  = ', ftempc(ji,jj)
303                  IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca(ji,jj)
304                  IF (lwp) write (numout,*) 'flat(',jk,')    = ',             &
305                                            abs(gphit(ji,jj))
306                  IF (lwp) write (numout,*) 'fcaco3(',jk,')  = ', fcaco3
307               endif
308# endif
309            ENDIF
310         ENDDO
311      ENDDO
312
313      !!----------------------------------------------------------
314      !! This version of MEDUSA offers a choice of three methods for
315      !! handling the remineralisation of fast detritus.  All three
316      !! do so in broadly the same way:
317      !!
318      !!   1.  Fast detritus is stored as a 2D array  [ ffastX  ]
319      !!   2.  Fast detritus is added level-by-level  [ ftempX  ]
320      !!   3.  Fast detritus is not remineralised in the top box
321      !!       [ freminX ]
322      !!   4.  Remaining fast detritus is remineralised in the
323      !!       bottom  [ fsedX   ] box
324      !!
325      !! The three remineralisation methods are:
326      !!   
327      !!   1.  Ballast model (i.e. that published in Yool et al.,
328      !!       2011)
329      !!  (1b. Ballast-sans-ballast model)
330      !!   2.  Martin et al. (1987)
331      !!   3.  Henson et al. (2011)
332      !!
333      !! The first of these couples C, N and Fe remineralisation to
334      !! the remineralisation of particulate Si and CaCO3, but the
335      !! latter two treat remineralisation of C, N, Fe, Si and CaCO3
336      !! completely separately.  At present a switch within the code
337      !! regulates which submodel is used, but this should be moved
338      !! to the namelist file.
339      !!
340      !! The ballast-sans-ballast submodel is an original development
341      !! feature of MEDUSA in which the ballast submodel's general
342      !! framework and parameterisation is used, but in which there
343      !! is no protection of organic material afforded by ballasting
344      !! minerals.  While similar, it is not the same as the Martin
345      !! et al. (1987) submodel.
346      !!
347      !! Since the three submodels behave the same in terms of
348      !! accumulating sinking material and remineralising it all at
349      !! the seafloor, these portions of the code below are common to
350      !! all three.
351      !!----------------------------------------------------------
352      if (jexport.eq.1) then
353         DO jj = 2,jpjm1
354            DO ji = 2,jpim1
355               if (tmask(ji,jj,jk) == 1) then
356                  !!=======================================================
357                  !! BALLAST SUBMODEL
358                  !!=======================================================
359                  !!
360                  !!-------------------------------------------------------
361                  !! Fast-sinking detritus fluxes, pt. 1: REMINERALISATION
362                  !! aside from explicitly modelled, slow-sinking detritus, the
363                  !! model includes an implicit representation of detrital
364                  !! particles that sink too quickly to be modelled with
365                  !! explicit state variables; this sinking flux is instead
366                  !! instantaneously remineralised down the water column using
367                  !! the version of Armstrong et al. (2002)'s ballast model
368                  !! used by Dunne et al. (2007); the version of this model
369                  !! here considers silicon and calcium carbonate ballast
370                  !! minerals; this section of the code redistributes the fast
371                  !! sinking material generated locally down the water column;
372                  !! this differs from Dunne et al. (2007) in that fast sinking
373                  !! material is distributed at *every* level below that it is
374                  !! generated, rather than at every level below some fixed
375                  !! depth; this scheme is also different in that sinking
376                  !! material generated in one level is aggregated with that
377                  !! generated by shallower levels; this should make the
378                  !! ballast model more self-consistent (famous last words)
379                  !!-------------------------------------------------------
380                  !!
381                  if (jk.eq.1) then
382                     !! this is the SURFACE OCEAN BOX (no remineralisation)
383                     !!
384                     freminc(ji,jj)  = 0.0
385                     freminn(ji,jj)  = 0.0
386                     freminfe(ji,jj) = 0.0
387                     freminsi(ji,jj) = 0.0
388                     freminca(ji,jj) = 0.0
389                  elseif (jk.le.mbathy(ji,jj)) then
390                     !! this is an OCEAN BOX (remineralise some material)
391                     !!
392                     !! set up CCD depth to be used depending on user choice
393                     if (jocalccd.eq.0) then
394                        !! use default CCD field
395                        fccd_dep(ji,jj) = ocal_ccd(ji,jj)
396                     elseif (jocalccd.eq.1) then
397                        !! use calculated CCD field
398                        fccd_dep(ji,jj) = f2_ccd_cal(ji,jj)
399                     endif
400                     !!
401                     !! === organic carbon ===
402                     !! how much organic C enters this box        (mol)
403                     fq0      = ffastc(ji,jj)
404                     if (iball.eq.1) then
405                        !! how much it weighs
406                        fq1      = (fq0 * xmassc)
407                        !! how much CaCO3 enters this box
408                        fq2      = (ffastca(ji,jj) * xmassca)
409                        !! how much  opal enters this box
410                        fq3      = (ffastsi(ji,jj) * xmasssi)
411                        !! total protected organic C
412                        fq4      = (fq2 * xprotca) + (fq3 * xprotsi)
413                        !! This next term is calculated for C but used for
414                        !! N and Fe as well
415                        !! It needs to be protected in case ALL C is protected
416                        if (fq4.lt.fq1) then
417                           !! protected fraction of total organic C (non-dim)
418                           fprotf   = (fq4 / (fq1 + tiny(fq1)))
419                        else
420                           !! all organic C is protected (non-dim)
421                           fprotf   = 1.0
422                        endif
423                        !! unprotected fraction of total organic C (non-dim)
424                        fq5      = (1.0 - fprotf)
425                        !! how much organic C is unprotected (mol)
426                        fq6      = (fq0 * fq5)
427                        !! how much unprotected C leaves this box (mol)
428                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))
429                        !! how much total C leaves this box (mol)
430                        fq8      = (fq7 + (fq0 * fprotf))
431                        !! C remineralisation in this box (mol)
432                        freminc(ji,jj)  = (fq0 - fq8) / fse3t(ji,jj,jk)
433                        ffastc(ji,jj) = fq8
434# if defined key_debug_medusa
435                        !! report in/out/remin fluxes of carbon for this level
436                           if (idf.eq.1.AND.idfval.eq.1) then
437                              IF (lwp) write (numout,*)                       &
438                                       '------------------------------'
439                              IF (lwp) write (numout,*) 'totalC(',jk,')  = ', &
440                                       fq1
441                              IF (lwp) write (numout,*) 'prtctC(',jk,')  = ', &
442                                       fq4
443                              IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', &
444                                       fprotf
445                              IF (lwp) write (numout,*)                       &
446                                       '------------------------------'
447                              IF (lwp) write (numout,*) 'IN   C(',jk,')  = ', &
448                                       fq0
449                              IF (lwp) write (numout,*) 'LOST C(',jk,')  = ', &
450                                       freminc(ji,jj) * fse3t(ji,jj,jk)
451                              IF (lwp) write (numout,*) 'OUT  C(',jk,')  = ', &
452                                       fq8
453                              IF (lwp) write (numout,*) 'NEW  C(',jk,')  = ', &
454                                       ftempc(ji,jj) * fse3t(ji,jj,jk)
455                           endif
456# endif
457                        else
458                        !! how much organic C leaves this box (mol)
459                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))
460                        !! C remineralisation in this box (mol)
461                        freminc(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)
462                        ffastc(ji,jj)  = fq1
463                     endif
464                     !!
465                     !! === organic nitrogen ===
466                     !! how much organic N enters this box (mol)
467                     fq0      = ffastn(ji,jj)
468                     if (iball.eq.1) then
469                        !! unprotected fraction of total organic N (non-dim)
470                        fq5      = (1.0 - fprotf)
471                        !! how much organic N is unprotected (mol)
472                        fq6      = (fq0 * fq5)
473                        !! how much unprotected N leaves this box (mol)
474                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))
475                        !! how much total N leaves this box (mol)
476                        fq8      = (fq7 + (fq0 * fprotf))
477                        !! N remineralisation in this box (mol)
478                        freminn(ji,jj)  = (fq0 - fq8) / fse3t(ji,jj,jk)
479                        ffastn(ji,jj)  = fq8
480# if defined key_debug_medusa
481                        !! report in/out/remin fluxes of carbon for this level
482                        if (idf.eq.1.AND.idfval.eq.1) then
483                           IF (lwp) write (numout,*)                          &
484                                    '------------------------------'
485                           IF (lwp) write (numout,*) 'totalN(',jk,')  = ', fq1
486                           IF (lwp) write (numout,*) 'prtctN(',jk,')  = ', fq4
487                           IF (lwp) write (numout,*) 'fprotf(',jk,')  = ',    &
488                                    fprotf
489                           IF (lwp) write (numout,*)                          &
490                                    '------------------------------'
491                           if (freminn(ji,jj) < 0.0) then
492                              IF (lwp) write (numout,*) '** FREMIN ERROR **'
493                           endif
494                           IF (lwp) write (numout,*) 'IN   N(',jk,')  = ', fq0
495                           IF (lwp) write (numout,*) 'LOST N(',jk,')  = ',    &
496                                    freminn(ji,jj) * fse3t(ji,jj,jk)
497                           IF (lwp) write (numout,*) 'OUT  N(',jk,')  = ', fq8
498                           IF (lwp) write (numout,*) 'NEW  N(',jk,')  = ',    &
499                                    ftempn(ji,jj) * fse3t(ji,jj,jk)
500                        endif
501# endif
502                     else
503                        !! how much organic N leaves this box (mol)
504                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))
505                        !! N remineralisation in this box (mol)
506                        freminn(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)
507                        ffastn(ji,jj)  = fq1
508                     endif
509                     !!
510                     !! === organic iron ===
511                     !! how much organic Fe enters this box (mol)
512                     fq0      = ffastfe(ji,jj)
513                     if (iball.eq.1) then
514                        !! unprotected fraction of total organic Fe (non-dim)
515                        fq5      = (1.0 - fprotf)
516                        !! how much organic Fe is unprotected (mol)
517                        fq6      = (fq0 * fq5)
518                        !! how much unprotected Fe leaves this box (mol)
519                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))
520                        !! how much total Fe leaves this box (mol)
521                        fq8      = (fq7 + (fq0 * fprotf))
522                        !! Fe remineralisation in this box (mol)
523                        freminfe(ji,jj) = (fq0 - fq8) / fse3t(ji,jj,jk)
524                        ffastfe(ji,jj) = fq8
525                     else
526                        !! how much total Fe leaves this box (mol)
527                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))
528                        !! Fe remineralisation in this box (mol)
529                        freminfe(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
530                        ffastfe(ji,jj) = fq1
531                     endif
532                     !!
533                     !! === biogenic silicon ===
534                     !! how much  opal centers this box (mol)
535                     fq0      = ffastsi(ji,jj)
536                     !! how much  opal leaves this box (mol)
537                     fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastsi))
538                     !! Si remineralisation in this box (mol)
539                     freminsi(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
540                     ffastsi(ji,jj) = fq1
541                     !!
542                     !! === biogenic calcium carbonate ===
543                     !! how much CaCO3 enters this box (mol)
544                     fq0      = ffastca(ji,jj)
545                     if (fsdepw(ji,jj,jk).le.fccd_dep(ji,jj)) then
546                        !! whole grid cell above CCD
547                        !! above lysocline, no Ca dissolves (mol)
548                        fq1      = fq0
549                        !! above lysocline, no Ca dissolves (mol)
550                        freminca(ji,jj) = 0.0
551                        !! which is the last level above the CCD?    (#)
552                        fccd(ji,jj) = real(jk)
553                     elseif (fsdepw(ji,jj,jk).ge.fccd_dep(ji,jj)) then
554                        !! whole grid cell below CCD
555                        !! how much CaCO3 leaves this box (mol)
556                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastca))
557                        !! Ca remineralisation in this box (mol)
558                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
559                     else
560                        !! partial grid cell below CCD
561                        !! amount of grid cell below CCD (m)
562                        fq2      = fdep1(ji,jj) - fccd_dep(ji,jj)
563                        !! how much CaCO3 leaves this box (mol)
564                        fq1      = fq0 * exp(-(fq2 / xfastca))
565                        !! Ca remineralisation in this box (mol)
566                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
567                     endif
568                     ffastca(ji,jj) = fq1 
569                  else
570                     !! this is BELOW THE LAST OCEAN BOX (do nothing)
571                     freminc(ji,jj)  = 0.0
572                     freminn(ji,jj)  = 0.0
573                     freminfe(ji,jj) = 0.0
574                     freminsi(ji,jj) = 0.0
575                     freminca(ji,jj) = 0.0             
576                  endif
577               ENDIF
578            ENDDO
579         ENDDO
580      elseif (jexport.eq.2.or.jexport.eq.3) then
581         DO jj = 2,jpjm1
582            DO ji = 2,jpim1
583               if (tmask(ji,jj,jk) == 1) then
584                  if (jexport.eq.2) then
585                     !!====================================================
586                     !! MARTIN ET AL. (1987) SUBMODEL
587                     !!====================================================
588                     !!
589                     !!----------------------------------------------------
590                     !! This submodel uses the classic Martin et al. (1987)
591                     !! curve to determine the attenuation of fast-sinking
592                     !! detritus down the water column.  All three organic
593                     !! elements, C, N and Fe, are handled identically, and
594                     !! their quantities in sinking particles attenuate
595                     !! according to a power relationship governed by
596                     !! parameter "b".  This is assigned a canonical value
597                     !! of -0.858.  Biogenic opal and calcium carbonate are
598                     !! attentuated using the same function as in the
599                     !! ballast submodel
600                     !!----------------------------------------------------
601                     !!
602                     fb_val = -0.858
603                  elseif (jexport.eq.3) then
604                     !!====================================================
605                     !! HENSON ET AL. (2011) SUBMODEL
606                     !!====================================================
607                     !!
608                     !!----------------------------------------------------
609                     !! This submodel reconfigures the Martin et al. (1987)
610                     !! curve by allowing the "b" value to vary
611                     !! geographically.  Its value is set, following Henson
612                     !! et al. (2011), as a function of local sea surface
613                     !! temperature:
614                     !!   b = -1.06 + (0.024 * SST)
615                     !! This means that remineralisation length scales are
616                     !! longer in warm, tropical areas and shorter in cold,
617                     !! polar areas.  This does seem back-to-front (i.e.
618                     !! one would expect GREATER remineralisation in warmer
619                     !! waters), but is an outcome of analysis of sediment
620                     !! trap data, and it may reflect details of ecosystem
621                     !! structure that pertain to particle production
622                     !! rather than simply Q10.
623                     !!----------------------------------------------------
624                     !!
625                     fl_sst = tsn(ji,jj,1,jp_tem)
626                     fb_val = -1.06 + (0.024 * fl_sst)
627                  endif
628                  !!   
629                  if (jk.eq.1) then
630                     !! this is the SURFACE OCEAN BOX (no remineralisation)
631                     !!
632                     freminc(ji,jj)  = 0.0
633                     freminn(ji,jj)  = 0.0
634                     freminfe(ji,jj) = 0.0
635                     freminsi(ji,jj) = 0.0
636                     freminca(ji,jj) = 0.0
637                  elseif (jk.le.mbathy(ji,jj)) then
638                     !! this is an OCEAN BOX (remineralise some material)
639                     !!
640                     !! === organic carbon ===
641                     !! how much organic C enters this box (mol)
642                     fq0      = ffastc(ji,jj)
643                     !! how much organic C leaves this box (mol)
644                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)
645                     !! C remineralisation in this box (mol)
646                     freminc(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)
647                     ffastc(ji,jj)  = fq1
648                     !!
649                     !! === organic nitrogen ===
650                     !! how much organic N enters this box (mol)
651                     fq0      = ffastn(ji,jj)
652                     !! how much organic N leaves this box (mol)
653                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)
654                     !! N remineralisation in this box (mol)
655                     freminn(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)
656                     ffastn(ji,jj)  = fq1
657                     !!
658                     !! === organic iron ===
659                     !! how much organic Fe enters this box (mol)
660                     fq0      = ffastfe(ji,jj)
661                     !! how much organic Fe leaves this box (mol)
662                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)
663                     !! Fe remineralisation in this box (mol)
664                     freminfe(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
665                     ffastfe(ji,jj) = fq1
666                     !!
667                     !! === biogenic silicon ===
668                     !! how much  opal centers this box (mol)
669                     fq0      = ffastsi(ji,jj)
670                     !! how much  opal leaves this box (mol)
671                     fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastsi))
672                     !! Si remineralisation in this box (mol)
673                     freminsi(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
674                     ffastsi(ji,jj) = fq1
675                     !!
676                     !! === biogenic calcium carbonate ===
677                     !! how much CaCO3 enters this box (mol)
678                     fq0      = ffastca(ji,jj)
679                     if (fsdepw(ji,jj,jk).le.ocal_ccd(ji,jj)) then
680                        !! whole grid cell above CCD
681                        !! above lysocline, no Ca dissolves (mol)
682                        fq1      = fq0
683                        !! above lysocline, no Ca dissolves (mol)
684                        freminca(ji,jj) = 0.0
685                        !! which is the last level above the CCD?    (#)
686                        fccd(ji,jj) = real(jk)
687                     elseif (fsdepw(ji,jj,jk).ge.ocal_ccd(ji,jj)) then
688                        !! whole grid cell below CCD
689                        !! how much CaCO3 leaves this box (mol)
690                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastca))
691                        !! Ca remineralisation in this box (mol)
692                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
693                     else
694                        !! partial grid cell below CCD
695                        !! amount of grid cell below CCD (m)
696                        fq2      = fdep1(ji,jj) - ocal_ccd(ji,jj)
697                        !! how much CaCO3 leaves this box (mol)
698                        fq1      = fq0 * exp(-(fq2 / xfastca))
699                        !! Ca remineralisation in this box (mol)
700                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)
701                     endif
702                     ffastca(ji,jj) = fq1 
703                  else
704                     !! this is BELOW THE LAST OCEAN BOX (do nothing)
705                     freminc(ji,jj)  = 0.0
706                     freminn(ji,jj)  = 0.0
707                     freminfe(ji,jj) = 0.0
708                     freminsi(ji,jj) = 0.0
709                     freminca(ji,jj) = 0.0             
710                  endif
711               ENDIF
712            ENDDO
713         ENDDO
714      endif
715
716      DO jj = 2,jpjm1
717         DO ji = 2,jpim1
718            if (tmask(ji,jj,jk) == 1) then
719               !!----------------------------------------------------------
720               !! Fast-sinking detritus fluxes, pt. 2: UPDATE FAST FLUXES
721               !! here locally calculated additions to the fast-sinking
722               !! flux are added to the total fast-sinking flux; this is
723               !! done here such that material produced in a particular
724               !! layer is only remineralised below this layer
725               !!----------------------------------------------------------
726               !!
727               !! add sinking material generated in this layer to running
728               !! totals
729               !!
730               !! === organic carbon ===
731               !! (diatom and mesozooplankton mortality)
732               ffastc(ji,jj)  = ffastc(ji,jj)  + (ftempc(ji,jj)  *           &
733                                                  fse3t(ji,jj,jk))
734               !!
735               !! === organic nitrogen ===
736               !! (diatom and mesozooplankton mortality)
737               ffastn(ji,jj)  = ffastn(ji,jj)  + (ftempn(ji,jj)  *           &
738                                                  fse3t(ji,jj,jk))
739               !!
740               !! === organic iron ===
741               !! (diatom and mesozooplankton mortality)
742               ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe(ji,jj) *          &
743                                                  fse3t(ji,jj,jk))
744               !!
745               !! === biogenic silicon ===
746               !! (diatom mortality and grazed diatoms)
747               ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi(ji,jj) *          &
748                                                  fse3t(ji,jj,jk))
749               !!
750               !! === biogenic calcium carbonate ===
751               !! (latitudinally-based fraction of total primary production)
752               ffastca(ji,jj) = ffastca(ji,jj) + (ftempca(ji,jj) *          &
753                                                  fse3t(ji,jj,jk))
754            ENDIF
755         ENDDO
756      ENDDO
757
758      DO jj = 2,jpjm1
759         DO ji = 2,jpim1
760            if (tmask(ji,jj,jk) == 1) then
761               !!----------------------------------------------------------
762               !! Fast-sinking detritus fluxes, pt. 3: SEAFLOOR
763               !! remineralise all remaining fast-sinking detritus to dissolved
764               !! nutrients; the sedimentation fluxes calculated here allow the
765               !! separation of what's remineralised sinking through the final
766               !! ocean box from that which is added to the final box by the
767               !! remineralisation of material that reaches the seafloor (i.e.
768               !! the model assumes that *all* material that hits the seafloor
769               !! is remineralised and that none is permanently buried; hey,
770               !! this is a giant GCM model that can't be run for long enough
771               !! to deal with burial fluxes!)
772               !!
773               !! in a change to this process, in part so that MEDUSA behaves
774               !! a little more like ERSEM et al., fast-sinking detritus (N, Fe
775               !! and C) is converted to slow sinking detritus at the seafloor
776               !! instead of being remineralised; the rationale is that in
777               !! shallower shelf regions (... that are not fully mixed!) this
778               !! allows the detrital material to return slowly to dissolved
779               !! nutrient rather than instantaneously as now; the alternative
780               !! would be to explicitly handle seafloor organic material - a
781               !! headache I don't wish to experience at this point; note that
782               !! fast-sinking Si and Ca detritus is just remineralised as
783               !! per usual
784               !!
785               !! AXY (13/01/12)
786               !! in a further change to this process, again so that MEDUSA is
787               !! a little more like ERSEM et al., material that reaches the
788               !! seafloor can now be added to sediment pools and stored for
789               !! slow release; there are new 2D arrays for organic nitrogen,
790               !! iron and carbon and inorganic silicon and carbon that allow
791               !! fast and slow detritus that reaches the seafloor to be held
792               !! and released back to the water column more slowly; these
793               !! arrays are transferred via the tracer restart files between
794               !! repeat submissions of the model
795               !!----------------------------------------------------------
796               !!
797               ffast2slowc(ji,jj)  = 0.0
798               ffast2slown(ji,jj)  = 0.0
799! I don't think this is used - marc 10/4/17
800!               ffast2slowfe(ji,jj) = 0.0
801               !!
802               if (jk.eq.mbathy(ji,jj)) then
803                  !! this is the BOTTOM OCEAN BOX (remineralise everything)
804                  !!
805                  !! AXY (17/01/12): tweaked to include benthos pools
806                  !!
807                  !! === organic carbon ===
808                  if (jfdfate.eq.0 .and. jorgben.eq.0) then
809                     !! C remineralisation in this box (mol/m3)
810                     freminc(ji,jj)  = freminc(ji,jj) + (ffastc(ji,jj) /     &
811                                                         fse3t(ji,jj,jk))
812                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
813                     !! fast C -> slow C (mol/m3)
814                     ffast2slowc(ji,jj) = ffastc(ji,jj) / fse3t(ji,jj,jk)
815                     fslowc(ji,jj)      = fslowc(ji,jj) + ffast2slowc(ji,jj)
816                  elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
817                     !! fast C -> benthic C (mol/m2)
818                     f_fbenin_c(ji,jj)  = ffastc(ji,jj)
819                  endif
820                  !! record seafloor C (mol/m2)
821                  fsedc(ji,jj)   = ffastc(ji,jj)
822                  ffastc(ji,jj)  = 0.0
823                  !!
824                  !! === organic nitrogen ===
825                  if (jfdfate.eq.0 .and. jorgben.eq.0) then
826                     !! N remineralisation in this box (mol/m3)
827                     freminn(ji,jj)  = freminn(ji,jj) + (ffastn(ji,jj) /     &
828                                                         fse3t(ji,jj,jk))
829                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
830                     !! fast N -> slow N (mol/m3)
831                     ffast2slown(ji,jj) = ffastn(ji,jj) / fse3t(ji,jj,jk)
832                     fslown(ji,jj)      = fslown(ji,jj) + ffast2slown(ji,jj)
833                  elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
834                     !! fast N -> benthic N (mol/m2)
835                     f_fbenin_n(ji,jj)  = ffastn(ji,jj)
836                  endif
837                  !! record seafloor N (mol/m2)
838                  fsedn(ji,jj)   = ffastn(ji,jj)
839                  ffastn(ji,jj)  = 0.0
840                  !!
841                  !! === organic iron ===
842                  if (jfdfate.eq.0 .and. jorgben.eq.0) then
843                     !! Fe remineralisation in this box (mol/m3)
844                     freminfe(ji,jj) = freminfe(ji,jj) + (ffastfe(ji,jj) /   &
845                                                          fse3t(ji,jj,jk))
846! I don't think ffast2slowfe is used - marc 10/4/17
847!                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
848!                     !! fast Fe -> slow Fe (mol/m3)
849!                     ffast2slowfe(ji,jj) = ffastn(ji,jj) / fse3t(ji,jj,jk)
850                  elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
851                     !! fast Fe -> benthic Fe (mol/m2)
852                     f_fbenin_fe(ji,jj) = ffastfe(ji,jj)
853                  endif
854                  !! record seafloor Fe (mol/m2)
855                  fsedfe(ji,jj)  = ffastfe(ji,jj)
856                  ffastfe(ji,jj) = 0.0
857                  !!
858                  !! === biogenic silicon ===
859                  if (jinorgben.eq.0) then
860                     !! Si remineralisation in this box (mol/m3)
861                     freminsi(ji,jj) = freminsi(ji,jj) + (ffastsi(ji,jj) /   &
862                                                          fse3t(ji,jj,jk))
863                  elseif (jinorgben.eq.1) then
864                     !! fast Si -> benthic Si
865                     f_fbenin_si(ji,jj) = ffastsi(ji,jj)
866                  endif
867                  !! record seafloor Si (mol/m2)
868                  fsedsi(ji,jj)   = ffastsi(ji,jj)
869                  ffastsi(ji,jj) = 0.0
870                  !!
871                  !! === biogenic calcium carbonate ===
872                  if (jinorgben.eq.0) then
873                     !! Ca remineralisation in this box (mol/m3)
874                     freminca(ji,jj) = freminca(ji,jj) + (ffastca(ji,jj) /   &
875                                                          fse3t(ji,jj,jk))
876                  elseif (jinorgben.eq.1) then
877                     !! fast Ca -> benthic Ca (mol/m2)
878                     f_fbenin_ca(ji,jj) = ffastca(ji,jj)
879                  endif
880                  !! record seafloor Ca (mol/m2)
881                  fsedca(ji,jj)   = ffastca(ji,jj)
882                  ffastca(ji,jj) = 0.0
883               endif
884
885# if defined key_debug_medusa
886               if (idf.eq.1) then
887                  !!-------------------------------------------------------
888                  !! Integrate total fast detritus remineralisation
889                  !!-------------------------------------------------------
890                  !!
891                  fofd_n(ji,jj)  = fofd_n(ji,jj)  + (freminn(ji,jj)  *       &
892                                                     fse3t(ji,jj,jk))
893                  fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi(ji,jj) *       &
894                                                     fse3t(ji,jj,jk))
895                  fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe(ji,jj) *       &
896                                                     fse3t(ji,jj,jk))
897#  if defined key_roam
898                  fofd_c(ji,jj)  = fofd_c(ji,jj)  + (freminc(ji,jj)  *       &
899                                                     fse3t(ji,jj,jk))
900#  endif
901               endif
902# endif
903            ENDIF
904         ENDDO
905      ENDDO
906
907      DO jj = 2,jpjm1
908         DO ji = 2,jpim1
909            if (tmask(ji,jj,jk) == 1) then
910               !!----------------------------------------------------------
911               !! Sort out remineralisation tally of fast-sinking detritus
912               !!----------------------------------------------------------
913               !!
914               !! update fast-sinking regeneration arrays
915               fregenfast(ji,jj)   = fregenfast(ji,jj)   +                  &
916                                     (freminn(ji,jj)  * fse3t(ji,jj,jk))
917               fregenfastsi(ji,jj) = fregenfastsi(ji,jj) +                  &
918                                     (freminsi(ji,jj) * fse3t(ji,jj,jk))
919# if defined key_roam
920               fregenfastc(ji,jj)  = fregenfastc(ji,jj)  +                  &
921                                     (freminc(ji,jj)  * fse3t(ji,jj,jk))
922# endif
923            ENDIF
924         ENDDO
925      ENDDO
926
927      DO jj = 2,jpjm1
928         DO ji = 2,jpim1
929            if (tmask(ji,jj,jk) == 1) then
930               !!----------------------------------------------------------
931               !! Benthic remineralisation fluxes
932               !!----------------------------------------------------------
933               !!
934               if (jk.eq.mbathy(ji,jj)) then
935                  !!
936                  !! organic components
937                  if (jorgben.eq.1) then
938                     f_benout_n(ji,jj)  = xsedn  * zn_sed_n(ji,jj)
939                     f_benout_fe(ji,jj) = xsedfe * zn_sed_fe(ji,jj)
940                     f_benout_c(ji,jj)  = xsedc  * zn_sed_c(ji,jj)
941                  endif
942                  !!
943                  !! inorganic components
944                  if (jinorgben.eq.1) then
945                     f_benout_si(ji,jj) = xsedsi * zn_sed_si(ji,jj)
946                     f_benout_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
947                     !!
948                     !! account for CaCO3 that dissolves when it shouldn't
949                     if ( fsdepw(ji,jj,jk) .le. fccd_dep(ji,jj) ) then
950                        f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
951                     endif
952                  endif
953               endif
954               CALL flush(numout)
955
956            ENDIF
957         ENDDO
958      ENDDO
959
960   END SUBROUTINE detritus_fast_sink
961
962#else
963   !!======================================================================
964   !!  Dummy module :                                   No MEDUSA bio-model
965   !!======================================================================
966CONTAINS
967   SUBROUTINE detritus_fast_sink( )                    ! Empty routine
968      WRITE(*,*) 'detritus_fast_sink: You should not have seen this print! error?'
969   END SUBROUTINE detritus_fast_sink
970#endif 
971
972   !!======================================================================
973END MODULE detritus_fast_sink_mod
Note: See TracBrowser for help on using the repository browser.