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

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

GMED ticket 406. CPP key fixes.

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