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.
trcbio_medusa.F90 in branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 7986

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

Pulled detritus processes out of trcbio_medusa.F90

File size: 119.9 KB
Line 
1MODULE trcbio_medusa
2   !!======================================================================
3   !!                         ***  MODULE trcbio  ***
4   !! TOP :   MEDUSA
5   !!======================================================================
6   !! History :
7   !!  -   !  1999-07  (M. Levy)              original code
8   !!  -   !  2000-12  (E. Kestenare)         assign parameters to name individual tracers
9   !!  -   !  2001-03  (M. Levy)              LNO3 + dia2d
10   !! 2.0  !  2007-12  (C. Deltel, G. Madec)  F90
11   !!  -   !  2008-08  (K. Popova)            adaptation for MEDUSA
12   !!  -   !  2008-11  (A. Yool)              continuing adaptation for MEDUSA
13   !!  -   !  2010-03  (A. Yool)              updated for branch inclusion
14   !!  -   !  2011-08  (A. Yool)              updated for ROAM (see below)
15   !!  -   !  2013-03  (A. Yool)              updated for iMARNET
16   !!  -   !  2013-05  (A. Yool)              updated for v3.5
17   !!  -   !  2014-08  (A. Yool, J. Palm)     Add DMS module for UKESM1 model
18   !!  -   !  2015-06  (A. Yool)              Update to include MOCSY
19   !!  -   !  2015-07  (A. Yool)              Update for rolling averages
20   !!  -   !  2015-10  (J. Palm)              Update for diag outputs through iom_use 
21   !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6
22   !!----------------------------------------------------------------------
23   !!
24#if defined key_roam
25   !!----------------------------------------------------------------------
26   !! Updates for the ROAM project include:
27   !!   - addition of DIC, alkalinity, detrital carbon and oxygen tracers
28   !!   - addition of air-sea fluxes of CO2 and oxygen
29   !!   - periodic (monthly) calculation of full 3D carbonate chemistry
30   !!   - detrital C:N ratio now free to evolve dynamically
31   !!   - benthic storage pools
32   !!
33   !! Opportunity also taken to add functionality:
34   !!   - switch for Liebig Law (= most-limiting) nutrient uptake
35   !!   - switch for accelerated seafloor detritus remineralisation
36   !!   - switch for fast -> slow detritus transfer at seafloor
37   !!   - switch for ballast vs. Martin vs. Henson fast detritus remin.
38   !!   - per GMD referee remarks, xfdfrac3 introduced for grazed PDS
39   !!----------------------------------------------------------------------
40#endif
41   !!
42#if defined key_mocsy
43   !!----------------------------------------------------------------------
44   !! Updates with the addition of MOCSY include:
45   !!   - option to use PML or MOCSY carbonate chemistry (the latter is
46   !!     preferred)
47   !!   - central calculation of gas transfer velocity, f_kw660; previously
48   !!     this was done separately for CO2 and O2 with predictable results
49   !!   - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux
50   !!     calculations and to those for O2 air-sea flux
51   !!   - extra diagnostics included for MOCSY
52   !!----------------------------------------------------------------------
53#endif
54   !!
55#if defined key_medusa
56   !!----------------------------------------------------------------------
57   !!                                        MEDUSA bio-model
58   !!----------------------------------------------------------------------
59   !!   trc_bio_medusa        : 
60   !!----------------------------------------------------------------------
61      USE oce_trc
62      USE trc
63      USE sms_medusa
64      USE lbclnk
65      USE prtctl_trc      ! Print control for debugging
66      USE trcsed_medusa
67      USE sbc_oce         ! surface forcing
68      USE sbcrnf          ! surface boundary condition: runoff variables
69      USE in_out_manager  ! I/O manager
70# if defined key_iomput
71      USE iom
72      USE trcnam_medusa         ! JPALM 13-11-2015 -- if iom_use for diag
73      !!USE trc_nam_iom_medusa  ! JPALM 13-11-2015 -- if iom_use for diag
74# endif
75# if defined key_roam
76      USE gastransfer
77#  if defined key_mocsy
78      USE mocsy_wrapper
79#  else
80      USE trcco2_medusa
81#  endif
82      USE trcoxy_medusa
83      !! Jpalm (08/08/2014)
84      USE trcdms_medusa
85# endif
86      !! AXY (18/01/12): brought in for benthic timestepping
87      USE trcnam_trp      ! AXY (24/05/2013)
88      USE trdmxl_trc
89      USE trdtrc_oce  ! AXY (24/05/2013)
90
91      !! AXY (30/01/14): necessary to find NaNs on HECTOR
92      USE, INTRINSIC :: ieee_arithmetic 
93
94      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
95      USE sbc_oce,                    ONLY: lk_oasis
96      USE oce,                        ONLY: CO2Flux_out_cpl, DMS_out_cpl,   &
97                                            PCO2a_in_cpl
98      USE bio_medusa_mod
99      USE bio_medusa_init_mod,        ONLY: bio_medusa_init
100      USE carb_chem_mod,              ONLY: carb_chem
101      USE air_sea_mod,                ONLY: air_sea
102      USE plankton_mod,               ONLY: plankton
103      USE iron_chem_scav_mod,         ONLY: iron_chem_scav
104      USE detritus_mod,               ONLY: detritus
105      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
106      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
107
108      IMPLICIT NONE
109      PRIVATE
110     
111      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
112
113   !!* Substitution
114#  include "domzgr_substitute.h90"
115   !!----------------------------------------------------------------------
116   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
117   !! $Id$
118   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
119   !!----------------------------------------------------------------------
120
121CONTAINS
122
123   SUBROUTINE trc_bio_medusa( kt )
124      !!---------------------------------------------------------------------
125      !!                     ***  ROUTINE trc_bio  ***
126      !!
127      !! ** Purpose :   compute the now trend due to biogeochemical processes
128      !!              and add it to the general trend of passive tracers equations
129      !!
130      !! ** Method  :   each now biological flux is calculated in function of now
131      !!              concentrations of tracers.
132      !!              depending on the tracer, these fluxes are sources or sinks.
133      !!              the total of the sources and sinks for each tracer
134      !!              is added to the general trend.
135      !!       
136      !!                      tra = tra + zf...tra - zftra...
137      !!                                     |         |
138      !!                                     |         |
139      !!                                  source      sink
140      !!       
141      !!              IF 'key_trc_diabio' defined , the biogeochemical trends
142      !!              for passive tracers are saved for futher diagnostics.
143      !!---------------------------------------------------------------------
144      !!
145      !!
146      !!----------------------------------------------------------------------           
147      !! Variable conventions
148      !!----------------------------------------------------------------------
149      !!
150      !! names: z*** - state variable
151      !!        f*** - function (or temporary variable used in part of a function)
152      !!        x*** - parameter
153      !!        b*** - right-hand part (sources and sinks)
154      !!        i*** - integer variable (usually used in yes/no flags)
155      !!
156      !! time (integer timestep)
157      INTEGER, INTENT( in ) ::    kt
158      !!
159      !! spatial array indices
160      INTEGER  ::    ji,jj,jk,jn
161      !!
162      !! AXY (27/07/10): add in indices for depth horizons (for sinking flux
163      !!                 and seafloor iron inputs)
164      !! INTEGER  ::    i0100, i0200, i0500, i1000, i1100
165      !!
166      !! model state variables
167!      REAL(wp), DIMENSION(jpi,jpj) ::    zchn,zchd,zphn,zphd,zpds,zzmi
168!      REAL(wp), DIMENSION(jpi,jpj) ::    zzme,zdet,zdtc,zdin,zsil,zfer
169! zage doesn't seem to be used - marc 19/4/17
170!      REAL(wp) ::    zage
171!# if defined key_roam
172!      REAL(wp), DIMENSION(jpi,jpj) ::    zdic, zalk, zoxy
173!      REAL(wp), DIMENSION(jpi,jpj) ::    ztmp, zsal
174!# endif
175!# if defined key_mocsy
176!      REAL(wp), DIMENSION(jpi,jpj) ::    zpho
177!# endif
178      !!
179      !! integrated source and sink terms
180!      REAL(wp) ::    b0
181      !! AXY (23/08/13): changed from individual variables for each flux to
182      !!                 an array that holds all fluxes
183      REAL(wp), DIMENSION(jpi,jpj,jp_medusa) ::    btra
184      !!
185      !! primary production and chl related quantities     
186!      REAL(wp), DIMENSION(jpi,jpj) ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
187!      REAL(wp), DIMENSION(jpi,jpj) ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
188      !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range)
189!      REAL(wp), DIMENSION(jpi,jpj) ::    fjlim_pn, fjlim_pd
190      !! AXY (03/02/11): add in Liebig terms
191!      REAL(wp), DIMENSION(jpi,jpj) ::    fpnlim, fpdlim
192      !! AXY (16/07/09): add in Eppley curve functionality
193!      REAL(wp), DIMENSION(jpi,jpj) ::    fun_T,xvpnT,xvpdT
194      INTEGER  ::    ieppley
195      !! AXY (16/05/11): per Katya's prompting, add in new T-dependence
196      !!                 for phytoplankton growth only (i.e. no change
197      !!                 for remineralisation)
198!      REAL(wp), DIMENSION(jpi,jpj) ::    fun_Q10
199      !! AXY (01/03/10): add in mixed layer PP diagnostics
200!      REAL(wp), DIMENSION(jpi,jpj) ::    fprn_ml,fprd_ml
201      !!
202      !! nutrient limiting factors
203!      REAL(wp), DIMENSION(jpi,jpj) ::    fnln,ffln2            !! N and Fe
204!      REAL(wp), DIMENSION(jpi,jpj) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si
205      !!
206      !! silicon cycle
207!      REAL(wp), DIMENSION(jpi,jpj) ::    fsin,fnsi,fprds,fsdiss
208      REAL(wp)                     ::    fsin1,fnsi1,fnsi2
209      !!
210      !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
211!      REAL(wp), DIMENSION(jpi,jpj) ::    ffetop,ffebot,ffescav
212      REAL(wp) ::    xLgF, xFeT, xFeF, xFeL         !! state variables for iron-ligand system
213!      REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system
214      REAL(wp) ::    xb_coef_tmp, xb2M4ac           !! iron-ligand parameters
215      REAL(wp) ::    xmaxFeF,fdeltaFe               !! max Fe' parameters
216      !!
217      !! local parameters for Moore et al. (2004) alternative scavenging scheme
218      REAL(wp) ::    fbase_scav,fscal_sink,fscal_part,fscal_scav
219      !!
220      !! local parameters for Moore et al. (2008) alternative scavenging scheme
221      REAL(wp) ::    fscal_csink,fscal_sisink,fscal_casink
222      !!
223      !! local parameters for Galbraith et al. (2010) alternative scavenging scheme
224      REAL(wp) ::    xCscav1, xCscav2, xk_org, xORGscav  !! organic portion of scavenging
225      REAL(wp) ::    xk_inorg, xINORGscav                !! inorganic portion of scavenging
226      !!
227      !! microzooplankton grazing
228!      REAL(wp), DIMENSION(jpi,jpj) ::    fmi1,fmi,fgmipn,fgmid,fgmidc
229!      REAL(wp), DIMENSION(jpi,jpj) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
230      !!
231      !! mesozooplankton grazing
232!      REAL(wp), DIMENSION(jpi,jpj) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
233!      REAL(wp), DIMENSION(jpi,jpj) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
234      !!
235      !! mortality/Remineralisation (defunct parameter "fz" removed)
236!      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
237# if defined key_roam
238!      REAL(wp), DIMENSION(jpi,jpj) ::    fddc
239# endif
240!      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
241!      REAL(wp), DIMENSION(jpi,jpj) ::    fslown, fslowc
242!      REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux
243      REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi
244!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi
245# if defined key_roam
246!! Doesn't look like this is used - marc 10/4/17
247!!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenc
248!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfastc
249# endif
250      !!
251      !! particle flux
252!      REAL(WP), DIMENSION(jpi,jpj) ::    fdep1,fcaco3
253!      REAL(WP), DIMENSION(jpi,jpj) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
254!      REAL(wp), DIMENSION(jpi,jpj) ::    freminn,freminsi,freminfe,freminc,freminca
255!      REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca
256!      REAL(wp), DIMENSION(jpi,jpj) ::    fprotf
257!      REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca
258!      REAL(wp), DIMENSION(jpi,jpj) ::    fccd
259!      REAL(wp), DIMENSION(jpi,jpj) ::    fccd_dep
260      !!
261      !! AXY (06/07/11): alternative fast detritus schemes
262      REAL(wp) ::    fb_val, fl_sst
263      !!
264      !! AXY (08/07/11): fate of fast detritus reaching the seafloor
265! I don't think ffast2slowfe is used - marc 10/4/17
266!      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowfe,ffast2slowc
267!      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowc
268      !!
269      !! conservation law
270      REAL(wp) ::    fnit0,fsil0,ffer0 
271# if defined key_roam
272      REAL(wp) ::    fcar0,falk0,foxy0 
273# endif     
274      !!
275      !! temporary variables
276      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8,fq9
277      !!
278      !! water column nutrient and flux integrals
279!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_n,ftot_si,ftot_fe
280!      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_n,fflx_si,fflx_fe
281!      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_n,fifd_si,fifd_fe
282!      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_n,fofd_si,fofd_fe
283# if defined key_roam
284!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_c,ftot_a,ftot_o2
285!      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_c,fflx_a,fflx_o2
286!      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_c,fifd_a,fifd_o2
287!      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_c,fofd_a,fofd_o2
288# endif
289      !!
290      !! zooplankton grazing integrals
291!      REAL(wp), DIMENSION(jpi,jpj) ::    fzmi_i,fzmi_o,fzme_i,fzme_o
292      !!
293      !! limitation term temporary variables
294!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_pn,ftot_pd
295!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_zmi,ftot_zme,ftot_det,ftot_dtc
296      !! use ballast scheme (1) or simple exponential scheme (0; a conservation test)
297      INTEGER  ::    iball
298      !! use biological fluxes (1) or not (0)
299      INTEGER  ::    ibio_switch
300      !!
301      !! diagnose fluxes (should only be used in 1D runs)
302!      INTEGER  ::    idf, idfval
303      !!
304      !! nitrogen and silicon production and consumption
305      REAL(wp) ::    fn_prod, fn_cons, fs_prod, fs_cons
306!      REAL(wp), DIMENSION(jpi,jpj) ::    fnit_prod, fnit_cons, fsil_prod, fsil_cons
307# if defined key_roam
308      !!
309      !! flags to help with calculating the position of the CCD
310! Moved into carb_chem.F90 - marc 20/4/17
311!      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
312      !!
313      !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
314!      REAL(wp)                     ::    f_xco2a
315!      REAL(wp), DIMENSION(jpi,jpj) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
316!      REAL(wp), DIMENSION(jpi,jpj) ::    f_TDIC, f_TALK, f_dcf, f_henry
317!      REAL(wp), DIMENSION(jpi,jpj) ::    f_pp0
318!      REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat
319      REAL(wp)                     ::    f_o2sat3
320!      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg
321      !!
322      !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
323!      REAL(wp), DIMENSION(jpi,jpj) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm
324!      REAL(wp), DIMENSION(jpi,jpj) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2
325      !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
326      REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol
327      REAL, PARAMETER :: secs_in_day    = 86400.0  !! s / d
328      REAL, PARAMETER :: CO2flux_conv   = (1.e-6 * weight_CO2_mol) / secs_in_day
329
330      !!
331!      INTEGER, DIMENSION(jpi,jpj)  ::    iters
332      REAL(wp) ::    f_year
333      INTEGER  ::    i_year
334      INTEGER  ::    iyr1, iyr2
335      !!
336      !! carbon, alkalinity production and consumption
337      REAL(wp) ::    fc_prod, fc_cons, fa_prod, fa_cons
338!      REAL(wp), DIMENSION(jpi,jpj) ::    fcomm_resp
339!      REAL(wp), DIMENSION(jpi,jpj) ::    fcar_prod, fcar_cons
340      !!
341      !! oxygen production and consumption (and non-consumption)
342      REAL(wp), DIMENSION(jpi,jpj) ::    fo2_prod, fo2_cons, fo2_ncons, fo2_ccons
343!      REAL(wp), DIMENSION(jpi,jpj) ::    foxy_prod, foxy_cons, foxy_anox
344      !! Jpalm (11-08-2014)
345      !! add DMS in MEDUSA for UKESM1 model
346!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf
347      !! AXY (13/03/15): add in other DMS calculations
348!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall
349
350# endif
351      !!
352      !! benthic fluxes
353!      INTEGER  ::    ibenthic
354!      REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe,              f_sbenin_c
355!      REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca
356!      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca
357      REAL(wp) ::    zfact
358      !!
359      !! benthic fluxes of CaCO3 that shouldn't happen because of lysocline
360!      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca
361      !!
362      !! riverine fluxes
363!      REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk
364      !! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface
365      REAL(wp), DIMENSION(jpi,jpj) :: f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk
366      !!---------------------------------------------------------------------
367
368# if defined key_debug_medusa
369      IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined'
370      CALL flush(numout)
371# endif 
372
373      !! AXY (20/11/14): alter this to report on first MEDUSA call
374      !! IF( kt == nit000 ) THEN
375      IF( kt == nittrc000 ) THEN
376         IF(lwp) WRITE(numout,*)
377         IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model'
378         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
379    IF(lwp) WRITE(numout,*) ' kt =',kt
380      ENDIF
381
382      !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
383      ibenthic = 1
384
385      !! not sure what this is for; it's not used anywhere; commenting out
386      !! fbodn(:,:) = 0.e0   
387
388      !!
389      IF( ln_diatrc ) THEN
390         !! blank 2D diagnostic array
391         trc2d(:,:,:) = 0.e0
392         !!
393         !! blank 3D diagnostic array
394         trc3d(:,:,:,:) = 0.e0
395      ENDIF
396
397      !!----------------------------------------------------------------------
398      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency
399      !! terms of all biological equations to 0.
400      !!----------------------------------------------------------------------
401      !!
402      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit
403      !!                 the bill for now; another item on the things-to-sort-
404      !!     out-in-the-future list ...
405# if defined key_kill_medusa
406      b0 = 0.
407# else
408      b0 = 1.
409# endif
410      !!----------------------------------------------------------------------
411      !! fast detritus ballast scheme (0 = no; 1 = yes)
412      !! alternative to ballast scheme is same scheme but with no ballast
413      !! protection (not dissimilar to Martin et al., 1987)
414      !!----------------------------------------------------------------------
415      !!
416      iball = 1
417
418      !!----------------------------------------------------------------------
419      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
420      !! these should *only* be used in 1D since they give comprehensive
421      !! output for ecological functions in the model; primarily used in
422      !! debugging
423      !!----------------------------------------------------------------------
424      !!
425      idf    = 0
426      !!
427      !! timer mechanism
428      if (kt/120*120.eq.kt) then
429         idfval = 1
430      else
431         idfval = 0
432      endif
433
434      !!----------------------------------------------------------------------
435      !! Initialise arrays to zero and set up arrays for diagnostics
436      !!----------------------------------------------------------------------
437! tmp - marc
438      write(numout,*) 'bbb13. before call to bio_medusa_init,kt=',kt
439      flush(numout)
440!
441      CALL bio_medusa_init( kt )
442! tmp - marc
443      write(numout,*) 'bbb14. after call to bio_medusa_init,kt=',kt
444      flush(numout)
445!
446       !!
447# if defined key_axy_nancheck
448       DO jn = 1,jptra
449          !! fq0 = MINVAL(trn(:,:,:,jn))
450          !! fq1 = MAXVAL(trn(:,:,:,jn))
451          fq2 = SUM(trn(:,:,:,jn))
452          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', &
453          !! &        kt, jn, fq0, fq1, fq2
454          !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR
455          !!                 and has been replaced here with a specialist routine
456          !! if (fq2 /= fq2 ) then
457          if ( ieee_is_nan( fq2 ) ) then
458             !! there's a NaN here
459             if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:'
460             DO jk = 1,jpk
461                DO jj = 1,jpj
462                   DO ji = 1,jpi
463                      !! AXY (30/01/14): "isnan" problem on HECTOR
464                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then
465                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
466                         if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', &
467                         &        tmask(ji,jj,jk), ji, jj, jk, jn
468                      endif
469                   enddo
470                enddo
471             enddo
472             CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
473          endif
474       ENDDO
475       CALL flush(numout)
476# endif
477
478# if defined key_debug_medusa
479      IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked'
480      CALL flush(numout)
481# endif 
482
483# if defined key_roam
484      !!----------------------------------------------------------------------
485      !! calculate atmospheric pCO2
486      !!----------------------------------------------------------------------
487      !!
488      !! what's atmospheric pCO2 doing? (data start in 1859)
489      iyr1  = nyear - 1859 + 1
490      iyr2  = iyr1 + 1
491      if (iyr1 .le. 1) then
492         !! before 1860
493         f_xco2a(:,:) = hist_pco2(1)
494      elseif (iyr2 .ge. 242) then
495         !! after 2099
496         f_xco2a(:,:) = hist_pco2(242)
497      else
498         !! just right
499         fq0 = hist_pco2(iyr1)
500         fq1 = hist_pco2(iyr2)
501         fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0)
502         !! AXY (14/06/12): tweaked to make more sense (and be correct)
503#  if defined key_bs_axy_yrlen
504         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  !! bugfix: for 360d year with HadGEM2-ES forcing
505#  else
506         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0  !! original use of 365 days (not accounting for leap year or 360d year)
507#  endif
508         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
509         f_xco2a(:,:) = fq4
510      endif
511#  if defined key_axy_pi_co2
512      f_xco2a(:,:) = 284.725          !! OCMIP pre-industrial pCO2
513#  endif
514      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
515      !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day)
516      !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)
517      !! AXY (29/01/14): remove surplus diagnostics
518      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0
519      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1
520      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
521      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
522      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
523# endif
524
525# if defined key_debug_medusa
526      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry'
527      IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt
528      IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000
529      CALL flush(numout)
530# endif 
531
532# if defined key_roam
533      !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every
534      !!                 month (this is hardwired as 960 timesteps but should
535      !!                 be calculated and done properly
536      !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN
537      !! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN
538      !!=============================
539      !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call :
540      !!          we don't want to call on the first time-step of all run submission,
541      !!          but only on the very first time-step, and then every month
542      !!          So we call on nittrc000 if not restarted run,
543      !!          else if one month after last call.
544      !!          assume one month is 30d --> 3600*24*30 : 2592000s
545      !!          try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt)   
546      !!          ++ need to pass carb-chem output var through restarts
547      If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN
548         !!----------------------------------------------------------------------
549         !! Calculate the carbonate chemistry for the whole ocean on the first
550         !! simulation timestep and every month subsequently; the resulting 3D
551         !! field of omega calcite is used to determine the depth of the CCD
552         !!----------------------------------------------------------------------
553         CALL carb_chem( kt )
554
555      ENDIF
556# endif
557
558# if defined key_debug_medusa
559      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations'
560      CALL flush(numout)
561# endif 
562
563      !!----------------------------------------------------------------------
564      !! MEDUSA has unified equation through the water column
565      !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)
566      !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1         
567      !!----------------------------------------------------------------------
568      !!
569      !! NOTE: the ordering of the loops below differs from that of some other
570      !! models; looping over the vertical dimension is the outermost loop and
571      !! this complicates some calculations (e.g. storage of vertical fluxes
572      !! that can otherwise be done via a singular variable require 2D fields
573      !! here); however, these issues are relatively easily resolved, but the
574      !! loops CANNOT be reordered without potentially causing code efficiency
575      !! problems (e.g. array indexing means that reordering the loops would
576      !! require skipping between widely-spaced memory location; potentially
577      !! outside those immediately cached)
578      !!
579      !! OPEN vertical loop
580      DO jk = 1,jpk
581         !! OPEN horizontal loops
582         DO jj = 2,jpjm1
583         DO ji = 2,jpim1
584            !! OPEN wet point IF..THEN loop
585            if (tmask(ji,jj,jk) == 1) then               
586               !!===========================================================
587               !! SETUP LOCAL GRID CELL
588               !!===========================================================
589               !!
590               !!-----------------------------------------------------------
591               !! Some notes on grid vertical structure
592               !! - fsdepw(ji,jj,jk) is the depth of the upper surface of
593               !!   level jk
594               !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of
595               !!   level jk
596               !! - fse3t(ji,jj,jk)  is the thickness of level jk
597               !!-----------------------------------------------------------
598               !!
599               !! AXY (01/03/10): set up level depth (bottom of level)
600               fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
601               !! AXY (28/11/16): local seafloor depth
602               !!                 previously mbathy(ji,jj) - 1, now
603               !!                 mbathy(ji,jj)
604               mbathy(ji,jj) = mbathy(ji,jj)
605               !!
606               !! set up model tracers
607               !! negative values of state variables are not allowed to
608               !! contribute to the calculated fluxes
609               !! non-diatom chlorophyll
610               zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
611               !! diatom chlorophyll
612               zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
613               !! non-diatoms
614               zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
615               !! diatoms
616               zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
617               !! diatom silicon
618               zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
619               !! AXY (28/01/10): probably need to take account of
620               !! chl/biomass connection
621               if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
622               if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
623               if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
624               if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
625          !! AXY (23/01/14): duh - why did I forget diatom silicon?
626          if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
627          if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
628               !! microzooplankton
629               zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
630               !! mesozooplankton
631               zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
632               !! detrital nitrogen
633               zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
634               !! dissolved inorganic nitrogen
635               zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
636               !! dissolved silicic acid
637               zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
638               !! dissolved "iron"
639               zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
640# if defined key_roam
641               !! detrital carbon
642               zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
643               !! dissolved inorganic carbon
644               zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
645               !! alkalinity
646               zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
647               !! oxygen
648               zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
649#  if defined key_axy_carbchem && defined key_mocsy
650               !! phosphate via DIN and Redfield
651               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
652#  endif
653               !!
654               !! also need physical parameters for gas exchange calculations
655               ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
656               zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
657               !!
658          !! AXY (28/02/14): check input fields
659               if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
660                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ',  &
661                     tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',     &
662                     ji, ',', jj, ',', jk, ') at time', kt
663        IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',&
664                     tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
665                  !! temperatur
666                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
667               endif
668               if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
669                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', &
670                     tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    &
671                     ji, ',', jj, ',', jk, ') at time', kt
672               endif
673# else
674               !! implicit detrital carbon
675               zdtc(ji,jj) = zdet(ji,jj) * xthetad
676# endif
677# if defined key_debug_medusa
678               if (idf.eq.1) then
679                  !! AXY (15/01/10)
680                  if (trn(ji,jj,jk,jpdin).lt.0.) then
681                     IF (lwp) write (numout,*) '------------------------------'
682                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',        &
683                        trn(ji,jj,jk,jpdin)
684                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',        &
685                        ji, jj, jk, kt
686                  endif
687                  if (trn(ji,jj,jk,jpsil).lt.0.) then
688                     IF (lwp) write (numout,*) '------------------------------'
689                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',        &
690                        trn(ji,jj,jk,jpsil)
691                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',        &
692                        ji, jj, jk, kt
693                  endif
694#  if defined key_roam
695                  if (trn(ji,jj,jk,jpdic).lt.0.) then
696                     IF (lwp) write (numout,*) '------------------------------'
697                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',        &
698                        trn(ji,jj,jk,jpdic)
699                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',        &
700                        ji, jj, jk, kt
701                  endif
702                  if (trn(ji,jj,jk,jpalk).lt.0.) then
703                     IF (lwp) write (numout,*) '------------------------------'
704                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',        &
705                        trn(ji,jj,jk,jpalk)
706                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',        &
707                        ji, jj, jk, kt
708                  endif
709                  if (trn(ji,jj,jk,jpoxy).lt.0.) then
710                     IF (lwp) write (numout,*) '------------------------------'
711                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',        &
712                        trn(ji,jj,jk,jpoxy)
713                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',        &
714                        ji, jj, jk, kt
715                  endif
716#  endif
717               endif
718# endif
719# if defined key_debug_medusa
720               !! report state variable values
721               if (idf.eq.1.AND.idfval.eq.1) then
722                  IF (lwp) write (numout,*) '------------------------------'
723                  IF (lwp) write (numout,*) 'fthk(',jk,') = ', fse3t(ji,jj,jk)
724                  IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
725                  IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
726                  IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
727                  IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
728                  IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
729                  IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
730                  IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
731                  IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
732                  IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
733#  if defined key_roam
734                  IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
735                  IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
736                  IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
737                  IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)                 
738#  endif
739               endif
740# endif
741
742# if defined key_debug_medusa
743               if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
744                  IF (lwp) write (numout,*) '------------------------------'
745                  IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
746               endif
747# endif
748
749               !! sum tracers for inventory checks
750               IF( lk_iomput ) THEN
751                  IF ( med_diag%INVTN%dgsave )   THEN
752                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        &
753                        (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) +     &
754                                             zzmi(ji,jj) + zzme(ji,jj) +     &
755                                             zdet(ji,jj) + zdin(ji,jj) ) )
756                  ENDIF
757                  IF ( med_diag%INVTSI%dgsave )  THEN
758                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       & 
759                        (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) )
760                  ENDIF
761                  IF ( med_diag%INVTFE%dgsave )  THEN
762                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       & 
763                        (fse3t(ji,jj,jk) * ( xrfn *                          &
764                                            ( zphn(ji,jj) + zphd(ji,jj) +    &
765                                              zzmi(ji,jj) + zzme(ji,jj) +    &
766                                              zdet(ji,jj) ) + zfer(ji,jj) ) )
767                  ENDIF
768# if defined key_roam
769                  IF ( med_diag%INVTC%dgsave )  THEN
770                     ftot_c(ji,jj)  = ftot_c(ji,jj) + & 
771                        (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) +      &
772                                             (xthetapd * zphd(ji,jj)) +      &
773                                             (xthetazmi * zzmi(ji,jj)) +     &
774                                             (xthetazme * zzme(ji,jj)) +     &
775                                             zdtc(ji,jj) + zdic(ji,jj) ) )
776                  ENDIF
777                  IF ( med_diag%INVTALK%dgsave ) THEN
778                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     &
779                                                       zalk(ji,jj))
780                  ENDIF
781                  IF ( med_diag%INVTO2%dgsave )  THEN
782                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    &
783                                                        zoxy(ji,jj))
784                  ENDIF
785                  !!
786                  !! AXY (10/11/16): CMIP6 diagnostics
787                  IF ( med_diag%INTDISSIC%dgsave ) THEN
788                     intdissic(ji,jj) = intdissic(ji,jj) +                   &
789                                        (fse3t(ji,jj,jk) * zdic(ji,jj))
790                  ENDIF
791                  IF ( med_diag%INTDISSIN%dgsave ) THEN
792                     intdissin(ji,jj) = intdissin(ji,jj) +                   &
793                                        (fse3t(ji,jj,jk) * zdin(ji,jj))
794                  ENDIF
795                  IF ( med_diag%INTDISSISI%dgsave ) THEN
796                     intdissisi(ji,jj) = intdissisi(ji,jj) +                 &
797                                         (fse3t(ji,jj,jk) * zsil(ji,jj))
798                  ENDIF
799                  IF ( med_diag%INTTALK%dgsave ) THEN
800                     inttalk(ji,jj) = inttalk(ji,jj) +                       &
801                                      (fse3t(ji,jj,jk) * zalk(ji,jj))
802                  ENDIF
803                  IF ( med_diag%O2min%dgsave ) THEN
804                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then
805                        o2min(ji,jj)  = zoxy(ji,jj)
806                        IF ( med_diag%ZO2min%dgsave ) THEN
807                           !! layer midpoint
808                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               &
809                                            fdep1(ji,jj)) / 2.0
810                        ENDIF
811                     endif
812                  ENDIF
813# endif
814               ENDIF
815
816               CALL flush(numout)
817
818            ENDIF
819         ENDDO
820         ENDDO
821
822         !!----------------------------------------------------------------
823         !! Calculate air-sea gas exchange and river inputs
824         !!----------------------------------------------------------------
825         IF ( jk == 1 ) THEN
826            call air_sea( kt )
827         END IF
828
829# if defined key_roam
830
831! Moved from above - marc 21/4/17
832! I think this should be moved into diagnostics at bottom - it
833! doesn't like it is used anyway else - marc 21/4/17
834         DO jj = 2,jpjm1
835         DO ji = 2,jpim1
836            !! OPEN wet point IF..THEN loop
837            if (tmask(ji,jj,jk) == 1) then
838
839               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
840               IF ( med_diag%O2SAT3%dgsave ) THEN
841                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 )
842                  o2sat3(ji, jj, jk) = f_o2sat3
843               ENDIF
844            ENDIF
845         ENDDO
846         ENDDO
847# endif
848
849         !!------------------------------------------------------------------
850         !! Phytoplankton growth, zooplankton grazing and miscellaneous
851         !! plankton losses.
852         !!------------------------------------------------------------------
853         CALL plankton( jk )
854
855         !!----------------------------------------------------------------
856         !! Iron chemistry and scavenging
857         !!----------------------------------------------------------------
858         CALL iron_chem_scav( jk )
859
860! Miscellaneous processes - marc
861
862         DO jj = 2,jpjm1
863         DO ji = 2,jpim1
864            !! OPEN wet point IF..THEN loop
865            if (tmask(ji,jj,jk) == 1) then
866
867               !!----------------------------------------------------------------------
868               !! Miscellaneous
869               !!----------------------------------------------------------------------
870               !!
871               !! diatom frustule dissolution
872               fsdiss(ji,jj)  = xsdiss * zpds(ji,jj)
873
874# if defined key_debug_medusa
875               !! report miscellaneous calculations
876               if (idf.eq.1.AND.idfval.eq.1) then
877                  IF (lwp) write (numout,*) '------------------------------'
878                  IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss(ji,jj)
879               endif
880# endif
881
882               !!----------------------------------------------------------------------
883               !! Slow detritus creation
884               !!----------------------------------------------------------------------
885               !! this variable integrates the creation of slow sinking detritus
886               !! to allow the split between fast and slow detritus to be
887               !! diagnosed
888               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj)) + &
889               ((1.0 - xfdfrac2) * fdzme(ji,jj)) + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))
890               !!
891               !! this variable records the slow detrital sinking flux at this
892               !! particular depth; it is used in the output of this flux at
893               !! standard depths in the diagnostic outputs; needs to be
894               !! adjusted from per second to per day because of parameter vsed
895               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400.
896# if defined key_roam
897               !!
898               !! and the same for detrital carbon
899               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) + (xthetazmi * fdzmi(ji,jj)) + &
900               (xthetapd * (1.0 - xfdfrac1) * fdpd(ji,jj)) + &
901               (xthetazme * (1.0 - xfdfrac2) * fdzme(ji,jj)) + &
902               ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))
903               !!
904               !! this variable records the slow detrital sinking flux at this
905               !! particular depth; it is used in the output of this flux at
906               !! standard depths in the diagnostic outputs; needs to be
907               !! adjusted from per second to per day because of parameter vsed
908               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400.
909# endif
910
911               !!----------------------------------------------------------------------
912               !! Nutrient regeneration
913               !! this variable integrates total nitrogen regeneration down the
914               !! watercolumn; its value is stored and output as a 2D diagnostic;
915               !! the corresponding dissolution flux of silicon (from sources
916               !! other than fast detritus) is also integrated; note that,
917               !! confusingly, the linear loss terms from plankton compartments
918               !! are labelled as fdX2 when one might have expected fdX or fdX1
919               !!----------------------------------------------------------------------
920               !!
921               !! nitrogen
922               fregen(ji,jj)   = (( (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) +                &  ! messy feeding
923               (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) +           &  ! messy feeding
924               fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) +                                &  ! excretion + D remin.
925               fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)) * fse3t(ji,jj,jk))                    ! linear mortality
926               !!
927               !! silicon
928               fregensi(ji,jj) = (( fsdiss(ji,jj) + ((1.0 - xfdfrac1) * fdpds(ji,jj)) +      &  ! dissolution + non-lin. mortality
929               ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +                           &  ! egestion by zooplankton
930               fdpds2(ji,jj)) * fse3t(ji,jj,jk))                                             ! linear mortality
931# if defined key_roam
932               !!
933               !! carbon
934! Doesn't look this is used - marc 10/4/17
935!               fregenc(ji,jj)  = (( (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))) +  &  ! messy feeding
936!               (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +     &  ! messy feeding
937!               (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))) +                       &  ! messy feeding
938!               fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) +                               &  ! respiration + D remin.
939!               (xthetapn * fdpn2(ji,jj)) + (xthetapd * fdpd2(ji,jj)) +                &  ! linear mortality
940!               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality
941# endif
942
943! MAYBE BUT A BREAK IN HERE, FAST-SINKINIG DETRITUS - marc 20/4/17
944! (miscellaneous processes is 79 lines)
945            ENDIF
946         ENDDO
947         ENDDO
948
949         !!------------------------------------------------------------------
950         !! Detritus process
951         !!------------------------------------------------------------------
952         CALL detritus( jk, iball )
953
954! Updating coming next - marc 28/4/17
955
956         DO jj = 2,jpjm1
957         DO ji = 2,jpim1
958            !! OPEN wet point IF..THEN loop
959            if (tmask(ji,jj,jk) == 1) then
960
961               !!======================================================================
962               !! LOCAL GRID CELL TRENDS
963               !!======================================================================
964               !!
965               !!----------------------------------------------------------------------
966               !! Determination of trends
967               !!----------------------------------------------------------------------
968               !!
969               !!----------------------------------------------------------------------
970               !! chlorophyll
971               btra(ji,jj,jpchn) = b0 * ( &
972                 + ((frn(ji,jj) * fprn(ji,jj) * zphn(ji,jj)) - fgmipn(ji,jj) - fgmepn(ji,jj) - fdpn(ji,jj) - fdpn2(ji,jj)) * (fthetan(ji,jj) / xxi) )
973               btra(ji,jj,jpchd) = b0 * ( &
974                 + ((frd(ji,jj) * fprd(ji,jj) * zphd(ji,jj)) - fgmepd(ji,jj) - fdpd(ji,jj) - fdpd2(ji,jj)) * (fthetad(ji,jj) / xxi) )
975               !!
976               !!----------------------------------------------------------------------
977               !! phytoplankton
978               btra(ji,jj,jpphn) = b0 * ( &
979                 + (fprn(ji,jj) * zphn(ji,jj)) - fgmipn(ji,jj) - fgmepn(ji,jj) - fdpn(ji,jj) - fdpn2(ji,jj) )
980               btra(ji,jj,jpphd) = b0 * ( &
981                 + (fprd(ji,jj) * zphd(ji,jj)) - fgmepd(ji,jj) - fdpd(ji,jj) - fdpd2(ji,jj) )
982               btra(ji,jj,jppds) = b0 * ( &
983                 + (fprds(ji,jj) * zpds(ji,jj)) - fgmepds(ji,jj) - fdpds(ji,jj) - fsdiss(ji,jj) - fdpds2(ji,jj) )
984               !!
985               !!----------------------------------------------------------------------
986               !! zooplankton
987               btra(ji,jj,jpzmi) = b0 * ( &
988                 + fmigrow(ji,jj) - fgmezmi(ji,jj) - fdzmi(ji,jj) - fdzmi2(ji,jj) )
989               btra(ji,jj,jpzme) = b0 * ( &
990                 + fmegrow(ji,jj) - fdzme(ji,jj) - fdzme2(ji,jj) )
991               !!
992               !!----------------------------------------------------------------------
993               !! detritus
994               btra(ji,jj,jpdet) = b0 * ( &
995                 + fdpn(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj))              &  ! mort. losses
996                 + fdzmi(ji,jj) + ((1.0 - xfdfrac2) * fdzme(ji,jj))            &  ! mort. losses
997                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))            &  ! assim. inefficiency
998                 - fgmid(ji,jj) - fgmed(ji,jj) - fdd(ji,jj)                           &  ! grazing and remin.
999                 + ffast2slown(ji,jj) )                                    ! seafloor fast->slow
1000               !!
1001               !!----------------------------------------------------------------------
1002               !! dissolved inorganic nitrogen nutrient
1003               fn_cons = 0.0  &
1004                 - (fprn(ji,jj) * zphn(ji,jj)) - (fprd(ji,jj) * zphd(ji,jj))                    ! primary production
1005               fn_prod = 0.0  &
1006                 + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)))                     &  ! messy feeding remin.
1007                 + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)))  &  ! messy feeding remin.
1008                 + fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) + freminn(ji,jj)             &  ! excretion and remin.
1009                 + fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)                  ! metab. losses
1010               !!
1011               !! riverine flux
1012               if ( jriver_n .gt. 0 ) then
1013                  f_riv_loc_n(ji,jj) = f_riv_n(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1014                  fn_prod = fn_prod + f_riv_loc_n(ji,jj)
1015               endif
1016               !! 
1017               !! benthic remineralisation
1018               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1019                  fn_prod = fn_prod + (f_benout_n(ji,jj) / fse3t(ji,jj,jk))
1020               endif
1021               !!
1022               btra(ji,jj,jpdin) = b0 * ( &
1023                 fn_prod + fn_cons )
1024               !!
1025               fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved nitrogen
1026                 fn_cons ) )
1027               fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved nitrogen
1028                 fn_prod ) )
1029               !!
1030               !!----------------------------------------------------------------------
1031               !! dissolved silicic acid nutrient
1032               fs_cons = 0.0  &
1033                 - (fprds(ji,jj) * zpds(ji,jj))                                   ! opal production
1034               fs_prod = 0.0  &
1035                 + fsdiss(ji,jj)                                        &  ! opal dissolution
1036                 + ((1.0 - xfdfrac1) * fdpds(ji,jj))                    &  ! mort. loss
1037                 + ((1.0 - xfdfrac3) * fgmepds(ji,jj))                  &  ! egestion of grazed Si
1038                 + freminsi(ji,jj) + fdpds2(ji,jj)                                ! fast diss. and metab. losses
1039               !!
1040               !! riverine flux
1041               if ( jriver_si .gt. 0 ) then
1042                  f_riv_loc_si(ji,jj) = f_riv_si(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1043                  fs_prod = fs_prod + f_riv_loc_si(ji,jj)
1044               endif
1045               !! 
1046               !! benthic remineralisation
1047               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
1048                  fs_prod = fs_prod + (f_benout_si(ji,jj) / fse3t(ji,jj,jk))
1049               endif
1050               !!
1051               btra(ji,jj,jpsil) = b0 * ( &
1052                 fs_prod + fs_cons )
1053               !!
1054               fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved silicon
1055                 fs_cons ) )
1056               fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved silicon
1057                 fs_prod ) )
1058               !!
1059               !!----------------------------------------------------------------------
1060               !! dissolved "iron" nutrient
1061               btra(ji,jj,jpfer) = b0 * ( &
1062               + (xrfn * btra(ji,jj,jpdin)) + ffetop(ji,jj) + ffebot(ji,jj) - ffescav(ji,jj) )
1063
1064# if defined key_roam
1065               !!
1066               !!----------------------------------------------------------------------
1067               !! AXY (26/11/08): implicit detrital carbon change
1068               btra(ji,jj,jpdtc) = b0 * ( &
1069                 + (xthetapn * fdpn(ji,jj)) + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))      &  ! mort. losses
1070                 + (xthetazmi * fdzmi(ji,jj)) + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))  &  ! mort. losses
1071                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))                              &  ! assim. inefficiency
1072                 - fgmidc(ji,jj) - fgmedc(ji,jj) - fddc(ji,jj)                                          &  ! grazing and remin.
1073                 + ffast2slowc(ji,jj) )                                                      ! seafloor fast->slow
1074               !!
1075               !!----------------------------------------------------------------------
1076               !! dissolved inorganic carbon
1077               fc_cons = 0.0  &
1078                 - (xthetapn * fprn(ji,jj) * zphn(ji,jj)) - (xthetapd * fprd(ji,jj) * zphd(ji,jj))                ! primary production
1079               fc_prod = 0.0  &
1080                 + (xthetapn * xphi * fgmipn(ji,jj)) + (xphi * fgmidc(ji,jj))                    &  ! messy feeding remin
1081                 + (xthetapn * xphi * fgmepn(ji,jj)) + (xthetapd * xphi * fgmepd(ji,jj))         &  ! messy feeding remin
1082                 + (xthetazmi * xphi * fgmezmi(ji,jj)) + (xphi * fgmedc(ji,jj))                  &  ! messy feeding remin
1083                 + fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) + freminc(ji,jj) + (xthetapn * fdpn2(ji,jj))         &  ! resp., remin., losses
1084                 + (xthetapd * fdpd2(ji,jj)) + (xthetazmi * fdzmi2(ji,jj))                       &  ! losses
1085                 + (xthetazme * fdzme2(ji,jj))                                               ! losses
1086               !!
1087               !! riverine flux
1088               if ( jriver_c .gt. 0 ) then
1089                  f_riv_loc_c(ji,jj) = f_riv_c(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1090                  fc_prod = fc_prod + f_riv_loc_c(ji,jj)
1091               endif
1092               !! 
1093               !! benthic remineralisation
1094               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1095                  fc_prod = fc_prod + (f_benout_c(ji,jj) / fse3t(ji,jj,jk))
1096               endif
1097               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
1098                  fc_prod = fc_prod + (f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
1099               endif
1100               !!
1101               !! community respiration (does not include CaCO3 terms - obviously!)
1102               fcomm_resp(ji,jj) = fcomm_resp(ji,jj) + fc_prod
1103               !!
1104               !! CaCO3
1105               fc_prod = fc_prod - ftempca(ji,jj) + freminca(ji,jj)
1106               !!
1107               !! riverine flux
1108               if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then
1109                  fc_prod = fc_prod + f_riv_c(ji,jj)
1110               endif
1111               !!
1112               btra(ji,jj,jpdic) = b0 * ( &
1113                 fc_prod + fc_cons )
1114               !!
1115               fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved carbon
1116                 fc_cons ) )
1117               fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved carbon
1118                 fc_prod ) )
1119               !!
1120               !!----------------------------------------------------------------------
1121               !! alkalinity
1122               fa_prod = 0.0  &
1123                 + (2.0 * freminca(ji,jj))                                                   ! CaCO3 dissolution
1124               fa_cons = 0.0  &
1125                 - (2.0 * ftempca(ji,jj))                                                    ! CaCO3 production
1126               !!
1127               !! riverine flux
1128               if ( jriver_alk .gt. 0 ) then
1129                  f_riv_loc_alk(ji,jj) = f_riv_alk(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1130                  fa_prod = fa_prod + f_riv_loc_alk(ji,jj)
1131               endif
1132               !! 
1133               !! benthic remineralisation
1134               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
1135                  fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
1136               endif
1137               !!
1138               btra(ji,jj,jpalk) = b0 * ( &
1139                 fa_prod + fa_cons )
1140               !!
1141               !!----------------------------------------------------------------------
1142               !! oxygen (has protection at low O2 concentrations; OCMIP-2 style)
1143               fo2_prod(ji,jj) = 0.0 &
1144                 + (xthetanit * fprn(ji,jj) * zphn(ji,jj))                                      & ! Pn primary production, N
1145                 + (xthetanit * fprd(ji,jj) * zphd(ji,jj))                                      & ! Pd primary production, N
1146                 + (xthetarem * xthetapn * fprn(ji,jj) * zphn(ji,jj))                           & ! Pn primary production, C
1147                 + (xthetarem * xthetapd * fprd(ji,jj) * zphd(ji,jj))                             ! Pd primary production, C
1148               fo2_ncons(ji,jj) = 0.0 &
1149                 - (xthetanit * xphi * fgmipn(ji,jj))                                    & ! Pn messy feeding remin., N
1150                 - (xthetanit * xphi * fgmid(ji,jj))                                     & ! D  messy feeding remin., N
1151                 - (xthetanit * xphi * fgmepn(ji,jj))                                    & ! Pn messy feeding remin., N
1152                 - (xthetanit * xphi * fgmepd(ji,jj))                                    & ! Pd messy feeding remin., N
1153                 - (xthetanit * xphi * fgmezmi(ji,jj))                                   & ! Zi messy feeding remin., N
1154                 - (xthetanit * xphi * fgmed(ji,jj))                                     & ! D  messy feeding remin., N
1155                 - (xthetanit * fmiexcr(ji,jj))                                          & ! microzoo excretion, N
1156                 - (xthetanit * fmeexcr(ji,jj))                                          & ! mesozoo  excretion, N
1157                 - (xthetanit * fdd(ji,jj))                                              & ! slow detritus remin., N
1158                 - (xthetanit * freminn(ji,jj))                                          & ! fast detritus remin., N
1159                 - (xthetanit * fdpn2(ji,jj))                                            & ! Pn  losses, N
1160                 - (xthetanit * fdpd2(ji,jj))                                            & ! Pd  losses, N
1161                 - (xthetanit * fdzmi2(ji,jj))                                           & ! Zmi losses, N
1162                 - (xthetanit * fdzme2(ji,jj))                                             ! Zme losses, N
1163               !! 
1164               !! benthic remineralisation
1165               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1166                  fo2_ncons(ji,jj) = fo2_ncons(ji,jj) - (xthetanit * f_benout_n(ji,jj) / fse3t(ji,jj,jk))
1167               endif
1168               fo2_ccons(ji,jj) = 0.0 &
1169                 - (xthetarem * xthetapn * xphi * fgmipn(ji,jj))                         & ! Pn messy feeding remin., C
1170                 - (xthetarem * xphi * fgmidc(ji,jj))                                    & ! D  messy feeding remin., C
1171                 - (xthetarem * xthetapn * xphi * fgmepn(ji,jj))                         & ! Pn messy feeding remin., C
1172                 - (xthetarem * xthetapd * xphi * fgmepd(ji,jj))                         & ! Pd messy feeding remin., C
1173                 - (xthetarem * xthetazmi * xphi * fgmezmi(ji,jj))                       & ! Zi messy feeding remin., C
1174                 - (xthetarem * xphi * fgmedc(ji,jj))                                    & ! D  messy feeding remin., C
1175                 - (xthetarem * fmiresp(ji,jj))                                          & ! microzoo respiration, C
1176                 - (xthetarem * fmeresp(ji,jj))                                          & ! mesozoo  respiration, C
1177                 - (xthetarem * fddc(ji,jj))                                             & ! slow detritus remin., C
1178                 - (xthetarem * freminc(ji,jj))                                          & ! fast detritus remin., C
1179                 - (xthetarem * xthetapn * fdpn2(ji,jj))                                 & ! Pn  losses, C
1180                 - (xthetarem * xthetapd * fdpd2(ji,jj))                                 & ! Pd  losses, C
1181                 - (xthetarem * xthetazmi * fdzmi2(ji,jj))                               & ! Zmi losses, C
1182                 - (xthetarem * xthetazme * fdzme2(ji,jj))                                 ! Zme losses, C
1183               !! 
1184               !! benthic remineralisation
1185               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1186                  fo2_ccons(ji,jj) = fo2_ccons(ji,jj) - (xthetarem * f_benout_c(ji,jj) / fse3t(ji,jj,jk))
1187               endif
1188               fo2_cons(ji,jj) = fo2_ncons(ji,jj) + fo2_ccons(ji,jj)
1189               !!
1190               !! is this a suboxic zone?
1191               if (zoxy(ji,jj).lt.xo2min) then  ! deficient O2; production fluxes only
1192                  btra(ji,jj,jpoxy) = b0 * ( &
1193                    fo2_prod(ji,jj) )
1194                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
1195                  foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
1196               else                      ! sufficient O2; production + consumption fluxes
1197                  btra(ji,jj,jpoxy) = b0 * ( &
1198                    fo2_prod(ji,jj) + fo2_cons(ji,jj) )
1199                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
1200                  foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
1201               endif
1202               !!
1203               !! air-sea fluxes (if this is the surface box)
1204               if (jk.eq.1) then
1205                  !!
1206                  !! CO2 flux
1207                  btra(ji,jj,jpdic) = btra(ji,jj,jpdic) + (b0 * f_co2flux(ji,jj))
1208                  !!
1209                  !! O2 flux (mol/m3/s -> mmol/m3/d)
1210                  btra(ji,jj,jpoxy) = btra(ji,jj,jpoxy) + (b0 * f_o2flux(ji,jj))
1211               endif
1212# endif
1213
1214# if defined key_debug_medusa
1215               !! report state variable fluxes (not including fast-sinking detritus)
1216               if (idf.eq.1.AND.idfval.eq.1) then
1217                  IF (lwp) write (numout,*) '------------------------------'
1218                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchn)(',jk,')  = ', btra(ji,jj,jpchn)
1219                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchd)(',jk,')  = ', btra(ji,jj,jpchd)
1220                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphn)(',jk,')  = ', btra(ji,jj,jpphn)
1221                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphd)(',jk,')  = ', btra(ji,jj,jpphd)
1222                  IF (lwp) write (numout,*) 'btra(ji,jj,jppds)(',jk,')  = ', btra(ji,jj,jppds)
1223                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzmi)(',jk,')  = ', btra(ji,jj,jpzmi)
1224                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzme)(',jk,')  = ', btra(ji,jj,jpzme)
1225                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdet)(',jk,')  = ', btra(ji,jj,jpdet)
1226                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdin)(',jk,')  = ', btra(ji,jj,jpdin)
1227                  IF (lwp) write (numout,*) 'btra(ji,jj,jpsil)(',jk,')  = ', btra(ji,jj,jpsil)
1228                  IF (lwp) write (numout,*) 'btra(ji,jj,jpfer)(',jk,')  = ', btra(ji,jj,jpfer)
1229#  if defined key_roam
1230                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdtc)(',jk,')  = ', btra(ji,jj,jpdtc)
1231                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdic)(',jk,')  = ', btra(ji,jj,jpdic)
1232                  IF (lwp) write (numout,*) 'btra(ji,jj,jpalk)(',jk,')  = ', btra(ji,jj,jpalk)
1233                  IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy)(',jk,')  = ', btra(ji,jj,jpoxy)
1234#  endif
1235               endif
1236# endif
1237
1238               !!----------------------------------------------------------------------
1239               !! Integrate calculated fluxes for mass balance
1240               !!----------------------------------------------------------------------
1241               !!
1242               !! === nitrogen ===
1243               fflx_n(ji,jj)  = fflx_n(ji,jj)  + &
1244                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpphn) + btra(ji,jj,jpphd) + btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) + btra(ji,jj,jpdet) + btra(ji,jj,jpdin) )
1245               !! === silicon ===
1246               fflx_si(ji,jj) = fflx_si(ji,jj) + &
1247                  fse3t(ji,jj,jk) * ( btra(ji,jj,jppds) + btra(ji,jj,jpsil) )
1248               !! === iron ===
1249               fflx_fe(ji,jj) = fflx_fe(ji,jj) + &
1250                  fse3t(ji,jj,jk) * ( ( xrfn * ( btra(ji,jj,jpphn) + btra(ji,jj,jpphd) + btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) + btra(ji,jj,jpdet)) ) + btra(ji,jj,jpfer) )
1251# if defined key_roam
1252               !! === carbon ===
1253               fflx_c(ji,jj)  = fflx_c(ji,jj)  + &
1254                  fse3t(ji,jj,jk) * ( (xthetapn * btra(ji,jj,jpphn)) + (xthetapd * btra(ji,jj,jpphd)) + &
1255                  (xthetazmi * btra(ji,jj,jpzmi)) + (xthetazme * btra(ji,jj,jpzme)) + btra(ji,jj,jpdtc) + btra(ji,jj,jpdic) )
1256               !! === alkalinity ===
1257               fflx_a(ji,jj)  = fflx_a(ji,jj)  + &
1258                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpalk) )
1259               !! === oxygen ===
1260               fflx_o2(ji,jj) = fflx_o2(ji,jj) + &
1261                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpoxy) )
1262# endif
1263
1264               !!----------------------------------------------------------------------
1265               !! Apply calculated tracer fluxes
1266               !!----------------------------------------------------------------------
1267               !!
1268               !! units: [unit of tracer] per second (fluxes are calculated above per day)
1269               !!
1270               ibio_switch = 1
1271# if defined key_gulf_finland
1272               !! AXY (17/05/13): fudge in a Gulf of Finland correction; uses longitude-
1273               !!                 latitude range to establish if this is a Gulf of Finland
1274               !!                 grid cell; if so, then BGC fluxes are ignored (though
1275               !!                 still calculated); for reference, this is meant to be a
1276               !!                 temporary fix to see if all of my problems can be done
1277               !!                 away with if I switch off BGC fluxes in the Gulf of
1278               !!                 Finland, which currently appears the source of trouble
1279               if ( glamt(ji,jj).gt.24.7 .and. glamt(ji,jj).lt.27.8 .and. &
1280                  &   gphit(ji,jj).gt.59.2 .and. gphit(ji,jj).lt.60.2 ) then
1281                  ibio_switch = 0
1282               endif
1283# endif               
1284               if (ibio_switch.eq.1) then
1285                  tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(ji,jj,jpchn) / 86400.)
1286                  tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(ji,jj,jpchd) / 86400.)
1287                  tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(ji,jj,jpphn) / 86400.)
1288                  tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(ji,jj,jpphd) / 86400.)
1289                  tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(ji,jj,jppds) / 86400.)
1290                  tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(ji,jj,jpzmi) / 86400.)
1291                  tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(ji,jj,jpzme) / 86400.)
1292                  tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(ji,jj,jpdet) / 86400.)
1293                  tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(ji,jj,jpdin) / 86400.)
1294                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(ji,jj,jpsil) / 86400.)
1295                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(ji,jj,jpfer) / 86400.)
1296# if defined key_roam
1297                  tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(ji,jj,jpdtc) / 86400.)
1298                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(ji,jj,jpdic) / 86400.)
1299                  tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(ji,jj,jpalk) / 86400.)
1300                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(ji,jj,jpoxy) / 86400.)
1301# endif
1302               endif               
1303
1304               !! AXY (18/11/16): CMIP6 diagnostics
1305               IF( med_diag%FBDDTALK%dgsave )  THEN
1306                  fbddtalk(ji,jj)  =  fbddtalk(ji,jj)  + (btra(ji,jj,jpalk) * fse3t(ji,jj,jk))
1307               ENDIF
1308               IF( med_diag%FBDDTDIC%dgsave )  THEN
1309                  fbddtdic(ji,jj)  =  fbddtdic(ji,jj)  + (btra(ji,jj,jpdic) * fse3t(ji,jj,jk))
1310               ENDIF
1311               IF( med_diag%FBDDTDIFE%dgsave ) THEN
1312                  fbddtdife(ji,jj) =  fbddtdife(ji,jj) + (btra(ji,jj,jpfer) * fse3t(ji,jj,jk))
1313               ENDIF
1314               IF( med_diag%FBDDTDIN%dgsave )  THEN
1315                  fbddtdin(ji,jj)  =  fbddtdin(ji,jj)  + (btra(ji,jj,jpdin) * fse3t(ji,jj,jk))
1316               ENDIF
1317               IF( med_diag%FBDDTDISI%dgsave ) THEN
1318                  fbddtdisi(ji,jj) =  fbddtdisi(ji,jj) + (btra(ji,jj,jpsil) * fse3t(ji,jj,jk))
1319               ENDIF
1320          !!
1321               IF( med_diag%BDDTALK3%dgsave )  THEN
1322                  bddtalk3(ji,jj,jk)  =  btra(ji,jj,jpalk)
1323               ENDIF
1324               IF( med_diag%BDDTDIC3%dgsave )  THEN
1325                  bddtdic3(ji,jj,jk)  =  btra(ji,jj,jpdic)
1326               ENDIF
1327               IF( med_diag%BDDTDIFE3%dgsave ) THEN
1328                  bddtdife3(ji,jj,jk) =  btra(ji,jj,jpfer)
1329               ENDIF
1330               IF( med_diag%BDDTDIN3%dgsave )  THEN
1331                  bddtdin3(ji,jj,jk)  =  btra(ji,jj,jpdin)
1332               ENDIF
1333               IF( med_diag%BDDTDISI3%dgsave ) THEN
1334                  bddtdisi3(ji,jj,jk) =  btra(ji,jj,jpsil)
1335               ENDIF
1336
1337#   if defined key_debug_medusa
1338               IF (lwp) write (numout,*) '------'
1339               IF (lwp) write (numout,*) 'trc_bio_medusa: end all calculations'
1340               IF (lwp) write (numout,*) 'trc_bio_medusa: now outputs'
1341                     CALL flush(numout)
1342#   endif
1343
1344# if defined key_axy_nancheck
1345               !!----------------------------------------------------------------------
1346               !! Check calculated tracer fluxes
1347               !!----------------------------------------------------------------------
1348               !!
1349               DO jn = 1,jptra
1350                  fq0 = btra(ji,jj,jn)
1351                  !! AXY (30/01/14): "isnan" problem on HECTOR
1352                  !! if (fq0 /= fq0 ) then
1353                  if ( ieee_is_nan( fq0 ) ) then
1354                     !! there's a NaN here
1355                     if (lwp) write(numout,*) 'NAN detected in btra(ji,jj,', ji, ',', &
1356                     & jj, ',', jk, ',', jn, ') at time', kt
1357           CALL ctl_stop( 'trcbio_medusa, NAN in btra field' )
1358                  endif
1359               ENDDO
1360               DO jn = 1,jptra
1361                  fq0 = tra(ji,jj,jk,jn)
1362                  !! AXY (30/01/14): "isnan" problem on HECTOR
1363                  !! if (fq0 /= fq0 ) then
1364                  if ( ieee_is_nan( fq0 ) ) then
1365                     !! there's a NaN here
1366                     if (lwp) write(numout,*) 'NAN detected in tra(', ji, ',', &
1367                     & jj, ',', jk, ',', jn, ') at time', kt
1368              CALL ctl_stop( 'trcbio_medusa, NAN in tra field' )
1369                  endif
1370               ENDDO
1371               CALL flush(numout)
1372# endif
1373
1374               !!----------------------------------------------------------------------
1375               !! Check model conservation
1376               !! these terms merely sum up the tendency terms of the relevant
1377               !! state variables, which should sum to zero; the iron cycle is
1378               !! complicated by fluxes that add (aeolian deposition and seafloor
1379               !! remineralisation) and remove (scavenging) dissolved iron from
1380               !! the model (i.e. the sum of iron fluxes is unlikely to be zero)
1381               !!----------------------------------------------------------------------
1382               !!
1383               !! fnit0 = btra(ji,jj,jpphn) + btra(ji,jj,jpphd) + btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) + btra(ji,jj,jpdet) + btra(ji,jj,jpdin)  ! + ftempn(ji,jj)
1384               !! fsil0 = btra(ji,jj,jppds) + btra(ji,jj,jpsil)                              ! + ftempsi(ji,jj)
1385               !! ffer0 = (xrfn * fnit0) + btra(ji,jj,jpfer)
1386# if defined key_roam
1387               !! fcar0 = 0.
1388               !! falk0 = 0.
1389               !! foxy0 = 0.
1390# endif
1391               !!
1392               !! if (kt/240*240.eq.kt) then
1393               !!    if (ji.eq.2.and.jj.eq.2.and.jk.eq.1) then
1394               !!       IF (lwp) write (*,*) '*******!MEDUSA Conservation!*******',kt
1395# if defined key_roam
1396               !!       IF (lwp) write (*,*) fnit0,fsil0,ffer0,fcar0,falk0,foxy0
1397# else
1398               !!       IF (lwp) write (*,*) fnit0,fsil0,ffer0
1399# endif
1400               !!    endif
1401               !! endif     
1402
1403! MAYBE BUT A BREAK IN HERE - marc 20/4/17
1404! (this would make the previous section about 470 lines and the one below
1405! about 700 lines)
1406            ENDIF
1407         ENDDO
1408         ENDDO
1409
1410         DO jj = 2,jpjm1
1411         DO ji = 2,jpim1
1412            !! OPEN wet point IF..THEN loop
1413            if (tmask(ji,jj,jk) == 1) then
1414
1415# if defined key_trc_diabio
1416               !!======================================================================
1417               !! LOCAL GRID CELL DIAGNOSTICS
1418               !!======================================================================
1419               !!
1420               !!----------------------------------------------------------------------
1421               !! Full diagnostics key_trc_diabio:
1422               !! LOBSTER and PISCES support full diagnistics option key_trc_diabio   
1423               !! which gives an option of FULL output of biological sourses and sinks.
1424               !! I cannot see any reason for doing this. If needed, it can be done
1425               !! as shown below.
1426               !!----------------------------------------------------------------------
1427               !!
1428               IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio'
1429               !!               trbio(ji,jj,jk, 1) = fprn(ji,jj)
1430# endif
1431
1432               IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
1433         !!----------------------------------------------------------------------
1434         !! Add in XML diagnostics stuff
1435         !!----------------------------------------------------------------------
1436         !!
1437         !! ** 2D diagnostics
1438#   if defined key_debug_medusa
1439                  IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk loop'
1440                  CALL flush(numout)
1441#   endif
1442                  IF ( med_diag%PRN%dgsave ) THEN
1443                      fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
1444                  ENDIF
1445                  IF ( med_diag%MPN%dgsave ) THEN
1446                      fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn(ji,jj)         * fse3t(ji,jj,jk))
1447                  ENDIF
1448                  IF ( med_diag%PRD%dgsave ) THEN
1449                      fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
1450                  ENDIF
1451                  IF( med_diag%MPD%dgsave ) THEN
1452                      fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd(ji,jj)         * fse3t(ji,jj,jk)) 
1453                  ENDIF
1454                  !  IF( med_diag%DSED%dgsave ) THEN
1455                  !      CALL iom_put( "DSED"  , ftot_n )
1456                  !  ENDIF
1457                  IF( med_diag%OPAL%dgsave ) THEN
1458                      fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds(ji,jj) * zpds(ji,jj) * fse3t(ji,jj,jk)) 
1459                  ENDIF
1460                  IF( med_diag%OPALDISS%dgsave ) THEN
1461                      fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss(ji,jj)  * fse3t(ji,jj,jk)) 
1462                  ENDIF
1463                  IF( med_diag%GMIPn%dgsave ) THEN
1464                      fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn(ji,jj)  * fse3t(ji,jj,jk)) 
1465                  ENDIF
1466                  IF( med_diag%GMID%dgsave ) THEN
1467                      fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid(ji,jj)   * fse3t(ji,jj,jk)) 
1468                  ENDIF
1469                  IF( med_diag%MZMI%dgsave ) THEN
1470                      fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi(ji,jj)   * fse3t(ji,jj,jk)) 
1471                  ENDIF
1472                  IF( med_diag%GMEPN%dgsave ) THEN
1473                      fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn(ji,jj)  * fse3t(ji,jj,jk))
1474                  ENDIF
1475                  IF( med_diag%GMEPD%dgsave ) THEN
1476                      fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd(ji,jj)  * fse3t(ji,jj,jk)) 
1477                  ENDIF
1478                  IF( med_diag%GMEZMI%dgsave ) THEN
1479                      fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi(ji,jj) * fse3t(ji,jj,jk)) 
1480                  ENDIF
1481                  IF( med_diag%GMED%dgsave ) THEN
1482                      fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed(ji,jj)   * fse3t(ji,jj,jk)) 
1483                  ENDIF
1484                  IF( med_diag%MZME%dgsave ) THEN
1485                      fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme(ji,jj)   * fse3t(ji,jj,jk)) 
1486                  ENDIF
1487                  !  IF( med_diag%DEXP%dgsave ) THEN
1488                  !      CALL iom_put( "DEXP"  , ftot_n )
1489                  !  ENDIF
1490                  IF( med_diag%DETN%dgsave ) THEN
1491                      fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown(ji,jj)  * fse3t(ji,jj,jk)) 
1492                  ENDIF
1493                  IF( med_diag%MDET%dgsave ) THEN
1494                      fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd(ji,jj)     * fse3t(ji,jj,jk)) 
1495                  ENDIF
1496                  IF( med_diag%AEOLIAN%dgsave ) THEN
1497                      ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop(ji,jj)  * fse3t(ji,jj,jk)) 
1498                  ENDIF
1499                  IF( med_diag%BENTHIC%dgsave ) THEN
1500                      ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot(ji,jj)  * fse3t(ji,jj,jk)) 
1501                  ENDIF
1502                  IF( med_diag%SCAVENGE%dgsave ) THEN
1503                      ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav(ji,jj) * fse3t(ji,jj,jk)) 
1504                  ENDIF
1505                  IF( med_diag%PN_JLIM%dgsave ) THEN
1506                      ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))
1507                      fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) 
1508                  ENDIF
1509                  IF( med_diag%PN_NLIM%dgsave ) THEN
1510                      fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
1511                  ENDIF
1512                  IF( med_diag%PN_FELIM%dgsave ) THEN
1513                      ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln2(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
1514                  ENDIF
1515                  IF( med_diag%PD_JLIM%dgsave ) THEN
1516                      ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
1517                      fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) 
1518                  ENDIF
1519                  IF( med_diag%PD_NLIM%dgsave ) THEN
1520                      fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
1521                  ENDIF
1522                  IF( med_diag%PD_FELIM%dgsave ) THEN
1523                      ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
1524                  ENDIF
1525                  IF( med_diag%PD_SILIM%dgsave ) THEN
1526                      fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) 
1527                  ENDIF
1528                  IF( med_diag%PDSILIM2%dgsave ) THEN
1529                      fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
1530                  ENDIF
1531                  !!
1532                  IF( med_diag%TOTREG_N%dgsave ) THEN
1533                      fregen2d(ji,jj) = fregen2d(ji,jj) + fregen(ji,jj)
1534                  ENDIF
1535                  IF( med_diag%TOTRG_SI%dgsave ) THEN
1536                      fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi(ji,jj)
1537                  ENDIF
1538                  !!
1539                  IF( med_diag%FASTN%dgsave ) THEN
1540                      ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
1541                  ENDIF
1542                  IF( med_diag%FASTSI%dgsave ) THEN
1543                      ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
1544                  ENDIF
1545                  IF( med_diag%FASTFE%dgsave ) THEN
1546                      ftempfe2d(ji,jj) =ftempfe2d(ji,jj)  + (ftempfe(ji,jj) * fse3t(ji,jj,jk)) 
1547                  ENDIF
1548                  IF( med_diag%FASTC%dgsave ) THEN
1549                      ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
1550                  ENDIF
1551                  IF( med_diag%FASTCA%dgsave ) THEN
1552                      ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca(ji,jj) * fse3t(ji,jj,jk))
1553                  ENDIF
1554                  !!
1555                  IF( med_diag%REMINN%dgsave ) THEN
1556                      freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn(ji,jj)  * fse3t(ji,jj,jk))
1557                  ENDIF
1558                  IF( med_diag%REMINSI%dgsave ) THEN
1559                      freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
1560                  ENDIF
1561                  IF( med_diag%REMINFE%dgsave ) THEN
1562                      freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe(ji,jj) * fse3t(ji,jj,jk)) 
1563                  ENDIF
1564                  IF( med_diag%REMINC%dgsave ) THEN
1565                      freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc(ji,jj)  * fse3t(ji,jj,jk)) 
1566                  ENDIF
1567                  IF( med_diag%REMINCA%dgsave ) THEN
1568                      freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca(ji,jj) * fse3t(ji,jj,jk)) 
1569                  ENDIF
1570                  !!
1571# if defined key_roam
1572                  !!
1573                  !! AXY (09/11/16): CMIP6 diagnostics
1574                  IF( med_diag%FD_NIT3%dgsave ) THEN
1575                     fd_nit3(ji,jj,jk) = ffastn(ji,jj)
1576                  ENDIF
1577                  IF( med_diag%FD_SIL3%dgsave ) THEN
1578                     fd_sil3(ji,jj,jk) = ffastsi(ji,jj)
1579                  ENDIF
1580                  IF( med_diag%FD_CAR3%dgsave ) THEN
1581                     fd_car3(ji,jj,jk) = ffastc(ji,jj)
1582                  ENDIF
1583                  IF( med_diag%FD_CAL3%dgsave ) THEN
1584                     fd_cal3(ji,jj,jk) = ffastca(ji,jj)
1585                  ENDIF
1586                  !!
1587                  IF (jk.eq.i0100) THEN
1588                     IF( med_diag%RR_0100%dgsave ) THEN
1589                        ffastca2d(ji,jj) =   &
1590                        ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
1591                     ENDIF                     
1592                  ELSE IF (jk.eq.i0500) THEN
1593                     IF( med_diag%RR_0500%dgsave ) THEN
1594                        ffastca2d(ji,jj) =   &
1595                        ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
1596                     ENDIF                       
1597                  ELSE IF (jk.eq.i1000) THEN
1598                     IF( med_diag%RR_1000%dgsave ) THEN
1599                        ffastca2d(ji,jj) =   &
1600                        ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
1601                     ENDIF
1602                  ELSE IF (jk.eq.mbathy(ji,jj)) THEN
1603                     IF( med_diag%IBEN_N%dgsave ) THEN
1604                        iben_n2d(ji,jj) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
1605                     ENDIF
1606                     IF( med_diag%IBEN_FE%dgsave ) THEN
1607                        iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
1608                     ENDIF
1609                     IF( med_diag%IBEN_C%dgsave ) THEN
1610                        iben_c2d(ji,jj) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
1611                     ENDIF
1612                     IF( med_diag%IBEN_SI%dgsave ) THEN
1613                        iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
1614                     ENDIF
1615                     IF( med_diag%IBEN_CA%dgsave ) THEN
1616                        iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
1617                     ENDIF
1618                     IF( med_diag%OBEN_N%dgsave ) THEN
1619                        oben_n2d(ji,jj) = f_benout_n(ji,jj)
1620                     ENDIF
1621                     IF( med_diag%OBEN_FE%dgsave ) THEN
1622                        oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
1623                     ENDIF
1624                     IF( med_diag%OBEN_C%dgsave ) THEN
1625                        oben_c2d(ji,jj) = f_benout_c(ji,jj)
1626                     ENDIF
1627                     IF( med_diag%OBEN_SI%dgsave ) THEN
1628                        oben_si2d(ji,jj) = f_benout_si(ji,jj)
1629                     ENDIF
1630                     IF( med_diag%OBEN_CA%dgsave ) THEN
1631                        oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
1632                     ENDIF
1633                     IF( med_diag%SFR_OCAL%dgsave ) THEN
1634                        sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
1635                     ENDIF
1636                     IF( med_diag%SFR_OARG%dgsave ) THEN
1637                        sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
1638                     ENDIF
1639                     IF( med_diag%LYSO_CA%dgsave ) THEN
1640                        lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
1641                     ENDIF
1642                  ENDIF
1643                  !! end bathy-1 diags
1644                  !!
1645                  IF( med_diag%RIV_N%dgsave ) THEN
1646                     rivn2d(ji,jj) = rivn2d(ji,jj) +  (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
1647                  ENDIF
1648                  IF( med_diag%RIV_SI%dgsave ) THEN
1649                     rivsi2d(ji,jj) = rivsi2d(ji,jj) +  (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
1650                  ENDIF
1651                  IF( med_diag%RIV_C%dgsave ) THEN
1652                     rivc2d(ji,jj) = rivc2d(ji,jj) +  (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
1653                  ENDIF
1654                  IF( med_diag%RIV_ALK%dgsave ) THEN
1655                     rivalk2d(ji,jj) = rivalk2d(ji,jj) +  (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
1656                  ENDIF
1657                  IF( med_diag%DETC%dgsave ) THEN
1658                     fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc(ji,jj)  * fse3t(ji,jj,jk))   
1659                  ENDIF
1660                  !!
1661                  !!             
1662                  !!
1663                  IF( med_diag%PN_LLOSS%dgsave ) THEN
1664                     fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
1665                  ENDIF
1666                  IF( med_diag%PD_LLOSS%dgsave ) THEN
1667                     fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
1668                  ENDIF
1669                  IF( med_diag%ZI_LLOSS%dgsave ) THEN
1670                     fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
1671                  ENDIF
1672                  IF( med_diag%ZE_LLOSS%dgsave ) THEN
1673                     fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2(ji,jj) * fse3t(ji,jj,jk))
1674                  ENDIF
1675                  IF( med_diag%ZI_MES_N%dgsave ) THEN
1676                     zimesn2d(ji,jj) = zimesn2d(ji,jj) +  &
1677                     (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)) * fse3t(ji,jj,jk))
1678                  ENDIF
1679                  IF( med_diag%ZI_MES_D%dgsave ) THEN
1680                     zimesd2d(ji,jj) = zimesd2d(ji,jj) + & 
1681                     ((1. - xbetan) * finmi(ji,jj) * fse3t(ji,jj,jk))
1682                  ENDIF
1683                  IF( med_diag%ZI_MES_C%dgsave ) THEN
1684                     zimesc2d(ji,jj) = zimesc2d(ji,jj) + &
1685                     (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj)) * fse3t(ji,jj,jk))
1686                  ENDIF
1687                  IF( med_diag%ZI_MESDC%dgsave ) THEN
1688                     zimesdc2d(ji,jj) = zimesdc2d(ji,jj) + &
1689                     ((1. - xbetac) * ficmi(ji,jj) * fse3t(ji,jj,jk))
1690                  ENDIF
1691                  IF( med_diag%ZI_EXCR%dgsave ) THEN
1692                     ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +  (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
1693                  ENDIF
1694                  IF( med_diag%ZI_RESP%dgsave ) THEN
1695                     ziresp2d(ji,jj) = ziresp2d(ji,jj) +  (fmiresp(ji,jj) * fse3t(ji,jj,jk))
1696                  ENDIF
1697                  IF( med_diag%ZI_GROW%dgsave ) THEN
1698                     zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow(ji,jj) * fse3t(ji,jj,jk))
1699                  ENDIF
1700                  IF( med_diag%ZE_MES_N%dgsave ) THEN
1701                     zemesn2d(ji,jj) = zemesn2d(ji,jj) + &
1702                     (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)) * fse3t(ji,jj,jk))
1703                  ENDIF
1704                  IF( med_diag%ZE_MES_D%dgsave ) THEN
1705                     zemesd2d(ji,jj) = zemesd2d(ji,jj) + &
1706                     ((1. - xbetan) * finme(ji,jj) * fse3t(ji,jj,jk))
1707                  ENDIF
1708                  IF( med_diag%ZE_MES_C%dgsave ) THEN
1709                     zemesc2d(ji,jj) = zemesc2d(ji,jj) +                         & 
1710                     (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +  &
1711                     (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj)) * fse3t(ji,jj,jk))
1712                  ENDIF
1713                  IF( med_diag%ZE_MESDC%dgsave ) THEN
1714                     zemesdc2d(ji,jj) = zemesdc2d(ji,jj) +  &
1715                     ((1. - xbetac) * ficme(ji,jj) * fse3t(ji,jj,jk))
1716                  ENDIF
1717                  IF( med_diag%ZE_EXCR%dgsave ) THEN
1718                     zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
1719                  ENDIF
1720                  IF( med_diag%ZE_RESP%dgsave ) THEN
1721                     zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp(ji,jj) * fse3t(ji,jj,jk))
1722                  ENDIF
1723                  IF( med_diag%ZE_GROW%dgsave ) THEN
1724                     zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow(ji,jj) * fse3t(ji,jj,jk))
1725                  ENDIF
1726                  IF( med_diag%MDETC%dgsave ) THEN
1727                     mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc(ji,jj) * fse3t(ji,jj,jk))
1728                  ENDIF
1729                  IF( med_diag%GMIDC%dgsave ) THEN
1730                     gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc(ji,jj) * fse3t(ji,jj,jk))
1731                  ENDIF
1732                  IF( med_diag%GMEDC%dgsave ) THEN
1733                     gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc(ji,jj)  * fse3t(ji,jj,jk))
1734                  ENDIF
1735                  !!
1736# endif                   
1737                  !!
1738                  !! ** 3D diagnostics
1739                  IF( med_diag%TPP3%dgsave ) THEN
1740                     tpp3d(ji,jj,jk) =  (fprn(ji,jj) * zphn(ji,jj)) + (fprd(ji,jj) * zphd(ji,jj))
1741                     !CALL iom_put( "TPP3"  , tpp3d )
1742                  ENDIF
1743                  IF( med_diag%TPPD3%dgsave ) THEN
1744                     tppd3(ji,jj,jk) =  (fprd(ji,jj) * zphd(ji,jj))
1745                  ENDIF
1746                 
1747                  IF( med_diag%REMIN3N%dgsave ) THEN
1748                     remin3dn(ji,jj,jk) = fregen(ji,jj) + (freminn(ji,jj) * fse3t(ji,jj,jk)) !! remineralisation
1749                     !CALL iom_put( "REMIN3N"  , remin3dn )
1750                  ENDIF
1751                  !! IF( med_diag%PH3%dgsave ) THEN
1752                  !!     CALL iom_put( "PH3"  , f3_pH )
1753                  !! ENDIF
1754                  !! IF( med_diag%OM_CAL3%dgsave ) THEN
1755                  !!     CALL iom_put( "OM_CAL3"  , f3_omcal )
1756                  !! ENDIF
1757        !!
1758        !! AXY (09/11/16): CMIP6 diagnostics
1759        IF ( med_diag%DCALC3%dgsave   ) THEN
1760                     dcalc3(ji,jj,jk) = freminca(ji,jj)
1761                  ENDIF
1762        IF ( med_diag%FEDISS3%dgsave  ) THEN
1763                     fediss3(ji,jj,jk) = ffetop(ji,jj)
1764                  ENDIF
1765        IF ( med_diag%FESCAV3%dgsave  ) THEN
1766                     fescav3(ji,jj,jk) = ffescav(ji,jj)
1767                  ENDIF
1768        IF ( med_diag%MIGRAZP3%dgsave ) THEN
1769                     migrazp3(ji,jj,jk) = fgmipn(ji,jj) * xthetapn
1770                  ENDIF
1771        IF ( med_diag%MIGRAZD3%dgsave ) THEN
1772                     migrazd3(ji,jj,jk) = fgmidc(ji,jj)
1773                  ENDIF
1774        IF ( med_diag%MEGRAZP3%dgsave ) THEN
1775                     megrazp3(ji,jj,jk) = (fgmepn(ji,jj) * xthetapn) + (fgmepd(ji,jj) * xthetapd)
1776                  ENDIF
1777        IF ( med_diag%MEGRAZD3%dgsave ) THEN
1778                     megrazd3(ji,jj,jk) = fgmedc(ji,jj)
1779                  ENDIF
1780        IF ( med_diag%MEGRAZZ3%dgsave ) THEN
1781                     megrazz3(ji,jj,jk) = (fgmezmi(ji,jj) * xthetazmi)
1782                  ENDIF
1783        IF ( med_diag%PBSI3%dgsave    ) THEN
1784                     pbsi3(ji,jj,jk)    = (fprds(ji,jj) * zpds(ji,jj))
1785                  ENDIF
1786        IF ( med_diag%PCAL3%dgsave    ) THEN
1787                     pcal3(ji,jj,jk)    = ftempca(ji,jj)
1788                  ENDIF
1789        IF ( med_diag%REMOC3%dgsave   ) THEN
1790                     remoc3(ji,jj,jk)   = freminc(ji,jj)
1791                  ENDIF
1792        IF ( med_diag%PNLIMJ3%dgsave  ) THEN
1793                     ! pnlimj3(ji,jj,jk)  = fjln(ji,jj)
1794                     pnlimj3(ji,jj,jk)  = fjlim_pn(ji,jj)
1795                  ENDIF
1796        IF ( med_diag%PNLIMN3%dgsave  ) THEN
1797                     pnlimn3(ji,jj,jk)  = fnln(ji,jj)
1798                  ENDIF
1799        IF ( med_diag%PNLIMFE3%dgsave ) THEN
1800                     pnlimfe3(ji,jj,jk) = ffln2(ji,jj)
1801                  ENDIF
1802        IF ( med_diag%PDLIMJ3%dgsave  ) THEN
1803                     ! pdlimj3(ji,jj,jk)  = fjld(ji,jj)
1804                     pdlimj3(ji,jj,jk)  = fjlim_pd(ji,jj)
1805                  ENDIF
1806        IF ( med_diag%PDLIMN3%dgsave  ) THEN
1807                     pdlimn3(ji,jj,jk)  = fnld(ji,jj)
1808                  ENDIF
1809        IF ( med_diag%PDLIMFE3%dgsave ) THEN
1810                     pdlimfe3(ji,jj,jk) = ffld(ji,jj)
1811                  ENDIF
1812        IF ( med_diag%PDLIMSI3%dgsave ) THEN
1813                     pdlimsi3(ji,jj,jk) = fsld2(ji,jj)
1814                  ENDIF
1815                  !!
1816                  !! ** Without using iom_use
1817               ELSE IF( ln_diatrc ) THEN
1818#   if defined key_debug_medusa
1819                  IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk ln_diatrc'
1820                  CALL flush(numout)
1821#   endif
1822                  !!----------------------------------------------------------------------
1823                  !! Prepare 2D diagnostics
1824                  !!----------------------------------------------------------------------
1825                  !!
1826                  !! if ((kt / 240*240).eq.kt) then
1827                  !!    IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt
1828                  !! endif     
1829                  trc2d(ji,jj,1)  =  ftot_n(ji,jj)                             !! nitrogen inventory
1830                  trc2d(ji,jj,2)  =  ftot_si(ji,jj)                            !! silicon  inventory
1831                  trc2d(ji,jj,3)  =  ftot_fe(ji,jj)                            !! iron     inventory
1832                  trc2d(ji,jj,4)  = trc2d(ji,jj,4)  + (fprn(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom production
1833                  trc2d(ji,jj,5)  = trc2d(ji,jj,5)  + (fdpn(ji,jj)         * fse3t(ji,jj,jk))    !! non-diatom non-grazing losses
1834                  trc2d(ji,jj,6)  = trc2d(ji,jj,6)  + (fprd(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom production
1835                  trc2d(ji,jj,7)  = trc2d(ji,jj,7)  + (fdpd(ji,jj)         * fse3t(ji,jj,jk))    !! diatom non-grazing losses
1836                  !! diagnostic field  8 is (ostensibly) supplied by trcsed.F           
1837                  trc2d(ji,jj,9)  = trc2d(ji,jj,9)  + (fprds(ji,jj) * zpds(ji,jj) * fse3t(ji,jj,jk))    !! diatom silicon production
1838                  trc2d(ji,jj,10) = trc2d(ji,jj,10) + (fsdiss(ji,jj)  * fse3t(ji,jj,jk))         !! diatom silicon dissolution
1839                  trc2d(ji,jj,11) = trc2d(ji,jj,11) + (fgmipn(ji,jj)  * fse3t(ji,jj,jk))         !! microzoo grazing on non-diatoms
1840                  trc2d(ji,jj,12) = trc2d(ji,jj,12) + (fgmid(ji,jj)   * fse3t(ji,jj,jk))         !! microzoo grazing on detrital nitrogen
1841                  trc2d(ji,jj,13) = trc2d(ji,jj,13) + (fdzmi(ji,jj)   * fse3t(ji,jj,jk))         !! microzoo non-grazing losses
1842                  trc2d(ji,jj,14) = trc2d(ji,jj,14) + (fgmepn(ji,jj)  * fse3t(ji,jj,jk))         !! mesozoo  grazing on non-diatoms
1843                  trc2d(ji,jj,15) = trc2d(ji,jj,15) + (fgmepd(ji,jj)  * fse3t(ji,jj,jk))         !! mesozoo  grazing on diatoms
1844                  trc2d(ji,jj,16) = trc2d(ji,jj,16) + (fgmezmi(ji,jj) * fse3t(ji,jj,jk))         !! mesozoo  grazing on microzoo
1845                  trc2d(ji,jj,17) = trc2d(ji,jj,17) + (fgmed(ji,jj)   * fse3t(ji,jj,jk))         !! mesozoo  grazing on detrital nitrogen
1846                  trc2d(ji,jj,18) = trc2d(ji,jj,18) + (fdzme(ji,jj)   * fse3t(ji,jj,jk))         !! mesozoo  non-grazing losses
1847                  !! diagnostic field 19 is (ostensibly) supplied by trcexp.F
1848                  trc2d(ji,jj,20) = trc2d(ji,jj,20) + (fslown(ji,jj)  * fse3t(ji,jj,jk))         !! slow sinking detritus N production
1849                  trc2d(ji,jj,21) = trc2d(ji,jj,21) + (fdd(ji,jj)     * fse3t(ji,jj,jk))         !! detrital remineralisation
1850                  trc2d(ji,jj,22) = trc2d(ji,jj,22) + (ffetop(ji,jj)  * fse3t(ji,jj,jk))         !! aeolian  iron addition
1851                  trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot(ji,jj)  * fse3t(ji,jj,jk))         !! seafloor iron addition
1852                  trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav(ji,jj) * fse3t(ji,jj,jk))         !! "free" iron scavenging
1853                  trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) !! non-diatom J  limitation term
1854                  trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom N  limitation term
1855                  trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln2(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom Fe limitation term
1856                  trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) !! diatom     J  limitation term
1857                  trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     N  limitation term
1858                  trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Fe limitation term
1859                  trc2d(ji,jj,31) = trc2d(ji,jj,31) + (fsld2(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Si limitation term
1860                  trc2d(ji,jj,32) = trc2d(ji,jj,32) + (fsld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Si uptake limitation term
1861                  if (jk.eq.i0100) trc2d(ji,jj,33) = fslownflux(ji,jj)         !! slow detritus flux at  100 m
1862                  if (jk.eq.i0200) trc2d(ji,jj,34) = fslownflux(ji,jj)         !! slow detritus flux at  200 m
1863                  if (jk.eq.i0500) trc2d(ji,jj,35) = fslownflux(ji,jj)         !! slow detritus flux at  500 m
1864                  if (jk.eq.i1000) trc2d(ji,jj,36) = fslownflux(ji,jj)         !! slow detritus flux at 1000 m
1865                  trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen(ji,jj)                   !! non-fast N  full column regeneration
1866                  trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi(ji,jj)                 !! non-fast Si full column regeneration
1867                  if (jk.eq.i0100) trc2d(ji,jj,39) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  100 m
1868                  if (jk.eq.i0200) trc2d(ji,jj,40) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  200 m
1869                  if (jk.eq.i0500) trc2d(ji,jj,41) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  500 m
1870                  if (jk.eq.i1000) trc2d(ji,jj,42) = trc2d(ji,jj,37)           !! non-fast N  regeneration to 1000 m
1871                  trc2d(ji,jj,43) = trc2d(ji,jj,43) + (ftempn(ji,jj)  * fse3t(ji,jj,jk))         !! fast sinking detritus N production
1872                  trc2d(ji,jj,44) = trc2d(ji,jj,44) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus Si production
1873                  trc2d(ji,jj,45) = trc2d(ji,jj,45) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus Fe production
1874                  trc2d(ji,jj,46) = trc2d(ji,jj,46) + (ftempc(ji,jj)  * fse3t(ji,jj,jk))         !! fast sinking detritus C production
1875                  trc2d(ji,jj,47) = trc2d(ji,jj,47) + (ftempca(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus CaCO3 production
1876                  if (jk.eq.i0100) trc2d(ji,jj,48) = ffastn(ji,jj)             !! fast detritus N  flux at  100 m
1877                  if (jk.eq.i0200) trc2d(ji,jj,49) = ffastn(ji,jj)             !! fast detritus N  flux at  200 m
1878                  if (jk.eq.i0500) trc2d(ji,jj,50) = ffastn(ji,jj)             !! fast detritus N  flux at  500 m
1879                  if (jk.eq.i1000) trc2d(ji,jj,51) = ffastn(ji,jj)             !! fast detritus N  flux at 1000 m
1880                  if (jk.eq.i0100) trc2d(ji,jj,52) = fregenfast(ji,jj)         !! N  regeneration to  100 m
1881                  if (jk.eq.i0200) trc2d(ji,jj,53) = fregenfast(ji,jj)         !! N  regeneration to  200 m
1882                  if (jk.eq.i0500) trc2d(ji,jj,54) = fregenfast(ji,jj)         !! N  regeneration to  500 m
1883                  if (jk.eq.i1000) trc2d(ji,jj,55) = fregenfast(ji,jj)         !! N  regeneration to 1000 m
1884                  if (jk.eq.i0100) trc2d(ji,jj,56) = ffastsi(ji,jj)            !! fast detritus Si flux at  100 m
1885                  if (jk.eq.i0200) trc2d(ji,jj,57) = ffastsi(ji,jj)            !! fast detritus Si flux at  200 m
1886                  if (jk.eq.i0500) trc2d(ji,jj,58) = ffastsi(ji,jj)            !! fast detritus Si flux at  500 m
1887                  if (jk.eq.i1000) trc2d(ji,jj,59) = ffastsi(ji,jj)            !! fast detritus Si flux at 1000 m
1888                  if (jk.eq.i0100) trc2d(ji,jj,60) = fregenfastsi(ji,jj)       !! Si regeneration to  100 m
1889                  if (jk.eq.i0200) trc2d(ji,jj,61) = fregenfastsi(ji,jj)       !! Si regeneration to  200 m
1890                  if (jk.eq.i0500) trc2d(ji,jj,62) = fregenfastsi(ji,jj)       !! Si regeneration to  500 m
1891                  if (jk.eq.i1000) trc2d(ji,jj,63) = fregenfastsi(ji,jj)       !! Si regeneration to 1000 m
1892                  trc2d(ji,jj,64) = trc2d(ji,jj,64) + (freminn(ji,jj)  * fse3t(ji,jj,jk))        !! sum of fast-sinking N  fluxes
1893                  trc2d(ji,jj,65) = trc2d(ji,jj,65) + (freminsi(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Si fluxes
1894                  trc2d(ji,jj,66) = trc2d(ji,jj,66) + (freminfe(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Fe fluxes
1895                  trc2d(ji,jj,67) = trc2d(ji,jj,67) + (freminc(ji,jj)  * fse3t(ji,jj,jk))        !! sum of fast-sinking C  fluxes
1896                  trc2d(ji,jj,68) = trc2d(ji,jj,68) + (freminca(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Ca fluxes
1897                  if (jk.eq.mbathy(ji,jj)) then
1898                     trc2d(ji,jj,69) = fsedn(ji,jj)                                   !! N  sedimentation flux                                 
1899                     trc2d(ji,jj,70) = fsedsi(ji,jj)                                  !! Si sedimentation flux
1900                     trc2d(ji,jj,71) = fsedfe(ji,jj)                                  !! Fe sedimentation flux
1901                     trc2d(ji,jj,72) = fsedc(ji,jj)                                   !! C  sedimentation flux
1902                     trc2d(ji,jj,73) = fsedca(ji,jj)                                  !! Ca sedimentation flux
1903                  endif
1904                  if (jk.eq.1)  trc2d(ji,jj,74) = qsr(ji,jj)
1905                  if (jk.eq.1)  trc2d(ji,jj,75) = xpar(ji,jj,jk)
1906                  !! if (jk.eq.1)  trc2d(ji,jj,75) = real(iters(ji,jj))
1907                  !! diagnostic fields 76 to 80 calculated below
1908                  trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj)           !! mixed layer non-diatom production
1909                  trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj)           !! mixed layer     diatom production
1910# if defined key_gulf_finland
1911                  if (jk.eq.1)  trc2d(ji,jj,83) = real(ibio_switch)            !! Gulf of Finland check
1912# else
1913                  trc2d(ji,jj,83) = ocal_ccd(ji,jj)                            !! calcite CCD depth
1914# endif
1915                  trc2d(ji,jj,84) = fccd(ji,jj)                                !! last model level above calcite CCD depth
1916                  if (jk.eq.1)     trc2d(ji,jj,85) = xFree(ji,jj)              !! surface "free" iron
1917                  if (jk.eq.i0200) trc2d(ji,jj,86) = xFree(ji,jj)              !! "free" iron at  100 m
1918                  if (jk.eq.i0200) trc2d(ji,jj,87) = xFree(ji,jj)              !! "free" iron at  200 m
1919                  if (jk.eq.i0500) trc2d(ji,jj,88) = xFree(ji,jj)              !! "free" iron at  500 m
1920                  if (jk.eq.i1000) trc2d(ji,jj,89) = xFree(ji,jj)              !! "free" iron at 1000 m
1921                  !! AXY (27/06/12): extract "euphotic depth"
1922                  if (jk.eq.1)     trc2d(ji,jj,90) = xze(ji,jj)
1923                  !!
1924# if defined key_roam
1925                  !! ROAM provisionally has access to a further 20 2D diagnostics
1926                  if (jk .eq. 1) then
1927                     trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + wndm(ji,jj)              !! surface wind
1928                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm(ji,jj)           !! atmospheric pCO2
1929                     trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph(ji,jj)                !! ocean pH
1930                     trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w(ji,jj)             !! ocean pCO2
1931                     trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3(ji,jj)             !! ocean H2CO3 conc.
1932                     trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3(ji,jj)              !! ocean HCO3 conc.
1933                     trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3(ji,jj)               !! ocean CO3 conc.
1934                     trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux(ji,jj)           !! air-sea CO2 flux
1935                     trc2d(ji,jj,99)  = trc2d(ji,jj,99)  + f_omcal(ji,jj)      !! ocean omega calcite
1936                     trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj)      !! ocean omega aragonite
1937                     trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC(ji,jj)              !! ocean TDIC
1938                     trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK(ji,jj)              !! ocean TALK
1939                     trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660(ji,jj)             !! surface kw660
1940                     trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0(ji,jj)               !! surface pressure
1941                     trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux(ji,jj)            !! air-sea O2 flux
1942                     trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat(ji,jj)             !! ocean O2 saturation
1943                     trc2d(ji,jj,107) = f2_ccd_cal(ji,jj)                      !! depth calcite CCD
1944                     trc2d(ji,jj,108) = f2_ccd_arg(ji,jj)                      !! depth aragonite CCD
1945                  endif
1946                  if (jk .eq. mbathy(ji,jj)) then
1947                     trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)                     !! seafloor omega calcite
1948                     trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)                     !! seafloor omega aragonite
1949                  endif
1950                  !! diagnostic fields 111 to 117 calculated below
1951                  if (jk.eq.i0100) trc2d(ji,jj,118) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at  100 m
1952                  if (jk.eq.i0500) trc2d(ji,jj,119) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at  500 m
1953                  if (jk.eq.i1000) trc2d(ji,jj,120) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at 1000 m
1954                  !! AXY (18/01/12): benthic flux diagnostics
1955                  if (jk.eq.mbathy(ji,jj)) then
1956                     trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
1957                     trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
1958                     trc2d(ji,jj,123) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
1959                     trc2d(ji,jj,124) = f_fbenin_si(ji,jj)
1960                     trc2d(ji,jj,125) = f_fbenin_ca(ji,jj)
1961                     trc2d(ji,jj,126) = f_benout_n(ji,jj)
1962                     trc2d(ji,jj,127) = f_benout_fe(ji,jj)
1963                     trc2d(ji,jj,128) = f_benout_c(ji,jj)
1964                     trc2d(ji,jj,129) = f_benout_si(ji,jj)
1965                     trc2d(ji,jj,130) = f_benout_ca(ji,jj)
1966                  endif
1967                  !! diagnostics fields 131 to 135 calculated below
1968                  trc2d(ji,jj,136) = f_runoff(ji,jj)
1969                  !! AXY (19/07/12): amended to allow for riverine nutrient addition below surface
1970                  trc2d(ji,jj,137) = trc2d(ji,jj,137) + (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
1971                  trc2d(ji,jj,138) = trc2d(ji,jj,138) + (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
1972                  trc2d(ji,jj,139) = trc2d(ji,jj,139) + (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
1973                  trc2d(ji,jj,140) = trc2d(ji,jj,140) + (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
1974                  trc2d(ji,jj,141) = trc2d(ji,jj,141) + (fslowc(ji,jj)  * fse3t(ji,jj,jk))       !! slow sinking detritus C production
1975                  if (jk.eq.i0100) trc2d(ji,jj,142) = fslowcflux(ji,jj)        !! slow detritus flux at  100 m
1976                  if (jk.eq.i0200) trc2d(ji,jj,143) = fslowcflux(ji,jj)        !! slow detritus flux at  200 m
1977                  if (jk.eq.i0500) trc2d(ji,jj,144) = fslowcflux(ji,jj)        !! slow detritus flux at  500 m
1978                  if (jk.eq.i1000) trc2d(ji,jj,145) = fslowcflux(ji,jj)        !! slow detritus flux at 1000 m
1979                  trc2d(ji,jj,146)  = trc2d(ji,jj,146)  + ftot_c(ji,jj)        !! carbon     inventory
1980                  trc2d(ji,jj,147)  = trc2d(ji,jj,147)  + ftot_a(ji,jj)        !! alkalinity inventory
1981                  trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)       !! oxygen     inventory
1982                  if (jk.eq.mbathy(ji,jj)) then
1983                     trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
1984                  endif
1985                  trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)                  !! community respiration
1986        !!
1987        !! AXY (14/02/14): a Valentines Day gift to BASIN - a shedload of new
1988                  !!                 diagnostics that they'll most likely never need!
1989                  !!                 (actually, as with all such gifts, I'm giving them
1990                  !!                 some things I'd like myself!)
1991                  !!
1992                  !! ----------------------------------------------------------------------
1993                  !! linear losses
1994                  !! non-diatom
1995                  trc2d(ji,jj,151) = trc2d(ji,jj,151) + (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
1996                  !! diatom
1997                  trc2d(ji,jj,152) = trc2d(ji,jj,152) + (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
1998                  !! microzooplankton
1999                  trc2d(ji,jj,153) = trc2d(ji,jj,153) + (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
2000                  !! mesozooplankton
2001                  trc2d(ji,jj,154) = trc2d(ji,jj,154) + (fdzme2(ji,jj) * fse3t(ji,jj,jk))
2002                  !! ----------------------------------------------------------------------
2003                  !! microzooplankton grazing
2004                  !! microzooplankton messy -> N
2005                  trc2d(ji,jj,155) = trc2d(ji,jj,155) + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)) * fse3t(ji,jj,jk))
2006                  !! microzooplankton messy -> D
2007                  trc2d(ji,jj,156) = trc2d(ji,jj,156) + ((1. - xbetan) * finmi(ji,jj) * fse3t(ji,jj,jk))
2008                  !! microzooplankton messy -> DIC
2009                  trc2d(ji,jj,157) = trc2d(ji,jj,157) + (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj)) * fse3t(ji,jj,jk))
2010                  !! microzooplankton messy -> Dc
2011                  trc2d(ji,jj,158) = trc2d(ji,jj,158) + ((1. - xbetac) * ficmi(ji,jj) * fse3t(ji,jj,jk))
2012                  !! microzooplankton excretion
2013                  trc2d(ji,jj,159) = trc2d(ji,jj,159) + (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
2014                  !! microzooplankton respiration
2015                  trc2d(ji,jj,160) = trc2d(ji,jj,160) + (fmiresp(ji,jj) * fse3t(ji,jj,jk))
2016                  !! microzooplankton growth
2017                  trc2d(ji,jj,161) = trc2d(ji,jj,161) + (fmigrow(ji,jj) * fse3t(ji,jj,jk))
2018                  !! ----------------------------------------------------------------------
2019                  !! mesozooplankton grazing
2020                  !! mesozooplankton messy -> N
2021                  trc2d(ji,jj,162) = trc2d(ji,jj,162) + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)) * fse3t(ji,jj,jk))
2022                  !! mesozooplankton messy -> D
2023                  trc2d(ji,jj,163) = trc2d(ji,jj,163) + ((1. - xbetan) * finme(ji,jj) * fse3t(ji,jj,jk))
2024                  !! mesozooplankton messy -> DIC
2025                  trc2d(ji,jj,164) = trc2d(ji,jj,164) + (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) + &
2026                  &                  (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj)) * fse3t(ji,jj,jk))
2027                  !! mesozooplankton messy -> Dc
2028                  trc2d(ji,jj,165) = trc2d(ji,jj,165) + ((1. - xbetac) * ficme(ji,jj) * fse3t(ji,jj,jk))
2029                  !! mesozooplankton excretion
2030                  trc2d(ji,jj,166) = trc2d(ji,jj,166) + (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
2031                  !! mesozooplankton respiration
2032                  trc2d(ji,jj,167) = trc2d(ji,jj,167) + (fmeresp(ji,jj) * fse3t(ji,jj,jk))
2033                  !! mesozooplankton growth
2034                  trc2d(ji,jj,168) = trc2d(ji,jj,168) + (fmegrow(ji,jj) * fse3t(ji,jj,jk))
2035                  !! ----------------------------------------------------------------------
2036                  !! miscellaneous
2037                  trc2d(ji,jj,169) = trc2d(ji,jj,169) + (fddc(ji,jj)    * fse3t(ji,jj,jk)) !! detrital C remineralisation
2038                  trc2d(ji,jj,170) = trc2d(ji,jj,170) + (fgmidc(ji,jj)  * fse3t(ji,jj,jk)) !! microzoo grazing on detrital carbon
2039                  trc2d(ji,jj,171) = trc2d(ji,jj,171) + (fgmedc(ji,jj)  * fse3t(ji,jj,jk)) !! mesozoo  grazing on detrital carbon
2040                  !!
2041                  !! ----------------------------------------------------------------------
2042        !!
2043        !! AXY (23/10/14): extract primary production related surface fields to
2044        !!                 deal with diel cycle issues; hijacking BASIN 150m
2045        !!                 diagnostics to do so (see commented out diagnostics
2046        !!                 below this section)
2047        !!
2048                  !! extract relevant BASIN fields at 150m
2049                  if (jk .eq. i0150) then
2050                     trc2d(ji,jj,172) = trc2d(ji,jj,4)    !! Pn PP
2051                     trc2d(ji,jj,173) = trc2d(ji,jj,151)  !! Pn linear loss
2052                     trc2d(ji,jj,174) = trc2d(ji,jj,5)    !! Pn non-linear loss
2053                     trc2d(ji,jj,175) = trc2d(ji,jj,11)   !! Pn grazing to Zmi
2054                     trc2d(ji,jj,176) = trc2d(ji,jj,14)   !! Pn grazing to Zme
2055                     trc2d(ji,jj,177) = trc2d(ji,jj,6)    !! Pd PP
2056                     trc2d(ji,jj,178) = trc2d(ji,jj,152)  !! Pd linear loss
2057                     trc2d(ji,jj,179) = trc2d(ji,jj,7)    !! Pd non-linear loss
2058                     trc2d(ji,jj,180) = trc2d(ji,jj,15)   !! Pd grazing to Zme
2059                     trc2d(ji,jj,181) = trc2d(ji,jj,12)   !! Zmi grazing on D
2060                     trc2d(ji,jj,182) = trc2d(ji,jj,170)  !! Zmi grazing on Dc
2061                     trc2d(ji,jj,183) = trc2d(ji,jj,155)  !! Zmi messy feeding loss to N
2062                     trc2d(ji,jj,184) = trc2d(ji,jj,156)  !! Zmi messy feeding loss to D
2063                     trc2d(ji,jj,185) = trc2d(ji,jj,157)  !! Zmi messy feeding loss to DIC
2064                     trc2d(ji,jj,186) = trc2d(ji,jj,158)  !! Zmi messy feeding loss to Dc
2065                     trc2d(ji,jj,187) = trc2d(ji,jj,159)  !! Zmi excretion
2066                     trc2d(ji,jj,188) = trc2d(ji,jj,160)  !! Zmi respiration
2067                     trc2d(ji,jj,189) = trc2d(ji,jj,161)  !! Zmi growth
2068                     trc2d(ji,jj,190) = trc2d(ji,jj,153)  !! Zmi linear loss
2069                     trc2d(ji,jj,191) = trc2d(ji,jj,13)   !! Zmi non-linear loss
2070                     trc2d(ji,jj,192) = trc2d(ji,jj,16)   !! Zmi grazing to Zme
2071                     trc2d(ji,jj,193) = trc2d(ji,jj,17)   !! Zme grazing on D
2072                     trc2d(ji,jj,194) = trc2d(ji,jj,171)  !! Zme grazing on Dc
2073                     trc2d(ji,jj,195) = trc2d(ji,jj,162)  !! Zme messy feeding loss to N
2074                     trc2d(ji,jj,196) = trc2d(ji,jj,163)  !! Zme messy feeding loss to D
2075                     trc2d(ji,jj,197) = trc2d(ji,jj,164)  !! Zme messy feeding loss to DIC
2076                     trc2d(ji,jj,198) = trc2d(ji,jj,165)  !! Zme messy feeding loss to Dc
2077                     trc2d(ji,jj,199) = trc2d(ji,jj,166)  !! Zme excretion
2078                     trc2d(ji,jj,200) = trc2d(ji,jj,167)  !! Zme respiration
2079                     trc2d(ji,jj,201) = trc2d(ji,jj,168)  !! Zme growth
2080                     trc2d(ji,jj,202) = trc2d(ji,jj,154)  !! Zme linear loss
2081                     trc2d(ji,jj,203) = trc2d(ji,jj,18)   !! Zme non-linear loss
2082                     trc2d(ji,jj,204) = trc2d(ji,jj,20)   !! Slow detritus production, N
2083                     trc2d(ji,jj,205) = trc2d(ji,jj,21)   !! Slow detritus remineralisation, N
2084                     trc2d(ji,jj,206) = trc2d(ji,jj,141)  !! Slow detritus production, C
2085                     trc2d(ji,jj,207) = trc2d(ji,jj,169)  !! Slow detritus remineralisation, C
2086                     trc2d(ji,jj,208) = trc2d(ji,jj,43)   !! Fast detritus production, N
2087                     trc2d(ji,jj,209) = trc2d(ji,jj,21)   !! Fast detritus remineralisation, N
2088                     trc2d(ji,jj,210) = trc2d(ji,jj,64)   !! Fast detritus production, C
2089                     trc2d(ji,jj,211) = trc2d(ji,jj,67)   !! Fast detritus remineralisation, C
2090                     trc2d(ji,jj,212) = trc2d(ji,jj,150)  !! Community respiration
2091                     trc2d(ji,jj,213) = fslownflux(ji,jj) !! Slow detritus N flux at 150 m
2092                     trc2d(ji,jj,214) = fslowcflux(ji,jj) !! Slow detritus C flux at 150 m
2093                     trc2d(ji,jj,215) = ffastn(ji,jj)     !! Fast detritus N flux at 150 m
2094                     trc2d(ji,jj,216) = ffastc(ji,jj)     !! Fast detritus C flux at 150 m
2095                  endif
2096                  !!
2097                  !! Jpalm (11-08-2014)
2098                  !! Add UKESM1 diagnoatics
2099                  !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2100                  if ((jk .eq. 1) .and.( jdms.eq.1)) then
2101                     trc2d(ji,jj,221) = dms_surf(ji,jj)          !! DMS surface concentration
2102                     !! AXY (13/03/15): add in other DMS estimates
2103                     trc2d(ji,jj,222) = dms_andr(ji,jj)          !! DMS surface concentration
2104                     trc2d(ji,jj,223) = dms_simo(ji,jj)          !! DMS surface concentration
2105                     trc2d(ji,jj,224) = dms_aran(ji,jj)          !! DMS surface concentration
2106                     trc2d(ji,jj,225) = dms_hall(ji,jj)          !! DMS surface concentration
2107                  endif
2108# endif
2109                  !! other possible future diagnostics include:
2110                  !!   - integrated tracer values (esp. biological)
2111                  !!   - mixed layer tracer values
2112                  !!   - sub-surface chlorophyll maxima (plus depth)
2113                  !!   - different mixed layer depth criteria (T, sigma, var. sigma)
2114
2115                  !!----------------------------------------------------------------------
2116                  !! Prepare 3D diagnostics
2117                  !!----------------------------------------------------------------------
2118                  !!
2119                  trc3d(ji,jj,jk,1)  = ((fprn(ji,jj) + fprd(ji,jj)) * zphn(ji,jj))     !! primary production 
2120                  trc3d(ji,jj,jk,2)  = fslownflux(ji,jj) + ffastn(ji,jj) !! detrital flux
2121                  trc3d(ji,jj,jk,3)  = fregen(ji,jj) + (freminn(ji,jj) * fse3t(ji,jj,jk))  !! remineralisation
2122# if defined key_roam
2123                  trc3d(ji,jj,jk,4)  = f3_pH(ji,jj,jk)            !! pH
2124                  trc3d(ji,jj,jk,5)  = f3_omcal(ji,jj,jk)         !! omega calcite
2125# else
2126                  trc3d(ji,jj,jk,4)  = ffastsi(ji,jj)             !! fast Si flux
2127# endif
2128             ENDIF   ! end of ln_diatrc option
2129             !! CLOSE wet point IF..THEN loop
2130            endif
2131         !! CLOSE horizontal loops
2132         ENDDO
2133         ENDDO
2134         !!
2135             IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
2136
2137                !!-------------------------------------------------------
2138                !! 2d specific k level diags
2139                !!-------------------------------------------------------
2140                CALL bio_medusa_diag_slice( jk )
2141
2142              ENDIF
2143      !! CLOSE vertical loop
2144      ENDDO
2145
2146      !!----------------------------------------------------------------------
2147      !! Final calculations for diagnostics
2148      !!----------------------------------------------------------------------
2149      CALL bio_medusa_fin( kt )
2150
2151# if defined key_trc_diabio
2152       !! Lateral boundary conditions on trcbio
2153       DO jn=1,jp_medusa_trd
2154          CALL lbc_lnk(trbio(:,:,1,jn),'T',1. )
2155       ENDDO 
2156# endif
2157
2158# if defined key_debug_medusa
2159       IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
2160       CALL flush(numout)
2161# endif
2162
2163   END SUBROUTINE trc_bio_medusa
2164
2165#else
2166   !!======================================================================
2167   !!  Dummy module :                                   No MEDUSA bio-model
2168   !!======================================================================
2169CONTAINS
2170   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
2171      INTEGER, INTENT( in ) ::   kt
2172      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
2173   END SUBROUTINE trc_bio_medusa
2174#endif 
2175
2176   !!======================================================================
2177END MODULE  trcbio_medusa
Note: See TracBrowser for help on using the repository browser.