source: branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90 @ 8419

Last change on this file since 8419 was 8419, checked in by jpalmier, 3 years ago

JPALM — 08-08-2017 — clean trcbio_medusa split only - running this time

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