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

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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.

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