New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bio_medusa_update.F90 in branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90 @ 7996

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

Pulled the updating of tracers out of trcbio_medusa.F90

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