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_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90 @ 8012

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

Pulled further code from trcbio_medusa.F90 into other files

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