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

Context Navigation

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

Last change on this file since 7958 was 7958, checked in by marc, 7 years ago
Pulled air-sea gas exchange and river inputs from trcbio_medusa.F90 into air_sea.F90
File size: 206.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: