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 @ 7978

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

Moved the iron chemistry and scavenging out of trcbio_medusa.F90

File size: 160.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 bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
105      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
106
107      IMPLICIT NONE
108      PRIVATE
109     
110      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
111
112   !!* Substitution
113#  include "domzgr_substitute.h90"
114   !!----------------------------------------------------------------------
115   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
116   !! $Id$
117   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
118   !!----------------------------------------------------------------------
119
120CONTAINS
121
122   SUBROUTINE trc_bio_medusa( kt )
123      !!---------------------------------------------------------------------
124      !!                     ***  ROUTINE trc_bio  ***
125      !!
126      !! ** Purpose :   compute the now trend due to biogeochemical processes
127      !!              and add it to the general trend of passive tracers equations
128      !!
129      !! ** Method  :   each now biological flux is calculated in function of now
130      !!              concentrations of tracers.
131      !!              depending on the tracer, these fluxes are sources or sinks.
132      !!              the total of the sources and sinks for each tracer
133      !!              is added to the general trend.
134      !!       
135      !!                      tra = tra + zf...tra - zftra...
136      !!                                     |         |
137      !!                                     |         |
138      !!                                  source      sink
139      !!       
140      !!              IF 'key_trc_diabio' defined , the biogeochemical trends
141      !!              for passive tracers are saved for futher diagnostics.
142      !!---------------------------------------------------------------------
143      !!
144      !!
145      !!----------------------------------------------------------------------           
146      !! Variable conventions
147      !!----------------------------------------------------------------------
148      !!
149      !! names: z*** - state variable
150      !!        f*** - function (or temporary variable used in part of a function)
151      !!        x*** - parameter
152      !!        b*** - right-hand part (sources and sinks)
153      !!        i*** - integer variable (usually used in yes/no flags)
154      !!
155      !! time (integer timestep)
156      INTEGER, INTENT( in ) ::    kt
157      !!
158      !! spatial array indices
159      INTEGER  ::    ji,jj,jk,jn
160      !!
161      !! AXY (27/07/10): add in indices for depth horizons (for sinking flux
162      !!                 and seafloor iron inputs)
163      !! INTEGER  ::    i0100, i0200, i0500, i1000, i1100
164      !!
165      !! model state variables
166!      REAL(wp), DIMENSION(jpi,jpj) ::    zchn,zchd,zphn,zphd,zpds,zzmi
167!      REAL(wp), DIMENSION(jpi,jpj) ::    zzme,zdet,zdtc,zdin,zsil,zfer
168! zage doesn't seem to be used - marc 19/4/17
169!      REAL(wp) ::    zage
170!# if defined key_roam
171!      REAL(wp), DIMENSION(jpi,jpj) ::    zdic, zalk, zoxy
172!      REAL(wp), DIMENSION(jpi,jpj) ::    ztmp, zsal
173!# endif
174!# if defined key_mocsy
175!      REAL(wp), DIMENSION(jpi,jpj) ::    zpho
176!# endif
177      !!
178      !! integrated source and sink terms
179      REAL(wp) ::    b0
180      !! AXY (23/08/13): changed from individual variables for each flux to
181      !!                 an array that holds all fluxes
182      REAL(wp), DIMENSION(jpi,jpj,jp_medusa) ::    btra
183      !!
184      !! primary production and chl related quantities     
185!      REAL(wp), DIMENSION(jpi,jpj) ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
186!      REAL(wp), DIMENSION(jpi,jpj) ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
187      !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range)
188!      REAL(wp), DIMENSION(jpi,jpj) ::    fjlim_pn, fjlim_pd
189      !! AXY (03/02/11): add in Liebig terms
190!      REAL(wp), DIMENSION(jpi,jpj) ::    fpnlim, fpdlim
191      !! AXY (16/07/09): add in Eppley curve functionality
192!      REAL(wp), DIMENSION(jpi,jpj) ::    fun_T,xvpnT,xvpdT
193      INTEGER  ::    ieppley
194      !! AXY (16/05/11): per Katya's prompting, add in new T-dependence
195      !!                 for phytoplankton growth only (i.e. no change
196      !!                 for remineralisation)
197!      REAL(wp), DIMENSION(jpi,jpj) ::    fun_Q10
198      !! AXY (01/03/10): add in mixed layer PP diagnostics
199!      REAL(wp), DIMENSION(jpi,jpj) ::    fprn_ml,fprd_ml
200      !!
201      !! nutrient limiting factors
202!      REAL(wp), DIMENSION(jpi,jpj) ::    fnln,ffln2            !! N and Fe
203!      REAL(wp), DIMENSION(jpi,jpj) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si
204      !!
205      !! silicon cycle
206!      REAL(wp), DIMENSION(jpi,jpj) ::    fsin,fnsi,fprds,fsdiss
207      REAL(wp)                     ::    fsin1,fnsi1,fnsi2
208      !!
209      !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
210!      REAL(wp), DIMENSION(jpi,jpj) ::    ffetop,ffebot,ffescav
211      REAL(wp) ::    xLgF, xFeT, xFeF, xFeL         !! state variables for iron-ligand system
212!      REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system
213      REAL(wp) ::    xb_coef_tmp, xb2M4ac           !! iron-ligand parameters
214      REAL(wp) ::    xmaxFeF,fdeltaFe               !! max Fe' parameters
215      !!
216      !! local parameters for Moore et al. (2004) alternative scavenging scheme
217      REAL(wp) ::    fbase_scav,fscal_sink,fscal_part,fscal_scav
218      !!
219      !! local parameters for Moore et al. (2008) alternative scavenging scheme
220      REAL(wp) ::    fscal_csink,fscal_sisink,fscal_casink
221      !!
222      !! local parameters for Galbraith et al. (2010) alternative scavenging scheme
223      REAL(wp) ::    xCscav1, xCscav2, xk_org, xORGscav  !! organic portion of scavenging
224      REAL(wp) ::    xk_inorg, xINORGscav                !! inorganic portion of scavenging
225      !!
226      !! microzooplankton grazing
227!      REAL(wp), DIMENSION(jpi,jpj) ::    fmi1,fmi,fgmipn,fgmid,fgmidc
228!      REAL(wp), DIMENSION(jpi,jpj) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
229      !!
230      !! mesozooplankton grazing
231!      REAL(wp), DIMENSION(jpi,jpj) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
232!      REAL(wp), DIMENSION(jpi,jpj) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
233      !!
234      !! mortality/Remineralisation (defunct parameter "fz" removed)
235!      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
236# if defined key_roam
237!      REAL(wp), DIMENSION(jpi,jpj) ::    fddc
238# endif
239!      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
240      REAL(wp), DIMENSION(jpi,jpj) ::    fslown, fslowc
241!      REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux
242      REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi
243!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi
244# if defined key_roam
245!! Doesn't look like this is used - marc 10/4/17
246!!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenc
247!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfastc
248# endif
249      !!
250      !! particle flux
251!      REAL(WP), DIMENSION(jpi,jpj) ::    fdep1,fcaco3
252!      REAL(WP), DIMENSION(jpi,jpj) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
253!      REAL(wp), DIMENSION(jpi,jpj) ::    freminn,freminsi,freminfe,freminc,freminca
254!      REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca
255!      REAL(wp), DIMENSION(jpi,jpj) ::    fprotf
256!      REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca
257!      REAL(wp), DIMENSION(jpi,jpj) ::    fccd
258!      REAL(wp), DIMENSION(jpi,jpj) ::    fccd_dep
259      !!
260      !! AXY (06/07/11): alternative fast detritus schemes
261      REAL(wp) ::    fb_val, fl_sst
262      !!
263      !! AXY (08/07/11): fate of fast detritus reaching the seafloor
264! I don't think ffast2slowfe is used - marc 10/4/17
265!      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowfe,ffast2slowc
266      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowc
267      !!
268      !! conservation law
269      REAL(wp) ::    fnit0,fsil0,ffer0 
270# if defined key_roam
271      REAL(wp) ::    fcar0,falk0,foxy0 
272# endif     
273      !!
274      !! temporary variables
275      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8,fq9
276      !!
277      !! water column nutrient and flux integrals
278!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_n,ftot_si,ftot_fe
279!      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_n,fflx_si,fflx_fe
280!      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_n,fifd_si,fifd_fe
281!      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_n,fofd_si,fofd_fe
282# if defined key_roam
283!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_c,ftot_a,ftot_o2
284!      REAL(wp), DIMENSION(jpi,jpj) ::    fflx_c,fflx_a,fflx_o2
285!      REAL(wp), DIMENSION(jpi,jpj) ::    fifd_c,fifd_a,fifd_o2
286!      REAL(wp), DIMENSION(jpi,jpj) ::    fofd_c,fofd_a,fofd_o2
287# endif
288      !!
289      !! zooplankton grazing integrals
290!      REAL(wp), DIMENSION(jpi,jpj) ::    fzmi_i,fzmi_o,fzme_i,fzme_o
291      !!
292      !! limitation term temporary variables
293!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_pn,ftot_pd
294!      REAL(wp), DIMENSION(jpi,jpj) ::    ftot_zmi,ftot_zme,ftot_det,ftot_dtc
295      !! use ballast scheme (1) or simple exponential scheme (0; a conservation test)
296      INTEGER  ::    iball
297      !! use biological fluxes (1) or not (0)
298      INTEGER  ::    ibio_switch
299      !!
300      !! diagnose fluxes (should only be used in 1D runs)
301!      INTEGER  ::    idf, idfval
302      !!
303      !! nitrogen and silicon production and consumption
304      REAL(wp) ::    fn_prod, fn_cons, fs_prod, fs_cons
305!      REAL(wp), DIMENSION(jpi,jpj) ::    fnit_prod, fnit_cons, fsil_prod, fsil_cons
306# if defined key_roam
307      !!
308      !! flags to help with calculating the position of the CCD
309! Moved into carb_chem.F90 - marc 20/4/17
310!      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
311      !!
312      !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
313!      REAL(wp)                     ::    f_xco2a
314!      REAL(wp), DIMENSION(jpi,jpj) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
315!      REAL(wp), DIMENSION(jpi,jpj) ::    f_TDIC, f_TALK, f_dcf, f_henry
316!      REAL(wp), DIMENSION(jpi,jpj) ::    f_pp0
317!      REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat
318      REAL(wp)                     ::    f_o2sat3
319!      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg
320      !!
321      !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
322!      REAL(wp), DIMENSION(jpi,jpj) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm
323!      REAL(wp), DIMENSION(jpi,jpj) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2
324      !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
325      REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol
326      REAL, PARAMETER :: secs_in_day    = 86400.0  !! s / d
327      REAL, PARAMETER :: CO2flux_conv   = (1.e-6 * weight_CO2_mol) / secs_in_day
328
329      !!
330!      INTEGER, DIMENSION(jpi,jpj)  ::    iters
331      REAL(wp) ::    f_year
332      INTEGER  ::    i_year
333      INTEGER  ::    iyr1, iyr2
334      !!
335      !! carbon, alkalinity production and consumption
336      REAL(wp) ::    fc_prod, fc_cons, fa_prod, fa_cons
337!      REAL(wp), DIMENSION(jpi,jpj) ::    fcomm_resp
338!      REAL(wp), DIMENSION(jpi,jpj) ::    fcar_prod, fcar_cons
339      !!
340      !! oxygen production and consumption (and non-consumption)
341      REAL(wp), DIMENSION(jpi,jpj) ::    fo2_prod, fo2_cons, fo2_ncons, fo2_ccons
342!      REAL(wp), DIMENSION(jpi,jpj) ::    foxy_prod, foxy_cons, foxy_anox
343      !! Jpalm (11-08-2014)
344      !! add DMS in MEDUSA for UKESM1 model
345!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf
346      !! AXY (13/03/15): add in other DMS calculations
347!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall
348
349# endif
350      !!
351      !! benthic fluxes
352!      INTEGER  ::    ibenthic
353!      REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe,              f_sbenin_c
354!      REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca
355!      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca
356      REAL(wp) ::    zfact
357      !!
358      !! benthic fluxes of CaCO3 that shouldn't happen because of lysocline
359!      REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca
360      !!
361      !! riverine fluxes
362!      REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk
363      !! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface
364      REAL(wp), DIMENSION(jpi,jpj) :: f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk
365      !!---------------------------------------------------------------------
366
367# if defined key_debug_medusa
368      IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined'
369      CALL flush(numout)
370# endif 
371
372      !! AXY (20/11/14): alter this to report on first MEDUSA call
373      !! IF( kt == nit000 ) THEN
374      IF( kt == nittrc000 ) THEN
375         IF(lwp) WRITE(numout,*)
376         IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model'
377         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
378    IF(lwp) WRITE(numout,*) ' kt =',kt
379      ENDIF
380
381      !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
382      ibenthic = 1
383
384      !! not sure what this is for; it's not used anywhere; commenting out
385      !! fbodn(:,:) = 0.e0   
386
387      !!
388      IF( ln_diatrc ) THEN
389         !! blank 2D diagnostic array
390         trc2d(:,:,:) = 0.e0
391         !!
392         !! blank 3D diagnostic array
393         trc3d(:,:,:,:) = 0.e0
394      ENDIF
395
396      !!----------------------------------------------------------------------
397      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency
398      !! terms of all biological equations to 0.
399      !!----------------------------------------------------------------------
400      !!
401      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit
402      !!                 the bill for now; another item on the things-to-sort-
403      !!     out-in-the-future list ...
404# if defined key_kill_medusa
405      b0 = 0.
406# else
407      b0 = 1.
408# endif
409      !!----------------------------------------------------------------------
410      !! fast detritus ballast scheme (0 = no; 1 = yes)
411      !! alternative to ballast scheme is same scheme but with no ballast
412      !! protection (not dissimilar to Martin et al., 1987)
413      !!----------------------------------------------------------------------
414      !!
415      iball = 1
416
417      !!----------------------------------------------------------------------
418      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
419      !! these should *only* be used in 1D since they give comprehensive
420      !! output for ecological functions in the model; primarily used in
421      !! debugging
422      !!----------------------------------------------------------------------
423      !!
424      idf    = 0
425      !!
426      !! timer mechanism
427      if (kt/120*120.eq.kt) then
428         idfval = 1
429      else
430         idfval = 0
431      endif
432
433      !!----------------------------------------------------------------------
434      !! Initialise arrays to zero and set up arrays for diagnostics
435      !!----------------------------------------------------------------------
436! tmp - marc
437      write(numout,*) 'bbb13. before call to bio_medusa_init,kt=',kt
438      flush(numout)
439!
440      CALL bio_medusa_init( kt )
441! tmp - marc
442      write(numout,*) 'bbb14. after call to bio_medusa_init,kt=',kt
443      flush(numout)
444!
445       !!
446# if defined key_axy_nancheck
447       DO jn = 1,jptra
448          !! fq0 = MINVAL(trn(:,:,:,jn))
449          !! fq1 = MAXVAL(trn(:,:,:,jn))
450          fq2 = SUM(trn(:,:,:,jn))
451          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', &
452          !! &        kt, jn, fq0, fq1, fq2
453          !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR
454          !!                 and has been replaced here with a specialist routine
455          !! if (fq2 /= fq2 ) then
456          if ( ieee_is_nan( fq2 ) ) then
457             !! there's a NaN here
458             if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:'
459             DO jk = 1,jpk
460                DO jj = 1,jpj
461                   DO ji = 1,jpi
462                      !! AXY (30/01/14): "isnan" problem on HECTOR
463                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then
464                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
465                         if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', &
466                         &        tmask(ji,jj,jk), ji, jj, jk, jn
467                      endif
468                   enddo
469                enddo
470             enddo
471             CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
472          endif
473       ENDDO
474       CALL flush(numout)
475# endif
476
477# if defined key_debug_medusa
478      IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked'
479      CALL flush(numout)
480# endif 
481
482# if defined key_roam
483      !!----------------------------------------------------------------------
484      !! calculate atmospheric pCO2
485      !!----------------------------------------------------------------------
486      !!
487      !! what's atmospheric pCO2 doing? (data start in 1859)
488      iyr1  = nyear - 1859 + 1
489      iyr2  = iyr1 + 1
490      if (iyr1 .le. 1) then
491         !! before 1860
492         f_xco2a(:,:) = hist_pco2(1)
493      elseif (iyr2 .ge. 242) then
494         !! after 2099
495         f_xco2a(:,:) = hist_pco2(242)
496      else
497         !! just right
498         fq0 = hist_pco2(iyr1)
499         fq1 = hist_pco2(iyr2)
500         fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0)
501         !! AXY (14/06/12): tweaked to make more sense (and be correct)
502#  if defined key_bs_axy_yrlen
503         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  !! bugfix: for 360d year with HadGEM2-ES forcing
504#  else
505         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0  !! original use of 365 days (not accounting for leap year or 360d year)
506#  endif
507         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
508         f_xco2a(:,:) = fq4
509      endif
510#  if defined key_axy_pi_co2
511      f_xco2a(:,:) = 284.725          !! OCMIP pre-industrial pCO2
512#  endif
513      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
514      !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day)
515      !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)
516      !! AXY (29/01/14): remove surplus diagnostics
517      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0
518      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1
519      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
520      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
521      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
522# endif
523
524# if defined key_debug_medusa
525      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry'
526      IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt
527      IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000
528      CALL flush(numout)
529# endif 
530
531# if defined key_roam
532      !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every
533      !!                 month (this is hardwired as 960 timesteps but should
534      !!                 be calculated and done properly
535      !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN
536      !! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN
537      !!=============================
538      !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call :
539      !!          we don't want to call on the first time-step of all run submission,
540      !!          but only on the very first time-step, and then every month
541      !!          So we call on nittrc000 if not restarted run,
542      !!          else if one month after last call.
543      !!          assume one month is 30d --> 3600*24*30 : 2592000s
544      !!          try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt)   
545      !!          ++ need to pass carb-chem output var through restarts
546      If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN
547         !!----------------------------------------------------------------------
548         !! Calculate the carbonate chemistry for the whole ocean on the first
549         !! simulation timestep and every month subsequently; the resulting 3D
550         !! field of omega calcite is used to determine the depth of the CCD
551         !!----------------------------------------------------------------------
552         CALL carb_chem( kt )
553
554      ENDIF
555# endif
556
557# if defined key_debug_medusa
558      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations'
559      CALL flush(numout)
560# endif 
561
562      !!----------------------------------------------------------------------
563      !! MEDUSA has unified equation through the water column
564      !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)
565      !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1         
566      !!----------------------------------------------------------------------
567      !!
568      !! NOTE: the ordering of the loops below differs from that of some other
569      !! models; looping over the vertical dimension is the outermost loop and
570      !! this complicates some calculations (e.g. storage of vertical fluxes
571      !! that can otherwise be done via a singular variable require 2D fields
572      !! here); however, these issues are relatively easily resolved, but the
573      !! loops CANNOT be reordered without potentially causing code efficiency
574      !! problems (e.g. array indexing means that reordering the loops would
575      !! require skipping between widely-spaced memory location; potentially
576      !! outside those immediately cached)
577      !!
578      !! OPEN vertical loop
579      DO jk = 1,jpk
580         !! OPEN horizontal loops
581         DO jj = 2,jpjm1
582         DO ji = 2,jpim1
583            !! OPEN wet point IF..THEN loop
584            if (tmask(ji,jj,jk) == 1) then               
585               !!===========================================================
586               !! SETUP LOCAL GRID CELL
587               !!===========================================================
588               !!
589               !!-----------------------------------------------------------
590               !! Some notes on grid vertical structure
591               !! - fsdepw(ji,jj,jk) is the depth of the upper surface of
592               !!   level jk
593               !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of
594               !!   level jk
595               !! - fse3t(ji,jj,jk)  is the thickness of level jk
596               !!-----------------------------------------------------------
597               !!
598               !! AXY (01/03/10): set up level depth (bottom of level)
599               fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
600               !! AXY (28/11/16): local seafloor depth
601               !!                 previously mbathy(ji,jj) - 1, now
602               !!                 mbathy(ji,jj)
603               mbathy(ji,jj) = mbathy(ji,jj)
604               !!
605               !! set up model tracers
606               !! negative values of state variables are not allowed to
607               !! contribute to the calculated fluxes
608               !! non-diatom chlorophyll
609               zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
610               !! diatom chlorophyll
611               zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
612               !! non-diatoms
613               zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
614               !! diatoms
615               zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
616               !! diatom silicon
617               zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
618               !! AXY (28/01/10): probably need to take account of
619               !! chl/biomass connection
620               if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
621               if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
622               if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
623               if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
624          !! AXY (23/01/14): duh - why did I forget diatom silicon?
625          if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
626          if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
627               !! microzooplankton
628               zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
629               !! mesozooplankton
630               zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
631               !! detrital nitrogen
632               zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
633               !! dissolved inorganic nitrogen
634               zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
635               !! dissolved silicic acid
636               zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
637               !! dissolved "iron"
638               zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
639# if defined key_roam
640               !! detrital carbon
641               zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
642               !! dissolved inorganic carbon
643               zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
644               !! alkalinity
645               zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
646               !! oxygen
647               zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
648#  if defined key_axy_carbchem && defined key_mocsy
649               !! phosphate via DIN and Redfield
650               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
651#  endif
652               !!
653               !! also need physical parameters for gas exchange calculations
654               ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
655               zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
656               !!
657          !! AXY (28/02/14): check input fields
658               if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
659                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ',  &
660                     tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',     &
661                     ji, ',', jj, ',', jk, ') at time', kt
662        IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',&
663                     tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
664                  !! temperatur
665                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
666               endif
667               if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
668                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', &
669                     tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    &
670                     ji, ',', jj, ',', jk, ') at time', kt
671               endif
672# else
673               !! implicit detrital carbon
674               zdtc(ji,jj) = zdet(ji,jj) * xthetad
675# endif
676# if defined key_debug_medusa
677               if (idf.eq.1) then
678                  !! AXY (15/01/10)
679                  if (trn(ji,jj,jk,jpdin).lt.0.) then
680                     IF (lwp) write (numout,*) '------------------------------'
681                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',        &
682                        trn(ji,jj,jk,jpdin)
683                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',        &
684                        ji, jj, jk, kt
685                  endif
686                  if (trn(ji,jj,jk,jpsil).lt.0.) then
687                     IF (lwp) write (numout,*) '------------------------------'
688                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',        &
689                        trn(ji,jj,jk,jpsil)
690                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',        &
691                        ji, jj, jk, kt
692                  endif
693#  if defined key_roam
694                  if (trn(ji,jj,jk,jpdic).lt.0.) then
695                     IF (lwp) write (numout,*) '------------------------------'
696                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',        &
697                        trn(ji,jj,jk,jpdic)
698                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',        &
699                        ji, jj, jk, kt
700                  endif
701                  if (trn(ji,jj,jk,jpalk).lt.0.) then
702                     IF (lwp) write (numout,*) '------------------------------'
703                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',        &
704                        trn(ji,jj,jk,jpalk)
705                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',        &
706                        ji, jj, jk, kt
707                  endif
708                  if (trn(ji,jj,jk,jpoxy).lt.0.) then
709                     IF (lwp) write (numout,*) '------------------------------'
710                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',        &
711                        trn(ji,jj,jk,jpoxy)
712                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',        &
713                        ji, jj, jk, kt
714                  endif
715#  endif
716               endif
717# endif
718# if defined key_debug_medusa
719               !! report state variable values
720               if (idf.eq.1.AND.idfval.eq.1) then
721                  IF (lwp) write (numout,*) '------------------------------'
722                  IF (lwp) write (numout,*) 'fthk(',jk,') = ', fse3t(ji,jj,jk)
723                  IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
724                  IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
725                  IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
726                  IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
727                  IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
728                  IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
729                  IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
730                  IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
731                  IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
732#  if defined key_roam
733                  IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
734                  IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
735                  IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
736                  IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)                 
737#  endif
738               endif
739# endif
740
741# if defined key_debug_medusa
742               if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
743                  IF (lwp) write (numout,*) '------------------------------'
744                  IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
745               endif
746# endif
747
748               !! sum tracers for inventory checks
749               IF( lk_iomput ) THEN
750                  IF ( med_diag%INVTN%dgsave )   THEN
751                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        &
752                        (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) +     &
753                                             zzmi(ji,jj) + zzme(ji,jj) +     &
754                                             zdet(ji,jj) + zdin(ji,jj) ) )
755                  ENDIF
756                  IF ( med_diag%INVTSI%dgsave )  THEN
757                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       & 
758                        (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) )
759                  ENDIF
760                  IF ( med_diag%INVTFE%dgsave )  THEN
761                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       & 
762                        (fse3t(ji,jj,jk) * ( xrfn *                          &
763                                            ( zphn(ji,jj) + zphd(ji,jj) +    &
764                                              zzmi(ji,jj) + zzme(ji,jj) +    &
765                                              zdet(ji,jj) ) + zfer(ji,jj) ) )
766                  ENDIF
767# if defined key_roam
768                  IF ( med_diag%INVTC%dgsave )  THEN
769                     ftot_c(ji,jj)  = ftot_c(ji,jj) + & 
770                        (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) +      &
771                                             (xthetapd * zphd(ji,jj)) +      &
772                                             (xthetazmi * zzmi(ji,jj)) +     &
773                                             (xthetazme * zzme(ji,jj)) +     &
774                                             zdtc(ji,jj) + zdic(ji,jj) ) )
775                  ENDIF
776                  IF ( med_diag%INVTALK%dgsave ) THEN
777                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     &
778                                                       zalk(ji,jj))
779                  ENDIF
780                  IF ( med_diag%INVTO2%dgsave )  THEN
781                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    &
782                                                        zoxy(ji,jj))
783                  ENDIF
784                  !!
785                  !! AXY (10/11/16): CMIP6 diagnostics
786                  IF ( med_diag%INTDISSIC%dgsave ) THEN
787                     intdissic(ji,jj) = intdissic(ji,jj) +                   &
788                                        (fse3t(ji,jj,jk) * zdic(ji,jj))
789                  ENDIF
790                  IF ( med_diag%INTDISSIN%dgsave ) THEN
791                     intdissin(ji,jj) = intdissin(ji,jj) +                   &
792                                        (fse3t(ji,jj,jk) * zdin(ji,jj))
793                  ENDIF
794                  IF ( med_diag%INTDISSISI%dgsave ) THEN
795                     intdissisi(ji,jj) = intdissisi(ji,jj) +                 &
796                                         (fse3t(ji,jj,jk) * zsil(ji,jj))
797                  ENDIF
798                  IF ( med_diag%INTTALK%dgsave ) THEN
799                     inttalk(ji,jj) = inttalk(ji,jj) +                       &
800                                      (fse3t(ji,jj,jk) * zalk(ji,jj))
801                  ENDIF
802                  IF ( med_diag%O2min%dgsave ) THEN
803                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then
804                        o2min(ji,jj)  = zoxy(ji,jj)
805                        IF ( med_diag%ZO2min%dgsave ) THEN
806                           !! layer midpoint
807                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               &
808                                            fdep1(ji,jj)) / 2.0
809                        ENDIF
810                     endif
811                  ENDIF
812# endif
813               ENDIF
814
815               CALL flush(numout)
816
817            ENDIF
818         ENDDO
819         ENDDO
820
821         !!----------------------------------------------------------------
822         !! Calculate air-sea gas exchange and river inputs
823         !!----------------------------------------------------------------
824         IF ( jk == 1 ) THEN
825            call air_sea( kt )
826         END IF
827
828# if defined key_roam
829
830! Moved from above - marc 21/4/17
831! I think this should be moved into diagnostics at bottom - it
832! doesn't like it is used anyway else - marc 21/4/17
833         DO jj = 2,jpjm1
834         DO ji = 2,jpim1
835            !! OPEN wet point IF..THEN loop
836            if (tmask(ji,jj,jk) == 1) then
837
838               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
839               IF ( med_diag%O2SAT3%dgsave ) THEN
840                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 )
841                  o2sat3(ji, jj, jk) = f_o2sat3
842               ENDIF
843            ENDIF
844         ENDDO
845         ENDDO
846# endif
847
848         !!------------------------------------------------------------------
849         !! Phytoplankton growth, zooplankton grazing and miscellaneous
850         !! plankton losses.
851         !!------------------------------------------------------------------
852         CALL plankton( jk )
853
854! Detritus process - marc
855
856         DO jj = 2,jpjm1
857         DO ji = 2,jpim1
858            !! OPEN wet point IF..THEN loop
859            if (tmask(ji,jj,jk) == 1) then
860
861               !!----------------------------------------------------------------------
862               !! Detritus remineralisation
863               !! Constant or temperature-dependent
864               !!----------------------------------------------------------------------
865               !!
866               if (jmd.eq.1) then
867                  !! temperature-dependent
868                  fdd(ji,jj)  = xmd  * fun_T(ji,jj) * zdet(ji,jj)
869# if defined key_roam
870                  fddc(ji,jj) = xmdc * fun_T(ji,jj) * zdtc(ji,jj)
871# endif
872               elseif (jmd.eq.2) then
873                  !! AXY (16/05/13): add in Q10-based parameterisation (def in nmlst)
874                  !! temperature-dependent
875                  fdd(ji,jj)  = xmd  * fun_Q10(ji,jj) * zdet(ji,jj)
876#if defined key_roam
877                  fddc(ji,jj) = xmdc * fun_Q10(ji,jj) * zdtc(ji,jj)
878#endif
879               else
880                  !! temperature-independent
881                  fdd(ji,jj)  = xmd  * zdet(ji,jj)
882# if defined key_roam
883                  fddc(ji,jj) = xmdc * zdtc(ji,jj)
884# endif
885               endif
886               !!
887               !! AXY (22/07/09): accelerate detrital remineralisation in the bottom box
888               if ((jk.eq.mbathy(ji,jj)) .and. jsfd.eq.1) then
889                  fdd(ji,jj)  = 1.0  * zdet(ji,jj)
890# if defined key_roam
891                  fddc(ji,jj) = 1.0  * zdtc(ji,jj)
892# endif
893               endif
894               
895# if defined key_debug_medusa
896               !! report plankton mortality and remineralisation
897               if (idf.eq.1.AND.idfval.eq.1) then
898                  IF (lwp) write (numout,*) '------------------------------'
899                  IF (lwp) write (numout,*) 'fdpn2(',jk,') = ', fdpn2(ji,jj)
900                  IF (lwp) write (numout,*) 'fdpd2(',jk,') = ', fdpd2(ji,jj)
901                  IF (lwp) write (numout,*) 'fdpds2(',jk,')= ', fdpds2(ji,jj)
902                  IF (lwp) write (numout,*) 'fdzmi2(',jk,')= ', fdzmi2(ji,jj)
903                  IF (lwp) write (numout,*) 'fdzme2(',jk,')= ', fdzme2(ji,jj)
904                  IF (lwp) write (numout,*) 'fdpn(',jk,')  = ', fdpn(ji,jj)
905                  IF (lwp) write (numout,*) 'fdpd(',jk,')  = ', fdpd(ji,jj)
906                  IF (lwp) write (numout,*) 'fdpds(',jk,') = ', fdpds(ji,jj)
907                  IF (lwp) write (numout,*) 'fdzmi(',jk,') = ', fdzmi(ji,jj)
908                  IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme(ji,jj)
909                  IF (lwp) write (numout,*) 'fdd(',jk,')   = ', fdd(ji,jj)
910#  if defined key_roam
911                  IF (lwp) write (numout,*) 'fddc(',jk,')  = ', fddc(ji,jj)
912#  endif
913               endif
914# endif
915
916               !!----------------------------------------------------------------------
917               !! Detritus addition to benthos
918               !! If activated, slow detritus in the bottom box will enter the
919               !! benthic pool
920               !!----------------------------------------------------------------------
921               !!
922               if ((jk.eq.mbathy(ji,jj)) .and. jorgben.eq.1) then
923                  !! this is the BOTTOM OCEAN BOX -> into the benthic pool!
924                  !!
925                  f_sbenin_n(ji,jj)  = (zdet(ji,jj) * vsed * 86400.)
926                  f_sbenin_fe(ji,jj) = (zdet(ji,jj) * vsed * 86400. * xrfn)
927# if defined key_roam
928                  f_sbenin_c(ji,jj)  = (zdtc(ji,jj) * vsed * 86400.)
929# else
930                  f_sbenin_c(ji,jj)  = (zdet(ji,jj) * vsed * 86400. * xthetad)
931# endif
932               endif
933
934! MAYBE BUT A BREAK IN HERE, IRON CHEMISTRY AND SCAVENGING - marc 20/4/17
935! (detritus processes is 74 lines)
936            ENDIF
937         ENDDO
938         ENDDO
939
940         !!----------------------------------------------------------------
941         !! Iron chemistry and scavenging
942         !!----------------------------------------------------------------
943         CALL iron_chem_scav( jk )
944
945! Miscellaneous processes - marc
946
947         DO jj = 2,jpjm1
948         DO ji = 2,jpim1
949            !! OPEN wet point IF..THEN loop
950            if (tmask(ji,jj,jk) == 1) then
951
952               !!----------------------------------------------------------------------
953               !! Miscellaneous
954               !!----------------------------------------------------------------------
955               !!
956               !! diatom frustule dissolution
957               fsdiss(ji,jj)  = xsdiss * zpds(ji,jj)
958
959# if defined key_debug_medusa
960               !! report miscellaneous calculations
961               if (idf.eq.1.AND.idfval.eq.1) then
962                  IF (lwp) write (numout,*) '------------------------------'
963                  IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss(ji,jj)
964               endif
965# endif
966
967               !!----------------------------------------------------------------------
968               !! Slow detritus creation
969               !!----------------------------------------------------------------------
970               !! this variable integrates the creation of slow sinking detritus
971               !! to allow the split between fast and slow detritus to be
972               !! diagnosed
973               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj)) + &
974               ((1.0 - xfdfrac2) * fdzme(ji,jj)) + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))
975               !!
976               !! this variable records the slow detrital sinking flux at this
977               !! particular depth; it is used in the output of this flux at
978               !! standard depths in the diagnostic outputs; needs to be
979               !! adjusted from per second to per day because of parameter vsed
980               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400.
981# if defined key_roam
982               !!
983               !! and the same for detrital carbon
984               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) + (xthetazmi * fdzmi(ji,jj)) + &
985               (xthetapd * (1.0 - xfdfrac1) * fdpd(ji,jj)) + &
986               (xthetazme * (1.0 - xfdfrac2) * fdzme(ji,jj)) + &
987               ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))
988               !!
989               !! this variable records the slow detrital sinking flux at this
990               !! particular depth; it is used in the output of this flux at
991               !! standard depths in the diagnostic outputs; needs to be
992               !! adjusted from per second to per day because of parameter vsed
993               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400.
994# endif
995
996               !!----------------------------------------------------------------------
997               !! Nutrient regeneration
998               !! this variable integrates total nitrogen regeneration down the
999               !! watercolumn; its value is stored and output as a 2D diagnostic;
1000               !! the corresponding dissolution flux of silicon (from sources
1001               !! other than fast detritus) is also integrated; note that,
1002               !! confusingly, the linear loss terms from plankton compartments
1003               !! are labelled as fdX2 when one might have expected fdX or fdX1
1004               !!----------------------------------------------------------------------
1005               !!
1006               !! nitrogen
1007               fregen(ji,jj)   = (( (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) +                &  ! messy feeding
1008               (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) +           &  ! messy feeding
1009               fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) +                                &  ! excretion + D remin.
1010               fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)) * fse3t(ji,jj,jk))                    ! linear mortality
1011               !!
1012               !! silicon
1013               fregensi(ji,jj) = (( fsdiss(ji,jj) + ((1.0 - xfdfrac1) * fdpds(ji,jj)) +      &  ! dissolution + non-lin. mortality
1014               ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +                           &  ! egestion by zooplankton
1015               fdpds2(ji,jj)) * fse3t(ji,jj,jk))                                             ! linear mortality
1016# if defined key_roam
1017               !!
1018               !! carbon
1019! Doesn't look this is used - marc 10/4/17
1020!               fregenc(ji,jj)  = (( (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))) +  &  ! messy feeding
1021!               (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +     &  ! messy feeding
1022!               (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))) +                       &  ! messy feeding
1023!               fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) +                               &  ! respiration + D remin.
1024!               (xthetapn * fdpn2(ji,jj)) + (xthetapd * fdpd2(ji,jj)) +                &  ! linear mortality
1025!               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality
1026# endif
1027
1028! MAYBE BUT A BREAK IN HERE, FAST-SINKINIG DETRITUS - marc 20/4/17
1029! (miscellaneous processes is 79 lines)
1030            ENDIF
1031         ENDDO
1032         ENDDO
1033
1034         DO jj = 2,jpjm1
1035         DO ji = 2,jpim1
1036            !! OPEN wet point IF..THEN loop
1037            if (tmask(ji,jj,jk) == 1) then
1038
1039               !!----------------------------------------------------------------------
1040               !! Fast-sinking detritus terms
1041               !! "local" variables declared so that conservation can be checked;
1042               !! the calculated terms are added to the fast-sinking flux later on
1043               !! only after the flux entering this level has experienced some
1044               !! remineralisation
1045               !! note: these fluxes need to be scaled by the level thickness
1046               !!----------------------------------------------------------------------
1047               !!
1048               !! nitrogen:   diatom and mesozooplankton mortality
1049               ftempn(ji,jj)         = b0 * ((xfdfrac1 * fdpd(ji,jj))  + (xfdfrac2 * fdzme(ji,jj)))
1050               !!
1051               !! silicon:    diatom mortality and grazed diatoms
1052               ftempsi(ji,jj)        = b0 * ((xfdfrac1 * fdpds(ji,jj)) + (xfdfrac3 * fgmepds(ji,jj)))
1053               !!
1054               !! iron:       diatom and mesozooplankton mortality
1055               ftempfe(ji,jj)        = b0 * (((xfdfrac1 * fdpd(ji,jj)) + (xfdfrac2 * fdzme(ji,jj))) * xrfn)
1056               !!
1057               !! carbon:     diatom and mesozooplankton mortality
1058               ftempc(ji,jj)         = b0 * ((xfdfrac1 * xthetapd * fdpd(ji,jj)) + & 
1059                                (xfdfrac2 * xthetazme * fdzme(ji,jj)))
1060               !!
1061# if defined key_roam
1062               if (jrratio.eq.0) then
1063                  !! CaCO3:      latitudinally-based fraction of total primary production
1064                  !!               0.10 at equator; 0.02 at pole
1065                  fcaco3(ji,jj)         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - abs(gphit(ji,jj))) / 90.0))
1066               elseif (jrratio.eq.1) then
1067                  !! CaCO3:      Ridgwell et al. (2007) submodel, version 1
1068                  !!             this uses SURFACE omega calcite to regulate rain ratio
1069                  if (f_omcal(ji,jj).ge.1.0) then
1070                     fq1 = (f_omcal(ji,jj) - 1.0)**0.81
1071                  else
1072                     fq1 = 0.
1073                  endif
1074                  fcaco3(ji,jj) = xridg_r0 * fq1
1075               elseif (jrratio.eq.2) then
1076                  !! CaCO3:      Ridgwell et al. (2007) submodel, version 2
1077                  !!             this uses FULL 3D omega calcite to regulate rain ratio
1078                  if (f3_omcal(ji,jj,jk).ge.1.0) then
1079                     fq1 = (f3_omcal(ji,jj,jk) - 1.0)**0.81
1080                  else
1081                     fq1 = 0.
1082                  endif
1083                  fcaco3(ji,jj) = xridg_r0 * fq1
1084               endif
1085# else
1086               !! CaCO3:      latitudinally-based fraction of total primary production
1087               !!               0.10 at equator; 0.02 at pole
1088               fcaco3(ji,jj)         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - abs(gphit(ji,jj))) / 90.0))
1089# endif
1090               !! AXY (09/03/09): convert CaCO3 production from function of
1091               !! primary production into a function of fast-sinking material;
1092               !! technically, this is what Dunne et al. (2007) do anyway; they
1093               !! convert total primary production estimated from surface
1094               !! chlorophyll to an export flux for which they apply conversion
1095               !! factors to estimate the various elemental fractions (Si, Ca)
1096               ftempca(ji,jj)        = ftempc(ji,jj) * fcaco3(ji,jj)
1097
1098# if defined key_debug_medusa
1099               !! integrate total fast detritus production
1100               if (idf.eq.1) then
1101                  fifd_n(ji,jj)  = fifd_n(ji,jj)  + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
1102                  fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
1103                  fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))
1104#  if defined key_roam
1105                  fifd_c(ji,jj)  = fifd_c(ji,jj)  + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
1106#  endif
1107               endif
1108
1109               !! report quantities of fast-sinking detritus for each component
1110               if (idf.eq.1.AND.idfval.eq.1) then
1111                  IF (lwp) write (numout,*) '------------------------------'
1112                  IF (lwp) write (numout,*) 'fdpd(',jk,')    = ', fdpd(ji,jj)
1113                  IF (lwp) write (numout,*) 'fdzme(',jk,')   = ', fdzme(ji,jj)
1114                  IF (lwp) write (numout,*) 'ftempn(',jk,')  = ', ftempn(ji,jj)
1115                  IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi(ji,jj)
1116                  IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe(ji,jj)
1117                  IF (lwp) write (numout,*) 'ftempc(',jk,')  = ', ftempc(ji,jj)
1118                  IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca(ji,jj)
1119                  IF (lwp) write (numout,*) 'flat(',jk,')    = ', abs(gphit(ji,jj))
1120                  IF (lwp) write (numout,*) 'fcaco3(',jk,')  = ', fcaco3(ji,jj)
1121               endif
1122# endif
1123
1124               !!----------------------------------------------------------------------
1125               !! This version of MEDUSA offers a choice of three methods for
1126               !! handling the remineralisation of fast detritus.  All three
1127               !! do so in broadly the same way:
1128               !!
1129               !!   1.  Fast detritus is stored as a 2D array                   [ ffastX  ]
1130               !!   2.  Fast detritus is added level-by-level                   [ ftempX  ]
1131               !!   3.  Fast detritus is not remineralised in the top box       [ freminX ]
1132               !!   4.  Remaining fast detritus is remineralised in the bottom  [ fsedX   ]
1133               !!       box
1134               !!
1135               !! The three remineralisation methods are:
1136               !!   
1137               !!   1.  Ballast model (i.e. that published in Yool et al., 2011)
1138               !!  (1b. Ballast-sans-ballast model)
1139               !!   2.  Martin et al. (1987)
1140               !!   3.  Henson et al. (2011)
1141               !!
1142               !! The first of these couples C, N and Fe remineralisation to
1143               !! the remineralisation of particulate Si and CaCO3, but the
1144               !! latter two treat remineralisation of C, N, Fe, Si and CaCO3
1145               !! completely separately.  At present a switch within the code
1146               !! regulates which submodel is used, but this should be moved
1147               !! to the namelist file.
1148               !!
1149               !! The ballast-sans-ballast submodel is an original development
1150               !! feature of MEDUSA in which the ballast submodel's general
1151               !! framework and parameterisation is used, but in which there
1152               !! is no protection of organic material afforded by ballasting
1153               !! minerals.  While similar, it is not the same as the Martin
1154               !! et al. (1987) submodel.
1155               !!
1156               !! Since the three submodels behave the same in terms of
1157               !! accumulating sinking material and remineralising it all at
1158               !! the seafloor, these portions of the code below are common to
1159               !! all three.
1160               !!----------------------------------------------------------------------
1161
1162               if (jexport.eq.1) then
1163                  !!======================================================================
1164                  !! BALLAST SUBMODEL
1165                  !!======================================================================
1166                  !!
1167                  !!----------------------------------------------------------------------
1168                  !! Fast-sinking detritus fluxes, pt. 1: REMINERALISATION
1169                  !! aside from explicitly modelled, slow-sinking detritus, the
1170                  !! model includes an implicit representation of detrital
1171                  !! particles that sink too quickly to be modelled with
1172                  !! explicit state variables; this sinking flux is instead
1173                  !! instantaneously remineralised down the water column using
1174                  !! the version of Armstrong et al. (2002)'s ballast model
1175                  !! used by Dunne et al. (2007); the version of this model
1176                  !! here considers silicon and calcium carbonate ballast
1177                  !! minerals; this section of the code redistributes the fast
1178                  !! sinking material generated locally down the water column;
1179                  !! this differs from Dunne et al. (2007) in that fast sinking
1180                  !! material is distributed at *every* level below that it is
1181                  !! generated, rather than at every level below some fixed
1182                  !! depth; this scheme is also different in that sinking material
1183                  !! generated in one level is aggregated with that generated by
1184                  !! shallower levels; this should make the ballast model more
1185                  !! self-consistent (famous last words)
1186                  !!----------------------------------------------------------------------
1187                  !!
1188                  if (jk.eq.1) then
1189                     !! this is the SURFACE OCEAN BOX (no remineralisation)
1190                     !!
1191                     freminc(ji,jj)  = 0.0
1192                     freminn(ji,jj)  = 0.0
1193                     freminfe(ji,jj) = 0.0
1194                     freminsi(ji,jj) = 0.0
1195                     freminca(ji,jj) = 0.0
1196                  elseif (jk.le.mbathy(ji,jj)) then
1197                     !! this is an OCEAN BOX (remineralise some material)
1198                     !!
1199                     !! set up CCD depth to be used depending on user choice
1200                     if (jocalccd.eq.0) then
1201                        !! use default CCD field
1202                        fccd_dep(ji,jj) = ocal_ccd(ji,jj)
1203                     elseif (jocalccd.eq.1) then
1204                        !! use calculated CCD field
1205                        fccd_dep(ji,jj) = f2_ccd_cal(ji,jj)
1206                     endif
1207                     !!
1208                     !! === organic carbon ===
1209                     fq0      = ffastc(ji,jj)                        !! how much organic C enters this box        (mol)
1210                     if (iball.eq.1) then
1211                        fq1      = (fq0 * xmassc)                    !! how much it weighs                        (mass)
1212                        fq2      = (ffastca(ji,jj) * xmassca)        !! how much CaCO3 enters this box            (mass)
1213                        fq3      = (ffastsi(ji,jj) * xmasssi)        !! how much  opal enters this box            (mass)
1214                        fq4      = (fq2 * xprotca) + (fq3 * xprotsi) !! total protected organic C                 (mass)
1215                        !! this next term is calculated for C but used for N and Fe as well
1216                        !! it needs to be protected in case ALL C is protected
1217                        if (fq4.lt.fq1) then
1218                           fprotf(ji,jj)   = (fq4 / (fq1 + tiny(fq1)))      !! protected fraction of total organic C     (non-dim)
1219                        else
1220                           fprotf(ji,jj)   = 1.0                            !! all organic C is protected                (non-dim)
1221                        endif
1222                        fq5      = (1.0 - fprotf(ji,jj))                    !! unprotected fraction of total organic C   (non-dim)
1223                        fq6      = (fq0 * fq5)                       !! how much organic C is unprotected         (mol)
1224                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))     !! how much unprotected C leaves this box    (mol)
1225                        fq8      = (fq7 + (fq0 * fprotf(ji,jj)))            !! how much total C leaves this box          (mol)
1226                        freminc(ji,jj)  = (fq0 - fq8) / fse3t(ji,jj,jk)                !! C remineralisation in this box            (mol)
1227                        ffastc(ji,jj) = fq8                         
1228# if defined key_debug_medusa
1229                        !! report in/out/remin fluxes of carbon for this level
1230                           if (idf.eq.1.AND.idfval.eq.1) then
1231                              IF (lwp) write (numout,*) '------------------------------'
1232                              IF (lwp) write (numout,*) 'totalC(',jk,')  = ', fq1
1233                              IF (lwp) write (numout,*) 'prtctC(',jk,')  = ', fq4
1234                              IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf(ji,jj)
1235                              IF (lwp) write (numout,*) '------------------------------'
1236                              IF (lwp) write (numout,*) 'IN   C(',jk,')  = ', fq0
1237                              IF (lwp) write (numout,*) 'LOST C(',jk,')  = ', freminc(ji,jj) * fse3t(ji,jj,jk)
1238                              IF (lwp) write (numout,*) 'OUT  C(',jk,')  = ', fq8
1239                              IF (lwp) write (numout,*) 'NEW  C(',jk,')  = ', ftempc(ji,jj) * fse3t(ji,jj,jk)
1240                           endif
1241# endif
1242                        else
1243                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))       !! how much organic C leaves this box        (mol)
1244                        freminc(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                !! C remineralisation in this box            (mol)
1245                        ffastc(ji,jj)  = fq1
1246                     endif
1247                     !!
1248                     !! === organic nitrogen ===
1249                     fq0      = ffastn(ji,jj)                        !! how much organic N enters this box        (mol)
1250                     if (iball.eq.1) then
1251                        fq5      = (1.0 - fprotf(ji,jj))                    !! unprotected fraction of total organic N   (non-dim)
1252                        fq6      = (fq0 * fq5)                       !! how much organic N is unprotected         (mol)
1253                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))     !! how much unprotected N leaves this box    (mol)
1254                        fq8      = (fq7 + (fq0 * fprotf(ji,jj)))            !! how much total N leaves this box          (mol)
1255                        freminn(ji,jj)  = (fq0 - fq8) / fse3t(ji,jj,jk)                !! N remineralisation in this box            (mol)
1256                        ffastn(ji,jj)  = fq8
1257# if defined key_debug_medusa
1258                        !! report in/out/remin fluxes of carbon for this level
1259                        if (idf.eq.1.AND.idfval.eq.1) then
1260                           IF (lwp) write (numout,*) '------------------------------'
1261                           IF (lwp) write (numout,*) 'totalN(',jk,')  = ', fq1
1262                           IF (lwp) write (numout,*) 'prtctN(',jk,')  = ', fq4
1263                           IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf(ji,jj)
1264                           IF (lwp) write (numout,*) '------------------------------'
1265                           if (freminn(ji,jj) < 0.0) then
1266                              IF (lwp) write (numout,*) '** FREMIN ERROR **'
1267                           endif
1268                           IF (lwp) write (numout,*) 'IN   N(',jk,')  = ', fq0
1269                           IF (lwp) write (numout,*) 'LOST N(',jk,')  = ', freminn(ji,jj) * fse3t(ji,jj,jk)
1270                           IF (lwp) write (numout,*) 'OUT  N(',jk,')  = ', fq8
1271                           IF (lwp) write (numout,*) 'NEW  N(',jk,')  = ', ftempn(ji,jj) * fse3t(ji,jj,jk)
1272                        endif
1273# endif
1274                     else
1275                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))       !! how much organic N leaves this box        (mol)
1276                        freminn(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                !! N remineralisation in this box            (mol)
1277                        ffastn(ji,jj)  = fq1
1278                     endif
1279                     !!
1280                     !! === organic iron ===
1281                     fq0      = ffastfe(ji,jj)                       !! how much organic Fe enters this box       (mol)
1282                     if (iball.eq.1) then
1283                        fq5      = (1.0 - fprotf(ji,jj))                    !! unprotected fraction of total organic Fe  (non-dim)
1284                        fq6      = (fq0 * fq5)                       !! how much organic Fe is unprotected        (mol)
1285                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))     !! how much unprotected Fe leaves this box   (mol)
1286                        fq8      = (fq7 + (fq0 * fprotf(ji,jj)))            !! how much total Fe leaves this box         (mol)           
1287                        freminfe(ji,jj) = (fq0 - fq8) / fse3t(ji,jj,jk)                !! Fe remineralisation in this box           (mol)
1288                        ffastfe(ji,jj) = fq8
1289                     else
1290                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))       !! how much total Fe leaves this box         (mol)     
1291                        freminfe(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Fe remineralisation in this box           (mol)
1292                        ffastfe(ji,jj) = fq1
1293                     endif
1294                     !!
1295                     !! === biogenic silicon ===
1296                     fq0      = ffastsi(ji,jj)                       !! how much  opal centers this box           (mol)
1297                     fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastsi))         !! how much  opal leaves this box            (mol)
1298                     freminsi(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! Si remineralisation in this box           (mol)
1299                     ffastsi(ji,jj) = fq1
1300                     !!
1301                     !! === biogenic calcium carbonate ===
1302                     fq0      = ffastca(ji,jj)                       !! how much CaCO3 enters this box            (mol)
1303                     if (fsdepw(ji,jj,jk).le.fccd_dep(ji,jj)) then
1304                        !! whole grid cell above CCD
1305                        fq1      = fq0                               !! above lysocline, no Ca dissolves          (mol)
1306                        freminca(ji,jj) = 0.0                               !! above lysocline, no Ca dissolves          (mol)
1307                        fccd(ji,jj) = real(jk)                       !! which is the last level above the CCD?    (#)
1308                     elseif (fsdepw(ji,jj,jk).ge.fccd_dep(ji,jj)) then
1309                        !! whole grid cell below CCD
1310                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastca))      !! how much CaCO3 leaves this box            (mol)
1311                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
1312                     else
1313                        !! partial grid cell below CCD
1314                        fq2      = fdep1(ji,jj) - fccd_dep(ji,jj)                  !! amount of grid cell below CCD             (m)
1315                        fq1      = fq0 * exp(-(fq2 / xfastca))       !! how much CaCO3 leaves this box            (mol)
1316                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
1317                     endif
1318                     ffastca(ji,jj) = fq1 
1319                  else
1320                     !! this is BELOW THE LAST OCEAN BOX (do nothing)
1321                     freminc(ji,jj)  = 0.0
1322                     freminn(ji,jj)  = 0.0
1323                     freminfe(ji,jj) = 0.0
1324                     freminsi(ji,jj) = 0.0
1325                     freminca(ji,jj) = 0.0             
1326                  endif
1327
1328               elseif (jexport.eq.2.or.jexport.eq.3) then
1329                  if (jexport.eq.2) then
1330                     !!======================================================================
1331                     !! MARTIN ET AL. (1987) SUBMODEL
1332                     !!======================================================================
1333                     !!
1334                     !!----------------------------------------------------------------------
1335                     !! This submodel uses the classic Martin et al. (1987) curve
1336                     !! to determine the attenuation of fast-sinking detritus down
1337                     !! the water column.  All three organic elements, C, N and Fe,
1338                     !! are handled identically, and their quantities in sinking
1339                     !! particles attenuate according to a power relationship
1340                     !! governed by parameter "b".  This is assigned a canonical
1341                     !! value of -0.858.  Biogenic opal and calcium carbonate are
1342                     !! attentuated using the same function as in the ballast
1343                     !! submodel
1344                     !!----------------------------------------------------------------------
1345                     !!
1346                     fb_val = -0.858
1347                  elseif (jexport.eq.3) then
1348                     !!======================================================================
1349                     !! HENSON ET AL. (2011) SUBMODEL
1350                     !!======================================================================
1351                     !!
1352                     !!----------------------------------------------------------------------
1353                     !! This submodel reconfigures the Martin et al. (1987) curve by
1354                     !! allowing the "b" value to vary geographically.  Its value is
1355                     !! set, following Henson et al. (2011), as a function of local
1356                     !! sea surface temperature:
1357                     !!   b = -1.06 + (0.024 * SST)
1358                     !! This means that remineralisation length scales are longer in
1359                     !! warm, tropical areas and shorter in cold, polar areas.  This
1360                     !! does seem back-to-front (i.e. one would expect GREATER
1361                     !! remineralisation in warmer waters), but is an outcome of
1362                     !! analysis of sediment trap data, and it may reflect details
1363                     !! of ecosystem structure that pertain to particle production
1364                     !! rather than simply Q10.
1365                     !!----------------------------------------------------------------------
1366                     !!
1367                     fl_sst = tsn(ji,jj,1,jp_tem)
1368                     fb_val = -1.06 + (0.024 * fl_sst)
1369                  endif
1370                  !!   
1371                  if (jk.eq.1) then
1372                     !! this is the SURFACE OCEAN BOX (no remineralisation)
1373                     !!
1374                     freminc(ji,jj)  = 0.0
1375                     freminn(ji,jj)  = 0.0
1376                     freminfe(ji,jj) = 0.0
1377                     freminsi(ji,jj) = 0.0
1378                     freminca(ji,jj) = 0.0
1379                  elseif (jk.le.mbathy(ji,jj)) then
1380                     !! this is an OCEAN BOX (remineralise some material)
1381                     !!
1382                     !! === organic carbon ===
1383                     fq0      = ffastc(ji,jj)                        !! how much organic C enters this box        (mol)
1384                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)         !! how much organic C leaves this box        (mol)
1385                     freminc(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! C remineralisation in this box            (mol)
1386                     ffastc(ji,jj)  = fq1
1387                     !!
1388                     !! === organic nitrogen ===
1389                     fq0      = ffastn(ji,jj)                        !! how much organic N enters this box        (mol)
1390                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)         !! how much organic N leaves this box        (mol)
1391                     freminn(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! N remineralisation in this box            (mol)
1392                     ffastn(ji,jj)  = fq1
1393                     !!
1394                     !! === organic iron ===
1395                     fq0      = ffastfe(ji,jj)                       !! how much organic Fe enters this box       (mol)
1396                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)         !! how much organic Fe leaves this box       (mol)
1397                     freminfe(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! Fe remineralisation in this box           (mol)
1398                     ffastfe(ji,jj) = fq1
1399                     !!
1400                     !! === biogenic silicon ===
1401                     fq0      = ffastsi(ji,jj)                       !! how much  opal centers this box           (mol)
1402                     fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastsi))         !! how much  opal leaves this box            (mol)
1403                     freminsi(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! Si remineralisation in this box           (mol)
1404                     ffastsi(ji,jj) = fq1
1405                     !!
1406                     !! === biogenic calcium carbonate ===
1407                     fq0      = ffastca(ji,jj)                       !! how much CaCO3 enters this box            (mol)
1408                     if (fsdepw(ji,jj,jk).le.ocal_ccd(ji,jj)) then
1409                        !! whole grid cell above CCD
1410                        fq1      = fq0                               !! above lysocline, no Ca dissolves          (mol)
1411                        freminca(ji,jj) = 0.0                               !! above lysocline, no Ca dissolves          (mol)
1412                        fccd(ji,jj) = real(jk)                       !! which is the last level above the CCD?    (#)
1413                     elseif (fsdepw(ji,jj,jk).ge.ocal_ccd(ji,jj)) then
1414                        !! whole grid cell below CCD
1415                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastca))      !! how much CaCO3 leaves this box            (mol)
1416                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
1417                     else
1418                        !! partial grid cell below CCD
1419                        fq2      = fdep1(ji,jj) - ocal_ccd(ji,jj)           !! amount of grid cell below CCD             (m)
1420                        fq1      = fq0 * exp(-(fq2 / xfastca))       !! how much CaCO3 leaves this box            (mol)
1421                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
1422                     endif
1423                     ffastca(ji,jj) = fq1 
1424                  else
1425                     !! this is BELOW THE LAST OCEAN BOX (do nothing)
1426                     freminc(ji,jj)  = 0.0
1427                     freminn(ji,jj)  = 0.0
1428                     freminfe(ji,jj) = 0.0
1429                     freminsi(ji,jj) = 0.0
1430                     freminca(ji,jj) = 0.0             
1431                  endif
1432
1433               endif
1434
1435               !!----------------------------------------------------------------------
1436               !! Fast-sinking detritus fluxes, pt. 2: UPDATE FAST FLUXES
1437               !! here locally calculated additions to the fast-sinking flux are added
1438               !! to the total fast-sinking flux; this is done here such that material
1439               !! produced in a particular layer is only remineralised below this
1440               !! layer
1441               !!----------------------------------------------------------------------
1442               !!
1443               !! add sinking material generated in this layer to running totals
1444               !!
1445               !! === organic carbon ===                          (diatom and mesozooplankton mortality)
1446               ffastc(ji,jj)  = ffastc(ji,jj)  + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
1447               !!
1448               !! === organic nitrogen ===                        (diatom and mesozooplankton mortality)
1449               ffastn(ji,jj)  = ffastn(ji,jj)  + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
1450               !!
1451               !! === organic iron ===                            (diatom and mesozooplankton mortality)
1452               ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))
1453               !!
1454               !! === biogenic silicon ===                        (diatom mortality and grazed diatoms)
1455               ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
1456               !!
1457               !! === biogenic calcium carbonate ===              (latitudinally-based fraction of total primary production)
1458               ffastca(ji,jj) = ffastca(ji,jj) + (ftempca(ji,jj) * fse3t(ji,jj,jk))
1459
1460               !!----------------------------------------------------------------------
1461               !! Fast-sinking detritus fluxes, pt. 3: SEAFLOOR
1462               !! remineralise all remaining fast-sinking detritus to dissolved
1463               !! nutrients; the sedimentation fluxes calculated here allow the
1464               !! separation of what's remineralised sinking through the final
1465               !! ocean box from that which is added to the final box by the
1466               !! remineralisation of material that reaches the seafloor (i.e.
1467               !! the model assumes that *all* material that hits the seafloor
1468               !! is remineralised and that none is permanently buried; hey,
1469               !! this is a giant GCM model that can't be run for long enough
1470               !! to deal with burial fluxes!)
1471               !!
1472               !! in a change to this process, in part so that MEDUSA behaves
1473               !! a little more like ERSEM et al., fast-sinking detritus (N, Fe
1474               !! and C) is converted to slow sinking detritus at the seafloor
1475               !! instead of being remineralised; the rationale is that in
1476               !! shallower shelf regions (... that are not fully mixed!) this
1477               !! allows the detrital material to return slowly to dissolved
1478               !! nutrient rather than instantaneously as now; the alternative
1479               !! would be to explicitly handle seafloor organic material - a
1480               !! headache I don't wish to experience at this point; note that
1481               !! fast-sinking Si and Ca detritus is just remineralised as
1482               !! per usual
1483               !!
1484               !! AXY (13/01/12)
1485               !! in a further change to this process, again so that MEDUSA is
1486               !! a little more like ERSEM et al., material that reaches the
1487               !! seafloor can now be added to sediment pools and stored for
1488               !! slow release; there are new 2D arrays for organic nitrogen,
1489               !! iron and carbon and inorganic silicon and carbon that allow
1490               !! fast and slow detritus that reaches the seafloor to be held
1491               !! and released back to the water column more slowly; these arrays
1492               !! are transferred via the tracer restart files between repeat
1493               !! submissions of the model
1494               !!----------------------------------------------------------------------
1495               !!
1496               ffast2slowc(ji,jj)  = 0.0
1497               ffast2slown(ji,jj)  = 0.0
1498! I don't think this is used - marc 10/4/17
1499!               ffast2slowfe(ji,jj) = 0.0
1500               !!
1501               if (jk.eq.mbathy(ji,jj)) then
1502                  !! this is the BOTTOM OCEAN BOX (remineralise everything)
1503                  !!
1504                  !! AXY (17/01/12): tweaked to include benthos pools
1505                  !!
1506                  !! === organic carbon ===
1507                  if (jfdfate.eq.0 .and. jorgben.eq.0) then
1508                     freminc(ji,jj)  = freminc(ji,jj) + (ffastc(ji,jj) / fse3t(ji,jj,jk))    !! C remineralisation in this box            (mol/m3)
1509                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
1510                     ffast2slowc(ji,jj) = ffastc(ji,jj) / fse3t(ji,jj,jk)             !! fast C -> slow C                          (mol/m3)
1511                     fslowc(ji,jj)      = fslowc(ji,jj) + ffast2slowc(ji,jj)
1512                  elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
1513                     f_fbenin_c(ji,jj)  = ffastc(ji,jj)             !! fast C -> benthic C                       (mol/m2)
1514                  endif
1515                  fsedc(ji,jj)   = ffastc(ji,jj)                          !! record seafloor C                         (mol/m2)
1516                  ffastc(ji,jj)  = 0.0
1517                  !!
1518                  !! === organic nitrogen ===
1519                  if (jfdfate.eq.0 .and. jorgben.eq.0) then
1520                     freminn(ji,jj)  = freminn(ji,jj) + (ffastn(ji,jj) / fse3t(ji,jj,jk))    !! N remineralisation in this box            (mol/m3)
1521                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
1522                     ffast2slown(ji,jj) = ffastn(ji,jj) / fse3t(ji,jj,jk)             !! fast N -> slow N                          (mol/m3)
1523                     fslown(ji,jj)      = fslown(ji,jj) + ffast2slown(ji,jj)
1524                  elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
1525                     f_fbenin_n(ji,jj)  = ffastn(ji,jj)             !! fast N -> benthic N                       (mol/m2)
1526                  endif
1527                  fsedn(ji,jj)   = ffastn(ji,jj)                          !! record seafloor N                         (mol/m2)
1528                  ffastn(ji,jj)  = 0.0
1529                  !!
1530                  !! === organic iron ===
1531                  if (jfdfate.eq.0 .and. jorgben.eq.0) then
1532                     freminfe(ji,jj) = freminfe(ji,jj) + (ffastfe(ji,jj) / fse3t(ji,jj,jk))  !! Fe remineralisation in this box           (mol/m3)
1533! I don't think ffast2slowfe is used - marc 10/4/17
1534!                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
1535!                     ffast2slowfe(ji,jj) = ffastn(ji,jj) / fse3t(ji,jj,jk)            !! fast Fe -> slow Fe                        (mol/m3)
1536                  elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
1537                     f_fbenin_fe(ji,jj) = ffastfe(ji,jj)            !! fast Fe -> benthic Fe                     (mol/m2)
1538                  endif
1539                  fsedfe(ji,jj)  = ffastfe(ji,jj)                         !! record seafloor Fe                        (mol/m2)
1540                  ffastfe(ji,jj) = 0.0
1541                  !!
1542                  !! === biogenic silicon ===
1543                  if (jinorgben.eq.0) then
1544                     freminsi(ji,jj) = freminsi(ji,jj) + (ffastsi(ji,jj) / fse3t(ji,jj,jk))  !! Si remineralisation in this box           (mol/m3)
1545                  elseif (jinorgben.eq.1) then
1546                     f_fbenin_si(ji,jj) = ffastsi(ji,jj)            !! fast Si -> benthic Si                     (mol/m2)
1547                  endif
1548                  fsedsi(ji,jj)   = ffastsi(ji,jj)                         !! record seafloor Si                        (mol/m2)
1549                  ffastsi(ji,jj) = 0.0
1550                  !!
1551                  !! === biogenic calcium carbonate ===
1552                  if (jinorgben.eq.0) then
1553                     freminca(ji,jj) = freminca(ji,jj) + (ffastca(ji,jj) / fse3t(ji,jj,jk))  !! Ca remineralisation in this box           (mol/m3)
1554                  elseif (jinorgben.eq.1) then
1555                     f_fbenin_ca(ji,jj) = ffastca(ji,jj)            !! fast Ca -> benthic Ca                     (mol/m2)
1556                  endif
1557                  fsedca(ji,jj)   = ffastca(ji,jj)                         !! record seafloor Ca                        (mol/m2)
1558                  ffastca(ji,jj) = 0.0
1559               endif
1560
1561# if defined key_debug_medusa
1562               if (idf.eq.1) then
1563                  !!----------------------------------------------------------------------
1564                  !! Integrate total fast detritus remineralisation
1565                  !!----------------------------------------------------------------------
1566                  !!
1567                  fofd_n(ji,jj)  = fofd_n(ji,jj)  + (freminn(ji,jj)  * fse3t(ji,jj,jk))
1568                  fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
1569                  fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe(ji,jj) * fse3t(ji,jj,jk))
1570#  if defined key_roam
1571                  fofd_c(ji,jj)  = fofd_c(ji,jj)  + (freminc(ji,jj)  * fse3t(ji,jj,jk))
1572#  endif
1573               endif
1574# endif
1575
1576               !!----------------------------------------------------------------------
1577               !! Sort out remineralisation tally of fast-sinking detritus
1578               !!----------------------------------------------------------------------
1579               !!
1580               !! update fast-sinking regeneration arrays
1581               fregenfast(ji,jj)   = fregenfast(ji,jj)   + (freminn(ji,jj)  * fse3t(ji,jj,jk))
1582               fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
1583# if defined key_roam
1584               fregenfastc(ji,jj)  = fregenfastc(ji,jj)  + (freminc(ji,jj)  * fse3t(ji,jj,jk))
1585# endif
1586
1587               !!----------------------------------------------------------------------
1588               !! Benthic remineralisation fluxes
1589               !!----------------------------------------------------------------------
1590               !!
1591               if (jk.eq.mbathy(ji,jj)) then
1592                  !!
1593                  !! organic components
1594                  if (jorgben.eq.1) then
1595                     f_benout_n(ji,jj)  = xsedn  * zn_sed_n(ji,jj)
1596                     f_benout_fe(ji,jj) = xsedfe * zn_sed_fe(ji,jj)
1597                     f_benout_c(ji,jj)  = xsedc  * zn_sed_c(ji,jj)
1598                  endif
1599                  !!
1600                  !! inorganic components
1601                  if (jinorgben.eq.1) then
1602                     f_benout_si(ji,jj) = xsedsi * zn_sed_si(ji,jj)
1603                     f_benout_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
1604                     !!
1605                     !! account for CaCO3 that dissolves when it shouldn't
1606                     if ( fsdepw(ji,jj,jk) .le. fccd_dep(ji,jj) ) then
1607                        f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
1608                     endif
1609                  endif
1610               endif
1611               CALL flush(numout)
1612
1613! MAYBE BUT A BREAK IN HERE, BUSINESS AND UPDATING - marc 20/4/17
1614! (fast sinking detritus is 576 lines)
1615            ENDIF
1616         ENDDO
1617         ENDDO
1618
1619         DO jj = 2,jpjm1
1620         DO ji = 2,jpim1
1621            !! OPEN wet point IF..THEN loop
1622            if (tmask(ji,jj,jk) == 1) then
1623
1624               !!======================================================================
1625               !! LOCAL GRID CELL TRENDS
1626               !!======================================================================
1627               !!
1628               !!----------------------------------------------------------------------
1629               !! Determination of trends
1630               !!----------------------------------------------------------------------
1631               !!
1632               !!----------------------------------------------------------------------
1633               !! chlorophyll
1634               btra(ji,jj,jpchn) = b0 * ( &
1635                 + ((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) )
1636               btra(ji,jj,jpchd) = b0 * ( &
1637                 + ((frd(ji,jj) * fprd(ji,jj) * zphd(ji,jj)) - fgmepd(ji,jj) - fdpd(ji,jj) - fdpd2(ji,jj)) * (fthetad(ji,jj) / xxi) )
1638               !!
1639               !!----------------------------------------------------------------------
1640               !! phytoplankton
1641               btra(ji,jj,jpphn) = b0 * ( &
1642                 + (fprn(ji,jj) * zphn(ji,jj)) - fgmipn(ji,jj) - fgmepn(ji,jj) - fdpn(ji,jj) - fdpn2(ji,jj) )
1643               btra(ji,jj,jpphd) = b0 * ( &
1644                 + (fprd(ji,jj) * zphd(ji,jj)) - fgmepd(ji,jj) - fdpd(ji,jj) - fdpd2(ji,jj) )
1645               btra(ji,jj,jppds) = b0 * ( &
1646                 + (fprds(ji,jj) * zpds(ji,jj)) - fgmepds(ji,jj) - fdpds(ji,jj) - fsdiss(ji,jj) - fdpds2(ji,jj) )
1647               !!
1648               !!----------------------------------------------------------------------
1649               !! zooplankton
1650               btra(ji,jj,jpzmi) = b0 * ( &
1651                 + fmigrow(ji,jj) - fgmezmi(ji,jj) - fdzmi(ji,jj) - fdzmi2(ji,jj) )
1652               btra(ji,jj,jpzme) = b0 * ( &
1653                 + fmegrow(ji,jj) - fdzme(ji,jj) - fdzme2(ji,jj) )
1654               !!
1655               !!----------------------------------------------------------------------
1656               !! detritus
1657               btra(ji,jj,jpdet) = b0 * ( &
1658                 + fdpn(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj))              &  ! mort. losses
1659                 + fdzmi(ji,jj) + ((1.0 - xfdfrac2) * fdzme(ji,jj))            &  ! mort. losses
1660                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))            &  ! assim. inefficiency
1661                 - fgmid(ji,jj) - fgmed(ji,jj) - fdd(ji,jj)                           &  ! grazing and remin.
1662                 + ffast2slown(ji,jj) )                                    ! seafloor fast->slow
1663               !!
1664               !!----------------------------------------------------------------------
1665               !! dissolved inorganic nitrogen nutrient
1666               fn_cons = 0.0  &
1667                 - (fprn(ji,jj) * zphn(ji,jj)) - (fprd(ji,jj) * zphd(ji,jj))                    ! primary production
1668               fn_prod = 0.0  &
1669                 + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)))                     &  ! messy feeding remin.
1670                 + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)))  &  ! messy feeding remin.
1671                 + fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) + freminn(ji,jj)             &  ! excretion and remin.
1672                 + fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)                  ! metab. losses
1673               !!
1674               !! riverine flux
1675               if ( jriver_n .gt. 0 ) then
1676                  f_riv_loc_n(ji,jj) = f_riv_n(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1677                  fn_prod = fn_prod + f_riv_loc_n(ji,jj)
1678               endif
1679               !! 
1680               !! benthic remineralisation
1681               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1682                  fn_prod = fn_prod + (f_benout_n(ji,jj) / fse3t(ji,jj,jk))
1683               endif
1684               !!
1685               btra(ji,jj,jpdin) = b0 * ( &
1686                 fn_prod + fn_cons )
1687               !!
1688               fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved nitrogen
1689                 fn_cons ) )
1690               fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved nitrogen
1691                 fn_prod ) )
1692               !!
1693               !!----------------------------------------------------------------------
1694               !! dissolved silicic acid nutrient
1695               fs_cons = 0.0  &
1696                 - (fprds(ji,jj) * zpds(ji,jj))                                   ! opal production
1697               fs_prod = 0.0  &
1698                 + fsdiss(ji,jj)                                        &  ! opal dissolution
1699                 + ((1.0 - xfdfrac1) * fdpds(ji,jj))                    &  ! mort. loss
1700                 + ((1.0 - xfdfrac3) * fgmepds(ji,jj))                  &  ! egestion of grazed Si
1701                 + freminsi(ji,jj) + fdpds2(ji,jj)                                ! fast diss. and metab. losses
1702               !!
1703               !! riverine flux
1704               if ( jriver_si .gt. 0 ) then
1705                  f_riv_loc_si(ji,jj) = f_riv_si(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1706                  fs_prod = fs_prod + f_riv_loc_si(ji,jj)
1707               endif
1708               !! 
1709               !! benthic remineralisation
1710               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
1711                  fs_prod = fs_prod + (f_benout_si(ji,jj) / fse3t(ji,jj,jk))
1712               endif
1713               !!
1714               btra(ji,jj,jpsil) = b0 * ( &
1715                 fs_prod + fs_cons )
1716               !!
1717               fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved silicon
1718                 fs_cons ) )
1719               fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved silicon
1720                 fs_prod ) )
1721               !!
1722               !!----------------------------------------------------------------------
1723               !! dissolved "iron" nutrient
1724               btra(ji,jj,jpfer) = b0 * ( &
1725               + (xrfn * btra(ji,jj,jpdin)) + ffetop(ji,jj) + ffebot(ji,jj) - ffescav(ji,jj) )
1726
1727# if defined key_roam
1728               !!
1729               !!----------------------------------------------------------------------
1730               !! AXY (26/11/08): implicit detrital carbon change
1731               btra(ji,jj,jpdtc) = b0 * ( &
1732                 + (xthetapn * fdpn(ji,jj)) + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))      &  ! mort. losses
1733                 + (xthetazmi * fdzmi(ji,jj)) + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))  &  ! mort. losses
1734                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))                              &  ! assim. inefficiency
1735                 - fgmidc(ji,jj) - fgmedc(ji,jj) - fddc(ji,jj)                                          &  ! grazing and remin.
1736                 + ffast2slowc(ji,jj) )                                                      ! seafloor fast->slow
1737               !!
1738               !!----------------------------------------------------------------------
1739               !! dissolved inorganic carbon
1740               fc_cons = 0.0  &
1741                 - (xthetapn * fprn(ji,jj) * zphn(ji,jj)) - (xthetapd * fprd(ji,jj) * zphd(ji,jj))                ! primary production
1742               fc_prod = 0.0  &
1743                 + (xthetapn * xphi * fgmipn(ji,jj)) + (xphi * fgmidc(ji,jj))                    &  ! messy feeding remin
1744                 + (xthetapn * xphi * fgmepn(ji,jj)) + (xthetapd * xphi * fgmepd(ji,jj))         &  ! messy feeding remin
1745                 + (xthetazmi * xphi * fgmezmi(ji,jj)) + (xphi * fgmedc(ji,jj))                  &  ! messy feeding remin
1746                 + fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) + freminc(ji,jj) + (xthetapn * fdpn2(ji,jj))         &  ! resp., remin., losses
1747                 + (xthetapd * fdpd2(ji,jj)) + (xthetazmi * fdzmi2(ji,jj))                       &  ! losses
1748                 + (xthetazme * fdzme2(ji,jj))                                               ! losses
1749               !!
1750               !! riverine flux
1751               if ( jriver_c .gt. 0 ) then
1752                  f_riv_loc_c(ji,jj) = f_riv_c(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1753                  fc_prod = fc_prod + f_riv_loc_c(ji,jj)
1754               endif
1755               !! 
1756               !! benthic remineralisation
1757               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1758                  fc_prod = fc_prod + (f_benout_c(ji,jj) / fse3t(ji,jj,jk))
1759               endif
1760               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
1761                  fc_prod = fc_prod + (f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
1762               endif
1763               !!
1764               !! community respiration (does not include CaCO3 terms - obviously!)
1765               fcomm_resp(ji,jj) = fcomm_resp(ji,jj) + fc_prod
1766               !!
1767               !! CaCO3
1768               fc_prod = fc_prod - ftempca(ji,jj) + freminca(ji,jj)
1769               !!
1770               !! riverine flux
1771               if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then
1772                  fc_prod = fc_prod + f_riv_c(ji,jj)
1773               endif
1774               !!
1775               btra(ji,jj,jpdic) = b0 * ( &
1776                 fc_prod + fc_cons )
1777               !!
1778               fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved carbon
1779                 fc_cons ) )
1780               fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved carbon
1781                 fc_prod ) )
1782               !!
1783               !!----------------------------------------------------------------------
1784               !! alkalinity
1785               fa_prod = 0.0  &
1786                 + (2.0 * freminca(ji,jj))                                                   ! CaCO3 dissolution
1787               fa_cons = 0.0  &
1788                 - (2.0 * ftempca(ji,jj))                                                    ! CaCO3 production
1789               !!
1790               !! riverine flux
1791               if ( jriver_alk .gt. 0 ) then
1792                  f_riv_loc_alk(ji,jj) = f_riv_alk(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
1793                  fa_prod = fa_prod + f_riv_loc_alk(ji,jj)
1794               endif
1795               !! 
1796               !! benthic remineralisation
1797               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
1798                  fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
1799               endif
1800               !!
1801               btra(ji,jj,jpalk) = b0 * ( &
1802                 fa_prod + fa_cons )
1803               !!
1804               !!----------------------------------------------------------------------
1805               !! oxygen (has protection at low O2 concentrations; OCMIP-2 style)
1806               fo2_prod(ji,jj) = 0.0 &
1807                 + (xthetanit * fprn(ji,jj) * zphn(ji,jj))                                      & ! Pn primary production, N
1808                 + (xthetanit * fprd(ji,jj) * zphd(ji,jj))                                      & ! Pd primary production, N
1809                 + (xthetarem * xthetapn * fprn(ji,jj) * zphn(ji,jj))                           & ! Pn primary production, C
1810                 + (xthetarem * xthetapd * fprd(ji,jj) * zphd(ji,jj))                             ! Pd primary production, C
1811               fo2_ncons(ji,jj) = 0.0 &
1812                 - (xthetanit * xphi * fgmipn(ji,jj))                                    & ! Pn messy feeding remin., N
1813                 - (xthetanit * xphi * fgmid(ji,jj))                                     & ! D  messy feeding remin., N
1814                 - (xthetanit * xphi * fgmepn(ji,jj))                                    & ! Pn messy feeding remin., N
1815                 - (xthetanit * xphi * fgmepd(ji,jj))                                    & ! Pd messy feeding remin., N
1816                 - (xthetanit * xphi * fgmezmi(ji,jj))                                   & ! Zi messy feeding remin., N
1817                 - (xthetanit * xphi * fgmed(ji,jj))                                     & ! D  messy feeding remin., N
1818                 - (xthetanit * fmiexcr(ji,jj))                                          & ! microzoo excretion, N
1819                 - (xthetanit * fmeexcr(ji,jj))                                          & ! mesozoo  excretion, N
1820                 - (xthetanit * fdd(ji,jj))                                              & ! slow detritus remin., N
1821                 - (xthetanit * freminn(ji,jj))                                          & ! fast detritus remin., N
1822                 - (xthetanit * fdpn2(ji,jj))                                            & ! Pn  losses, N
1823                 - (xthetanit * fdpd2(ji,jj))                                            & ! Pd  losses, N
1824                 - (xthetanit * fdzmi2(ji,jj))                                           & ! Zmi losses, N
1825                 - (xthetanit * fdzme2(ji,jj))                                             ! Zme losses, N
1826               !! 
1827               !! benthic remineralisation
1828               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1829                  fo2_ncons(ji,jj) = fo2_ncons(ji,jj) - (xthetanit * f_benout_n(ji,jj) / fse3t(ji,jj,jk))
1830               endif
1831               fo2_ccons(ji,jj) = 0.0 &
1832                 - (xthetarem * xthetapn * xphi * fgmipn(ji,jj))                         & ! Pn messy feeding remin., C
1833                 - (xthetarem * xphi * fgmidc(ji,jj))                                    & ! D  messy feeding remin., C
1834                 - (xthetarem * xthetapn * xphi * fgmepn(ji,jj))                         & ! Pn messy feeding remin., C
1835                 - (xthetarem * xthetapd * xphi * fgmepd(ji,jj))                         & ! Pd messy feeding remin., C
1836                 - (xthetarem * xthetazmi * xphi * fgmezmi(ji,jj))                       & ! Zi messy feeding remin., C
1837                 - (xthetarem * xphi * fgmedc(ji,jj))                                    & ! D  messy feeding remin., C
1838                 - (xthetarem * fmiresp(ji,jj))                                          & ! microzoo respiration, C
1839                 - (xthetarem * fmeresp(ji,jj))                                          & ! mesozoo  respiration, C
1840                 - (xthetarem * fddc(ji,jj))                                             & ! slow detritus remin., C
1841                 - (xthetarem * freminc(ji,jj))                                          & ! fast detritus remin., C
1842                 - (xthetarem * xthetapn * fdpn2(ji,jj))                                 & ! Pn  losses, C
1843                 - (xthetarem * xthetapd * fdpd2(ji,jj))                                 & ! Pd  losses, C
1844                 - (xthetarem * xthetazmi * fdzmi2(ji,jj))                               & ! Zmi losses, C
1845                 - (xthetarem * xthetazme * fdzme2(ji,jj))                                 ! Zme losses, C
1846               !! 
1847               !! benthic remineralisation
1848               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
1849                  fo2_ccons(ji,jj) = fo2_ccons(ji,jj) - (xthetarem * f_benout_c(ji,jj) / fse3t(ji,jj,jk))
1850               endif
1851               fo2_cons(ji,jj) = fo2_ncons(ji,jj) + fo2_ccons(ji,jj)
1852               !!
1853               !! is this a suboxic zone?
1854               if (zoxy(ji,jj).lt.xo2min) then  ! deficient O2; production fluxes only
1855                  btra(ji,jj,jpoxy) = b0 * ( &
1856                    fo2_prod(ji,jj) )
1857                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
1858                  foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
1859               else                      ! sufficient O2; production + consumption fluxes
1860                  btra(ji,jj,jpoxy) = b0 * ( &
1861                    fo2_prod(ji,jj) + fo2_cons(ji,jj) )
1862                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
1863                  foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
1864               endif
1865               !!
1866               !! air-sea fluxes (if this is the surface box)
1867               if (jk.eq.1) then
1868                  !!
1869                  !! CO2 flux
1870                  btra(ji,jj,jpdic) = btra(ji,jj,jpdic) + (b0 * f_co2flux(ji,jj))
1871                  !!
1872                  !! O2 flux (mol/m3/s -> mmol/m3/d)
1873                  btra(ji,jj,jpoxy) = btra(ji,jj,jpoxy) + (b0 * f_o2flux(ji,jj))
1874               endif
1875# endif
1876
1877# if defined key_debug_medusa
1878               !! report state variable fluxes (not including fast-sinking detritus)
1879               if (idf.eq.1.AND.idfval.eq.1) then
1880                  IF (lwp) write (numout,*) '------------------------------'
1881                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchn)(',jk,')  = ', btra(ji,jj,jpchn)
1882                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchd)(',jk,')  = ', btra(ji,jj,jpchd)
1883                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphn)(',jk,')  = ', btra(ji,jj,jpphn)
1884                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphd)(',jk,')  = ', btra(ji,jj,jpphd)
1885                  IF (lwp) write (numout,*) 'btra(ji,jj,jppds)(',jk,')  = ', btra(ji,jj,jppds)
1886                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzmi)(',jk,')  = ', btra(ji,jj,jpzmi)
1887                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzme)(',jk,')  = ', btra(ji,jj,jpzme)
1888                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdet)(',jk,')  = ', btra(ji,jj,jpdet)
1889                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdin)(',jk,')  = ', btra(ji,jj,jpdin)
1890                  IF (lwp) write (numout,*) 'btra(ji,jj,jpsil)(',jk,')  = ', btra(ji,jj,jpsil)
1891                  IF (lwp) write (numout,*) 'btra(ji,jj,jpfer)(',jk,')  = ', btra(ji,jj,jpfer)
1892#  if defined key_roam
1893                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdtc)(',jk,')  = ', btra(ji,jj,jpdtc)
1894                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdic)(',jk,')  = ', btra(ji,jj,jpdic)
1895                  IF (lwp) write (numout,*) 'btra(ji,jj,jpalk)(',jk,')  = ', btra(ji,jj,jpalk)
1896                  IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy)(',jk,')  = ', btra(ji,jj,jpoxy)
1897#  endif
1898               endif
1899# endif
1900
1901               !!----------------------------------------------------------------------
1902               !! Integrate calculated fluxes for mass balance
1903               !!----------------------------------------------------------------------
1904               !!
1905               !! === nitrogen ===
1906               fflx_n(ji,jj)  = fflx_n(ji,jj)  + &
1907                  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) )
1908               !! === silicon ===
1909               fflx_si(ji,jj) = fflx_si(ji,jj) + &
1910                  fse3t(ji,jj,jk) * ( btra(ji,jj,jppds) + btra(ji,jj,jpsil) )
1911               !! === iron ===
1912               fflx_fe(ji,jj) = fflx_fe(ji,jj) + &
1913                  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) )
1914# if defined key_roam
1915               !! === carbon ===
1916               fflx_c(ji,jj)  = fflx_c(ji,jj)  + &
1917                  fse3t(ji,jj,jk) * ( (xthetapn * btra(ji,jj,jpphn)) + (xthetapd * btra(ji,jj,jpphd)) + &
1918                  (xthetazmi * btra(ji,jj,jpzmi)) + (xthetazme * btra(ji,jj,jpzme)) + btra(ji,jj,jpdtc) + btra(ji,jj,jpdic) )
1919               !! === alkalinity ===
1920               fflx_a(ji,jj)  = fflx_a(ji,jj)  + &
1921                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpalk) )
1922               !! === oxygen ===
1923               fflx_o2(ji,jj) = fflx_o2(ji,jj) + &
1924                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpoxy) )
1925# endif
1926
1927               !!----------------------------------------------------------------------
1928               !! Apply calculated tracer fluxes
1929               !!----------------------------------------------------------------------
1930               !!
1931               !! units: [unit of tracer] per second (fluxes are calculated above per day)
1932               !!
1933               ibio_switch = 1
1934# if defined key_gulf_finland
1935               !! AXY (17/05/13): fudge in a Gulf of Finland correction; uses longitude-
1936               !!                 latitude range to establish if this is a Gulf of Finland
1937               !!                 grid cell; if so, then BGC fluxes are ignored (though
1938               !!                 still calculated); for reference, this is meant to be a
1939               !!                 temporary fix to see if all of my problems can be done
1940               !!                 away with if I switch off BGC fluxes in the Gulf of
1941               !!                 Finland, which currently appears the source of trouble
1942               if ( glamt(ji,jj).gt.24.7 .and. glamt(ji,jj).lt.27.8 .and. &
1943                  &   gphit(ji,jj).gt.59.2 .and. gphit(ji,jj).lt.60.2 ) then
1944                  ibio_switch = 0
1945               endif
1946# endif               
1947               if (ibio_switch.eq.1) then
1948                  tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(ji,jj,jpchn) / 86400.)
1949                  tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(ji,jj,jpchd) / 86400.)
1950                  tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(ji,jj,jpphn) / 86400.)
1951                  tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(ji,jj,jpphd) / 86400.)
1952                  tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(ji,jj,jppds) / 86400.)
1953                  tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(ji,jj,jpzmi) / 86400.)
1954                  tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(ji,jj,jpzme) / 86400.)
1955                  tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(ji,jj,jpdet) / 86400.)
1956                  tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(ji,jj,jpdin) / 86400.)
1957                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(ji,jj,jpsil) / 86400.)
1958                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(ji,jj,jpfer) / 86400.)
1959# if defined key_roam
1960                  tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(ji,jj,jpdtc) / 86400.)
1961                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(ji,jj,jpdic) / 86400.)
1962                  tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(ji,jj,jpalk) / 86400.)
1963                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(ji,jj,jpoxy) / 86400.)
1964# endif
1965               endif               
1966
1967               !! AXY (18/11/16): CMIP6 diagnostics
1968               IF( med_diag%FBDDTALK%dgsave )  THEN
1969                  fbddtalk(ji,jj)  =  fbddtalk(ji,jj)  + (btra(ji,jj,jpalk) * fse3t(ji,jj,jk))
1970               ENDIF
1971               IF( med_diag%FBDDTDIC%dgsave )  THEN
1972                  fbddtdic(ji,jj)  =  fbddtdic(ji,jj)  + (btra(ji,jj,jpdic) * fse3t(ji,jj,jk))
1973               ENDIF
1974               IF( med_diag%FBDDTDIFE%dgsave ) THEN
1975                  fbddtdife(ji,jj) =  fbddtdife(ji,jj) + (btra(ji,jj,jpfer) * fse3t(ji,jj,jk))
1976               ENDIF
1977               IF( med_diag%FBDDTDIN%dgsave )  THEN
1978                  fbddtdin(ji,jj)  =  fbddtdin(ji,jj)  + (btra(ji,jj,jpdin) * fse3t(ji,jj,jk))
1979               ENDIF
1980               IF( med_diag%FBDDTDISI%dgsave ) THEN
1981                  fbddtdisi(ji,jj) =  fbddtdisi(ji,jj) + (btra(ji,jj,jpsil) * fse3t(ji,jj,jk))
1982               ENDIF
1983          !!
1984               IF( med_diag%BDDTALK3%dgsave )  THEN
1985                  bddtalk3(ji,jj,jk)  =  btra(ji,jj,jpalk)
1986               ENDIF
1987               IF( med_diag%BDDTDIC3%dgsave )  THEN
1988                  bddtdic3(ji,jj,jk)  =  btra(ji,jj,jpdic)
1989               ENDIF
1990               IF( med_diag%BDDTDIFE3%dgsave ) THEN
1991                  bddtdife3(ji,jj,jk) =  btra(ji,jj,jpfer)
1992               ENDIF
1993               IF( med_diag%BDDTDIN3%dgsave )  THEN
1994                  bddtdin3(ji,jj,jk)  =  btra(ji,jj,jpdin)
1995               ENDIF
1996               IF( med_diag%BDDTDISI3%dgsave ) THEN
1997                  bddtdisi3(ji,jj,jk) =  btra(ji,jj,jpsil)
1998               ENDIF
1999
2000#   if defined key_debug_medusa
2001               IF (lwp) write (numout,*) '------'
2002               IF (lwp) write (numout,*) 'trc_bio_medusa: end all calculations'
2003               IF (lwp) write (numout,*) 'trc_bio_medusa: now outputs'
2004                     CALL flush(numout)
2005#   endif
2006
2007# if defined key_axy_nancheck
2008               !!----------------------------------------------------------------------
2009               !! Check calculated tracer fluxes
2010               !!----------------------------------------------------------------------
2011               !!
2012               DO jn = 1,jptra
2013                  fq0 = btra(ji,jj,jn)
2014                  !! AXY (30/01/14): "isnan" problem on HECTOR
2015                  !! if (fq0 /= fq0 ) then
2016                  if ( ieee_is_nan( fq0 ) ) then
2017                     !! there's a NaN here
2018                     if (lwp) write(numout,*) 'NAN detected in btra(ji,jj,', ji, ',', &
2019                     & jj, ',', jk, ',', jn, ') at time', kt
2020           CALL ctl_stop( 'trcbio_medusa, NAN in btra field' )
2021                  endif
2022               ENDDO
2023               DO jn = 1,jptra
2024                  fq0 = tra(ji,jj,jk,jn)
2025                  !! AXY (30/01/14): "isnan" problem on HECTOR
2026                  !! if (fq0 /= fq0 ) then
2027                  if ( ieee_is_nan( fq0 ) ) then
2028                     !! there's a NaN here
2029                     if (lwp) write(numout,*) 'NAN detected in tra(', ji, ',', &
2030                     & jj, ',', jk, ',', jn, ') at time', kt
2031              CALL ctl_stop( 'trcbio_medusa, NAN in tra field' )
2032                  endif
2033               ENDDO
2034               CALL flush(numout)
2035# endif
2036
2037               !!----------------------------------------------------------------------
2038               !! Check model conservation
2039               !! these terms merely sum up the tendency terms of the relevant
2040               !! state variables, which should sum to zero; the iron cycle is
2041               !! complicated by fluxes that add (aeolian deposition and seafloor
2042               !! remineralisation) and remove (scavenging) dissolved iron from
2043               !! the model (i.e. the sum of iron fluxes is unlikely to be zero)
2044               !!----------------------------------------------------------------------
2045               !!
2046               !! 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)
2047               !! fsil0 = btra(ji,jj,jppds) + btra(ji,jj,jpsil)                              ! + ftempsi(ji,jj)
2048               !! ffer0 = (xrfn * fnit0) + btra(ji,jj,jpfer)
2049# if defined key_roam
2050               !! fcar0 = 0.
2051               !! falk0 = 0.
2052               !! foxy0 = 0.
2053# endif
2054               !!
2055               !! if (kt/240*240.eq.kt) then
2056               !!    if (ji.eq.2.and.jj.eq.2.and.jk.eq.1) then
2057               !!       IF (lwp) write (*,*) '*******!MEDUSA Conservation!*******',kt
2058# if defined key_roam
2059               !!       IF (lwp) write (*,*) fnit0,fsil0,ffer0,fcar0,falk0,foxy0
2060# else
2061               !!       IF (lwp) write (*,*) fnit0,fsil0,ffer0
2062# endif
2063               !!    endif
2064               !! endif     
2065
2066! MAYBE BUT A BREAK IN HERE - marc 20/4/17
2067! (this would make the previous section about 470 lines and the one below
2068! about 700 lines)
2069            ENDIF
2070         ENDDO
2071         ENDDO
2072
2073         DO jj = 2,jpjm1
2074         DO ji = 2,jpim1
2075            !! OPEN wet point IF..THEN loop
2076            if (tmask(ji,jj,jk) == 1) then
2077
2078# if defined key_trc_diabio
2079               !!======================================================================
2080               !! LOCAL GRID CELL DIAGNOSTICS
2081               !!======================================================================
2082               !!
2083               !!----------------------------------------------------------------------
2084               !! Full diagnostics key_trc_diabio:
2085               !! LOBSTER and PISCES support full diagnistics option key_trc_diabio   
2086               !! which gives an option of FULL output of biological sourses and sinks.
2087               !! I cannot see any reason for doing this. If needed, it can be done
2088               !! as shown below.
2089               !!----------------------------------------------------------------------
2090               !!
2091               IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio'
2092               !!               trbio(ji,jj,jk, 1) = fprn(ji,jj)
2093# endif
2094
2095               IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
2096         !!----------------------------------------------------------------------
2097         !! Add in XML diagnostics stuff
2098         !!----------------------------------------------------------------------
2099         !!
2100         !! ** 2D diagnostics
2101#   if defined key_debug_medusa
2102                  IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk loop'
2103                  CALL flush(numout)
2104#   endif
2105                  IF ( med_diag%PRN%dgsave ) THEN
2106                      fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
2107                  ENDIF
2108                  IF ( med_diag%MPN%dgsave ) THEN
2109                      fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn(ji,jj)         * fse3t(ji,jj,jk))
2110                  ENDIF
2111                  IF ( med_diag%PRD%dgsave ) THEN
2112                      fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
2113                  ENDIF
2114                  IF( med_diag%MPD%dgsave ) THEN
2115                      fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd(ji,jj)         * fse3t(ji,jj,jk)) 
2116                  ENDIF
2117                  !  IF( med_diag%DSED%dgsave ) THEN
2118                  !      CALL iom_put( "DSED"  , ftot_n )
2119                  !  ENDIF
2120                  IF( med_diag%OPAL%dgsave ) THEN
2121                      fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds(ji,jj) * zpds(ji,jj) * fse3t(ji,jj,jk)) 
2122                  ENDIF
2123                  IF( med_diag%OPALDISS%dgsave ) THEN
2124                      fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss(ji,jj)  * fse3t(ji,jj,jk)) 
2125                  ENDIF
2126                  IF( med_diag%GMIPn%dgsave ) THEN
2127                      fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn(ji,jj)  * fse3t(ji,jj,jk)) 
2128                  ENDIF
2129                  IF( med_diag%GMID%dgsave ) THEN
2130                      fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid(ji,jj)   * fse3t(ji,jj,jk)) 
2131                  ENDIF
2132                  IF( med_diag%MZMI%dgsave ) THEN
2133                      fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi(ji,jj)   * fse3t(ji,jj,jk)) 
2134                  ENDIF
2135                  IF( med_diag%GMEPN%dgsave ) THEN
2136                      fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn(ji,jj)  * fse3t(ji,jj,jk))
2137                  ENDIF
2138                  IF( med_diag%GMEPD%dgsave ) THEN
2139                      fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd(ji,jj)  * fse3t(ji,jj,jk)) 
2140                  ENDIF
2141                  IF( med_diag%GMEZMI%dgsave ) THEN
2142                      fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi(ji,jj) * fse3t(ji,jj,jk)) 
2143                  ENDIF
2144                  IF( med_diag%GMED%dgsave ) THEN
2145                      fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed(ji,jj)   * fse3t(ji,jj,jk)) 
2146                  ENDIF
2147                  IF( med_diag%MZME%dgsave ) THEN
2148                      fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme(ji,jj)   * fse3t(ji,jj,jk)) 
2149                  ENDIF
2150                  !  IF( med_diag%DEXP%dgsave ) THEN
2151                  !      CALL iom_put( "DEXP"  , ftot_n )
2152                  !  ENDIF
2153                  IF( med_diag%DETN%dgsave ) THEN
2154                      fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown(ji,jj)  * fse3t(ji,jj,jk)) 
2155                  ENDIF
2156                  IF( med_diag%MDET%dgsave ) THEN
2157                      fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd(ji,jj)     * fse3t(ji,jj,jk)) 
2158                  ENDIF
2159                  IF( med_diag%AEOLIAN%dgsave ) THEN
2160                      ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop(ji,jj)  * fse3t(ji,jj,jk)) 
2161                  ENDIF
2162                  IF( med_diag%BENTHIC%dgsave ) THEN
2163                      ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot(ji,jj)  * fse3t(ji,jj,jk)) 
2164                  ENDIF
2165                  IF( med_diag%SCAVENGE%dgsave ) THEN
2166                      ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav(ji,jj) * fse3t(ji,jj,jk)) 
2167                  ENDIF
2168                  IF( med_diag%PN_JLIM%dgsave ) THEN
2169                      ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))
2170                      fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) 
2171                  ENDIF
2172                  IF( med_diag%PN_NLIM%dgsave ) THEN
2173                      fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
2174                  ENDIF
2175                  IF( med_diag%PN_FELIM%dgsave ) THEN
2176                      ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln2(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
2177                  ENDIF
2178                  IF( med_diag%PD_JLIM%dgsave ) THEN
2179                      ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
2180                      fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) 
2181                  ENDIF
2182                  IF( med_diag%PD_NLIM%dgsave ) THEN
2183                      fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
2184                  ENDIF
2185                  IF( med_diag%PD_FELIM%dgsave ) THEN
2186                      ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
2187                  ENDIF
2188                  IF( med_diag%PD_SILIM%dgsave ) THEN
2189                      fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) 
2190                  ENDIF
2191                  IF( med_diag%PDSILIM2%dgsave ) THEN
2192                      fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
2193                  ENDIF
2194                  !!
2195                  IF( med_diag%TOTREG_N%dgsave ) THEN
2196                      fregen2d(ji,jj) = fregen2d(ji,jj) + fregen(ji,jj)
2197                  ENDIF
2198                  IF( med_diag%TOTRG_SI%dgsave ) THEN
2199                      fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi(ji,jj)
2200                  ENDIF
2201                  !!
2202                  IF( med_diag%FASTN%dgsave ) THEN
2203                      ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
2204                  ENDIF
2205                  IF( med_diag%FASTSI%dgsave ) THEN
2206                      ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
2207                  ENDIF
2208                  IF( med_diag%FASTFE%dgsave ) THEN
2209                      ftempfe2d(ji,jj) =ftempfe2d(ji,jj)  + (ftempfe(ji,jj) * fse3t(ji,jj,jk)) 
2210                  ENDIF
2211                  IF( med_diag%FASTC%dgsave ) THEN
2212                      ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
2213                  ENDIF
2214                  IF( med_diag%FASTCA%dgsave ) THEN
2215                      ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca(ji,jj) * fse3t(ji,jj,jk))
2216                  ENDIF
2217                  !!
2218                  IF( med_diag%REMINN%dgsave ) THEN
2219                      freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn(ji,jj)  * fse3t(ji,jj,jk))
2220                  ENDIF
2221                  IF( med_diag%REMINSI%dgsave ) THEN
2222                      freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
2223                  ENDIF
2224                  IF( med_diag%REMINFE%dgsave ) THEN
2225                      freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe(ji,jj) * fse3t(ji,jj,jk)) 
2226                  ENDIF
2227                  IF( med_diag%REMINC%dgsave ) THEN
2228                      freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc(ji,jj)  * fse3t(ji,jj,jk)) 
2229                  ENDIF
2230                  IF( med_diag%REMINCA%dgsave ) THEN
2231                      freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca(ji,jj) * fse3t(ji,jj,jk)) 
2232                  ENDIF
2233                  !!
2234# if defined key_roam
2235                  !!
2236                  !! AXY (09/11/16): CMIP6 diagnostics
2237                  IF( med_diag%FD_NIT3%dgsave ) THEN
2238                     fd_nit3(ji,jj,jk) = ffastn(ji,jj)
2239                  ENDIF
2240                  IF( med_diag%FD_SIL3%dgsave ) THEN
2241                     fd_sil3(ji,jj,jk) = ffastsi(ji,jj)
2242                  ENDIF
2243                  IF( med_diag%FD_CAR3%dgsave ) THEN
2244                     fd_car3(ji,jj,jk) = ffastc(ji,jj)
2245                  ENDIF
2246                  IF( med_diag%FD_CAL3%dgsave ) THEN
2247                     fd_cal3(ji,jj,jk) = ffastca(ji,jj)
2248                  ENDIF
2249                  !!
2250                  IF (jk.eq.i0100) THEN
2251                     IF( med_diag%RR_0100%dgsave ) THEN
2252                        ffastca2d(ji,jj) =   &
2253                        ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
2254                     ENDIF                     
2255                  ELSE IF (jk.eq.i0500) THEN
2256                     IF( med_diag%RR_0500%dgsave ) THEN
2257                        ffastca2d(ji,jj) =   &
2258                        ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
2259                     ENDIF                       
2260                  ELSE IF (jk.eq.i1000) THEN
2261                     IF( med_diag%RR_1000%dgsave ) THEN
2262                        ffastca2d(ji,jj) =   &
2263                        ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
2264                     ENDIF
2265                  ELSE IF (jk.eq.mbathy(ji,jj)) THEN
2266                     IF( med_diag%IBEN_N%dgsave ) THEN
2267                        iben_n2d(ji,jj) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
2268                     ENDIF
2269                     IF( med_diag%IBEN_FE%dgsave ) THEN
2270                        iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
2271                     ENDIF
2272                     IF( med_diag%IBEN_C%dgsave ) THEN
2273                        iben_c2d(ji,jj) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
2274                     ENDIF
2275                     IF( med_diag%IBEN_SI%dgsave ) THEN
2276                        iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
2277                     ENDIF
2278                     IF( med_diag%IBEN_CA%dgsave ) THEN
2279                        iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
2280                     ENDIF
2281                     IF( med_diag%OBEN_N%dgsave ) THEN
2282                        oben_n2d(ji,jj) = f_benout_n(ji,jj)
2283                     ENDIF
2284                     IF( med_diag%OBEN_FE%dgsave ) THEN
2285                        oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
2286                     ENDIF
2287                     IF( med_diag%OBEN_C%dgsave ) THEN
2288                        oben_c2d(ji,jj) = f_benout_c(ji,jj)
2289                     ENDIF
2290                     IF( med_diag%OBEN_SI%dgsave ) THEN
2291                        oben_si2d(ji,jj) = f_benout_si(ji,jj)
2292                     ENDIF
2293                     IF( med_diag%OBEN_CA%dgsave ) THEN
2294                        oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
2295                     ENDIF
2296                     IF( med_diag%SFR_OCAL%dgsave ) THEN
2297                        sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
2298                     ENDIF
2299                     IF( med_diag%SFR_OARG%dgsave ) THEN
2300                        sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
2301                     ENDIF
2302                     IF( med_diag%LYSO_CA%dgsave ) THEN
2303                        lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
2304                     ENDIF
2305                  ENDIF
2306                  !! end bathy-1 diags
2307                  !!
2308                  IF( med_diag%RIV_N%dgsave ) THEN
2309                     rivn2d(ji,jj) = rivn2d(ji,jj) +  (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
2310                  ENDIF
2311                  IF( med_diag%RIV_SI%dgsave ) THEN
2312                     rivsi2d(ji,jj) = rivsi2d(ji,jj) +  (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
2313                  ENDIF
2314                  IF( med_diag%RIV_C%dgsave ) THEN
2315                     rivc2d(ji,jj) = rivc2d(ji,jj) +  (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
2316                  ENDIF
2317                  IF( med_diag%RIV_ALK%dgsave ) THEN
2318                     rivalk2d(ji,jj) = rivalk2d(ji,jj) +  (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
2319                  ENDIF
2320                  IF( med_diag%DETC%dgsave ) THEN
2321                     fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc(ji,jj)  * fse3t(ji,jj,jk))   
2322                  ENDIF
2323                  !!
2324                  !!             
2325                  !!
2326                  IF( med_diag%PN_LLOSS%dgsave ) THEN
2327                     fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
2328                  ENDIF
2329                  IF( med_diag%PD_LLOSS%dgsave ) THEN
2330                     fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
2331                  ENDIF
2332                  IF( med_diag%ZI_LLOSS%dgsave ) THEN
2333                     fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
2334                  ENDIF
2335                  IF( med_diag%ZE_LLOSS%dgsave ) THEN
2336                     fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2(ji,jj) * fse3t(ji,jj,jk))
2337                  ENDIF
2338                  IF( med_diag%ZI_MES_N%dgsave ) THEN
2339                     zimesn2d(ji,jj) = zimesn2d(ji,jj) +  &
2340                     (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)) * fse3t(ji,jj,jk))
2341                  ENDIF
2342                  IF( med_diag%ZI_MES_D%dgsave ) THEN
2343                     zimesd2d(ji,jj) = zimesd2d(ji,jj) + & 
2344                     ((1. - xbetan) * finmi(ji,jj) * fse3t(ji,jj,jk))
2345                  ENDIF
2346                  IF( med_diag%ZI_MES_C%dgsave ) THEN
2347                     zimesc2d(ji,jj) = zimesc2d(ji,jj) + &
2348                     (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj)) * fse3t(ji,jj,jk))
2349                  ENDIF
2350                  IF( med_diag%ZI_MESDC%dgsave ) THEN
2351                     zimesdc2d(ji,jj) = zimesdc2d(ji,jj) + &
2352                     ((1. - xbetac) * ficmi(ji,jj) * fse3t(ji,jj,jk))
2353                  ENDIF
2354                  IF( med_diag%ZI_EXCR%dgsave ) THEN
2355                     ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +  (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
2356                  ENDIF
2357                  IF( med_diag%ZI_RESP%dgsave ) THEN
2358                     ziresp2d(ji,jj) = ziresp2d(ji,jj) +  (fmiresp(ji,jj) * fse3t(ji,jj,jk))
2359                  ENDIF
2360                  IF( med_diag%ZI_GROW%dgsave ) THEN
2361                     zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow(ji,jj) * fse3t(ji,jj,jk))
2362                  ENDIF
2363                  IF( med_diag%ZE_MES_N%dgsave ) THEN
2364                     zemesn2d(ji,jj) = zemesn2d(ji,jj) + &
2365                     (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)) * fse3t(ji,jj,jk))
2366                  ENDIF
2367                  IF( med_diag%ZE_MES_D%dgsave ) THEN
2368                     zemesd2d(ji,jj) = zemesd2d(ji,jj) + &
2369                     ((1. - xbetan) * finme(ji,jj) * fse3t(ji,jj,jk))
2370                  ENDIF
2371                  IF( med_diag%ZE_MES_C%dgsave ) THEN
2372                     zemesc2d(ji,jj) = zemesc2d(ji,jj) +                         & 
2373                     (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +  &
2374                     (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj)) * fse3t(ji,jj,jk))
2375                  ENDIF
2376                  IF( med_diag%ZE_MESDC%dgsave ) THEN
2377                     zemesdc2d(ji,jj) = zemesdc2d(ji,jj) +  &
2378                     ((1. - xbetac) * ficme(ji,jj) * fse3t(ji,jj,jk))
2379                  ENDIF
2380                  IF( med_diag%ZE_EXCR%dgsave ) THEN
2381                     zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
2382                  ENDIF
2383                  IF( med_diag%ZE_RESP%dgsave ) THEN
2384                     zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp(ji,jj) * fse3t(ji,jj,jk))
2385                  ENDIF
2386                  IF( med_diag%ZE_GROW%dgsave ) THEN
2387                     zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow(ji,jj) * fse3t(ji,jj,jk))
2388                  ENDIF
2389                  IF( med_diag%MDETC%dgsave ) THEN
2390                     mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc(ji,jj) * fse3t(ji,jj,jk))
2391                  ENDIF
2392                  IF( med_diag%GMIDC%dgsave ) THEN
2393                     gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc(ji,jj) * fse3t(ji,jj,jk))
2394                  ENDIF
2395                  IF( med_diag%GMEDC%dgsave ) THEN
2396                     gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc(ji,jj)  * fse3t(ji,jj,jk))
2397                  ENDIF
2398                  !!
2399# endif                   
2400                  !!
2401                  !! ** 3D diagnostics
2402                  IF( med_diag%TPP3%dgsave ) THEN
2403                     tpp3d(ji,jj,jk) =  (fprn(ji,jj) * zphn(ji,jj)) + (fprd(ji,jj) * zphd(ji,jj))
2404                     !CALL iom_put( "TPP3"  , tpp3d )
2405                  ENDIF
2406                  IF( med_diag%TPPD3%dgsave ) THEN
2407                     tppd3(ji,jj,jk) =  (fprd(ji,jj) * zphd(ji,jj))
2408                  ENDIF
2409                 
2410                  IF( med_diag%REMIN3N%dgsave ) THEN
2411                     remin3dn(ji,jj,jk) = fregen(ji,jj) + (freminn(ji,jj) * fse3t(ji,jj,jk)) !! remineralisation
2412                     !CALL iom_put( "REMIN3N"  , remin3dn )
2413                  ENDIF
2414                  !! IF( med_diag%PH3%dgsave ) THEN
2415                  !!     CALL iom_put( "PH3"  , f3_pH )
2416                  !! ENDIF
2417                  !! IF( med_diag%OM_CAL3%dgsave ) THEN
2418                  !!     CALL iom_put( "OM_CAL3"  , f3_omcal )
2419                  !! ENDIF
2420        !!
2421        !! AXY (09/11/16): CMIP6 diagnostics
2422        IF ( med_diag%DCALC3%dgsave   ) THEN
2423                     dcalc3(ji,jj,jk) = freminca(ji,jj)
2424                  ENDIF
2425        IF ( med_diag%FEDISS3%dgsave  ) THEN
2426                     fediss3(ji,jj,jk) = ffetop(ji,jj)
2427                  ENDIF
2428        IF ( med_diag%FESCAV3%dgsave  ) THEN
2429                     fescav3(ji,jj,jk) = ffescav(ji,jj)
2430                  ENDIF
2431        IF ( med_diag%MIGRAZP3%dgsave ) THEN
2432                     migrazp3(ji,jj,jk) = fgmipn(ji,jj) * xthetapn
2433                  ENDIF
2434        IF ( med_diag%MIGRAZD3%dgsave ) THEN
2435                     migrazd3(ji,jj,jk) = fgmidc(ji,jj)
2436                  ENDIF
2437        IF ( med_diag%MEGRAZP3%dgsave ) THEN
2438                     megrazp3(ji,jj,jk) = (fgmepn(ji,jj) * xthetapn) + (fgmepd(ji,jj) * xthetapd)
2439                  ENDIF
2440        IF ( med_diag%MEGRAZD3%dgsave ) THEN
2441                     megrazd3(ji,jj,jk) = fgmedc(ji,jj)
2442                  ENDIF
2443        IF ( med_diag%MEGRAZZ3%dgsave ) THEN
2444                     megrazz3(ji,jj,jk) = (fgmezmi(ji,jj) * xthetazmi)
2445                  ENDIF
2446        IF ( med_diag%PBSI3%dgsave    ) THEN
2447                     pbsi3(ji,jj,jk)    = (fprds(ji,jj) * zpds(ji,jj))
2448                  ENDIF
2449        IF ( med_diag%PCAL3%dgsave    ) THEN
2450                     pcal3(ji,jj,jk)    = ftempca(ji,jj)
2451                  ENDIF
2452        IF ( med_diag%REMOC3%dgsave   ) THEN
2453                     remoc3(ji,jj,jk)   = freminc(ji,jj)
2454                  ENDIF
2455        IF ( med_diag%PNLIMJ3%dgsave  ) THEN
2456                     ! pnlimj3(ji,jj,jk)  = fjln(ji,jj)
2457                     pnlimj3(ji,jj,jk)  = fjlim_pn(ji,jj)
2458                  ENDIF
2459        IF ( med_diag%PNLIMN3%dgsave  ) THEN
2460                     pnlimn3(ji,jj,jk)  = fnln(ji,jj)
2461                  ENDIF
2462        IF ( med_diag%PNLIMFE3%dgsave ) THEN
2463                     pnlimfe3(ji,jj,jk) = ffln2(ji,jj)
2464                  ENDIF
2465        IF ( med_diag%PDLIMJ3%dgsave  ) THEN
2466                     ! pdlimj3(ji,jj,jk)  = fjld(ji,jj)
2467                     pdlimj3(ji,jj,jk)  = fjlim_pd(ji,jj)
2468                  ENDIF
2469        IF ( med_diag%PDLIMN3%dgsave  ) THEN
2470                     pdlimn3(ji,jj,jk)  = fnld(ji,jj)
2471                  ENDIF
2472        IF ( med_diag%PDLIMFE3%dgsave ) THEN
2473                     pdlimfe3(ji,jj,jk) = ffld(ji,jj)
2474                  ENDIF
2475        IF ( med_diag%PDLIMSI3%dgsave ) THEN
2476                     pdlimsi3(ji,jj,jk) = fsld2(ji,jj)
2477                  ENDIF
2478                  !!
2479                  !! ** Without using iom_use
2480               ELSE IF( ln_diatrc ) THEN
2481#   if defined key_debug_medusa
2482                  IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk ln_diatrc'
2483                  CALL flush(numout)
2484#   endif
2485                  !!----------------------------------------------------------------------
2486                  !! Prepare 2D diagnostics
2487                  !!----------------------------------------------------------------------
2488                  !!
2489                  !! if ((kt / 240*240).eq.kt) then
2490                  !!    IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt
2491                  !! endif     
2492                  trc2d(ji,jj,1)  =  ftot_n(ji,jj)                             !! nitrogen inventory
2493                  trc2d(ji,jj,2)  =  ftot_si(ji,jj)                            !! silicon  inventory
2494                  trc2d(ji,jj,3)  =  ftot_fe(ji,jj)                            !! iron     inventory
2495                  trc2d(ji,jj,4)  = trc2d(ji,jj,4)  + (fprn(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom production
2496                  trc2d(ji,jj,5)  = trc2d(ji,jj,5)  + (fdpn(ji,jj)         * fse3t(ji,jj,jk))    !! non-diatom non-grazing losses
2497                  trc2d(ji,jj,6)  = trc2d(ji,jj,6)  + (fprd(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom production
2498                  trc2d(ji,jj,7)  = trc2d(ji,jj,7)  + (fdpd(ji,jj)         * fse3t(ji,jj,jk))    !! diatom non-grazing losses
2499                  !! diagnostic field  8 is (ostensibly) supplied by trcsed.F           
2500                  trc2d(ji,jj,9)  = trc2d(ji,jj,9)  + (fprds(ji,jj) * zpds(ji,jj) * fse3t(ji,jj,jk))    !! diatom silicon production
2501                  trc2d(ji,jj,10) = trc2d(ji,jj,10) + (fsdiss(ji,jj)  * fse3t(ji,jj,jk))         !! diatom silicon dissolution
2502                  trc2d(ji,jj,11) = trc2d(ji,jj,11) + (fgmipn(ji,jj)  * fse3t(ji,jj,jk))         !! microzoo grazing on non-diatoms
2503                  trc2d(ji,jj,12) = trc2d(ji,jj,12) + (fgmid(ji,jj)   * fse3t(ji,jj,jk))         !! microzoo grazing on detrital nitrogen
2504                  trc2d(ji,jj,13) = trc2d(ji,jj,13) + (fdzmi(ji,jj)   * fse3t(ji,jj,jk))         !! microzoo non-grazing losses
2505                  trc2d(ji,jj,14) = trc2d(ji,jj,14) + (fgmepn(ji,jj)  * fse3t(ji,jj,jk))         !! mesozoo  grazing on non-diatoms
2506                  trc2d(ji,jj,15) = trc2d(ji,jj,15) + (fgmepd(ji,jj)  * fse3t(ji,jj,jk))         !! mesozoo  grazing on diatoms
2507                  trc2d(ji,jj,16) = trc2d(ji,jj,16) + (fgmezmi(ji,jj) * fse3t(ji,jj,jk))         !! mesozoo  grazing on microzoo
2508                  trc2d(ji,jj,17) = trc2d(ji,jj,17) + (fgmed(ji,jj)   * fse3t(ji,jj,jk))         !! mesozoo  grazing on detrital nitrogen
2509                  trc2d(ji,jj,18) = trc2d(ji,jj,18) + (fdzme(ji,jj)   * fse3t(ji,jj,jk))         !! mesozoo  non-grazing losses
2510                  !! diagnostic field 19 is (ostensibly) supplied by trcexp.F
2511                  trc2d(ji,jj,20) = trc2d(ji,jj,20) + (fslown(ji,jj)  * fse3t(ji,jj,jk))         !! slow sinking detritus N production
2512                  trc2d(ji,jj,21) = trc2d(ji,jj,21) + (fdd(ji,jj)     * fse3t(ji,jj,jk))         !! detrital remineralisation
2513                  trc2d(ji,jj,22) = trc2d(ji,jj,22) + (ffetop(ji,jj)  * fse3t(ji,jj,jk))         !! aeolian  iron addition
2514                  trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot(ji,jj)  * fse3t(ji,jj,jk))         !! seafloor iron addition
2515                  trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav(ji,jj) * fse3t(ji,jj,jk))         !! "free" iron scavenging
2516                  trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) !! non-diatom J  limitation term
2517                  trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom N  limitation term
2518                  trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln2(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom Fe limitation term
2519                  trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) !! diatom     J  limitation term
2520                  trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     N  limitation term
2521                  trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Fe limitation term
2522                  trc2d(ji,jj,31) = trc2d(ji,jj,31) + (fsld2(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Si limitation term
2523                  trc2d(ji,jj,32) = trc2d(ji,jj,32) + (fsld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Si uptake limitation term
2524                  if (jk.eq.i0100) trc2d(ji,jj,33) = fslownflux(ji,jj)         !! slow detritus flux at  100 m
2525                  if (jk.eq.i0200) trc2d(ji,jj,34) = fslownflux(ji,jj)         !! slow detritus flux at  200 m
2526                  if (jk.eq.i0500) trc2d(ji,jj,35) = fslownflux(ji,jj)         !! slow detritus flux at  500 m
2527                  if (jk.eq.i1000) trc2d(ji,jj,36) = fslownflux(ji,jj)         !! slow detritus flux at 1000 m
2528                  trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen(ji,jj)                   !! non-fast N  full column regeneration
2529                  trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi(ji,jj)                 !! non-fast Si full column regeneration
2530                  if (jk.eq.i0100) trc2d(ji,jj,39) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  100 m
2531                  if (jk.eq.i0200) trc2d(ji,jj,40) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  200 m
2532                  if (jk.eq.i0500) trc2d(ji,jj,41) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  500 m
2533                  if (jk.eq.i1000) trc2d(ji,jj,42) = trc2d(ji,jj,37)           !! non-fast N  regeneration to 1000 m
2534                  trc2d(ji,jj,43) = trc2d(ji,jj,43) + (ftempn(ji,jj)  * fse3t(ji,jj,jk))         !! fast sinking detritus N production
2535                  trc2d(ji,jj,44) = trc2d(ji,jj,44) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus Si production
2536                  trc2d(ji,jj,45) = trc2d(ji,jj,45) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus Fe production
2537                  trc2d(ji,jj,46) = trc2d(ji,jj,46) + (ftempc(ji,jj)  * fse3t(ji,jj,jk))         !! fast sinking detritus C production
2538                  trc2d(ji,jj,47) = trc2d(ji,jj,47) + (ftempca(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus CaCO3 production
2539                  if (jk.eq.i0100) trc2d(ji,jj,48) = ffastn(ji,jj)             !! fast detritus N  flux at  100 m
2540                  if (jk.eq.i0200) trc2d(ji,jj,49) = ffastn(ji,jj)             !! fast detritus N  flux at  200 m
2541                  if (jk.eq.i0500) trc2d(ji,jj,50) = ffastn(ji,jj)             !! fast detritus N  flux at  500 m
2542                  if (jk.eq.i1000) trc2d(ji,jj,51) = ffastn(ji,jj)             !! fast detritus N  flux at 1000 m
2543                  if (jk.eq.i0100) trc2d(ji,jj,52) = fregenfast(ji,jj)         !! N  regeneration to  100 m
2544                  if (jk.eq.i0200) trc2d(ji,jj,53) = fregenfast(ji,jj)         !! N  regeneration to  200 m
2545                  if (jk.eq.i0500) trc2d(ji,jj,54) = fregenfast(ji,jj)         !! N  regeneration to  500 m
2546                  if (jk.eq.i1000) trc2d(ji,jj,55) = fregenfast(ji,jj)         !! N  regeneration to 1000 m
2547                  if (jk.eq.i0100) trc2d(ji,jj,56) = ffastsi(ji,jj)            !! fast detritus Si flux at  100 m
2548                  if (jk.eq.i0200) trc2d(ji,jj,57) = ffastsi(ji,jj)            !! fast detritus Si flux at  200 m
2549                  if (jk.eq.i0500) trc2d(ji,jj,58) = ffastsi(ji,jj)            !! fast detritus Si flux at  500 m
2550                  if (jk.eq.i1000) trc2d(ji,jj,59) = ffastsi(ji,jj)            !! fast detritus Si flux at 1000 m
2551                  if (jk.eq.i0100) trc2d(ji,jj,60) = fregenfastsi(ji,jj)       !! Si regeneration to  100 m
2552                  if (jk.eq.i0200) trc2d(ji,jj,61) = fregenfastsi(ji,jj)       !! Si regeneration to  200 m
2553                  if (jk.eq.i0500) trc2d(ji,jj,62) = fregenfastsi(ji,jj)       !! Si regeneration to  500 m
2554                  if (jk.eq.i1000) trc2d(ji,jj,63) = fregenfastsi(ji,jj)       !! Si regeneration to 1000 m
2555                  trc2d(ji,jj,64) = trc2d(ji,jj,64) + (freminn(ji,jj)  * fse3t(ji,jj,jk))        !! sum of fast-sinking N  fluxes
2556                  trc2d(ji,jj,65) = trc2d(ji,jj,65) + (freminsi(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Si fluxes
2557                  trc2d(ji,jj,66) = trc2d(ji,jj,66) + (freminfe(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Fe fluxes
2558                  trc2d(ji,jj,67) = trc2d(ji,jj,67) + (freminc(ji,jj)  * fse3t(ji,jj,jk))        !! sum of fast-sinking C  fluxes
2559                  trc2d(ji,jj,68) = trc2d(ji,jj,68) + (freminca(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Ca fluxes
2560                  if (jk.eq.mbathy(ji,jj)) then
2561                     trc2d(ji,jj,69) = fsedn(ji,jj)                                   !! N  sedimentation flux                                 
2562                     trc2d(ji,jj,70) = fsedsi(ji,jj)                                  !! Si sedimentation flux
2563                     trc2d(ji,jj,71) = fsedfe(ji,jj)                                  !! Fe sedimentation flux
2564                     trc2d(ji,jj,72) = fsedc(ji,jj)                                   !! C  sedimentation flux
2565                     trc2d(ji,jj,73) = fsedca(ji,jj)                                  !! Ca sedimentation flux
2566                  endif
2567                  if (jk.eq.1)  trc2d(ji,jj,74) = qsr(ji,jj)
2568                  if (jk.eq.1)  trc2d(ji,jj,75) = xpar(ji,jj,jk)
2569                  !! if (jk.eq.1)  trc2d(ji,jj,75) = real(iters(ji,jj))
2570                  !! diagnostic fields 76 to 80 calculated below
2571                  trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj)           !! mixed layer non-diatom production
2572                  trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj)           !! mixed layer     diatom production
2573# if defined key_gulf_finland
2574                  if (jk.eq.1)  trc2d(ji,jj,83) = real(ibio_switch)            !! Gulf of Finland check
2575# else
2576                  trc2d(ji,jj,83) = ocal_ccd(ji,jj)                            !! calcite CCD depth
2577# endif
2578                  trc2d(ji,jj,84) = fccd(ji,jj)                                !! last model level above calcite CCD depth
2579                  if (jk.eq.1)     trc2d(ji,jj,85) = xFree(ji,jj)              !! surface "free" iron
2580                  if (jk.eq.i0200) trc2d(ji,jj,86) = xFree(ji,jj)              !! "free" iron at  100 m
2581                  if (jk.eq.i0200) trc2d(ji,jj,87) = xFree(ji,jj)              !! "free" iron at  200 m
2582                  if (jk.eq.i0500) trc2d(ji,jj,88) = xFree(ji,jj)              !! "free" iron at  500 m
2583                  if (jk.eq.i1000) trc2d(ji,jj,89) = xFree(ji,jj)              !! "free" iron at 1000 m
2584                  !! AXY (27/06/12): extract "euphotic depth"
2585                  if (jk.eq.1)     trc2d(ji,jj,90) = xze(ji,jj)
2586                  !!
2587# if defined key_roam
2588                  !! ROAM provisionally has access to a further 20 2D diagnostics
2589                  if (jk .eq. 1) then
2590                     trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + wndm(ji,jj)              !! surface wind
2591                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm(ji,jj)           !! atmospheric pCO2
2592                     trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph(ji,jj)                !! ocean pH
2593                     trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w(ji,jj)             !! ocean pCO2
2594                     trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3(ji,jj)             !! ocean H2CO3 conc.
2595                     trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3(ji,jj)              !! ocean HCO3 conc.
2596                     trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3(ji,jj)               !! ocean CO3 conc.
2597                     trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux(ji,jj)           !! air-sea CO2 flux
2598                     trc2d(ji,jj,99)  = trc2d(ji,jj,99)  + f_omcal(ji,jj)      !! ocean omega calcite
2599                     trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj)      !! ocean omega aragonite
2600                     trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC(ji,jj)              !! ocean TDIC
2601                     trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK(ji,jj)              !! ocean TALK
2602                     trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660(ji,jj)             !! surface kw660
2603                     trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0(ji,jj)               !! surface pressure
2604                     trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux(ji,jj)            !! air-sea O2 flux
2605                     trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat(ji,jj)             !! ocean O2 saturation
2606                     trc2d(ji,jj,107) = f2_ccd_cal(ji,jj)                      !! depth calcite CCD
2607                     trc2d(ji,jj,108) = f2_ccd_arg(ji,jj)                      !! depth aragonite CCD
2608                  endif
2609                  if (jk .eq. mbathy(ji,jj)) then
2610                     trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)                     !! seafloor omega calcite
2611                     trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)                     !! seafloor omega aragonite
2612                  endif
2613                  !! diagnostic fields 111 to 117 calculated below
2614                  if (jk.eq.i0100) trc2d(ji,jj,118) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at  100 m
2615                  if (jk.eq.i0500) trc2d(ji,jj,119) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at  500 m
2616                  if (jk.eq.i1000) trc2d(ji,jj,120) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at 1000 m
2617                  !! AXY (18/01/12): benthic flux diagnostics
2618                  if (jk.eq.mbathy(ji,jj)) then
2619                     trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
2620                     trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
2621                     trc2d(ji,jj,123) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
2622                     trc2d(ji,jj,124) = f_fbenin_si(ji,jj)
2623                     trc2d(ji,jj,125) = f_fbenin_ca(ji,jj)
2624                     trc2d(ji,jj,126) = f_benout_n(ji,jj)
2625                     trc2d(ji,jj,127) = f_benout_fe(ji,jj)
2626                     trc2d(ji,jj,128) = f_benout_c(ji,jj)
2627                     trc2d(ji,jj,129) = f_benout_si(ji,jj)
2628                     trc2d(ji,jj,130) = f_benout_ca(ji,jj)
2629                  endif
2630                  !! diagnostics fields 131 to 135 calculated below
2631                  trc2d(ji,jj,136) = f_runoff(ji,jj)
2632                  !! AXY (19/07/12): amended to allow for riverine nutrient addition below surface
2633                  trc2d(ji,jj,137) = trc2d(ji,jj,137) + (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
2634                  trc2d(ji,jj,138) = trc2d(ji,jj,138) + (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
2635                  trc2d(ji,jj,139) = trc2d(ji,jj,139) + (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
2636                  trc2d(ji,jj,140) = trc2d(ji,jj,140) + (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
2637                  trc2d(ji,jj,141) = trc2d(ji,jj,141) + (fslowc(ji,jj)  * fse3t(ji,jj,jk))       !! slow sinking detritus C production
2638                  if (jk.eq.i0100) trc2d(ji,jj,142) = fslowcflux(ji,jj)        !! slow detritus flux at  100 m
2639                  if (jk.eq.i0200) trc2d(ji,jj,143) = fslowcflux(ji,jj)        !! slow detritus flux at  200 m
2640                  if (jk.eq.i0500) trc2d(ji,jj,144) = fslowcflux(ji,jj)        !! slow detritus flux at  500 m
2641                  if (jk.eq.i1000) trc2d(ji,jj,145) = fslowcflux(ji,jj)        !! slow detritus flux at 1000 m
2642                  trc2d(ji,jj,146)  = trc2d(ji,jj,146)  + ftot_c(ji,jj)        !! carbon     inventory
2643                  trc2d(ji,jj,147)  = trc2d(ji,jj,147)  + ftot_a(ji,jj)        !! alkalinity inventory
2644                  trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)       !! oxygen     inventory
2645                  if (jk.eq.mbathy(ji,jj)) then
2646                     trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
2647                  endif
2648                  trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)                  !! community respiration
2649        !!
2650        !! AXY (14/02/14): a Valentines Day gift to BASIN - a shedload of new
2651                  !!                 diagnostics that they'll most likely never need!
2652                  !!                 (actually, as with all such gifts, I'm giving them
2653                  !!                 some things I'd like myself!)
2654                  !!
2655                  !! ----------------------------------------------------------------------
2656                  !! linear losses
2657                  !! non-diatom
2658                  trc2d(ji,jj,151) = trc2d(ji,jj,151) + (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
2659                  !! diatom
2660                  trc2d(ji,jj,152) = trc2d(ji,jj,152) + (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
2661                  !! microzooplankton
2662                  trc2d(ji,jj,153) = trc2d(ji,jj,153) + (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
2663                  !! mesozooplankton
2664                  trc2d(ji,jj,154) = trc2d(ji,jj,154) + (fdzme2(ji,jj) * fse3t(ji,jj,jk))
2665                  !! ----------------------------------------------------------------------
2666                  !! microzooplankton grazing
2667                  !! microzooplankton messy -> N
2668                  trc2d(ji,jj,155) = trc2d(ji,jj,155) + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)) * fse3t(ji,jj,jk))
2669                  !! microzooplankton messy -> D
2670                  trc2d(ji,jj,156) = trc2d(ji,jj,156) + ((1. - xbetan) * finmi(ji,jj) * fse3t(ji,jj,jk))
2671                  !! microzooplankton messy -> DIC
2672                  trc2d(ji,jj,157) = trc2d(ji,jj,157) + (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj)) * fse3t(ji,jj,jk))
2673                  !! microzooplankton messy -> Dc
2674                  trc2d(ji,jj,158) = trc2d(ji,jj,158) + ((1. - xbetac) * ficmi(ji,jj) * fse3t(ji,jj,jk))
2675                  !! microzooplankton excretion
2676                  trc2d(ji,jj,159) = trc2d(ji,jj,159) + (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
2677                  !! microzooplankton respiration
2678                  trc2d(ji,jj,160) = trc2d(ji,jj,160) + (fmiresp(ji,jj) * fse3t(ji,jj,jk))
2679                  !! microzooplankton growth
2680                  trc2d(ji,jj,161) = trc2d(ji,jj,161) + (fmigrow(ji,jj) * fse3t(ji,jj,jk))
2681                  !! ----------------------------------------------------------------------
2682                  !! mesozooplankton grazing
2683                  !! mesozooplankton messy -> N
2684                  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))
2685                  !! mesozooplankton messy -> D
2686                  trc2d(ji,jj,163) = trc2d(ji,jj,163) + ((1. - xbetan) * finme(ji,jj) * fse3t(ji,jj,jk))
2687                  !! mesozooplankton messy -> DIC
2688                  trc2d(ji,jj,164) = trc2d(ji,jj,164) + (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) + &
2689                  &                  (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj)) * fse3t(ji,jj,jk))
2690                  !! mesozooplankton messy -> Dc
2691                  trc2d(ji,jj,165) = trc2d(ji,jj,165) + ((1. - xbetac) * ficme(ji,jj) * fse3t(ji,jj,jk))
2692                  !! mesozooplankton excretion
2693                  trc2d(ji,jj,166) = trc2d(ji,jj,166) + (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
2694                  !! mesozooplankton respiration
2695                  trc2d(ji,jj,167) = trc2d(ji,jj,167) + (fmeresp(ji,jj) * fse3t(ji,jj,jk))
2696                  !! mesozooplankton growth
2697                  trc2d(ji,jj,168) = trc2d(ji,jj,168) + (fmegrow(ji,jj) * fse3t(ji,jj,jk))
2698                  !! ----------------------------------------------------------------------
2699                  !! miscellaneous
2700                  trc2d(ji,jj,169) = trc2d(ji,jj,169) + (fddc(ji,jj)    * fse3t(ji,jj,jk)) !! detrital C remineralisation
2701                  trc2d(ji,jj,170) = trc2d(ji,jj,170) + (fgmidc(ji,jj)  * fse3t(ji,jj,jk)) !! microzoo grazing on detrital carbon
2702                  trc2d(ji,jj,171) = trc2d(ji,jj,171) + (fgmedc(ji,jj)  * fse3t(ji,jj,jk)) !! mesozoo  grazing on detrital carbon
2703                  !!
2704                  !! ----------------------------------------------------------------------
2705        !!
2706        !! AXY (23/10/14): extract primary production related surface fields to
2707        !!                 deal with diel cycle issues; hijacking BASIN 150m
2708        !!                 diagnostics to do so (see commented out diagnostics
2709        !!                 below this section)
2710        !!
2711                  !! extract relevant BASIN fields at 150m
2712                  if (jk .eq. i0150) then
2713                     trc2d(ji,jj,172) = trc2d(ji,jj,4)    !! Pn PP
2714                     trc2d(ji,jj,173) = trc2d(ji,jj,151)  !! Pn linear loss
2715                     trc2d(ji,jj,174) = trc2d(ji,jj,5)    !! Pn non-linear loss
2716                     trc2d(ji,jj,175) = trc2d(ji,jj,11)   !! Pn grazing to Zmi
2717                     trc2d(ji,jj,176) = trc2d(ji,jj,14)   !! Pn grazing to Zme
2718                     trc2d(ji,jj,177) = trc2d(ji,jj,6)    !! Pd PP
2719                     trc2d(ji,jj,178) = trc2d(ji,jj,152)  !! Pd linear loss
2720                     trc2d(ji,jj,179) = trc2d(ji,jj,7)    !! Pd non-linear loss
2721                     trc2d(ji,jj,180) = trc2d(ji,jj,15)   !! Pd grazing to Zme
2722                     trc2d(ji,jj,181) = trc2d(ji,jj,12)   !! Zmi grazing on D
2723                     trc2d(ji,jj,182) = trc2d(ji,jj,170)  !! Zmi grazing on Dc
2724                     trc2d(ji,jj,183) = trc2d(ji,jj,155)  !! Zmi messy feeding loss to N
2725                     trc2d(ji,jj,184) = trc2d(ji,jj,156)  !! Zmi messy feeding loss to D
2726                     trc2d(ji,jj,185) = trc2d(ji,jj,157)  !! Zmi messy feeding loss to DIC
2727                     trc2d(ji,jj,186) = trc2d(ji,jj,158)  !! Zmi messy feeding loss to Dc
2728                     trc2d(ji,jj,187) = trc2d(ji,jj,159)  !! Zmi excretion
2729                     trc2d(ji,jj,188) = trc2d(ji,jj,160)  !! Zmi respiration
2730                     trc2d(ji,jj,189) = trc2d(ji,jj,161)  !! Zmi growth
2731                     trc2d(ji,jj,190) = trc2d(ji,jj,153)  !! Zmi linear loss
2732                     trc2d(ji,jj,191) = trc2d(ji,jj,13)   !! Zmi non-linear loss
2733                     trc2d(ji,jj,192) = trc2d(ji,jj,16)   !! Zmi grazing to Zme
2734                     trc2d(ji,jj,193) = trc2d(ji,jj,17)   !! Zme grazing on D
2735                     trc2d(ji,jj,194) = trc2d(ji,jj,171)  !! Zme grazing on Dc
2736                     trc2d(ji,jj,195) = trc2d(ji,jj,162)  !! Zme messy feeding loss to N
2737                     trc2d(ji,jj,196) = trc2d(ji,jj,163)  !! Zme messy feeding loss to D
2738                     trc2d(ji,jj,197) = trc2d(ji,jj,164)  !! Zme messy feeding loss to DIC
2739                     trc2d(ji,jj,198) = trc2d(ji,jj,165)  !! Zme messy feeding loss to Dc
2740                     trc2d(ji,jj,199) = trc2d(ji,jj,166)  !! Zme excretion
2741                     trc2d(ji,jj,200) = trc2d(ji,jj,167)  !! Zme respiration
2742                     trc2d(ji,jj,201) = trc2d(ji,jj,168)  !! Zme growth
2743                     trc2d(ji,jj,202) = trc2d(ji,jj,154)  !! Zme linear loss
2744                     trc2d(ji,jj,203) = trc2d(ji,jj,18)   !! Zme non-linear loss
2745                     trc2d(ji,jj,204) = trc2d(ji,jj,20)   !! Slow detritus production, N
2746                     trc2d(ji,jj,205) = trc2d(ji,jj,21)   !! Slow detritus remineralisation, N
2747                     trc2d(ji,jj,206) = trc2d(ji,jj,141)  !! Slow detritus production, C
2748                     trc2d(ji,jj,207) = trc2d(ji,jj,169)  !! Slow detritus remineralisation, C
2749                     trc2d(ji,jj,208) = trc2d(ji,jj,43)   !! Fast detritus production, N
2750                     trc2d(ji,jj,209) = trc2d(ji,jj,21)   !! Fast detritus remineralisation, N
2751                     trc2d(ji,jj,210) = trc2d(ji,jj,64)   !! Fast detritus production, C
2752                     trc2d(ji,jj,211) = trc2d(ji,jj,67)   !! Fast detritus remineralisation, C
2753                     trc2d(ji,jj,212) = trc2d(ji,jj,150)  !! Community respiration
2754                     trc2d(ji,jj,213) = fslownflux(ji,jj) !! Slow detritus N flux at 150 m
2755                     trc2d(ji,jj,214) = fslowcflux(ji,jj) !! Slow detritus C flux at 150 m
2756                     trc2d(ji,jj,215) = ffastn(ji,jj)     !! Fast detritus N flux at 150 m
2757                     trc2d(ji,jj,216) = ffastc(ji,jj)     !! Fast detritus C flux at 150 m
2758                  endif
2759                  !!
2760                  !! Jpalm (11-08-2014)
2761                  !! Add UKESM1 diagnoatics
2762                  !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2763                  if ((jk .eq. 1) .and.( jdms.eq.1)) then
2764                     trc2d(ji,jj,221) = dms_surf(ji,jj)          !! DMS surface concentration
2765                     !! AXY (13/03/15): add in other DMS estimates
2766                     trc2d(ji,jj,222) = dms_andr(ji,jj)          !! DMS surface concentration
2767                     trc2d(ji,jj,223) = dms_simo(ji,jj)          !! DMS surface concentration
2768                     trc2d(ji,jj,224) = dms_aran(ji,jj)          !! DMS surface concentration
2769                     trc2d(ji,jj,225) = dms_hall(ji,jj)          !! DMS surface concentration
2770                  endif
2771# endif
2772                  !! other possible future diagnostics include:
2773                  !!   - integrated tracer values (esp. biological)
2774                  !!   - mixed layer tracer values
2775                  !!   - sub-surface chlorophyll maxima (plus depth)
2776                  !!   - different mixed layer depth criteria (T, sigma, var. sigma)
2777
2778                  !!----------------------------------------------------------------------
2779                  !! Prepare 3D diagnostics
2780                  !!----------------------------------------------------------------------
2781                  !!
2782                  trc3d(ji,jj,jk,1)  = ((fprn(ji,jj) + fprd(ji,jj)) * zphn(ji,jj))     !! primary production 
2783                  trc3d(ji,jj,jk,2)  = fslownflux(ji,jj) + ffastn(ji,jj) !! detrital flux
2784                  trc3d(ji,jj,jk,3)  = fregen(ji,jj) + (freminn(ji,jj) * fse3t(ji,jj,jk))  !! remineralisation
2785# if defined key_roam
2786                  trc3d(ji,jj,jk,4)  = f3_pH(ji,jj,jk)            !! pH
2787                  trc3d(ji,jj,jk,5)  = f3_omcal(ji,jj,jk)         !! omega calcite
2788# else
2789                  trc3d(ji,jj,jk,4)  = ffastsi(ji,jj)             !! fast Si flux
2790# endif
2791             ENDIF   ! end of ln_diatrc option
2792             !! CLOSE wet point IF..THEN loop
2793            endif
2794         !! CLOSE horizontal loops
2795         ENDDO
2796         ENDDO
2797         !!
2798             IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
2799
2800                !!-------------------------------------------------------
2801                !! 2d specific k level diags
2802                !!-------------------------------------------------------
2803                CALL bio_medusa_diag_slice( jk )
2804
2805              ENDIF
2806      !! CLOSE vertical loop
2807      ENDDO
2808
2809      !!----------------------------------------------------------------------
2810      !! Final calculations for diagnostics
2811      !!----------------------------------------------------------------------
2812      CALL bio_medusa_fin( kt )
2813
2814# if defined key_trc_diabio
2815       !! Lateral boundary conditions on trcbio
2816       DO jn=1,jp_medusa_trd
2817          CALL lbc_lnk(trbio(:,:,1,jn),'T',1. )
2818       ENDDO 
2819# endif
2820
2821# if defined key_debug_medusa
2822       IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
2823       CALL flush(numout)
2824# endif
2825
2826   END SUBROUTINE trc_bio_medusa
2827
2828#else
2829   !!======================================================================
2830   !!  Dummy module :                                   No MEDUSA bio-model
2831   !!======================================================================
2832CONTAINS
2833   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
2834      INTEGER, INTENT( in ) ::   kt
2835      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
2836   END SUBROUTINE trc_bio_medusa
2837#endif 
2838
2839   !!======================================================================
2840END MODULE  trcbio_medusa
Note: See TracBrowser for help on using the repository browser.