Context Navigation

trcbio_medusa.F90 @ 8224

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