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

source: branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90 @ 11726

Last change on this file since 11726 was 11726, checked in by jpalmier, 5 years ago

add Changes to enable Alkalinity geo-ingeneering tests within MEDUSA

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