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

source: branches/UKMO/dev_r5518_GO6_fix_cpp_keys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90 @ 10004

Last change on this file since 10004 was 10004, checked in by frrh, 6 years ago

Various fixes

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