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

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

Pulled the updating of tracers out of trcbio_medusa.F90

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