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

Last change on this file since 8442 was 8442, checked in by frrh, 3 years ago

Commit changes relating to Met Office GMED ticket 340 for the
tidying of MEDUSA related code and debugging statements in the TOP code.

Only code introduced at revision 8434 of branch
http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
is included here, all previous revisions of that branch having been dealt with
under GMED ticket 339.

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.