source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90 @ 8521

Last change on this file since 8521 was 8521, checked in by marc, 3 years ago

gmed:#346 Merging branches/NERC/dev_r5518_GO6_MEDUSA_conserv

File size: 42.7 KB
Line 
1MODULE bio_medusa_update_mod
2   !!======================================================================
3   !!                         ***  MODULE bio_medusa_update_mod  ***
4   !! Update tracer fields
5   !!======================================================================
6   !! History :
7   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
8   !!   -   ! 2017-08 (A. Yool)            Amend slow-detritus bug
9   !!   -   ! 2017-08 (A. Yool)            Reformatting for clarity
10   !!----------------------------------------------------------------------
11#if defined key_medusa
12   !!----------------------------------------------------------------------
13   !!                                                   MEDUSA bio-model
14   !!----------------------------------------------------------------------
15
16   IMPLICIT NONE
17   PRIVATE
18     
19   PUBLIC   bio_medusa_update        ! Called in trcbio_medusa.F90
20
21   !!----------------------------------------------------------------------
22   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
23   !! $Id$
24   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
25   !!----------------------------------------------------------------------
26
27CONTAINS
28
29   SUBROUTINE bio_medusa_update( kt, jk )
30      !!---------------------------------------------------------------------
31      !!                     ***  ROUTINE bio_medusa_update  ***
32      !! This called from TRC_BIO_MEDUSA and updates the tracer fields
33      !!---------------------------------------------------------------------
34      USE bio_medusa_mod,    ONLY: b0, bddtalk3, bddtdic3, bddtdife3,        &
35                                   bddtdin3, bddtdisi3,                      &
36                                   ibenthic, ibio_switch, idf, idfval,       &
37                                   f_benout_c, f_benout_ca, f_benout_n,      &
38                                   f_benout_si,                              &
39                                   f_co2flux, f_o2flux,                      &
40                                   f_riv_loc_alk, f_riv_loc_c,               &
41                                   f_riv_loc_n, f_riv_loc_si,                & 
42                                   f_riv_alk, f_riv_c, f_riv_n, f_riv_si,    &
43                                   fbddtalk, fbddtdic, fbddtdife,            &
44                                   fbddtdin, fbddtdisi,                      & 
45                                   fdd, fdpd, fdpd2, fdpds, fdpds2,          &
46                                   fdpn, fdpn2,                              &
47                                   fdzme, fdzme2, fdzmi, fdzmi2,             &
48                                   ffast2slowc, ffast2slown,                 &
49                                   ffebot, ffetop, ffescav,                  &
50                                   fflx_fe, fflx_n, fflx_si,                 &
51                                   fgmed, fgmepd, fgmedc, fgmepd, fgmepds,   &
52                                   fgmepn, fgmezmi,                          &
53                                   fgmid, fgmidc, fgmipn,                    &
54                                   ficme, ficmi, finme, finmi,               &
55                                   fmeexcr, fmegrow, fmeresp,                &
56                                   fmiexcr, fmigrow, fmiresp,                &
57                                   fnit_cons, fnit_prod,                     &
58                                   fprd, fprds, fprn,                        &
59                                   frd,                                      &
60                                   freminc, freminca, freminn, freminsi,     &
61                                   frn,                                      &
62                                   fsil_cons, fsil_prod, fsdiss,             &
63                                   ftempca, fthetad, fthetan,                &
64                                   fslowsink, fslowgain, fslowloss,          & ! AXY (22/08/17)
65                                   f_sbenin_n, f_sbenin_c,                   &
66# if defined key_roam
67                                   fslowsinkc, fslowgainc, fslowlossc,       & ! AXY (22/08/17)
68                                   fcar_cons, fcar_prod, fcomm_resp,         &
69                                   fddc, fflx_a, fflx_c, fflx_o2, zoxy,      &
70                                   foxy_anox, foxy_cons, foxy_prod,          &
71# endif
72                                   zpds, zphd, zphn
73      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, mbathy, tmask
74      USE in_out_manager,    ONLY: lwp, numout
75      USE lib_mpp,           ONLY: ctl_stop
76      USE par_kind,          ONLY: wp
77      USE par_medusa,        ONLY: jp_medusa,                                &
78                                   jpalk, jpchd, jpchn, jpdet, jpdic,        &
79                                   jpdin, jpdtc, jpfer, jpoxy, jppds,        &
80                                   jpphd, jpphn, jpsil, jpzme, jpzmi,        &
81                                   jpalk_lc, jpchd_lc, jpchn_lc, jpdet_lc,   & 
82                                   jpdic_lc, jpdin_lc, jpdtc_lc, jpfer_lc,   & 
83                                   jpoxy_lc, jppds_lc, jpphd_lc, jpphn_lc,   &
84                                   jpsil_lc, jpzme_lc, jpzmi_lc
85      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
86      USE par_trc,           ONLY: jptra
87      USE sms_medusa,        ONLY: friver_dep,                               &
88                                   jinorgben, jorgben,                       &
89                                   jriver_alk, jriver_c,                     &
90                                   jriver_n, jriver_si,                      &
91                                   xbetac, xbetan,                           &
92                                   xfdfrac1, xfdfrac2, xfdfrac3,             &
93                                   xo2min, xphi, xrfn,                       &
94                                   xthetanit, xthetapd, xthetapn,            &
95                                   xthetarem, xthetazme, xthetazmi,          &
96                                   xxi
97      USE trc,               ONLY: med_diag, tra
98
99   !!* Substitution
100#  include "domzgr_substitute.h90"
101
102      !! time (integer timestep)
103      INTEGER, INTENT( in ) :: kt
104      !! Level
105      INTEGER, INTENT( in ) :: jk
106
107      !! Loop variables
108      INTEGER :: ji, jj, jn
109
110      !! AXY (23/08/13): changed from individual variables for each flux to
111      !!                 an array that holds all fluxes
112      REAL(wp), DIMENSION(jpi,jpj,jp_medusa) :: btra
113
114      !! nitrogen and silicon production and consumption
115      REAL(wp) :: fn_prod, fn_cons, fs_prod, fs_cons
116
117      !! carbon, alkalinity production and consumption
118      REAL(wp) :: fc_prod, fc_cons, fa_prod, fa_cons
119
120      !! oxygen production and consumption (and non-consumption)
121      REAL(wp), DIMENSION(jpi,jpj) :: fo2_prod, fo2_cons
122      REAL(wp), DIMENSION(jpi,jpj) :: fo2_ncons, fo2_ccons
123
124      !! temporary variables
125      REAL(wp) :: fq0
126
127      !!==========================================================
128      !! LOCAL GRID CELL TRENDS
129      !!==========================================================
130      !!
131      !!----------------------------------------------------------
132      !! Determination of trends
133      !!----------------------------------------------------------
134      DO jj = 2,jpjm1
135         DO ji = 2,jpim1
136            !! OPEN wet point IF..THEN loop
137            if (tmask(ji,jj,jk) == 1) then
138               !!
139               !!----------------------------------------------------------
140               !! chlorophyll
141               btra(ji,jj,jpchn_lc) = b0 * ( ( (frn(ji,jj) * fprn(ji,jj) *      &
142                                             zphn(ji,jj) ) -                 &
143                                           fgmipn(ji,jj) - fgmepn(ji,jj) -   &
144                                           fdpn(ji,jj) - fdpn2(ji,jj) ) *    &
145                                          (fthetan(ji,jj) / xxi) )
146               btra(ji,jj,jpchd_lc) = b0 * ( ( (frd(ji,jj) * fprd(ji,jj) *      &
147                                             zphd(ji,jj) ) -                 &
148                                           fgmepd(ji,jj) - fdpd(ji,jj) -     &
149                                           fdpd2(ji,jj) ) *                  &
150                                          (fthetad(ji,jj) / xxi) )
151            ENDIF
152         ENDDO
153      ENDDO
154
155      DO jj = 2,jpjm1
156         DO ji = 2,jpim1
157            if (tmask(ji,jj,jk) == 1) then
158               !!
159               !!----------------------------------------------------------
160               !! phytoplankton
161               btra(ji,jj,jpphn_lc) = b0 * ( (fprn(ji,jj) * zphn(ji,jj)) -      &
162                                          fgmipn(ji,jj) - fgmepn(ji,jj) -    &
163                                          fdpn(ji,jj) - fdpn2(ji,jj) )
164               btra(ji,jj,jpphd_lc) = b0 * ( (fprd(ji,jj) * zphd(ji,jj)) -      &
165                                          fgmepd(ji,jj) - fdpd(ji,jj) -      &
166                                          fdpd2(ji,jj) )
167               btra(ji,jj,jppds_lc) = b0 * ( (fprds(ji,jj) * zpds(ji,jj)) -     &
168                                          fgmepds(ji,jj) - fdpds(ji,jj) -    &
169                                          fsdiss(ji,jj) - fdpds2(ji,jj) )
170            ENDIF
171         ENDDO
172      ENDDO
173
174      DO jj = 2,jpjm1
175         DO ji = 2,jpim1
176            if (tmask(ji,jj,jk) == 1) then
177               !!
178               !!----------------------------------------------------------
179               !! zooplankton
180               btra(ji,jj,jpzmi_lc) = b0 * (fmigrow(ji,jj) - fgmezmi(ji,jj) -   &
181                                         fdzmi(ji,jj) - fdzmi2(ji,jj))
182               btra(ji,jj,jpzme_lc) = b0 * (fmegrow(ji,jj) - fdzme(ji,jj) -     &
183                                         fdzme2(ji,jj))
184            ENDIF
185         ENDDO
186      ENDDO
187
188      !!----------------------------------------------------------
189      !! detritus
190      DO jj = 2,jpjm1
191         DO ji = 2,jpim1
192            if (tmask(ji,jj,jk) == 1) then
193               !!
194               btra(ji,jj,jpdet_lc) = b0 * (                           &
195                   fdpn(ji,jj)                                         & ! mort. losses
196                 + ((1.0 - xfdfrac1) * fdpd(ji,jj))                    & ! mort. losses
197                 + fdzmi(ji,jj)                                        & ! mort. losses
198                 + ((1.0 - xfdfrac2) * fdzme(ji,jj))                   & ! mort. losses
199                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))    & ! assim. inefficiency
200                 - fgmid(ji,jj) - fgmed(ji,jj)                         & ! grazing
201                 - fdd(ji,jj)                                          & ! remin.
202                 + fslowgain(ji,jj) - fslowloss(ji,jj)                 & ! slow-sinking
203                 - (f_sbenin_n(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor
204                 + ffast2slown(ji,jj) )                                  ! seafloor fast->slow
205            ENDIF
206         ENDDO
207      ENDDO
208
209      DO jj = 2,jpjm1
210         DO ji = 2,jpim1
211            if (tmask(ji,jj,jk) == 1) then
212               !!----------------------------------------------------------
213               !! dissolved inorganic nitrogen nutrient
214               !! primary production
215               fn_cons = - (fprn(ji,jj) * zphn(ji,jj)) -                     &
216                           (fprd(ji,jj) * zphd(ji,jj))
217               !!
218               fn_prod =                                                     &
219                                         ! messy feeding remin.
220                         (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) +           &
221                         (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) +            &
222                                  fgmezmi(ji,jj) + fgmed(ji,jj))) +          &
223                                         ! excretion and remin.
224                         fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) +      &
225                         freminn(ji,jj) +                                    &
226                                         ! metab. losses
227                         fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) +       &
228                         fdzme2(ji,jj)
229               !!
230               !! riverine flux
231               if ( jriver_n .gt. 0 ) then
232                  f_riv_loc_n(ji,jj) = f_riv_n(ji,jj) *                      &
233                     friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
234                  fn_prod = fn_prod + f_riv_loc_n(ji,jj)
235               endif
236               !! 
237               !! benthic remineralisation
238               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and.              &
239                   ibenthic.eq.1) then
240                  fn_prod = fn_prod + (f_benout_n(ji,jj) / fse3t(ji,jj,jk))
241               endif
242               !!
243               btra(ji,jj,jpdin_lc) = b0 * ( fn_prod + fn_cons )
244               !!
245               !! consumption of dissolved nitrogen
246               fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fse3t(ji,jj,jk) *     &
247                                                       fn_cons )
248               !! production of dissolved nitrogen
249               fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fse3t(ji,jj,jk) *     &
250                                                       fn_prod )
251            ENDIF
252         ENDDO
253      ENDDO
254
255      DO jj = 2,jpjm1
256         DO ji = 2,jpim1
257            if (tmask(ji,jj,jk) == 1) then
258               !!
259               !!----------------------------------------------------------
260               !! dissolved silicic acid nutrient
261               !! opal production
262               fs_cons = - (fprds(ji,jj) * zpds(ji,jj))
263               !!
264               fs_prod =                                                     &
265                             ! opal dissolution
266                         fsdiss(ji,jj) +                                     &
267                             ! mort. loss
268                         ((1.0 - xfdfrac1) * fdpds(ji,jj)) +                 &
269                             &  ! egestion of grazed Si
270                         ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +               &
271                             ! fast diss. and metab. losses
272                         freminsi(ji,jj) + fdpds2(ji,jj)
273               !!
274               !! riverine flux
275               if ( jriver_si .gt. 0 ) then
276                  f_riv_loc_si(ji,jj) = f_riv_si(ji,jj) *                    &
277                                        friver_dep(jk,mbathy(ji,jj)) /       &
278                                        fse3t(ji,jj,jk)
279                  fs_prod = fs_prod + f_riv_loc_si(ji,jj)
280               endif
281               !! 
282               !! benthic remineralisation
283               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and.            &
284                   ibenthic.eq.1) then
285                  fs_prod = fs_prod + (f_benout_si(ji,jj) / fse3t(ji,jj,jk))
286               endif
287               !!
288               btra(ji,jj,jpsil_lc) = b0 * ( &
289                 fs_prod + fs_cons )
290               !! consumption of dissolved silicon
291               fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fse3t(ji,jj,jk) *     &
292                                                       fs_cons )
293               !! production of dissolved silicon
294               fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fse3t(ji,jj,jk) *     &
295                                                       fs_prod )
296            ENDIF
297         ENDDO
298      ENDDO
299
300      DO jj = 2,jpjm1
301         DO ji = 2,jpim1
302            if (tmask(ji,jj,jk) == 1) then !!
303               !!----------------------------------------------------------
304               !! dissolved "iron" nutrient
305               btra(ji,jj,jpfer_lc) = b0 * ( (xrfn * btra(ji,jj,jpdin_lc)) +       &
306                                          ffetop(ji,jj) + ffebot(ji,jj) -    &
307                                          ffescav(ji,jj) )
308            ENDIF
309         ENDDO
310      ENDDO
311
312# if defined key_roam
313      !!----------------------------------------------------------
314      !! AXY (26/11/08): implicit detrital carbon change
315      DO jj = 2,jpjm1
316         DO ji = 2,jpim1
317            if (tmask(ji,jj,jk) == 1) then 
318               !!
319               btra(ji,jj,jpdtc_lc) = b0 * (                           &
320                   (xthetapn * fdpn(ji,jj))                            & ! mort. losses
321                 + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))       & ! mort. losses
322                 + (xthetazmi * fdzmi(ji,jj))                          & ! mort. losses
323                 + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))     & ! mort. losses
324                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))    & ! assim. inefficiency
325                 - fgmidc(ji,jj) - fgmedc(ji,jj)                       & ! grazing
326                 - fddc(ji,jj)                                         & ! remin.
327                 + fslowgainc(ji,jj) - fslowlossc(ji,jj)               & ! slow-sinking
328                 - (f_sbenin_c(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor
329                 + ffast2slowc(ji,jj) )                                  ! seafloor fast->slow
330            ENDIF
331         ENDDO
332      ENDDO
333
334      DO jj = 2,jpjm1
335         DO ji = 2,jpim1
336            if (tmask(ji,jj,jk) == 1) then
337               !!
338               !!----------------------------------------------------------
339               !! dissolved inorganic carbon
340               !! primary production
341               fc_cons = - (xthetapn * fprn(ji,jj) * zphn(ji,jj)) -          &
342                           (xthetapd * fprd(ji,jj) * zphd(ji,jj))
343               !!
344               fc_prod =                                                     &
345                            ! messy feeding remin
346                         (xthetapn * xphi * fgmipn(ji,jj)) +                 &
347                         (xphi * fgmidc(ji,jj)) +                            &
348                         (xthetapn * xphi * fgmepn(ji,jj)) +                 &
349                         (xthetapd * xphi * fgmepd(ji,jj)) +                 &
350                         (xthetazmi * xphi * fgmezmi(ji,jj)) +               &
351                         (xphi * fgmedc(ji,jj)) +                            &
352                            ! resp., remin., losses
353                         fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) +     &
354                         freminc(ji,jj) + (xthetapn * fdpn2(ji,jj)) +        &
355                            ! losses
356                         (xthetapd * fdpd2(ji,jj)) +                         &
357                         (xthetazmi * fdzmi2(ji,jj)) +                       &
358                         (xthetazme * fdzme2(ji,jj))
359               !!
360               !! riverine flux
361               if ( jriver_c .gt. 0 ) then
362                  f_riv_loc_c(ji,jj) = f_riv_c(ji,jj) *                      &
363                                       friver_dep(jk,mbathy(ji,jj)) /        &
364                                       fse3t(ji,jj,jk)
365                  fc_prod = fc_prod + f_riv_loc_c(ji,jj)
366               endif
367               !! 
368               !! benthic remineralisation
369               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and.              &
370                   ibenthic.eq.1) then
371                  fc_prod = fc_prod + (f_benout_c(ji,jj) / fse3t(ji,jj,jk))
372               endif
373               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and.            &
374                   ibenthic.eq.1) then
375                  fc_prod = fc_prod + (f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
376               endif
377               !!
378               !! community respiration (does not include CaCO3 terms -
379               !! obviously!)
380               fcomm_resp(ji,jj) = fcomm_resp(ji,jj) + fc_prod
381               !!
382               !! CaCO3
383               fc_prod = fc_prod - ftempca(ji,jj) + freminca(ji,jj)
384               !!
385               !! riverine flux
386               if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then
387                  fc_prod = fc_prod + f_riv_c(ji,jj)
388               endif
389               !!
390               btra(ji,jj,jpdic_lc) = b0 * ( fc_prod + fc_cons )
391               !! consumption of dissolved carbon
392               fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fse3t(ji,jj,jk) *     &
393                                                       fc_cons )
394               !! production of dissolved carbon
395               fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fse3t(ji,jj,jk) *     &
396                                                       fc_prod )
397            ENDIF
398         ENDDO
399      ENDDO
400
401      DO jj = 2,jpjm1
402         DO ji = 2,jpim1
403            if (tmask(ji,jj,jk) == 1) then
404               !!
405               !!----------------------------------------------------------
406               !! alkalinity
407               !! CaCO3 dissolution
408               fa_prod = 2.0 * freminca(ji,jj)
409               !! CaCO3 production
410               fa_cons = - 2.0 * ftempca(ji,jj)
411               !!
412               !! riverine flux
413               if ( jriver_alk .gt. 0 ) then
414                  f_riv_loc_alk(ji,jj) = f_riv_alk(ji,jj) *                  &
415                                         friver_dep(jk,mbathy(ji,jj)) /      &
416                                         fse3t(ji,jj,jk)
417                  fa_prod = fa_prod + f_riv_loc_alk(ji,jj)
418               endif
419               !! 
420               !! benthic remineralisation
421               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and.            &
422                   ibenthic.eq.1) then
423                  fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) /            &
424                                       fse3t(ji,jj,jk))
425               endif
426               !!
427               btra(ji,jj,jpalk_lc) = b0 * ( fa_prod + fa_cons )
428            ENDIF
429         ENDDO
430      ENDDO
431
432      DO jj = 2,jpjm1
433         DO ji = 2,jpim1
434            if (tmask(ji,jj,jk) == 1) then
435               !!
436               !!----------------------------------------------------------
437               !! oxygen (has protection at low O2 concentrations;
438               !! OCMIP-2 style)
439               fo2_prod(ji,jj) =                                             &
440                                     ! Pn primary production, N
441                                 (xthetanit * fprn(ji,jj) * zphn(ji,jj)) +   &
442                                     ! Pd primary production, N
443                                 (xthetanit * fprd(ji,jj) * zphd(ji,jj)) +   &
444                                     ! Pn primary production, C
445                                 (xthetarem * xthetapn * fprn(ji,jj) *       &
446                                  zphn(ji,jj)) +                             &
447                                     ! Pd primary production, C
448                                  (xthetarem * xthetapd * fprd(ji,jj) *      &
449                                   zphd(ji,jj))
450               fo2_ncons(ji,jj) =                                            &
451                                     ! Pn messy feeding remin., N
452                                   - (xthetanit * xphi * fgmipn(ji,jj))      &
453                                     ! D  messy feeding remin., N
454                                   - (xthetanit * xphi * fgmid(ji,jj))       &
455                                     ! Pn messy feeding remin., N
456                                   - (xthetanit * xphi * fgmepn(ji,jj))      &
457                                     ! Pd messy feeding remin., N
458                                   - (xthetanit * xphi * fgmepd(ji,jj))      &
459                                     ! Zi messy feeding remin., N
460                                   - (xthetanit * xphi * fgmezmi(ji,jj))     &
461                                     ! D  messy feeding remin., N
462                                   - (xthetanit * xphi * fgmed(ji,jj))       &
463                                     ! microzoo excretion, N
464                                   - (xthetanit * fmiexcr(ji,jj))            &
465                                     ! mesozoo  excretion, N
466                                   - (xthetanit * fmeexcr(ji,jj))            &
467                                     ! slow detritus remin., N
468                                   - (xthetanit * fdd(ji,jj))                &
469                                     ! fast detritus remin., N
470                                   - (xthetanit * freminn(ji,jj))            &
471                                     ! Pn  losses, N
472                                   - (xthetanit * fdpn2(ji,jj))              &
473                                     ! Pd  losses, N
474                                   - (xthetanit * fdpd2(ji,jj))              &
475                                     ! Zmi losses, N
476                                   - (xthetanit * fdzmi2(ji,jj))             &
477                                     ! Zme losses, N
478                                   - (xthetanit * fdzme2(ji,jj))
479               !! 
480               !! benthic remineralisation
481               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and.              &
482                   ibenthic.eq.1) then
483                  fo2_ncons(ji,jj) = fo2_ncons(ji,jj) -                      &
484                                     (xthetanit * f_benout_n(ji,jj) /        &
485                                      fse3t(ji,jj,jk))
486               endif
487            ENDIF
488         ENDDO
489      ENDDO
490
491      DO jj = 2,jpjm1
492         DO ji = 2,jpim1
493            if (tmask(ji,jj,jk) == 1) then
494               fo2_ccons(ji,jj) =                                            &
495                                     ! Pn messy feeding remin., C
496                                  - (xthetarem * xthetapn * xphi *           &
497                                     fgmipn(ji,jj))                          &
498                                     ! D  messy feeding remin., C
499                                  - (xthetarem * xphi * fgmidc(ji,jj))       &
500                                     ! Pn messy feeding remin., C
501                                  - (xthetarem * xthetapn * xphi *           &
502                                     fgmepn(ji,jj))                          &
503                                     ! Pd messy feeding remin., C
504                                  - (xthetarem * xthetapd * xphi *           &
505                                     fgmepd(ji,jj))                          &
506                                     ! Zi messy feeding remin., C
507                                  - (xthetarem * xthetazmi * xphi *          &
508                                     fgmezmi(ji,jj))                         &
509                                     ! D  messy feeding remin., C
510                                  - (xthetarem * xphi * fgmedc(ji,jj))       &
511                                     ! microzoo respiration, C
512                                  - (xthetarem * fmiresp(ji,jj))             &
513                                     ! mesozoo  respiration, C
514                                  - (xthetarem * fmeresp(ji,jj))             &
515                                     ! slow detritus remin., C
516                                  - (xthetarem * fddc(ji,jj))                &
517                                     ! fast detritus remin., C
518                                  - (xthetarem * freminc(ji,jj))             &
519                                     ! Pn  losses, C
520                                  - (xthetarem * xthetapn * fdpn2(ji,jj))    &
521                                     ! Pd  losses, C
522                                  - (xthetarem * xthetapd * fdpd2(ji,jj))    &
523                                     ! Zmi losses, C
524                                  - (xthetarem * xthetazmi * fdzmi2(ji,jj))  &
525                                     ! Zme losses, C
526                                  - (xthetarem * xthetazme * fdzme2(ji,jj))
527               !! 
528               !! benthic remineralisation
529               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and.              &
530                   ibenthic.eq.1) then
531                  fo2_ccons(ji,jj) = fo2_ccons(ji,jj) - (xthetarem *         &
532                                                         f_benout_c(ji,jj) / &
533                                                         fse3t(ji,jj,jk))
534               endif
535               fo2_cons(ji,jj) = fo2_ncons(ji,jj) + fo2_ccons(ji,jj)
536            ENDIF
537         ENDDO
538      ENDDO
539
540      DO jj = 2,jpjm1
541         DO ji = 2,jpim1
542            if (tmask(ji,jj,jk) == 1) then
543               !!
544               !! is this a suboxic zone?
545               !! deficient O2; production fluxes only
546               if (zoxy(ji,jj).lt.xo2min) then
547                  btra(ji,jj,jpoxy_lc) = b0 * fo2_prod(ji,jj)
548                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) *  &
549                                                          fo2_prod(ji,jj) )
550                  foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fse3t(ji,jj,jk) *  &
551                                                          fo2_cons(ji,jj) )
552               else
553                  !! sufficient O2; production + consumption fluxes
554                  btra(ji,jj,jpoxy_lc) = b0 * ( fo2_prod(ji,jj) +               &
555                                             fo2_cons(ji,jj) )
556                  foxy_prod(ji,jj) = foxy_prod(ji,jj) +                      &
557                                     ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
558                  foxy_cons(ji,jj) = foxy_cons(ji,jj) +                      &
559                                     ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
560               endif
561            ENDIF
562         ENDDO
563      ENDDO
564
565      DO jj = 2,jpjm1
566         DO ji = 2,jpim1
567            if (tmask(ji,jj,jk) == 1) then
568               !!
569               !! air-sea fluxes (if this is the surface box)
570               if (jk.eq.1) then
571                  !!
572                  !! CO2 flux
573                  btra(ji,jj,jpdic_lc) = btra(ji,jj,jpdic_lc) + (b0 *              &
574                                                           f_co2flux(ji,jj))
575                  !!
576                  !! O2 flux (mol/m3/s -> mmol/m3/d)
577                  btra(ji,jj,jpoxy_lc) = btra(ji,jj,jpoxy_lc) + (b0 *              &
578                                                           f_o2flux(ji,jj))
579               endif
580            ENDIF
581         ENDDO
582      ENDDO
583# endif
584
585# if defined key_debug_medusa
586! I DON'T THIS IS MUCH USE, NOW IT'S BEEN PULLED OUT OF THE MAIN DO LOOP
587! - marc 5/5/17
588      DO jj = 2,jpjm1
589         DO ji = 2,jpim1
590            if (tmask(ji,jj,jk) == 1) then
591               !! report state variable fluxes (not including
592               !! fast-sinking detritus)
593               if (idf.eq.1.AND.idfval.eq.1) then
594                  IF (lwp) write (numout,*) '------------------------------'
595                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchn_lc)(',jk,')  = ', &
596                                            btra(ji,jj,jpchn_lc)
597                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchd_lc)(',jk,')  = ', &
598                                            btra(ji,jj,jpchd_lc)
599                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphn_lc)(',jk,')  = ', &
600                                            btra(ji,jj,jpphn_lc)
601                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphd_lc)(',jk,')  = ', &
602                                            btra(ji,jj,jpphd_lc)
603                  IF (lwp) write (numout,*) 'btra(ji,jj,jppds_lc)(',jk,')  = ', &
604                                            btra(ji,jj,jppds_lc)
605                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzmi_lc)(',jk,')  = ', &
606                                            btra(ji,jj,jpzmi_lc)
607                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzme_lc)(',jk,')  = ', &
608                                            btra(ji,jj,jpzme_lc)
609                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdet_lc)(',jk,')  = ', &
610                                            btra(ji,jj,jpdet_lc)
611                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdin_lc)(',jk,')  = ', &
612                                            btra(ji,jj,jpdin_lc)
613                  IF (lwp) write (numout,*) 'btra(ji,jj,jpsil_lc)(',jk,')  = ', &
614                                            btra(ji,jj,jpsil_lc)
615                  IF (lwp) write (numout,*) 'btra(ji,jj,jpfer_lc)(',jk,')  = ', &
616                                            btra(ji,jj,jpfer_lc)
617#  if defined key_roam
618                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdtc_lc)(',jk,')  = ', &
619                                            btra(ji,jj,jpdtc_lc)
620                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdic_lc)(',jk,')  = ', &
621                                            btra(ji,jj,jpdic_lc)
622                  IF (lwp) write (numout,*) 'btra(ji,jj,jpalk_lc)(',jk,')  = ', &
623                                            btra(ji,jj,jpalk_lc)
624                  IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy_lc)(',jk,')  = ', &
625                                            btra(ji,jj,jpoxy_lc)
626#  endif
627               endif
628            ENDIF
629         ENDDO
630      ENDDO
631# endif
632
633      !!----------------------------------------------------------
634      !! Integrate calculated fluxes for mass balance
635      !!----------------------------------------------------------
636      DO jj = 2,jpjm1
637         DO ji = 2,jpim1
638            if (tmask(ji,jj,jk) == 1) then
639               !!
640               !! === nitrogen ===
641               fflx_n(ji,jj)  = fflx_n(ji,jj) + fse3t(ji,jj,jk) *            &
642                                ( btra(ji,jj,jpphn_lc) + btra(ji,jj,jpphd_lc) +    &
643                                  btra(ji,jj,jpzmi_lc) + btra(ji,jj,jpzme_lc) +    &
644                                  btra(ji,jj,jpdet_lc) + btra(ji,jj,jpdin_lc) )
645               !! === silicon ===
646               fflx_si(ji,jj) = fflx_si(ji,jj) + fse3t(ji,jj,jk) *           &
647                                ( btra(ji,jj,jppds_lc) + btra(ji,jj,jpsil_lc) )
648               !! === iron ===
649               fflx_fe(ji,jj) = fflx_fe(ji,jj) + fse3t(ji,jj,jk) *           &
650                                ( (xrfn *                                    &
651                                   (btra(ji,jj,jpphn_lc) + btra(ji,jj,jpphd_lc) +  &
652                                    btra(ji,jj,jpzmi_lc) + btra(ji,jj,jpzme_lc) +  &
653                                    btra(ji,jj,jpdet_lc))) + btra(ji,jj,jpfer_lc) )
654# if defined key_roam
655               !! === carbon ===
656               fflx_c(ji,jj)  = fflx_c(ji,jj) + fse3t(ji,jj,jk) *            &
657                                ( (xthetapn * btra(ji,jj,jpphn_lc)) +           &
658                                  (xthetapd * btra(ji,jj,jpphd_lc)) +           &
659                                  (xthetazmi * btra(ji,jj,jpzmi_lc)) +          &
660                                  (xthetazme * btra(ji,jj,jpzme_lc)) +          &
661                                  btra(ji,jj,jpdtc_lc) + btra(ji,jj,jpdic_lc) )
662               !! === alkalinity ===
663               fflx_a(ji,jj)  = fflx_a(ji,jj) + fse3t(ji,jj,jk) *            &
664                                btra(ji,jj,jpalk_lc)
665               !! === oxygen ===
666               fflx_o2(ji,jj) = fflx_o2(ji,jj) + fse3t(ji,jj,jk) *           &
667                                btra(ji,jj,jpoxy_lc)
668# endif
669            ENDIF
670         ENDDO
671      ENDDO
672
673      !!----------------------------------------------------------
674      !! Apply calculated tracer fluxes
675      !!----------------------------------------------------------
676      !!
677      !! units: [unit of tracer] per second (fluxes are calculated
678      !! above per day)
679      !!
680      DO jj = 2,jpjm1
681         DO ji = 2,jpim1
682            if (tmask(ji,jj,jk) == 1) then
683               ibio_switch = 1
684# if defined key_gulf_finland
685               !! AXY (17/05/13): fudge in a Gulf of Finland correction;
686               !!                 uses longitude-latitude range to
687               !!                 establish if this is a Gulf of Finland
688               !!                 grid cell; if so, then BGC fluxes are
689               !!                 ignored (though still calculated); for
690               !!                 reference, this is meant to be a
691               !!                 temporary fix to see if all of my
692               !!                 problems can be done away with if I
693               !!                 switch off BGC fluxes in the Gulf of
694               !!                 Finland, which currently appears the
695               !!                 source of trouble
696               if ( glamt(ji,jj).gt.24.7 .and. glamt(ji,jj).lt.27.8 .and.    &
697                    gphit(ji,jj).gt.59.2 .and. gphit(ji,jj).lt.60.2 ) then
698                  ibio_switch = 0
699               endif
700# endif               
701               if (ibio_switch.eq.1) then
702                  tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) +                &
703                                        (btra(ji,jj,jpchn_lc) / 86400.)
704                  tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) +                &
705                                        (btra(ji,jj,jpchd_lc) / 86400.)
706                  tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) +                &
707                                        (btra(ji,jj,jpphn_lc) / 86400.)
708                  tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) +                &
709                                        (btra(ji,jj,jpphd_lc) / 86400.)
710                  tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) +                &
711                                        (btra(ji,jj,jppds_lc) / 86400.)
712                  tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) +                &
713                                        (btra(ji,jj,jpzmi_lc) / 86400.)
714                  tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) +                &
715                                        (btra(ji,jj,jpzme_lc) / 86400.)
716                  tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) +                &
717                                        (btra(ji,jj,jpdet_lc) / 86400.)
718                  tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) +                &
719                                        (btra(ji,jj,jpdin_lc) / 86400.)
720                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +                &
721                                        (btra(ji,jj,jpsil_lc) / 86400.)
722                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +                &
723                                        (btra(ji,jj,jpfer_lc) / 86400.)
724# if defined key_roam
725                  tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) +                &
726                                        (btra(ji,jj,jpdtc_lc) / 86400.)
727                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +                &
728                                        (btra(ji,jj,jpdic_lc) / 86400.)
729                  tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) +                &
730                                        (btra(ji,jj,jpalk_lc) / 86400.)
731                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) +                &
732                                        (btra(ji,jj,jpoxy_lc) / 86400.)
733# endif
734               endif
735            ENDIF
736         ENDDO
737      ENDDO
738
739      DO jj = 2,jpjm1
740         DO ji = 2,jpim1
741            if (tmask(ji,jj,jk) == 1) then
742
743               !! AXY (18/11/16): CMIP6 diagnostics
744               IF( med_diag%FBDDTALK%dgsave )  THEN
745                  fbddtalk(ji,jj)  =  fbddtalk(ji,jj)  +                     &
746                                      (btra(ji,jj,jpalk_lc) * fse3t(ji,jj,jk))
747               ENDIF
748               IF( med_diag%FBDDTDIC%dgsave )  THEN
749                  fbddtdic(ji,jj)  =  fbddtdic(ji,jj)  +                     &
750                                      (btra(ji,jj,jpdic_lc) * fse3t(ji,jj,jk))
751               ENDIF
752               IF( med_diag%FBDDTDIFE%dgsave ) THEN
753                  fbddtdife(ji,jj) =  fbddtdife(ji,jj) +                     &
754                                      (btra(ji,jj,jpfer_lc) * fse3t(ji,jj,jk))
755               ENDIF
756               IF( med_diag%FBDDTDIN%dgsave )  THEN
757                  fbddtdin(ji,jj)  =  fbddtdin(ji,jj)  +                     &
758                                      (btra(ji,jj,jpdin_lc) * fse3t(ji,jj,jk))
759               ENDIF
760               IF( med_diag%FBDDTDISI%dgsave ) THEN
761                  fbddtdisi(ji,jj) =  fbddtdisi(ji,jj) +                     &
762                                      (btra(ji,jj,jpsil_lc) * fse3t(ji,jj,jk))
763               ENDIF
764          !!
765               IF( med_diag%BDDTALK3%dgsave )  THEN
766                  bddtalk3(ji,jj,jk)  =  btra(ji,jj,jpalk_lc)
767               ENDIF
768               IF( med_diag%BDDTDIC3%dgsave )  THEN
769                  bddtdic3(ji,jj,jk)  =  btra(ji,jj,jpdic_lc)
770               ENDIF
771               IF( med_diag%BDDTDIFE3%dgsave ) THEN
772                  bddtdife3(ji,jj,jk) =  btra(ji,jj,jpfer_lc)
773               ENDIF
774               IF( med_diag%BDDTDIN3%dgsave )  THEN
775                  bddtdin3(ji,jj,jk)  =  btra(ji,jj,jpdin_lc)
776               ENDIF
777               IF( med_diag%BDDTDISI3%dgsave ) THEN
778                  bddtdisi3(ji,jj,jk) =  btra(ji,jj,jpsil_lc)
779               ENDIF
780            ENDIF
781         ENDDO
782      ENDDO
783
784#   if defined key_debug_medusa
785      IF (lwp) write (numout,*) '------'
786      IF (lwp) write (numout,*) 'bio_medusa_update: end all calculations'
787      IF (lwp) write (numout,*) 'bio_medusa_update: now outputs kt = ', kt
788      CALL flush(numout)
789#   endif
790
791# if defined key_axy_nancheck
792      !!----------------------------------------------------------
793      !! Check calculated tracer fluxes
794      !!----------------------------------------------------------
795      DO jj = 2,jpjm1
796         DO ji = 2,jpim1
797            if (tmask(ji,jj,jk) == 1) then
798               !!
799               DO jn = 1,jp_medusa
800                  fq0 = btra(ji,jj,jn)
801                  !! AXY (30/01/14): "isnan" problem on HECTOR
802                  !! if (fq0 /= fq0 ) then
803                  if ( ieee_is_nan( fq0 ) ) then
804                     !! there's a NaN here
805                     if (lwp) write(numout,*) 'NAN detected in btra(ji,jj,',  &
806                        ji, ',', jj, ',', jk, ',', jn, ') at time', kt
807           CALL ctl_stop( 'trcbio_medusa, NAN in btra field' )
808                  endif
809               ENDDO
810               DO jn = jp_msa0,jp_msa1
811                  fq0 = tra(ji,jj,jk,jn)
812                  !! AXY (30/01/14): "isnan" problem on HECTOR
813                  !! if (fq0 /= fq0 ) then
814                  if ( ieee_is_nan( fq0 ) ) then
815                     !! there's a NaN here
816                     if (lwp) write(numout,*) 'NAN detected in tra(', ji, &
817                        ',', jj, ',', jk, ',', jn, ') at time', kt
818              CALL ctl_stop( 'trcbio_medusa, NAN in tra field' )
819                  endif
820               ENDDO
821               CALL flush(numout)
822            ENDIF
823         ENDDO
824      ENDDO
825# endif
826
827
828   END SUBROUTINE bio_medusa_update
829
830#else
831   !!======================================================================
832   !!  Dummy module :                                   No MEDUSA bio-model
833   !!======================================================================
834CONTAINS
835   SUBROUTINE bio_medusa_update( )                    ! Empty routine
836      WRITE(*,*) 'bio_medusa_update: You should not have seen this print! error?'
837   END SUBROUTINE bio_medusa_update
838#endif 
839
840   !!======================================================================
841END MODULE bio_medusa_update_mod
Note: See TracBrowser for help on using the repository browser.