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 @ 10149

Last change on this file since 10149 was 10020, checked in by marc, 6 years ago

GMED ticket 406. CPP key fixes.

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