source: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90 @ 9070

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

merge with GO6 8356:8447

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