New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

trcbio_medusa.F90 in branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

Context Navigation

source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 7894

Last change on this file since 7894 was 7894, checked in by jpalmier, 7 years ago
JPALM -- 11-04-2017 -- MEDUSA spring tidy-up refreshning session
File size: 242.0 KB

Line
1	MODULE trcbio_medusa
2	!!======================================================================
3	!! * MODULE trcbio *
4	!! TOP : MEDUSA
5	!!======================================================================
6	!! History :
7	!! - ! 1999-07 (M. Levy) original code
8	!! - ! 2000-12 (E. Kestenare) assign parameters to name individual tracers
9	!! - ! 2001-03 (M. Levy) LNO3 + dia2d
10	!! 2.0 ! 2007-12 (C. Deltel, G. Madec) F90
11	!! - ! 2008-08 (K. Popova) adaptation for MEDUSA
12	!! - ! 2008-11 (A. Yool) continuing adaptation for MEDUSA
13	!! - ! 2010-03 (A. Yool) updated for branch inclusion
14	!! - ! 2011-08 (A. Yool) updated for ROAM (see below)
15	!! - ! 2013-03 (A. Yool) updated for iMARNET
16	!! - ! 2013-05 (A. Yool) updated for v3.5
17	!! - ! 2014-08 (A. Yool, J. Palm) Add DMS module for UKESM1 model
18	!! - ! 2015-06 (A. Yool) Update to include MOCSY
19	!! - ! 2015-07 (A. Yool) Update for rolling averages
20	!! - ! 2015-10 (J. Palm) Update for diag outputs through iom_use
21	!! - ! 2016-11 (A. Yool) Updated diags for CMIP6
22	!! - ! 2017-03 (A. Yool) Updated DMS for DIN limitation
23	!! - ! 2017-04 (A. Yool) Simplify code to remove unused options, etc.
24	!! - remove ln_diatrc, etc. code
25	!! - remove PML carbonate chemistry code
26	!! - remove defunct iron scavenging code
27	!! - remove defunct debug diagnostic code
28	!!----------------------------------------------------------------------
29	!!
30	#if defined key_roam
31	!!----------------------------------------------------------------------
32	!! Updates for the ROAM project include:
33	!! - addition of DIC, alkalinity, detrital carbon and oxygen tracers
34	!! - addition of air-sea fluxes of CO2 and oxygen (updated with MOCSY)
35	!! - periodic (monthly) calculation of full 3D carbonate chemistry
36	!! - detrital C:N ratio now free to evolve dynamically
37	!! - benthic storage pools
38	!!
39	!! Opportunity also taken to add functionality:
40	!! - switch for Liebig Law (= most-limiting) nutrient uptake
41	!! - switch for accelerated seafloor detritus remineralisation
42	!! - switch for fast -> slow detritus transfer at seafloor
43	!! - switch for ballast vs. Martin vs. Henson fast detritus remin.
44	!! - per GMD referee remarks, xfdfrac3 introduced for grazed PDS
45	!!
46	!! Updates with the addition of MOCSY include:
47	!! - option to use PML or MOCSY carbonate chemistry (the latter is
48	!! preferred)
49	!! - central calculation of gas transfer velocity, f_kw660; previously
50	!! this was done separately for CO2 and O2 with predictable results
51	!! - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux
52	!! calculations and to those for O2 air-sea flux
53	!! - extra diagnostics included for MOCSY
54	!!----------------------------------------------------------------------
55	#endif
56	!!
57	#if defined key_medusa
58	!!----------------------------------------------------------------------
59	!! MEDUSA bio-model
60	!!----------------------------------------------------------------------
61	!! trc_bio_medusa :
62	!!----------------------------------------------------------------------
63	USE oce_trc
64	USE trc
65	USE sms_medusa
66	USE lbclnk
67	USE prtctl_trc ! Print control for debugging
68	USE trcsed_medusa
69	USE sbc_oce ! surface forcing
70	USE sbcrnf ! surface boundary condition: runoff variables
71	USE in_out_manager ! I/O manager
72	# if defined key_iomput
73	USE iom
74	USE trcnam_medusa ! JPALM 13-11-2015 -- if iom_use for diag
75	!!USE trc_nam_iom_medusa ! JPALM 13-11-2015 -- if iom_use for diag
76	# endif
77	# if defined key_roam
78	USE gastransfer
79	USE mocsy_wrapper
80	# endif
81	USE trcoxy_medusa
82	!! Jpalm (08/08/2014)
83	USE trcdms_medusa
84	# endif
85	!! AXY (18/01/12): brought in for benthic timestepping
86	USE trcnam_trp ! AXY (24/05/2013)
87	USE trdmxl_trc
88	USE trdtrc_oce ! AXY (24/05/2013)
89
90	!! AXY (30/01/14): necessary to find NaNs on HECTOR
91	USE, INTRINSIC :: ieee_arithmetic
92
93	!! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
94	USE sbc_oce, ONLY: lk_oasis
95	USE oce, ONLY: CO2Flux_out_cpl, DMS_out_cpl, PCO2a_in_cpl
96
97	IMPLICIT NONE
98	PRIVATE
99
100	PUBLIC trc_bio_medusa ! called in ???
101
102	!!* Substitution
103	# include "domzgr_substitute.h90"
104	!!----------------------------------------------------------------------
105	!! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
106	!! $Id$
107	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
108	!!----------------------------------------------------------------------
109
110	CONTAINS
111
112	SUBROUTINE trc_bio_medusa( kt )
113	!!---------------------------------------------------------------------
114	!! * ROUTINE trc_bio *
115	!!
116	!! ** Purpose : compute the now trend due to biogeochemical processes
117	!! and add it to the general trend of passive tracers equations
118	!!
119	!! ** Method : each now biological flux is calculated in function of now
120	!! concentrations of tracers.
121	!! depending on the tracer, these fluxes are sources or sinks.
122	!! the total of the sources and sinks for each tracer
123	!! is added to the general trend.
124	!!
125	!! tra = tra + zf...tra - zftra...
126	!! \| \|
127	!! \| \|
128	!! source sink
129	!!---------------------------------------------------------------------
130	!!
131	!!
132	!!----------------------------------------------------------------------
133	!! Variable conventions
134	!!----------------------------------------------------------------------
135	!!
136	!! names: z*** - state variable
137	!! f*** - function (or temporary variable used in part of a function)
138	!! x*** - parameter
139	!! b*** - right-hand part (sources and sinks)
140	!! i*** - integer variable (usually used in yes/no flags)
141	!!
142	!! time (integer timestep)
143	INTEGER, INTENT( in ) :: kt
144	!!
145	!! spatial array indices
146	INTEGER :: ji,jj,jk,jn
147	!!
148	!! AXY (27/07/10): add in indices for depth horizons (for sinking flux
149	!! and seafloor iron inputs)
150	!! INTEGER :: i0100, i0200, i0500, i1000, i1100
151	!!
152	!! model state variables
153	REAL(wp) :: zchn,zchd,zphn,zphd,zpds,zzmi
154	REAL(wp) :: zzme,zdet,zdtc,zdin,zsil,zfer
155	# if defined key_roam
156	REAL(wp) :: zdic, zalk, zoxy
157	REAL(wp) :: ztmp, zsal
158	REAL(wp) :: zpho
159	# endif
160	!!
161	!! integrated source and sink terms
162	REAL(wp) :: b0
163	!! AXY (23/08/13): changed from individual variables for each flux to
164	!! an array that holds all fluxes
165	REAL(wp), DIMENSION(jp_medusa) :: btra
166	!!
167	!! primary production and chl related quantities
168	REAL(wp) :: fthetan,faln,fchn1,fchn,fjln,fprn,frn
169	REAL(wp) :: fthetad,fald,fchd1,fchd,fjld,fprd,frd
170	!! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range)
171	REAL(wp) :: fjlim_pn, fjlim_pd
172	!! AXY (03/02/11): add in Liebig terms
173	REAL(wp) :: fpnlim, fpdlim
174	!! AXY (16/07/09): add in Eppley curve functionality
175	REAL(wp) :: loc_T,fun_T,xvpnT,xvpdT
176	INTEGER :: ieppley
177	!! AXY (16/05/11): per Katya's prompting, add in new T-dependence
178	!! for phytoplankton growth only (i.e. no change
179	!! for remineralisation)
180	REAL(wp) :: fun_Q10
181	!! AXY (01/03/10): add in mixed layer PP diagnostics
182	REAL(wp), DIMENSION(jpi,jpj) :: fprn_ml,fprd_ml
183	!!
184	!! nutrient limiting factors
185	REAL(wp) :: fnln,ffln !! N and Fe
186	REAL(wp) :: fnld,ffld,fsld,fsld2 !! N, Fe and Si
187	!!
188	!! silicon cycle
189	REAL(wp) :: fsin,fnsi,fsin1,fnsi1,fnsi2,fprds,fsdiss
190	!!
191	!! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
192	REAL(wp) :: ffetop,ffebot,ffescav
193	REAL(wp) :: xLgF, xFeT, xFeF, xFeL !! state variables for iron-ligand system
194	REAL(wp), DIMENSION(jpi,jpj) :: xFree !! state variables for iron-ligand system
195	REAL(wp) :: xb_coef_tmp, xb2M4ac !! iron-ligand parameters
196	REAL(wp) :: xmaxFeF,fdeltaFe !! max Fe' parameters
197	!!
198	!! local parameters for Moore et al. (2004) alternative scavenging scheme
199	REAL(wp) :: fbase_scav,fscal_sink,fscal_part,fscal_scav
200	!!
201	!! local parameters for Moore et al. (2008) alternative scavenging scheme
202	REAL(wp) :: fscal_csink,fscal_sisink,fscal_casink
203	!!
204	!! local parameters for Galbraith et al. (2010) alternative scavenging scheme
205	REAL(wp) :: xCscav1, xCscav2, xk_org, xORGscav !! organic portion of scavenging
206	REAL(wp) :: xk_inorg, xINORGscav !! inorganic portion of scavenging
207	!!
208	!! microzooplankton grazing
209	REAL(wp) :: fmi1,fmi,fgmipn,fgmid,fgmidc
210	REAL(wp) :: finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
211	!!
212	!! mesozooplankton grazing
213	REAL(wp) :: fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
214	REAL(wp) :: finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
215	!!
216	!! mortality/Remineralisation (defunct parameter "fz" removed)
217	REAL(wp) :: fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
218	# if defined key_roam
219	REAL(wp) :: fddc
220	# endif
221	REAL(wp) :: fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
222	REAL(wp) :: fslown, fslowc
223	REAL(wp), DIMENSION(jpi,jpj) :: fslownflux, fslowcflux
224	REAL(wp) :: fregen,fregensi
225	REAL(wp), DIMENSION(jpi,jpj) :: fregenfast,fregenfastsi
226	# if defined key_roam
227	REAL(wp) :: fregenc
228	REAL(wp), DIMENSION(jpi,jpj) :: fregenfastc
229	# endif
230	!!
231	!! particle flux
232	REAL(WP) :: fthk,fdep,fdep1,fdep2,flat,fcaco3
233	REAL(WP) :: ftempn,ftempsi,ftempfe,ftempc,ftempca
234	REAL(wp) :: freminn,freminsi,freminfe,freminc,freminca
235	REAL(wp), DIMENSION(jpi,jpj) :: ffastn,ffastsi,ffastfe,ffastc,ffastca
236	REAL(wp) :: fleftn,fleftsi,fleftfe,fleftc,fleftca
237	REAL(wp) :: fheren,fheresi,fherefe,fherec,fhereca
238	REAL(wp) :: fprotf
239	REAL(wp), DIMENSION(jpi,jpj) :: fsedn,fsedsi,fsedfe,fsedc,fsedca
240	REAL(wp), DIMENSION(jpi,jpj) :: fccd
241	REAL(wp) :: fccd_dep
242	!! AXY (28/11/16): fix mbathy bug
243	INTEGER :: jmbathy
244	!!
245	!! AXY (06/07/11): alternative fast detritus schemes
246	REAL(wp) :: fb_val, fl_sst
247	!!
248	!! AXY (08/07/11): fate of fast detritus reaching the seafloor
249	REAL(wp) :: ffast2slown,ffast2slowfe,ffast2slowc
250	!!
251	!! conservation law
252	REAL(wp) :: fnit0,fsil0,ffer0
253	# if defined key_roam
254	REAL(wp) :: fcar0,falk0,foxy0
255	# endif
256	!!
257	!! temporary variables
258	REAL(wp) :: fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8,fq9
259	!!
260	!! water column nutrient and flux integrals
261	REAL(wp), DIMENSION(jpi,jpj) :: ftot_n,ftot_si,ftot_fe
262	REAL(wp), DIMENSION(jpi,jpj) :: fflx_n,fflx_si,fflx_fe
263	REAL(wp), DIMENSION(jpi,jpj) :: fifd_n,fifd_si,fifd_fe
264	REAL(wp), DIMENSION(jpi,jpj) :: fofd_n,fofd_si,fofd_fe
265	# if defined key_roam
266	REAL(wp), DIMENSION(jpi,jpj) :: ftot_c,ftot_a,ftot_o2
267	REAL(wp), DIMENSION(jpi,jpj) :: fflx_c,fflx_a,fflx_o2
268	REAL(wp), DIMENSION(jpi,jpj) :: fifd_c,fifd_a,fifd_o2
269	REAL(wp), DIMENSION(jpi,jpj) :: fofd_c,fofd_a,fofd_o2
270	# endif
271	!!
272	!! zooplankton grazing integrals
273	REAL(wp), DIMENSION(jpi,jpj) :: fzmi_i,fzmi_o,fzme_i,fzme_o
274	!!
275	!! limitation term temporary variables
276	REAL(wp), DIMENSION(jpi,jpj) :: ftot_pn,ftot_pd
277	REAL(wp), DIMENSION(jpi,jpj) :: ftot_zmi,ftot_zme,ftot_det,ftot_dtc
278	!! use ballast scheme (1) or simple exponential scheme (0; a conservation test)
279	INTEGER :: iball
280	!! use biological fluxes (1) or not (0)
281	INTEGER :: ibio_switch
282	!!
283	!! diagnose fluxes (should only be used in 1D runs)
284	INTEGER :: idf, idfval
285	!!
286	!! nitrogen and silicon production and consumption
287	REAL(wp) :: fn_prod, fn_cons, fs_prod, fs_cons
288	REAL(wp), DIMENSION(jpi,jpj) :: fnit_prod, fnit_cons, fsil_prod, fsil_cons
289	# if defined key_roam
290	!!
291	!! flags to help with calculating the position of the CCD
292	INTEGER, DIMENSION(jpi,jpj) :: i2_omcal,i2_omarg
293	!!
294	!! ROAM air-sea flux and diagnostic parameters
295	REAL(wp) :: f_wind
296	!! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
297	REAL(wp) :: f_xco2a
298	REAL(wp) :: f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
299	REAL(wp) :: f_TDIC, f_TALK, f_dcf, f_henry
300	REAL(wp) :: f_uwind, f_vwind, f_pp0
301	REAL(wp) :: f_kw660, f_o2flux, f_o2sat, f_o2sat3
302	REAL(wp), DIMENSION(jpi,jpj) :: f_omcal, f_omarg
303	!!
304	!! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
305	REAL(wp) :: f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm
306	REAL(wp) :: f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2
307	!! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
308	REAL, PARAMETER :: weight_CO2_mol = 44.0095 !! g / mol
309	REAL, PARAMETER :: secs_in_day = 86400.0 !! s / d
310	REAL, PARAMETER :: CO2flux_conv = (1.e-6 * weight_CO2_mol) / secs_in_day
311	!!
312	INTEGER :: iters
313	REAL(wp) :: f_year
314	INTEGER :: i_year
315	INTEGER :: iyr1, iyr2
316	!!
317	!! carbon, alkalinity production and consumption
318	REAL(wp) :: fc_prod, fc_cons, fa_prod, fa_cons
319	REAL(wp), DIMENSION(jpi,jpj) :: fcomm_resp
320	REAL(wp), DIMENSION(jpi,jpj) :: fcar_prod, fcar_cons
321	!!
322	!! oxygen production and consumption (and non-consumption)
323	REAL(wp) :: fo2_prod, fo2_cons, fo2_ncons, fo2_ccons
324	REAL(wp), DIMENSION(jpi,jpj) :: foxy_prod, foxy_cons, foxy_anox
325	!! Jpalm (11-08-2014)
326	!! add DMS in MEDUSA for UKESM1 model
327	REAL(wp) :: dms_surf
328	!! AXY (13/03/15): add in other DMS calculations
329	REAL(wp) :: dms_andr, dms_simo, dms_aran, dms_hall, dms_nlim, dms_wtkn
330	# endif
331	!!
332	!! benthic fluxes
333	INTEGER :: ibenthic
334	REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe, f_sbenin_c
335	REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca
336	REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca
337	REAL(wp) :: zfact
338	!!
339	!! benthic fluxes of CaCO3 that shouldn't happen because of lysocline
340	REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca
341	!!
342	!! riverine fluxes
343	REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk
344	!! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface
345	REAL(wp) :: f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk
346	!!
347	!! horizontal grid location
348	REAL(wp) :: flatx, flonx
349	!!
350	!! Jpalm -- 11-10-2015 -- adapt diag to iom_use
351	!! 2D var for diagnostics.
352	REAL(wp), POINTER, DIMENSION(:,: ) :: fprn2d, fdpn2d, fprd2d, fdpd2d, fprds2d, fsdiss2d, fgmipn2d
353	REAL(wp), POINTER, DIMENSION(:,: ) :: fgmid2d, fdzmi2d, fgmepn2d, fgmepd2d, fgmezmi2d, fgmed2d
354	REAL(wp), POINTER, DIMENSION(:,: ) :: fdzme2d, fslown2d, fdd2d, ffetop2d, ffebot2d, ffescav2d
355	REAL(wp), POINTER, DIMENSION(:,: ) :: fjln2d, fnln2d, ffln2d, fjld2d, fnld2d, ffld2d, fsld2d2
356	REAL(wp), POINTER, DIMENSION(:,: ) :: fsld2d, fregen2d, fregensi2d, ftempn2d, ftempsi2d, ftempfe2d
357	REAL(wp), POINTER, DIMENSION(:,: ) :: ftempc2d, ftempca2d, freminn2d, freminsi2d, freminfe2d
358	REAL(wp), POINTER, DIMENSION(:,: ) :: freminc2d, freminca2d
359	REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d
360	# if defined key_roam
361	REAL(wp), POINTER, DIMENSION(:,: ) :: ffastca2d, rivn2d, rivsi2d, rivc2d, rivalk2d, fslowc2d
362	REAL(wp), POINTER, DIMENSION(:,: ) :: fdpn22d, fdpd22d, fdzmi22d, fdzme22d, zimesn2d, zimesd2d
363	REAL(wp), POINTER, DIMENSION(:,: ) :: zimesc2d, zimesdc2d, ziexcr2d, ziresp2d, zigrow2d, zemesn2d
364	REAL(wp), POINTER, DIMENSION(:,: ) :: zemesd2d, zemesc2d, zemesdc2d, zeexcr2d, zeresp2d, zegrow2d
365	REAL(wp), POINTER, DIMENSION(:,: ) :: mdetc2d, gmidc2d, gmedc2d, f_pco2a2d, f_pco2w2d, f_co2flux2d
366	REAL(wp), POINTER, DIMENSION(:,: ) :: f_TDIC2d, f_TALK2d, f_kw6602d, f_pp02d, f_o2flux2d, f_o2sat2d
367	REAL(wp), POINTER, DIMENSION(:,: ) :: dms_andr2d, dms_simo2d, dms_aran2d, dms_hall2d, dms_surf2d
368	REAL(wp), POINTER, DIMENSION(:,: ) :: iben_n2d, iben_fe2d, iben_c2d, iben_si2d, iben_ca2d, oben_n2d
369	REAL(wp), POINTER, DIMENSION(:,: ) :: oben_fe2d, oben_c2d, oben_si2d, oben_ca2d, sfr_ocal2d
370	REAL(wp), POINTER, DIMENSION(:,: ) :: sfr_oarg2d, lyso_ca2d
371	!! AXY (23/11/16): extra MOCSY diagnostics
372	REAL(wp), POINTER, DIMENSION(:,: ) :: f_xco2a_2d, f_fco2w_2d, f_fco2a_2d
373	REAL(wp), POINTER, DIMENSION(:,: ) :: f_ocnrhosw_2d, f_ocnschco2_2d, f_ocnkwco2_2d
374	REAL(wp), POINTER, DIMENSION(:,: ) :: f_ocnk0_2d, f_co2starair_2d, f_ocndpco2_2d
375	# endif
376	!!
377	!! 3D var for diagnostics.
378	REAL(wp), POINTER, DIMENSION(:,:,:) :: tpp3d, detflux3d, remin3dn
379	!!
380	# if defined key_roam
381	!! AXY (04/11/16)
382	!! 2D var for new CMIP6 diagnostics (behind a key_roam ifdef for simplicity)
383	REAL(wp), POINTER, DIMENSION(:,: ) :: fgco2, intdissic, intdissin, intdissisi, inttalk, o2min, zo2min
384	REAL(wp), POINTER, DIMENSION(:,: ) :: fbddtalk, fbddtdic, fbddtdife, fbddtdin, fbddtdisi
385	!!
386	!! 3D var for new CMIP6 diagnostics
387	REAL(wp), POINTER, DIMENSION(:,:,:) :: tppd3
388	REAL(wp), POINTER, DIMENSION(:,:,:) :: bddtalk3, bddtdic3, bddtdife3, bddtdin3, bddtdisi3
389	REAL(wp), POINTER, DIMENSION(:,:,:) :: fd_nit3, fd_sil3, fd_car3, fd_cal3
390	REAL(wp), POINTER, DIMENSION(:,:,:) :: co33, co3satarag3, co3satcalc3, dcalc3
391	REAL(wp), POINTER, DIMENSION(:,:,:) :: expc3, expn3
392	REAL(wp), POINTER, DIMENSION(:,:,:) :: fediss3, fescav3
393	REAL(wp), POINTER, DIMENSION(:,:,:) :: migrazp3, migrazd3, megrazp3, megrazd3, megrazz3
394	REAL(wp), POINTER, DIMENSION(:,:,:) :: o2sat3, pbsi3, pcal3, remoc3
395	REAL(wp), POINTER, DIMENSION(:,:,:) :: pnlimj3, pnlimn3, pnlimfe3, pdlimj3, pdlimn3, pdlimfe3, pdlimsi3
396	# endif
397	!!---------------------------------------------------------------------
398
399	# if defined key_debug_medusa
400	IF ( lwp ) write (numout,*) 'trc_bio_medusa: variables defined'
401	CALL flush(numout)
402	# endif
403
404	!! AXY (20/11/14): alter this to report on first MEDUSA call
405	!! IF( kt == nit000 ) THEN
406	IF( kt == nittrc000 ) THEN
407	IF(lwp) WRITE(numout,*)
408	IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model'
409	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
410	IF(lwp) WRITE(numout,*) ' kt =',kt
411	ENDIF
412
413	!! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
414	ibenthic = 1
415
416	!!----------------------------------------------------------------------
417	!! b0 is present for debugging purposes; using b0 = 0 sets the tendency
418	!! terms of all biological equations to 0.
419	!!----------------------------------------------------------------------
420	!!
421	!! AXY (03/09/14): probably not the smartest move ever, but it'll fit
422	!! the bill for now; another item on the things-to-sort-
423	!! out-in-the-future list ...
424	# if defined key_kill_medusa
425	b0 = 0.
426	# else
427	b0 = 1.
428	# endif
429	!!----------------------------------------------------------------------
430	!! fast detritus ballast scheme (0 = no; 1 = yes)
431	!! alternative to ballast scheme is same scheme but with no ballast
432	!! protection (not dissimilar to Martin et al., 1987)
433	!!----------------------------------------------------------------------
434	!!
435	iball = 1
436
437	!!----------------------------------------------------------------------
438	!! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
439	!! these should only be used in 1D since they give comprehensive
440	!! output for ecological functions in the model; primarily used in
441	!! debugging
442	!!----------------------------------------------------------------------
443	!!
444	idf = 0
445	!!
446	!! timer mechanism
447	if (kt/120*120.eq.kt) then
448	idfval = 1
449	else
450	idfval = 0
451	endif
452
453	!!----------------------------------------------------------------------
454	!! blank fast-sinking detritus 2D fields
455	!!----------------------------------------------------------------------
456	!!
457	ffastn(:,:) = 0.0 !! organic nitrogen
458	ffastsi(:,:) = 0.0 !! biogenic silicon
459	ffastfe(:,:) = 0.0 !! organic iron
460	ffastc(:,:) = 0.0 !! organic carbon
461	ffastca(:,:) = 0.0 !! biogenic calcium carbonate
462	!!
463	fsedn(:,:) = 0.0 !! Seafloor flux of N
464	fsedsi(:,:) = 0.0 !! Seafloor flux of Si
465	fsedfe(:,:) = 0.0 !! Seafloor flux of Fe
466	fsedc(:,:) = 0.0 !! Seafloor flux of C
467	fsedca(:,:) = 0.0 !! Seafloor flux of CaCO3
468	!!
469	fregenfast(:,:) = 0.0 !! integrated N regeneration (fast detritus)
470	fregenfastsi(:,:) = 0.0 !! integrated Si regeneration (fast detritus)
471	# if defined key_roam
472	fregenfastc(:,:) = 0.0 !! integrated C regeneration (fast detritus)
473	# endif
474	!!
475	fccd(:,:) = 0.0 !! last depth level before CCD
476
477	!!----------------------------------------------------------------------
478	!! blank nutrient/flux inventories
479	!!----------------------------------------------------------------------
480	!!
481	fflx_n(:,:) = 0.0 !! nitrogen flux total
482	fflx_si(:,:) = 0.0 !! silicon flux total
483	fflx_fe(:,:) = 0.0 !! iron flux total
484	fifd_n(:,:) = 0.0 !! nitrogen fast detritus production
485	fifd_si(:,:) = 0.0 !! silicon fast detritus production
486	fifd_fe(:,:) = 0.0 !! iron fast detritus production
487	fofd_n(:,:) = 0.0 !! nitrogen fast detritus remineralisation
488	fofd_si(:,:) = 0.0 !! silicon fast detritus remineralisation
489	fofd_fe(:,:) = 0.0 !! iron fast detritus remineralisation
490	# if defined key_roam
491	fflx_c(:,:) = 0.0 !! carbon flux total
492	fflx_a(:,:) = 0.0 !! alkalinity flux total
493	fflx_o2(:,:) = 0.0 !! oxygen flux total
494	ftot_c(:,:) = 0.0 !! carbon inventory
495	ftot_a(:,:) = 0.0 !! alkalinity inventory
496	ftot_o2(:,:) = 0.0 !! oxygen inventory
497	fifd_c(:,:) = 0.0 !! carbon fast detritus production
498	fifd_a(:,:) = 0.0 !! alkalinity fast detritus production
499	fifd_o2(:,:) = 0.0 !! oxygen fast detritus production
500	fofd_c(:,:) = 0.0 !! carbon fast detritus remineralisation
501	fofd_a(:,:) = 0.0 !! alkalinity fast detritus remineralisation
502	fofd_o2(:,:) = 0.0 !! oxygen fast detritus remineralisation
503	!!
504	fnit_prod(:,:) = 0.0 !! (organic) nitrogen production
505	fnit_cons(:,:) = 0.0 !! (organic) nitrogen consumption
506	fsil_prod(:,:) = 0.0 !! (inorganic) silicon production
507	fsil_cons(:,:) = 0.0 !! (inorganic) silicon consumption
508	fcar_prod(:,:) = 0.0 !! (organic) carbon production
509	fcar_cons(:,:) = 0.0 !! (organic) carbon consumption
510	!!
511	foxy_prod(:,:) = 0.0 !! oxygen production
512	foxy_cons(:,:) = 0.0 !! oxygen consumption
513	foxy_anox(:,:) = 0.0 !! unrealised oxygen consumption
514	!!
515	# endif
516	ftot_n(:,:) = 0.0 !! N inventory
517	ftot_si(:,:) = 0.0 !! Si inventory
518	ftot_fe(:,:) = 0.0 !! Fe inventory
519	ftot_pn(:,:) = 0.0 !! integrated non-diatom phytoplankton
520	ftot_pd(:,:) = 0.0 !! integrated diatom phytoplankton
521	ftot_zmi(:,:) = 0.0 !! integrated microzooplankton
522	ftot_zme(:,:) = 0.0 !! integrated mesozooplankton
523	ftot_det(:,:) = 0.0 !! integrated slow detritus, nitrogen
524	ftot_dtc(:,:) = 0.0 !! integrated slow detritus, carbon
525	!!
526	fzmi_i(:,:) = 0.0 !! material grazed by microzooplankton
527	fzmi_o(:,:) = 0.0 !! ... sum of fate of this material
528	fzme_i(:,:) = 0.0 !! material grazed by mesozooplankton
529	fzme_o(:,:) = 0.0 !! ... sum of fate of this material
530	!!
531	f_sbenin_n(:,:) = 0.0 !! slow detritus N -> benthic pool
532	f_sbenin_fe(:,:) = 0.0 !! slow detritus Fe -> benthic pool
533	f_sbenin_c(:,:) = 0.0 !! slow detritus C -> benthic pool
534	f_fbenin_n(:,:) = 0.0 !! fast detritus N -> benthic pool
535	f_fbenin_fe(:,:) = 0.0 !! fast detritus Fe -> benthic pool
536	f_fbenin_si(:,:) = 0.0 !! fast detritus Si -> benthic pool
537	f_fbenin_c(:,:) = 0.0 !! fast detritus C -> benthic pool
538	f_fbenin_ca(:,:) = 0.0 !! fast detritus Ca -> benthic pool
539	!!
540	f_benout_n(:,:) = 0.0 !! benthic N pool -> dissolved
541	f_benout_fe(:,:) = 0.0 !! benthic Fe pool -> dissolved
542	f_benout_si(:,:) = 0.0 !! benthic Si pool -> dissolved
543	f_benout_c(:,:) = 0.0 !! benthic C pool -> dissolved
544	f_benout_ca(:,:) = 0.0 !! benthic Ca pool -> dissolved
545	!!
546	f_benout_lyso_ca(:,:) = 0.0 !! benthic Ca pool -> dissolved (when it shouldn't!)
547	!!
548	f_runoff(:,:) = 0.0 !! riverine runoff
549	f_riv_n(:,:) = 0.0 !! riverine N input
550	f_riv_si(:,:) = 0.0 !! riverine Si input
551	f_riv_c(:,:) = 0.0 !! riverine C input
552	f_riv_alk(:,:) = 0.0 !! riverine alk input
553	!!
554	!! Jpalm -- 06-03-2017 -- Forgotten var to init
555	f_omarg(:,:) = 0.0 !!
556	f_omcal(:,:) = 0.0
557	xFree(:,:) = 0.0 !! state variables for iron-ligand system
558	fcomm_resp(:,:) = 0.0
559	fprn_ml(:,:) = 0.0 !! mixed layer PP diagnostics
560	fprd_ml(:,:) = 0.0 !! mixed layer PP diagnostics
561
562	!!----------------------------------------------------------------------
563	!! allocate and initiate 2D diag
564	!!----------------------------------------------------------------------
565	!!
566	IF ( lk_iomput .AND. .NOT. ln_diatrc ) THEN
567	!! Juju :: add kt condition !!
568	if ( kt == nittrc000 ) CALL trc_nam_iom_medusa !! initialise iom_use test
569	!!
570	CALL wrk_alloc( jpi, jpj, zw2d )
571	zw2d(:,:) = 0.0 !!
572	IF ( med_diag%PRN%dgsave ) THEN
573	CALL wrk_alloc( jpi, jpj, fprn2d )
574	fprn2d(:,:) = 0.0 !!
575	ENDIF
576	IF ( med_diag%MPN%dgsave ) THEN
577	CALL wrk_alloc( jpi, jpj, fdpn2d )
578	fdpn2d(:,:) = 0.0 !!
579	ENDIF
580	IF ( med_diag%PRD%dgsave ) THEN
581	CALL wrk_alloc( jpi, jpj, fprd2d )
582	fprd2d(:,:) = 0.0 !!
583	ENDIF
584	IF( med_diag%MPD%dgsave ) THEN
585	CALL wrk_alloc( jpi, jpj, fdpd2d )
586	fdpd2d(:,:) = 0.0 !!
587	ENDIF
588	IF( med_diag%OPAL%dgsave ) THEN
589	CALL wrk_alloc( jpi, jpj, fprds2d )
590	fprds2d(:,:) = 0.0 !!
591	ENDIF
592	IF( med_diag%OPALDISS%dgsave ) THEN
593	CALL wrk_alloc( jpi, jpj, fsdiss2d )
594	fsdiss2d(:,:) = 0.0 !!
595	ENDIF
596	IF( med_diag%GMIPn%dgsave ) THEN
597	CALL wrk_alloc( jpi, jpj, fgmipn2d )
598	fgmipn2d(:,:) = 0.0 !!
599	ENDIF
600	IF( med_diag%GMID%dgsave ) THEN
601	CALL wrk_alloc( jpi, jpj, fgmid2d )
602	fgmid2d(:,:) = 0.0 !!
603	ENDIF
604	IF( med_diag%MZMI%dgsave ) THEN
605	CALL wrk_alloc( jpi, jpj, fdzmi2d )
606	fdzmi2d(:,:) = 0.0 !!
607	ENDIF
608	IF( med_diag%GMEPN%dgsave ) THEN
609	CALL wrk_alloc( jpi, jpj, fgmepn2d )
610	fgmepn2d(:,:) = 0.0 !!
611	ENDIF
612	IF( med_diag%GMEPD%dgsave ) THEN
613	CALL wrk_alloc( jpi, jpj, fgmepd2d )
614	fgmepd2d(:,:) = 0.0 !!
615	ENDIF
616	IF( med_diag%GMEZMI%dgsave ) THEN
617	CALL wrk_alloc( jpi, jpj, fgmezmi2d )
618	fgmezmi2d(:,:) = 0.0 !!
619	ENDIF
620	IF( med_diag%GMED%dgsave ) THEN
621	CALL wrk_alloc( jpi, jpj, fgmed2d )
622	fgmed2d(:,:) = 0.0 !!
623	ENDIF
624	IF( med_diag%MZME%dgsave ) THEN
625	CALL wrk_alloc( jpi, jpj, fdzme2d )
626	fdzme2d(:,:) = 0.0 !!
627	ENDIF
628	IF( med_diag%DETN%dgsave ) THEN
629	CALL wrk_alloc( jpi, jpj, fslown2d )
630	fslown2d(:,:) = 0.0 !!
631	ENDIF
632	IF( med_diag%MDET%dgsave ) THEN
633	CALL wrk_alloc( jpi, jpj, fdd2d )
634	fdd2d(:,:) = 0.0 !!
635	ENDIF
636	IF( med_diag%AEOLIAN%dgsave ) THEN
637	CALL wrk_alloc( jpi, jpj, ffetop2d )
638	ffetop2d(:,:) = 0.0 !!
639	ENDIF
640	IF( med_diag%BENTHIC%dgsave ) THEN
641	CALL wrk_alloc( jpi, jpj, ffebot2d )
642	ffebot2d(:,:) = 0.0 !!
643	ENDIF
644	IF( med_diag%SCAVENGE%dgsave ) THEN
645	CALL wrk_alloc( jpi, jpj, ffescav2d )
646	ffescav2d(:,:) = 0.0 !!
647	ENDIF
648	IF( med_diag%PN_JLIM%dgsave ) THEN
649	CALL wrk_alloc( jpi, jpj, fjln2d )
650	fjln2d(:,:) = 0.0 !!
651	ENDIF
652	IF( med_diag%PN_NLIM%dgsave ) THEN
653	CALL wrk_alloc( jpi, jpj, fnln2d )
654	fnln2d(:,:) = 0.0 !!
655	ENDIF
656	IF( med_diag%PN_FELIM%dgsave ) THEN
657	CALL wrk_alloc( jpi, jpj, ffln2d )
658	ffln2d(:,:) = 0.0 !!
659	ENDIF
660	IF( med_diag%PD_JLIM%dgsave ) THEN
661	CALL wrk_alloc( jpi, jpj, fjld2d )
662	fjld2d(:,:) = 0.0 !!
663	ENDIF
664	IF( med_diag%PD_NLIM%dgsave ) THEN
665	CALL wrk_alloc( jpi, jpj, fnld2d )
666	fnld2d(:,:) = 0.0 !!
667	ENDIF
668	IF( med_diag%PD_FELIM%dgsave ) THEN
669	CALL wrk_alloc( jpi, jpj, ffld2d )
670	ffld2d(:,:) = 0.0 !!
671	ENDIF
672	IF( med_diag%PD_SILIM%dgsave ) THEN
673	CALL wrk_alloc( jpi, jpj, fsld2d2 )
674	fsld2d2(:,:) = 0.0 !!
675	ENDIF
676	IF( med_diag%PDSILIM2%dgsave ) THEN
677	CALL wrk_alloc( jpi, jpj, fsld2d )
678	fsld2d(:,:) = 0.0 !!
679	ENDIF
680	!!
681	!! skip SDT_XXXX diagnostics here
682	!!
683	IF( med_diag%TOTREG_N%dgsave ) THEN
684	CALL wrk_alloc( jpi, jpj, fregen2d )
685	fregen2d(:,:) = 0.0 !!
686	ENDIF
687	IF( med_diag%TOTRG_SI%dgsave ) THEN
688	CALL wrk_alloc( jpi, jpj, fregensi2d )
689	fregensi2d(:,:) = 0.0 !!
690	ENDIF
691	!!
692	!! skip REG_XXXX diagnostics here
693	!!
694	IF( med_diag%FASTN%dgsave ) THEN
695	CALL wrk_alloc( jpi, jpj, ftempn2d )
696	ftempn2d(:,:) = 0.0 !!
697	ENDIF
698	IF( med_diag%FASTSI%dgsave ) THEN
699	CALL wrk_alloc( jpi, jpj, ftempsi2d )
700	ftempsi2d(:,:) = 0.0 !!
701	ENDIF
702	IF( med_diag%FASTFE%dgsave ) THEN
703	CALL wrk_alloc( jpi, jpj, ftempfe2d )
704	ftempfe2d(:,:) = 0.0 !!
705	ENDIF
706	IF( med_diag%FASTC%dgsave ) THEN
707	CALL wrk_alloc( jpi, jpj, ftempc2d )
708	ftempc2d(:,:) = 0.0 !!
709	ENDIF
710	IF( med_diag%FASTCA%dgsave ) THEN
711	CALL wrk_alloc( jpi, jpj, ftempca2d )
712	ftempca2d(:,:) = 0.0 !!
713	ENDIF
714	!!
715	!! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here
716	!!
717	IF( med_diag%REMINN%dgsave ) THEN
718	CALL wrk_alloc( jpi, jpj, freminn2d )
719	freminn2d(:,:) = 0.0 !!
720	ENDIF
721	IF( med_diag%REMINSI%dgsave ) THEN
722	CALL wrk_alloc( jpi, jpj, freminsi2d )
723	freminsi2d(:,:) = 0.0 !!
724	ENDIF
725	IF( med_diag%REMINFE%dgsave ) THEN
726	CALL wrk_alloc( jpi, jpj, freminfe2d )
727	freminfe2d(:,:) = 0.0 !!
728	ENDIF
729	IF( med_diag%REMINC%dgsave ) THEN
730	CALL wrk_alloc( jpi, jpj, freminc2d )
731	freminc2d(:,:) = 0.0 !!
732	ENDIF
733	IF( med_diag%REMINCA%dgsave ) THEN
734	CALL wrk_alloc( jpi, jpj, freminca2d )
735	freminca2d(:,:) = 0.0 !!
736	ENDIF
737	# if defined key_roam
738	!!
739	!! skip SEAFLRXX, MED_XXXX, INTFLX_XX, INT_XX, ML_XXX, OCAL_XXX, FE_XXXX, MED_XZE, WIND diagnostics here
740	!!
741	IF( med_diag%RR_0100%dgsave ) THEN
742	CALL wrk_alloc( jpi, jpj, ffastca2d )
743	ffastca2d(:,:) = 0.0 !!
744	ENDIF
745
746	IF( med_diag%ATM_PCO2%dgsave ) THEN
747	CALL wrk_alloc( jpi, jpj, f_pco2a2d )
748	f_pco2a2d(:,:) = 0.0 !!
749	ENDIF
750	!!
751	!! skip OCN_PH diagnostic here
752	!!
753	IF( med_diag%OCN_PCO2%dgsave ) THEN
754	CALL wrk_alloc( jpi, jpj, f_pco2w2d )
755	f_pco2w2d(:,:) = 0.0 !!
756	ENDIF
757	!!
758	!! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here
759	!!
760	IF( med_diag%CO2FLUX%dgsave ) THEN
761	CALL wrk_alloc( jpi, jpj, f_co2flux2d )
762	f_co2flux2d(:,:) = 0.0 !!
763	ENDIF
764	!!
765	!! skip OM_XXX diagnostics here
766	!!
767	IF( med_diag%TCO2%dgsave ) THEN
768	CALL wrk_alloc( jpi, jpj, f_TDIC2d )
769	f_TDIC2d(:,:) = 0.0 !!
770	ENDIF
771	IF( med_diag%TALK%dgsave ) THEN
772	CALL wrk_alloc( jpi, jpj, f_TALK2d )
773	f_TALK2d(:,:) = 0.0 !!
774	ENDIF
775	IF( med_diag%KW660%dgsave ) THEN
776	CALL wrk_alloc( jpi, jpj, f_kw6602d )
777	f_kw6602d(:,:) = 0.0 !!
778	ENDIF
779	IF( med_diag%ATM_PP0%dgsave ) THEN
780	CALL wrk_alloc( jpi, jpj, f_pp02d )
781	f_pp02d(:,:) = 0.0 !!
782	ENDIF
783	IF( med_diag%O2FLUX%dgsave ) THEN
784	CALL wrk_alloc( jpi, jpj, f_o2flux2d )
785	f_o2flux2d(:,:) = 0.0 !!
786	ENDIF
787	IF( med_diag%O2SAT%dgsave ) THEN
788	CALL wrk_alloc( jpi, jpj, f_o2sat2d )
789	f_o2sat2d(:,:) = 0.0 !!
790	ENDIF
791	!!
792	!! skip XXX_CCD diagnostics here
793	!!
794	IF( med_diag%SFR_OCAL%dgsave ) THEN
795	CALL wrk_alloc( jpi, jpj, sfr_ocal2d )
796	sfr_ocal2d(:,:) = 0.0 !!
797	ENDIF
798	IF( med_diag%SFR_OARG%dgsave ) THEN
799	CALL wrk_alloc( jpi, jpj, sfr_oarg2d )
800	sfr_oarg2d(:,:) = 0.0 !!
801	ENDIF
802	!!
803	!! skip XX_PROD, XX_CONS, O2_ANOX, RR_XXXX diagnostics here
804	!!
805	IF( med_diag%IBEN_N%dgsave ) THEN
806	CALL wrk_alloc( jpi, jpj, iben_n2d )
807	iben_n2d(:,:) = 0.0 !!
808	ENDIF
809	IF( med_diag%IBEN_FE%dgsave ) THEN
810	CALL wrk_alloc( jpi, jpj, iben_fe2d )
811	iben_fe2d(:,:) = 0.0 !!
812	ENDIF
813	IF( med_diag%IBEN_C%dgsave ) THEN
814	CALL wrk_alloc( jpi, jpj, iben_c2d )
815	iben_c2d(:,:) = 0.0 !!
816	ENDIF
817	IF( med_diag%IBEN_SI%dgsave ) THEN
818	CALL wrk_alloc( jpi, jpj, iben_si2d )
819	iben_si2d(:,:) = 0.0 !!
820	ENDIF
821	IF( med_diag%IBEN_CA%dgsave ) THEN
822	CALL wrk_alloc( jpi, jpj, iben_ca2d )
823	iben_ca2d(:,:) = 0.0 !!
824	ENDIF
825	IF( med_diag%OBEN_N%dgsave ) THEN
826	CALL wrk_alloc( jpi, jpj, oben_n2d )
827	oben_n2d(:,:) = 0.0 !!
828	ENDIF
829	IF( med_diag%OBEN_FE%dgsave ) THEN
830	CALL wrk_alloc( jpi, jpj, oben_fe2d )
831	oben_fe2d(:,:) = 0.0 !!
832	ENDIF
833	IF( med_diag%OBEN_C%dgsave ) THEN
834	CALL wrk_alloc( jpi, jpj, oben_c2d )
835	oben_c2d(:,:) = 0.0 !!
836	ENDIF
837	IF( med_diag%OBEN_SI%dgsave ) THEN
838	CALL wrk_alloc( jpi, jpj, oben_si2d )
839	oben_si2d(:,:) = 0.0 !!
840	ENDIF
841	IF( med_diag%OBEN_CA%dgsave ) THEN
842	CALL wrk_alloc( jpi, jpj, oben_ca2d )
843	oben_ca2d(:,:) = 0.0 !!
844	ENDIF
845	!!
846	!! skip BEN_XX diagnostics here
847	!!
848	IF( med_diag%RIV_N%dgsave ) THEN
849	CALL wrk_alloc( jpi, jpj, rivn2d )
850	rivn2d(:,:) = 0.0 !!
851	ENDIF
852	IF( med_diag%RIV_SI%dgsave ) THEN
853	CALL wrk_alloc( jpi, jpj, rivsi2d )
854	rivsi2d(:,:) = 0.0 !!
855	ENDIF
856	IF( med_diag%RIV_C%dgsave ) THEN
857	CALL wrk_alloc( jpi, jpj, rivc2d )
858	rivc2d(:,:) = 0.0 !!
859	ENDIF
860	IF( med_diag%RIV_ALK%dgsave ) THEN
861	CALL wrk_alloc( jpi, jpj, rivalk2d )
862	rivalk2d(:,:) = 0.0 !!
863	ENDIF
864	IF( med_diag%DETC%dgsave ) THEN
865	CALL wrk_alloc( jpi, jpj, fslowc2d )
866	fslowc2d(:,:) = 0.0 !!
867	ENDIF
868	!!
869	!! skip SDC_XXXX, INVTXXX diagnostics here
870	!!
871	IF( med_diag%LYSO_CA%dgsave ) THEN
872	CALL wrk_alloc( jpi, jpj, lyso_ca2d )
873	lyso_ca2d(:,:) = 0.0 !!
874	ENDIF
875	!!
876	!! skip COM_RESP diagnostic here
877	!!
878	IF( med_diag%PN_LLOSS%dgsave ) THEN
879	CALL wrk_alloc( jpi, jpj, fdpn22d )
880	fdpn22d(:,:) = 0.0 !!
881	ENDIF
882	IF( med_diag%PD_LLOSS%dgsave ) THEN
883	CALL wrk_alloc( jpi, jpj, fdpd22d )
884	fdpd22d(:,:) = 0.0 !!
885	ENDIF
886	IF( med_diag%ZI_LLOSS%dgsave ) THEN
887	CALL wrk_alloc( jpi, jpj, fdzmi22d )
888	fdzmi22d(:,:) = 0.0 !!
889	ENDIF
890	IF( med_diag%ZE_LLOSS%dgsave ) THEN
891	CALL wrk_alloc( jpi, jpj, fdzme22d )
892	fdzme22d(:,:) = 0.0 !!
893	ENDIF
894	IF( med_diag%ZI_MES_N%dgsave ) THEN
895	CALL wrk_alloc( jpi, jpj, zimesn2d )
896	zimesn2d(:,:) = 0.0 !!
897	ENDIF
898	IF( med_diag%ZI_MES_D%dgsave ) THEN
899	CALL wrk_alloc( jpi, jpj, zimesd2d )
900	zimesd2d(:,:) = 0.0 !!
901	ENDIF
902	IF( med_diag%ZI_MES_C%dgsave ) THEN
903	CALL wrk_alloc( jpi, jpj, zimesc2d )
904	zimesc2d(:,:) = 0.0 !!
905	ENDIF
906	IF( med_diag%ZI_MESDC%dgsave ) THEN
907	CALL wrk_alloc( jpi, jpj, zimesdc2d )
908	zimesdc2d(:,:) = 0.0 !!
909	ENDIF
910	IF( med_diag%ZI_EXCR%dgsave ) THEN
911	CALL wrk_alloc( jpi, jpj, ziexcr2d )
912	ziexcr2d(:,:) = 0.0 !!
913	ENDIF
914	IF( med_diag%ZI_RESP%dgsave ) THEN
915	CALL wrk_alloc( jpi, jpj, ziresp2d )
916	ziresp2d(:,:) = 0.0 !!
917	ENDIF
918	IF( med_diag%ZI_GROW%dgsave ) THEN
919	CALL wrk_alloc( jpi, jpj, zigrow2d )
920	zigrow2d(:,:) = 0.0 !!
921	ENDIF
922	IF( med_diag%ZE_MES_N%dgsave ) THEN
923	CALL wrk_alloc( jpi, jpj, zemesn2d )
924	zemesn2d(:,:) = 0.0 !!
925	ENDIF
926	IF( med_diag%ZE_MES_D%dgsave ) THEN
927	CALL wrk_alloc( jpi, jpj, zemesd2d )
928	zemesd2d(:,:) = 0.0 !!
929	ENDIF
930	IF( med_diag%ZE_MES_C%dgsave ) THEN
931	CALL wrk_alloc( jpi, jpj, zemesc2d )
932	zemesc2d(:,:) = 0.0 !!
933	ENDIF
934	IF( med_diag%ZE_MESDC%dgsave ) THEN
935	CALL wrk_alloc( jpi, jpj, zemesdc2d )
936	zemesdc2d(:,:) = 0.0 !!
937	ENDIF
938	IF( med_diag%ZE_EXCR%dgsave ) THEN
939	CALL wrk_alloc( jpi, jpj, zeexcr2d )
940	zeexcr2d(:,:) = 0.0 !!
941	ENDIF
942	IF( med_diag%ZE_RESP%dgsave ) THEN
943	CALL wrk_alloc( jpi, jpj, zeresp2d )
944	zeresp2d(:,:) = 0.0 !!
945	ENDIF
946	IF( med_diag%ZE_GROW%dgsave ) THEN
947	CALL wrk_alloc( jpi, jpj, zegrow2d )
948	zegrow2d(:,:) = 0.0 !!
949	ENDIF
950	IF( med_diag%MDETC%dgsave ) THEN
951	CALL wrk_alloc( jpi, jpj, mdetc2d )
952	mdetc2d(:,:) = 0.0 !!
953	ENDIF
954	IF( med_diag%GMIDC%dgsave ) THEN
955	CALL wrk_alloc( jpi, jpj, gmidc2d )
956	gmidc2d(:,:) = 0.0 !!
957	ENDIF
958	IF( med_diag%GMEDC%dgsave ) THEN
959	CALL wrk_alloc( jpi, jpj, gmedc2d )
960	gmedc2d(:,:) = 0.0 !!
961	ENDIF
962	!!
963	!! skip INT_XXX diagnostics here
964	!!
965	IF (jdms .eq. 1) THEN
966	IF( med_diag%DMS_SURF%dgsave ) THEN
967	CALL wrk_alloc( jpi, jpj, dms_surf2d )
968	dms_surf2d(:,:) = 0.0 !!
969	ENDIF
970	IF( med_diag%DMS_ANDR%dgsave ) THEN
971	CALL wrk_alloc( jpi, jpj, dms_andr2d )
972	dms_andr2d(:,:) = 0.0 !!
973	ENDIF
974	IF( med_diag%DMS_SIMO%dgsave ) THEN
975	CALL wrk_alloc( jpi, jpj, dms_simo2d )
976	dms_simo2d(:,:) = 0.0 !!
977	ENDIF
978	IF( med_diag%DMS_ARAN%dgsave ) THEN
979	CALL wrk_alloc( jpi, jpj, dms_aran2d )
980	dms_aran2d(:,:) = 0.0 !!
981	ENDIF
982	IF( med_diag%DMS_HALL%dgsave ) THEN
983	CALL wrk_alloc( jpi, jpj, dms_hall2d )
984	dms_hall2d(:,:) = 0.0 !!
985	ENDIF
986	ENDIF
987	!!
988	!! AXY (24/11/16): extra MOCSY diagnostics, 2D
989	IF( med_diag%ATM_XCO2%dgsave ) THEN
990	CALL wrk_alloc( jpi, jpj, f_xco2a_2d )
991	f_xco2a_2d(:,:) = 0.0 !!
992	ENDIF
993	IF( med_diag%OCN_FCO2%dgsave ) THEN
994	CALL wrk_alloc( jpi, jpj, f_fco2w_2d )
995	f_fco2w_2d(:,:) = 0.0 !!
996	ENDIF
997	IF( med_diag%ATM_FCO2%dgsave ) THEN
998	CALL wrk_alloc( jpi, jpj, f_fco2a_2d )
999	f_fco2a_2d(:,:) = 0.0 !!
1000	ENDIF
1001	IF( med_diag%OCN_RHOSW%dgsave ) THEN
1002	CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d )
1003	f_ocnrhosw_2d(:,:) = 0.0 !!
1004	ENDIF
1005	IF( med_diag%OCN_SCHCO2%dgsave ) THEN
1006	CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d )
1007	f_ocnschco2_2d(:,:) = 0.0 !!
1008	ENDIF
1009	IF( med_diag%OCN_KWCO2%dgsave ) THEN
1010	CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d )
1011	f_ocnkwco2_2d(:,:) = 0.0 !!
1012	ENDIF
1013	IF( med_diag%OCN_K0%dgsave ) THEN
1014	CALL wrk_alloc( jpi, jpj, f_ocnk0_2d )
1015	f_ocnk0_2d(:,:) = 0.0 !!
1016	ENDIF
1017	IF( med_diag%CO2STARAIR%dgsave ) THEN
1018	CALL wrk_alloc( jpi, jpj, f_co2starair_2d )
1019	f_co2starair_2d(:,:) = 0.0 !!
1020	ENDIF
1021	IF( med_diag%OCN_DPCO2%dgsave ) THEN
1022	CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d )
1023	f_ocndpco2_2d(:,:) = 0.0 !!
1024	ENDIF
1025	# endif
1026	IF( med_diag%TPP3%dgsave ) THEN
1027	CALL wrk_alloc( jpi, jpj, jpk, tpp3d )
1028	tpp3d(:,:,:) = 0.0 !!
1029	ENDIF
1030	IF( med_diag%DETFLUX3%dgsave ) THEN
1031	CALL wrk_alloc( jpi, jpj, jpk, detflux3d )
1032	detflux3d(:,:,:) = 0.0 !!
1033	ENDIF
1034	IF( med_diag%REMIN3N%dgsave ) THEN
1035	CALL wrk_alloc( jpi, jpj, jpk, remin3dn )
1036	remin3dn(:,:,:) = 0.0 !!
1037	ENDIF
1038	!!
1039	!! AXY (10/11/16): CMIP6 diagnostics, 2D
1040	!! JPALM -- 17-11-16 -- put fgco2 alloc out of diag request
1041	!! needed for coupling/passed through restart
1042	!! IF( med_diag%FGCO2%dgsave ) THEN
1043	CALL wrk_alloc( jpi, jpj, fgco2 )
1044	fgco2(:,:) = 0.0 !!
1045	!! ENDIF
1046	IF( med_diag%INTDISSIC%dgsave ) THEN
1047	CALL wrk_alloc( jpi, jpj, intdissic )
1048	intdissic(:,:) = 0.0 !!
1049	ENDIF
1050	IF( med_diag%INTDISSIN%dgsave ) THEN
1051	CALL wrk_alloc( jpi, jpj, intdissin )
1052	intdissin(:,:) = 0.0 !!
1053	ENDIF
1054	IF( med_diag%INTDISSISI%dgsave ) THEN
1055	CALL wrk_alloc( jpi, jpj, intdissisi )
1056	intdissisi(:,:) = 0.0 !!
1057	ENDIF
1058	IF( med_diag%INTTALK%dgsave ) THEN
1059	CALL wrk_alloc( jpi, jpj, inttalk )
1060	inttalk(:,:) = 0.0 !!
1061	ENDIF
1062	IF( med_diag%O2min%dgsave ) THEN
1063	CALL wrk_alloc( jpi, jpj, o2min )
1064	o2min(:,:) = 1.e3 !! set to high value as we're looking for min(o2)
1065	ENDIF
1066	IF( med_diag%ZO2min%dgsave ) THEN
1067	CALL wrk_alloc( jpi, jpj, zo2min )
1068	zo2min(:,:) = 0.0 !!
1069	ENDIF
1070	IF( med_diag%FBDDTALK%dgsave ) THEN
1071	CALL wrk_alloc( jpi, jpj, fbddtalk )
1072	fbddtalk(:,:) = 0.0 !!
1073	ENDIF
1074	IF( med_diag%FBDDTDIC%dgsave ) THEN
1075	CALL wrk_alloc( jpi, jpj, fbddtdic )
1076	fbddtdic(:,:) = 0.0 !!
1077	ENDIF
1078	IF( med_diag%FBDDTDIFE%dgsave ) THEN
1079	CALL wrk_alloc( jpi, jpj, fbddtdife )
1080	fbddtdife(:,:) = 0.0 !!
1081	ENDIF
1082	IF( med_diag%FBDDTDIN%dgsave ) THEN
1083	CALL wrk_alloc( jpi, jpj, fbddtdin )
1084	fbddtdin(:,:) = 0.0 !!
1085	ENDIF
1086	IF( med_diag%FBDDTDISI%dgsave ) THEN
1087	CALL wrk_alloc( jpi, jpj, fbddtdisi )
1088	fbddtdisi(:,:) = 0.0 !!
1089	ENDIF
1090	!!
1091	!! AXY (10/11/16): CMIP6 diagnostics, 3D
1092	IF( med_diag%TPPD3%dgsave ) THEN
1093	CALL wrk_alloc( jpi, jpj, jpk, tppd3 )
1094	tppd3(:,:,:) = 0.0 !!
1095	ENDIF
1096	IF( med_diag%BDDTALK3%dgsave ) THEN
1097	CALL wrk_alloc( jpi, jpj, jpk, bddtalk3 )
1098	bddtalk3(:,:,:) = 0.0 !!
1099	ENDIF
1100	IF( med_diag%BDDTDIC3%dgsave ) THEN
1101	CALL wrk_alloc( jpi, jpj, jpk, bddtdic3 )
1102	bddtdic3(:,:,:) = 0.0 !!
1103	ENDIF
1104	IF( med_diag%BDDTDIFE3%dgsave ) THEN
1105	CALL wrk_alloc( jpi, jpj, jpk, bddtdife3 )
1106	bddtdife3(:,:,:) = 0.0 !!
1107	ENDIF
1108	IF( med_diag%BDDTDIN3%dgsave ) THEN
1109	CALL wrk_alloc( jpi, jpj, jpk, bddtdin3 )
1110	bddtdin3(:,:,:) = 0.0 !!
1111	ENDIF
1112	IF( med_diag%BDDTDISI3%dgsave ) THEN
1113	CALL wrk_alloc( jpi, jpj, jpk, bddtdisi3 )
1114	bddtdisi3(:,:,:) = 0.0 !!
1115	ENDIF
1116	IF( med_diag%FD_NIT3%dgsave ) THEN
1117	CALL wrk_alloc( jpi, jpj, jpk, fd_nit3 )
1118	fd_nit3(:,:,:) = 0.0 !!
1119	ENDIF
1120	IF( med_diag%FD_SIL3%dgsave ) THEN
1121	CALL wrk_alloc( jpi, jpj, jpk, fd_sil3 )
1122	fd_sil3(:,:,:) = 0.0 !!
1123	ENDIF
1124	IF( med_diag%FD_CAR3%dgsave ) THEN
1125	CALL wrk_alloc( jpi, jpj, jpk, fd_car3 )
1126	fd_car3(:,:,:) = 0.0 !!
1127	ENDIF
1128	IF( med_diag%FD_CAL3%dgsave ) THEN
1129	CALL wrk_alloc( jpi, jpj, jpk, fd_cal3 )
1130	fd_cal3(:,:,:) = 0.0 !!
1131	ENDIF
1132	IF( med_diag%DCALC3%dgsave ) THEN
1133	CALL wrk_alloc( jpi, jpj, jpk, dcalc3 )
1134	dcalc3(:,:,: ) = 0.0 !!
1135	ENDIF
1136	IF( med_diag%EXPC3%dgsave ) THEN
1137	CALL wrk_alloc( jpi, jpj, jpk, expc3 )
1138	expc3(:,:,: ) = 0.0 !!
1139	ENDIF
1140	IF( med_diag%EXPN3%dgsave ) THEN
1141	CALL wrk_alloc( jpi, jpj, jpk, expn3 )
1142	expn3(:,:,: ) = 0.0 !!
1143	ENDIF
1144	IF( med_diag%FEDISS3%dgsave ) THEN
1145	CALL wrk_alloc( jpi, jpj, jpk, fediss3 )
1146	fediss3(:,:,: ) = 0.0 !!
1147	ENDIF
1148	IF( med_diag%FESCAV3%dgsave ) THEN
1149	CALL wrk_alloc( jpi, jpj, jpk, fescav3 )
1150	fescav3(:,:,: ) = 0.0 !!
1151	ENDIF
1152	IF( med_diag%MIGRAZP3%dgsave ) THEN
1153	CALL wrk_alloc( jpi, jpj, jpk, migrazp3 )
1154	migrazp3(:,:,: ) = 0.0 !!
1155	ENDIF
1156	IF( med_diag%MIGRAZD3%dgsave ) THEN
1157	CALL wrk_alloc( jpi, jpj, jpk, migrazd3 )
1158	migrazd3(:,:,: ) = 0.0 !!
1159	ENDIF
1160	IF( med_diag%MEGRAZP3%dgsave ) THEN
1161	CALL wrk_alloc( jpi, jpj, jpk, megrazp3 )
1162	megrazp3(:,:,: ) = 0.0 !!
1163	ENDIF
1164	IF( med_diag%MEGRAZD3%dgsave ) THEN
1165	CALL wrk_alloc( jpi, jpj, jpk, megrazd3 )
1166	megrazd3(:,:,: ) = 0.0 !!
1167	ENDIF
1168	IF( med_diag%MEGRAZZ3%dgsave ) THEN
1169	CALL wrk_alloc( jpi, jpj, jpk, megrazz3 )
1170	megrazz3(:,:,: ) = 0.0 !!
1171	ENDIF
1172	IF( med_diag%O2SAT3%dgsave ) THEN
1173	CALL wrk_alloc( jpi, jpj, jpk, o2sat3 )
1174	o2sat3(:,:,: ) = 0.0 !!
1175	ENDIF
1176	IF( med_diag%PBSI3%dgsave ) THEN
1177	CALL wrk_alloc( jpi, jpj, jpk, pbsi3 )
1178	pbsi3(:,:,: ) = 0.0 !!
1179	ENDIF
1180	IF( med_diag%PCAL3%dgsave ) THEN
1181	CALL wrk_alloc( jpi, jpj, jpk, pcal3 )
1182	pcal3(:,:,: ) = 0.0 !!
1183	ENDIF
1184	IF( med_diag%REMOC3%dgsave ) THEN
1185	CALL wrk_alloc( jpi, jpj, jpk, remoc3 )
1186	remoc3(:,:,: ) = 0.0 !!
1187	ENDIF
1188	IF( med_diag%PNLIMJ3%dgsave ) THEN
1189	CALL wrk_alloc( jpi, jpj, jpk, pnlimj3 )
1190	pnlimj3(:,:,: ) = 0.0 !!
1191	ENDIF
1192	IF( med_diag%PNLIMN3%dgsave ) THEN
1193	CALL wrk_alloc( jpi, jpj, jpk, pnlimn3 )
1194	pnlimn3(:,:,: ) = 0.0 !!
1195	ENDIF
1196	IF( med_diag%PNLIMFE3%dgsave ) THEN
1197	CALL wrk_alloc( jpi, jpj, jpk, pnlimfe3 )
1198	pnlimfe3(:,:,: ) = 0.0 !!
1199	ENDIF
1200	IF( med_diag%PDLIMJ3%dgsave ) THEN
1201	CALL wrk_alloc( jpi, jpj, jpk, pdlimj3 )
1202	pdlimj3(:,:,: ) = 0.0 !!
1203	ENDIF
1204	IF( med_diag%PDLIMN3%dgsave ) THEN
1205	CALL wrk_alloc( jpi, jpj, jpk, pdlimn3 )
1206	pdlimn3(:,:,: ) = 0.0 !!
1207	ENDIF
1208	IF( med_diag%PDLIMFE3%dgsave ) THEN
1209	CALL wrk_alloc( jpi, jpj, jpk, pdlimfe3 )
1210	pdlimfe3(:,:,: ) = 0.0 !!
1211	ENDIF
1212	IF( med_diag%PDLIMSI3%dgsave ) THEN
1213	CALL wrk_alloc( jpi, jpj, jpk, pdlimsi3 )
1214	pdlimsi3(:,:,: ) = 0.0 !!
1215	ENDIF
1216	ENDIF !! lk_iomput
1217	!!
1218	# if defined key_axy_nancheck
1219	DO jn = 1,jptra
1220	fq2 = SUM(trn(:,:,:,jn))
1221	if ( ieee_is_nan( fq2 ) ) then
1222	!! there's a NaN here
1223	if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:'
1224	DO jk = 1,jpk
1225	DO jj = 1,jpj
1226	DO ji = 1,jpi
1227	if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
1228	if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', &
1229	& tmask(ji,jj,jk), ji, jj, jk, jn
1230	endif
1231	enddo
1232	enddo
1233	enddo
1234	CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
1235	endif
1236	ENDDO
1237	CALL flush(numout)
1238	# endif
1239
1240	# if defined key_debug_medusa
1241	IF ( lwp ) write (numout,*) 'trc_bio_medusa: variables initialised and checked'
1242	CALL flush(numout)
1243	# endif
1244
1245	# if defined key_roam
1246	!!----------------------------------------------------------------------
1247	!! calculate atmospheric pCO2
1248	!!----------------------------------------------------------------------
1249	!!
1250	# if defined key_axy_pi_co2
1251	f_xco2a = 284.725 !! OCMIP pre-industrial pCO2
1252	# else
1253	f_xco2a = 284.725 !! OCMIP pre-industrial pCO2
1254	# endif
1255	IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a
1256	# endif
1257
1258	# if defined key_debug_medusa
1259	IF ( lwp ) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry'
1260	IF ( lwp ) write (numout,*) 'trc_bio_medusa: kt = ', kt
1261	IF ( lwp ) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000
1262	CALL flush(numout)
1263	# endif
1264
1265	!!======================================================================
1266	!! AXY (07/04/17): possible subroutine block; ocean interior carbonate chemistry
1267	!!======================================================================
1268	# if defined key_roam
1269	!! AXY (20/11/14): alter to call on first MEDUSA timestep and then every
1270	!! month (this is hardwired as 960 timesteps but should
1271	!! be calculated and done properly
1272	!! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN
1273	!! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN
1274	!!=============================
1275	!! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call :
1276	!! we don't want to call on the first time-step of all run submission,
1277	!! but only on the very first time-step, and then every month
1278	!! So we call on nittrc000 if not restarted run,
1279	!! else if one month after last call.
1280	!! assume one month is 30d --> 36002430 : 2592000s
1281	!! try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt)
1282	!! ++ need to pass carb-chem output var through restarts
1283	IF ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN
1284	!!----------------------------------------------------------------------
1285	!! Calculate the carbonate chemistry for the whole ocean on the first
1286	!! simulation timestep and every month subsequently; the resulting 3D
1287	!! field of omega calcite is used to determine the depth of the CCD
1288	!!----------------------------------------------------------------------
1289	!!
1290	IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt
1291	CALL flush(numout)
1292	!! blank flags
1293	i2_omcal(:,:) = 0
1294	i2_omarg(:,:) = 0
1295	!! loop over 3D space
1296	DO jk = 1,jpk
1297	DO jj = 2,jpjm1
1298	DO ji = 2,jpim1
1299	!! OPEN wet point IF..THEN loop
1300	if (tmask(ji,jj,jk).eq.1) then
1301	IF ( lk_oasis ) THEN
1302	f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling
1303	ENDIF
1304	!! AXY (06/04/17): where am I?
1305	flatx = gphit(ji,jj)
1306	!! do carbonate chemistry
1307	!!
1308	fdep2 = fsdept(ji,jj,jk) !! set up level midpoint
1309	!! AXY (28/11/16): seafloor depth; previously mbathy(ji,jj) - 1, now mbathy(ji,jj)
1310	jmbathy = mbathy(ji,jj)
1311	!!
1312	!! set up required state variables
1313	zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon
1314	zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity
1315	ztmp = tsn(ji,jj,jk,jp_tem) !! temperature
1316	zsal = tsn(ji,jj,jk,jp_sal) !! salinity
1317	zsil = max(0.,trn(ji,jj,jk,jpsil)) !! silicic acid
1318	zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
1319	!!
1320	!! AXY (28/02/14): check input fields
1321	if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then
1322	IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 3D, ', &
1323	tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', &
1324	ji, ',', jj, ',', jk, ') at time', kt
1325	IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', &
1326	tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
1327	ztmp = tsb(ji,jj,jk,jp_tem) !! temperature
1328	endif
1329	if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then
1330	IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 3D, ', &
1331	tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', &
1332	ji, ',', jj, ',', jk, ') at time', kt
1333	endif
1334	!!
1335	!! blank input variables not used at this stage (they relate to air-sea flux)
1336	f_kw660 = 1.0
1337	f_pp0 = 1.0
1338	!!
1339	!! calculate carbonate chemistry at grid cell midpoint
1340	!! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
1341	!! chemistry package
1342	CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs
1343	f_pp0, fdep2, flatx, f_kw660, f_xco2a, 1, & ! inputs
1344	f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs
1345	f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs
1346	f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs
1347	f_co2starair, f_co2flux, f_dpco2 ) ! outputs
1348	!!
1349	f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg
1350	f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 -> ueq / kg
1351	f_dcf = f_rhosw
1352	!!
1353	!! store 3D outputs
1354	f3_pH(ji,jj,jk) = f_ph
1355	f3_h2co3(ji,jj,jk) = f_h2co3
1356	f3_hco3(ji,jj,jk) = f_hco3
1357	f3_co3(ji,jj,jk) = f_co3
1358	f3_omcal(ji,jj,jk) = f_omcal(ji,jj)
1359	f3_omarg(ji,jj,jk) = f_omarg(ji,jj)
1360	!!
1361	!! CCD calculation: calcite
1362	if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then
1363	if (jk .eq. 1) then
1364	f2_ccd_cal(ji,jj) = fdep2
1365	else
1366	fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj)
1367	fq1 = f3_omcal(ji,jj,jk-1) - 1.0
1368	fq2 = fq1 / (fq0 + tiny(fq0))
1369	fq3 = fdep2 - fsdept(ji,jj,jk-1)
1370	fq4 = fq2 * fq3
1371	f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
1372	endif
1373	i2_omcal(ji,jj) = 1
1374	endif
1375	if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then
1376	!! reached seafloor and still no dissolution; set to seafloor (W-point)
1377	f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1)
1378	i2_omcal(ji,jj) = 1
1379	endif
1380	!!
1381	!! CCD calculation: aragonite
1382	if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then
1383	if (jk .eq. 1) then
1384	f2_ccd_arg(ji,jj) = fdep2
1385	else
1386	fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj)
1387	fq1 = f3_omarg(ji,jj,jk-1) - 1.0
1388	fq2 = fq1 / (fq0 + tiny(fq0))
1389	fq3 = fdep2 - fsdept(ji,jj,jk-1)
1390	fq4 = fq2 * fq3
1391	f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
1392	endif
1393	i2_omarg(ji,jj) = 1
1394	endif
1395	if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then
1396	!! reached seafloor and still no dissolution; set to seafloor (W-point)
1397	f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1)
1398	i2_omarg(ji,jj) = 1
1399	endif
1400	endif
1401	ENDDO
1402	ENDDO
1403	ENDDO
1404	ENDIF
1405	# endif
1406
1407	# if defined key_debug_medusa
1408	IF ( lwp ) write (numout,*) 'trc_bio_medusa: ready for full domain calculations'
1409	CALL flush(numout)
1410	# endif
1411
1412	!!----------------------------------------------------------------------
1413	!! MEDUSA has unified equation through the water column
1414	!! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)
1415	!! Statement below in LOBSTER is different: DO jk = 1, jpkbm1
1416	!!----------------------------------------------------------------------
1417	!!
1418	!! NOTE: the ordering of the loops below differs from that of some other
1419	!! models; looping over the vertical dimension is the outermost loop and
1420	!! this complicates some calculations (e.g. storage of vertical fluxes
1421	!! that can otherwise be done via a singular variable require 2D fields
1422	!! here); however, these issues are relatively easily resolved, but the
1423	!! loops CANNOT be reordered without potentially causing code efficiency
1424	!! problems (e.g. array indexing means that reordering the loops would
1425	!! require skipping between widely-spaced memory location; potentially
1426	!! outside those immediately cached)
1427	!!
1428	!! OPEN vertical loop
1429	DO jk = 1,jpk
1430	!! OPEN horizontal loops
1431	DO jj = 2,jpjm1
1432	DO ji = 2,jpim1
1433	!! OPEN wet point IF..THEN loop
1434	if (tmask(ji,jj,jk).eq.1) then
1435	!!======================================================================
1436	!! SETUP LOCAL GRID CELL
1437	!!======================================================================
1438	!!
1439	!!---------------------------------------------------------------------
1440	!! Some notes on grid vertical structure
1441	!! - fsdepw(ji,jj,jk) is the depth of the upper surface of level jk
1442	!! - fsde3w(ji,jj,jk) is approximately the midpoint of level jk
1443	!! - fse3t(ji,jj,jk) is the thickness of level jk
1444	!!---------------------------------------------------------------------
1445	!!
1446	!! AXY (11/12/08): set up level thickness
1447	fthk = fse3t(ji,jj,jk)
1448	!! AXY (25/02/10): set up level depth (top of level)
1449	fdep = fsdepw(ji,jj,jk)
1450	!! AXY (01/03/10): set up level depth (bottom of level)
1451	fdep1 = fdep + fthk
1452	!! AXY (17/05/13): where am I?
1453	flatx = gphit(ji,jj)
1454	flonx = glamt(ji,jj)
1455	!! AXY (28/11/16): local seafloor depth
1456	!! previously mbathy(ji,jj) - 1, now mbathy(ji,jj)
1457	jmbathy = mbathy(ji,jj)
1458	!!
1459	!! set up model tracers
1460	!! negative values of state variables are not allowed to
1461	!! contribute to the calculated fluxes
1462	zchn = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll
1463	zchd = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll
1464	zphn = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms
1465	zphd = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms
1466	zpds = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon
1467	!! AXY (28/01/10): probably need to take account of chl/biomass connection
1468	if (zchn.eq.0.) zphn = 0.
1469	if (zchd.eq.0.) zphd = 0.
1470	if (zphn.eq.0.) zchn = 0.
1471	if (zphd.eq.0.) zchd = 0.
1472	!! AXY (23/01/14): duh - why did I forget diatom silicon?
1473	if (zpds.eq.0.) zphd = 0.
1474	if (zphd.eq.0.) zpds = 0.
1475	zzmi = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton
1476	zzme = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton
1477	zdet = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen
1478	zdin = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen
1479	zsil = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid
1480	zfer = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron"
1481	# if defined key_roam
1482	zdtc = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon
1483	zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon
1484	zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity
1485	zoxy = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen
1486	# if defined key_axy_carbchem
1487	zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
1488	# endif
1489	!!
1490	!! also need physical parameters for gas exchange calculations
1491	ztmp = tsn(ji,jj,jk,jp_tem)
1492	zsal = tsn(ji,jj,jk,jp_sal)
1493	!!
1494	!! AXY (28/02/14): check input fields
1495	if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then
1496	IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', &
1497	tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', &
1498	ji, ',', jj, ',', jk, ') at time', kt
1499	IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', &
1500	tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
1501	ztmp = tsb(ji,jj,jk,jp_tem) !! temperature
1502	endif
1503	if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then
1504	IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', &
1505	tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', &
1506	ji, ',', jj, ',', jk, ') at time', kt
1507	endif
1508	# else
1509	zdtc = zdet * xthetad !! implicit detrital carbon
1510	# endif
1511	# if defined key_debug_medusa
1512	if (idf.eq.1) then
1513	!! AXY (15/01/10)
1514	if (trn(ji,jj,jk,jpdin).lt.0.) then
1515	IF ( lwp ) write (numout,*) '------------------------------'
1516	IF ( lwp ) write (numout,*) 'NEGATIVE DIN ERROR =', trn(ji,jj,jk,jpdin)
1517	IF ( lwp ) write (numout,*) 'NEGATIVE DIN ERROR @', ji, jj, jk, kt
1518	endif
1519	if (trn(ji,jj,jk,jpsil).lt.0.) then
1520	IF ( lwp ) write (numout,*) '------------------------------'
1521	IF ( lwp ) write (numout,*) 'NEGATIVE SIL ERROR =', trn(ji,jj,jk,jpsil)
1522	IF ( lwp ) write (numout,*) 'NEGATIVE SIL ERROR @', ji, jj, jk, kt
1523	endif
1524	# if defined key_roam
1525	if (trn(ji,jj,jk,jpdic).lt.0.) then
1526	IF ( lwp ) write (numout,*) '------------------------------'
1527	IF ( lwp ) write (numout,*) 'NEGATIVE DIC ERROR =', trn(ji,jj,jk,jpdic)
1528	IF ( lwp ) write (numout,*) 'NEGATIVE DIC ERROR @', ji, jj, jk, kt
1529	endif
1530	if (trn(ji,jj,jk,jpalk).lt.0.) then
1531	IF ( lwp ) write (numout,*) '------------------------------'
1532	IF ( lwp ) write (numout,*) 'NEGATIVE ALK ERROR =', trn(ji,jj,jk,jpalk)
1533	IF ( lwp ) write (numout,*) 'NEGATIVE ALK ERROR @', ji, jj, jk, kt
1534	endif
1535	if (trn(ji,jj,jk,jpoxy).lt.0.) then
1536	IF ( lwp ) write (numout,*) '------------------------------'
1537	IF ( lwp ) write (numout,*) 'NEGATIVE OXY ERROR =', trn(ji,jj,jk,jpoxy)
1538	IF ( lwp ) write (numout,*) 'NEGATIVE OXY ERROR @', ji, jj, jk, kt
1539	endif
1540	# endif
1541	endif
1542	# endif
1543	!! sum tracers for inventory checks
1544	IF( lk_iomput ) THEN
1545	IF ( med_diag%INVTN%dgsave ) THEN
1546	ftot_n(ji,jj) = ftot_n(ji,jj) + &
1547	(fthk * ( zphn + zphd + zzmi + zzme + zdet + zdin ) )
1548	ENDIF
1549	IF ( med_diag%INVTSI%dgsave ) THEN
1550	ftot_si(ji,jj) = ftot_si(ji,jj) + &
1551	(fthk * ( zpds + zsil ) )
1552	ENDIF
1553	IF ( med_diag%INVTFE%dgsave ) THEN
1554	ftot_fe(ji,jj) = ftot_fe(ji,jj) + &
1555	(fthk * ( xrfn * ( zphn + zphd + zzmi + zzme + zdet ) + zfer ) )
1556	ENDIF
1557	# if defined key_roam
1558	IF ( med_diag%INVTC%dgsave ) THEN
1559	ftot_c(ji,jj) = ftot_c(ji,jj) + &
1560	(fthk * ( (xthetapn * zphn) + (xthetapd * zphd) + &
1561	(xthetazmi * zzmi) + (xthetazme * zzme) + zdtc + &
1562	zdic ) )
1563	ENDIF
1564	IF ( med_diag%INVTALK%dgsave ) THEN
1565	ftot_a(ji,jj) = ftot_a(ji,jj) + (fthk * ( zalk ) )
1566	ENDIF
1567	IF ( med_diag%INVTO2%dgsave ) THEN
1568	ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fthk * ( zoxy ) )
1569	ENDIF
1570	!!
1571	!! AXY (10/11/16): CMIP6 diagnostics
1572	IF ( med_diag%INTDISSIC%dgsave ) THEN
1573	intdissic(ji,jj) = intdissic(ji,jj) + (fthk * zdic)
1574	ENDIF
1575	IF ( med_diag%INTDISSIN%dgsave ) THEN
1576	intdissin(ji,jj) = intdissin(ji,jj) + (fthk * zdin)
1577	ENDIF
1578	IF ( med_diag%INTDISSISI%dgsave ) THEN
1579	intdissisi(ji,jj) = intdissisi(ji,jj) + (fthk * zsil)
1580	ENDIF
1581	IF ( med_diag%INTTALK%dgsave ) THEN
1582	inttalk(ji,jj) = inttalk(ji,jj) + (fthk * zalk)
1583	ENDIF
1584	IF ( med_diag%O2min%dgsave ) THEN
1585	if ( zoxy < o2min(ji,jj) ) then
1586	o2min(ji,jj) = zoxy
1587	IF ( med_diag%ZO2min%dgsave ) THEN
1588	zo2min(ji,jj) = (fdep + fdep1) / 2. !! layer midpoint
1589	ENDIF
1590	endif
1591	ENDIF
1592	# endif
1593	ENDIF
1594
1595	CALL flush(numout)
1596
1597	!!======================================================================
1598	!! LOCAL GRID CELL CALCULATIONS
1599	!!======================================================================
1600	!!
1601	!!======================================================================
1602	!! AXY (07/04/17): possible subroutine block; air-sea gas exchange
1603	!!======================================================================
1604	# if defined key_roam
1605	if ( jk .eq. 1 ) then
1606	!!----------------------------------------------------------------------
1607	!! Air-sea gas exchange
1608	!!----------------------------------------------------------------------
1609	!!
1610	!! AXY (17/07/14): zwind_i and zwind_j do not exist in this
1611	!! version of NEMO because it does not include
1612	!! the SBC changes that our local version has
1613	!! for accessing the HadGEM2 forcing; they
1614	!! could be added, but an alternative approach
1615	!! is to make use of wndm from oce_trc.F90
1616	!! which is wind speed at 10m (which is what
1617	!! is required here; this may need to be
1618	!! revisited when MEDUSA properly interacts
1619	!! with UKESM1 physics
1620	!!
1621	f_wind = wndm(ji,jj)
1622	IF ( lk_oasis ) THEN
1623	f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling
1624	ENDIF
1625	!!
1626	!! AXY (23/06/15): as part of an effort to update the carbonate chemistry
1627	!! in MEDUSA, the gas transfer velocity used in the carbon
1628	!! and oxygen cycles has been harmonised and is calculated
1629	!! by the same function here; this harmonisation includes
1630	!! changes to the PML carbonate chemistry scheme so that
1631	!! it too makes use of the same gas transfer velocity; the
1632	!! preferred parameterisation of this is Wanninkhof (2014),
1633	!! option 7
1634	!!
1635	# if defined key_debug_medusa
1636	IF ( lwp ) write (numout,*) 'trc_bio_medusa: entering gas_transfer'
1637	CALL flush(numout)
1638	# endif
1639	CALL gas_transfer( f_wind, 1, 7, & ! inputs
1640	f_kw660 ) ! outputs
1641	# if defined key_debug_medusa
1642	IF ( lwp ) write (numout,*) 'trc_bio_medusa: exiting gas_transfer'
1643	CALL flush(numout)
1644	# endif
1645	!!
1646	!! air pressure (atm); ultimately this will use air pressure at the base
1647	!! of the UKESM1 atmosphere
1648	!!
1649	f_pp0 = 1.0
1650	!!
1651	!! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp
1652	!! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj)
1653	!! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj)
1654	!! IF(lwp) WRITE(numout,*) ' MEDUSA f_wind =', f_wind
1655	!! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj)
1656	!!
1657	# if defined key_axy_carbchem
1658	!!
1659	!! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
1660	!! chemistry package; note that depth is set to
1661	!! zero in this call
1662	CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs
1663	f_pp0, 0.0, flatx, f_kw660, f_xco2a, 1, & ! inputs
1664	f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs
1665	f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs
1666	f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs
1667	f_co2starair, f_co2flux, f_dpco2 ) ! outputs
1668	!!
1669	f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg
1670	f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 -> ueq / kg
1671	f_dcf = f_rhosw
1672	# else
1673	!! AXY (18/04/13): switch off carbonate chemistry calculations; provide
1674	!! quasi-sensible alternatives
1675	f_ph = 8.1
1676	f_pco2w = f_xco2a
1677	f_h2co3 = 0.005 * zdic
1678	f_hco3 = 0.865 * zdic
1679	f_co3 = 0.130 * zdic
1680	f_omcal(ji,jj) = 4.
1681	f_omarg(ji,jj) = 2.
1682	f_co2flux = 0.
1683	f_TDIC = zdic
1684	f_TALK = zalk
1685	f_dcf = 1.026
1686	f_henry = 1.
1687	!! AXY (23/06/15): add in some extra MOCSY diagnostics
1688	f_fco2w = f_xco2a
1689	f_BetaD = 1.
1690	f_rhosw = 1.026
1691	f_opres = 0.
1692	f_insitut = ztmp
1693	f_pco2atm = f_xco2a
1694	f_fco2atm = f_xco2a
1695	f_schmidtco2 = 660.
1696	f_kwco2 = 0.
1697	f_K0 = 0.
1698	f_co2starair = f_xco2a
1699	f_dpco2 = 0.
1700	# endif
1701	!!
1702	!! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness
1703	f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux * 86400. / fthk
1704	!!
1705	!! oxygen (O2); OCMIP-2 code
1706	!! AXY (23/06/15): amend input list for oxygen to account for common gas
1707	!! transfer velocity
1708	!! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs
1709	!! f_kw660, f_o2flux, f_o2sat ) ! outputs
1710	CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy, & ! inputs
1711	f_kwo2, f_o2flux, f_o2sat ) ! outputs
1712	!!
1713	!! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness
1714	f_o2flux = (1. - fr_i(ji,jj)) * f_o2flux * 86400. / fthk
1715	!!
1716	!! Jpalm (08-2014)
1717	!! DMS surface concentration calculation; initialy added using MET-OFFICE subroutine
1718	!! air-sea flux calculated in atmospheric chemistry from atm and ocn concentrations
1719	!!
1720	!! AXY (13/03/15): this is amended to calculate all of the DMS
1721	!! estimates examined during UKESM1 (see comments
1722	!! in trcdms_medusa.F90)
1723	!!
1724	!! AXY (28/03/17): amended to pass DIN limitation instead of DIN concentration;
1725	!! accounts for differences in nutrient half-saturations; changes
1726	!! also made in trc_dms_medusa
1727	!!
1728	IF (jdms .eq. 1) THEN
1729	!!
1730	!! calculate weighted half-saturation for DIN uptake
1731	dms_wtkn = ((zphn * xnln) + (zphd * xnld)) / (zphn + zphd)
1732	!!
1733	!! feed in correct inputs
1734	if (jdms_input .eq. 0) then
1735	!! use instantaneous inputs
1736	dms_nlim = zdin / (zdin + dms_wtkn)
1737	!!
1738	CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), dms_nlim, & ! inputs
1739	dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs
1740	else
1741	!! use diel-average inputs
1742	dms_nlim = zn_dms_din(ji,jj) / (zn_dms_din(ji,jj) + dms_wtkn)
1743	!!
1744	CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs
1745	zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), dms_nlim, & ! inputs
1746	dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs
1747	endif
1748	!!
1749	!! assign correct output to variable passed to atmosphere
1750	if (jdms_model .eq. 1) then
1751	dms_surf = dms_andr
1752	elseif (jdms_model .eq. 2) then
1753	dms_surf = dms_simo
1754	elseif (jdms_model .eq. 3) then
1755	dms_surf = dms_aran
1756	elseif (jdms_model .eq. 4) then
1757	dms_surf = dms_hall
1758	endif
1759	!!
1760	!! 2D diag through iom_use
1761	IF( lk_iomput ) THEN
1762	IF( med_diag%DMS_SURF%dgsave ) THEN
1763	dms_surf2d(ji,jj) = dms_surf
1764	ENDIF
1765	IF( med_diag%DMS_ANDR%dgsave ) THEN
1766	dms_andr2d(ji,jj) = dms_andr
1767	ENDIF
1768	IF( med_diag%DMS_SIMO%dgsave ) THEN
1769	dms_simo2d(ji,jj) = dms_simo
1770	ENDIF
1771	IF( med_diag%DMS_ARAN%dgsave ) THEN
1772	dms_aran2d(ji,jj) = dms_aran
1773	ENDIF
1774	IF( med_diag%DMS_HALL%dgsave ) THEN
1775	dms_hall2d(ji,jj) = dms_hall
1776	ENDIF
1777	# if defined key_debug_medusa
1778	IF ( lwp ) write (numout,*) 'trc_bio_medusa: finish calculating dms'
1779	CALL flush(numout)
1780	# endif
1781	ENDIF
1782	!! End iom
1783	ENDIF
1784	!! End DMS Loop
1785	!!
1786	!! store 2D outputs
1787	!!
1788	!! JPALM -- 17-11-16 -- put fgco2 out of diag request
1789	!! is needed for coupling; pass through restart
1790	!! IF( med_diag%FGCO2%dgsave ) THEN
1791	!! convert from mol/m2/day to kg/m2/s
1792	fgco2(ji,jj) = f_co2flux * fthk * CO2flux_conv !! mmol-C/m3/d -> kg-CO2/m2/s
1793	!! ENDIF
1794	IF ( lk_iomput ) THEN
1795	IF( med_diag%ATM_PCO2%dgsave ) THEN
1796	f_pco2a2d(ji,jj) = f_pco2atm
1797	ENDIF
1798	IF( med_diag%OCN_PCO2%dgsave ) THEN
1799	f_pco2w2d(ji,jj) = f_pco2w
1800	ENDIF
1801	IF( med_diag%CO2FLUX%dgsave ) THEN
1802	f_co2flux2d(ji,jj) = f_co2flux * fthk !! mmol/m3/d -> mmol/m2/d
1803	ENDIF
1804	IF( med_diag%TCO2%dgsave ) THEN
1805	f_TDIC2d(ji,jj) = f_TDIC
1806	ENDIF
1807	IF( med_diag%TALK%dgsave ) THEN
1808	f_TALK2d(ji,jj) = f_TALK
1809	ENDIF
1810	IF( med_diag%KW660%dgsave ) THEN
1811	f_kw6602d(ji,jj) = f_kw660
1812	ENDIF
1813	IF( med_diag%ATM_PP0%dgsave ) THEN
1814	f_pp02d(ji,jj) = f_pp0
1815	ENDIF
1816	IF( med_diag%O2FLUX%dgsave ) THEN
1817	f_o2flux2d(ji,jj) = f_o2flux
1818	ENDIF
1819	IF( med_diag%O2SAT%dgsave ) THEN
1820	f_o2sat2d(ji,jj) = f_o2sat
1821	ENDIF
1822	!! AXY (24/11/16): add in extra MOCSY diagnostics
1823	IF( med_diag%ATM_XCO2%dgsave ) THEN
1824	f_xco2a_2d(ji,jj) = f_xco2a
1825	ENDIF
1826	IF( med_diag%OCN_FCO2%dgsave ) THEN
1827	f_fco2w_2d(ji,jj) = f_fco2w
1828	ENDIF
1829	IF( med_diag%ATM_FCO2%dgsave ) THEN
1830	f_fco2a_2d(ji,jj) = f_fco2atm
1831	ENDIF
1832	IF( med_diag%OCN_RHOSW%dgsave ) THEN
1833	f_ocnrhosw_2d(ji,jj) = f_rhosw
1834	ENDIF
1835	IF( med_diag%OCN_SCHCO2%dgsave ) THEN
1836	f_ocnschco2_2d(ji,jj) = f_schmidtco2
1837	ENDIF
1838	IF( med_diag%OCN_KWCO2%dgsave ) THEN
1839	f_ocnkwco2_2d(ji,jj) = f_kwco2
1840	ENDIF
1841	IF( med_diag%OCN_K0%dgsave ) THEN
1842	f_ocnk0_2d(ji,jj) = f_K0
1843	ENDIF
1844	IF( med_diag%CO2STARAIR%dgsave ) THEN
1845	f_co2starair_2d(ji,jj) = f_co2starair
1846	ENDIF
1847	IF( med_diag%OCN_DPCO2%dgsave ) THEN
1848	f_ocndpco2_2d(ji,jj) = f_dpco2
1849	ENDIF
1850	ENDIF
1851	!!
1852	endif
1853	!! End jk = 1 loop within ROAM key
1854
1855	!! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
1856	IF ( med_diag%O2SAT3%dgsave ) THEN
1857	call oxy_sato( ztmp, zsal, f_o2sat3 )
1858	o2sat3(ji, jj, jk) = f_o2sat3
1859	ENDIF
1860
1861	# endif
1862
1863	!!======================================================================
1864	!! AXY (07/04/17): possible subroutine block; riverine inputs (or delete; it's unused presently)
1865	!!======================================================================
1866	if ( jk .eq. 1 ) then
1867	!!----------------------------------------------------------------------
1868	!! River inputs
1869	!!----------------------------------------------------------------------
1870	!!
1871	!! runoff comes in as kg / m2 / s
1872	!! used and written out as m3 / m2 / d (= m / d)
1873	!! where 1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d
1874	!!
1875	!! AXY (17/07/14): the compiler doesn't like this line for some reason;
1876	!! as MEDUSA doesn't even use runoff for riverine inputs,
1877	!! a temporary solution is to switch off runoff entirely
1878	!! here; again, this change is one of several that will
1879	!! need revisiting once MEDUSA has bedded down in UKESM1;
1880	!! particularly so if the land scheme provides information
1881	!! concerning nutrient fluxes
1882	!!
1883	!! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24.
1884	f_runoff(ji,jj) = 0.0
1885	!!
1886	!! nutrients are added via rivers to the model in one of two ways:
1887	!! 1. via river concentration; i.e. the average nutrient concentration
1888	!! of a river water is described by a spatial file, and this is
1889	!! multiplied by runoff to give a nutrient flux
1890	!! 2. via direct river flux; i.e. the average nutrient flux due to
1891	!! rivers is described by a spatial file, and this is simply applied
1892	!! as a direct nutrient flux (i.e. it does not relate or respond to
1893	!! model runoff)
1894	!! nutrient fields are derived from the GlobalNEWS 2 database; carbon and
1895	!! alkalinity are derived from continent-scale DIC estimates (Huang et al.,
1896	!! 2012) and some Arctic river alkalinity estimates (Katya?)
1897	!!
1898	!! as of 19/07/12, riverine nutrients can now be spread vertically across
1899	!! several grid cells rather than just poured into the surface box; this
1900	!! block of code is still executed, however, to set up the total amounts
1901	!! of nutrient entering via rivers
1902	!!
1903	!! nitrogen
1904	if (jriver_n .eq. 1) then
1905	!! river concentration specified; use runoff to calculate input
1906	f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj)
1907	elseif (jriver_n .eq. 2) then
1908	!! river flux specified; independent of runoff
1909	f_riv_n(ji,jj) = riv_n(ji,jj)
1910	endif
1911	!!
1912	!! silicon
1913	if (jriver_si .eq. 1) then
1914	!! river concentration specified; use runoff to calculate input
1915	f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj)
1916	elseif (jriver_si .eq. 2) then
1917	!! river flux specified; independent of runoff
1918	f_riv_si(ji,jj) = riv_si(ji,jj)
1919	endif
1920	!!
1921	!! carbon
1922	if (jriver_c .eq. 1) then
1923	!! river concentration specified; use runoff to calculate input
1924	f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj)
1925	elseif (jriver_c .eq. 2) then
1926	!! river flux specified; independent of runoff
1927	f_riv_c(ji,jj) = riv_c(ji,jj)
1928	endif
1929	!!
1930	!! alkalinity
1931	if (jriver_alk .eq. 1) then
1932	!! river concentration specified; use runoff to calculate input
1933	f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj)
1934	elseif (jriver_alk .eq. 2) then
1935	!! river flux specified; independent of runoff
1936	f_riv_alk(ji,jj) = riv_alk(ji,jj)
1937	endif
1938
1939	endif
1940
1941	!!======================================================================
1942	!! AXY (07/04/17): possible subroutine block; phytoplankton growth
1943	!!======================================================================
1944
1945	!!----------------------------------------------------------------------
1946	!! Chlorophyll calculations
1947	!!----------------------------------------------------------------------
1948	!!
1949	!! non-diatoms
1950	if (zphn.GT.rsmall) then
1951	fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn)))
1952	faln = xaln * fthetan
1953	else
1954	fthetan = 0.
1955	faln = 0.
1956	endif
1957	!!
1958	!! diatoms
1959	if (zphd.GT.rsmall) then
1960	fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd)))
1961	fald = xald * fthetad
1962	else
1963	fthetad = 0.
1964	fald = 0.
1965	endif
1966
1967	!!----------------------------------------------------------------------
1968	!! Phytoplankton light limitation
1969	!!----------------------------------------------------------------------
1970	!!
1971	!! It is assumed xpar is the depth-averaged (vertical layer) PAR
1972	!! Light limitation (check self-shading) in W/m2
1973	!!
1974	!! Note that there is no temperature dependence in phytoplankton
1975	!! growth rate or any other function.
1976	!! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid
1977	!! NaNs in case of Phy==0.
1978	!!
1979	!! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and non-diat:
1980	!! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012
1981	!!
1982	!! AXY (16/07/09)
1983	!! temperature for new Eppley style phytoplankton growth
1984	loc_T = tsn(ji,jj,jk,jp_tem)
1985	fun_T = 1.066*(1.0 loc_T)
1986	!! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for
1987	!phytoplankton
1988	!! growth; remin. unaffected
1989	fun_Q10 = xq10**((loc_T - 0.0) / 10.0)
1990	if (jphy.eq.1) then
1991	xvpnT = xvpn * fun_T
1992	xvpdT = xvpd * fun_T
1993	elseif (jphy.eq.2) then
1994	xvpnT = xvpn * fun_Q10
1995	xvpdT = xvpd * fun_Q10
1996	else
1997	xvpnT = xvpn
1998	xvpdT = xvpd
1999	endif
2000	!!
2001	!! non-diatoms
2002	fchn1 = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk))
2003	if (fchn1.GT.rsmall) then
2004	fchn = xvpnT / (sqrt(fchn1) + tiny(fchn1))
2005	else
2006	fchn = 0.
2007	endif
2008	fjln = fchn * faln * xpar(ji,jj,jk) !! non-diatom J term
2009	fjlim_pn = fjln / xvpnT
2010	!!
2011	!! diatoms
2012	fchd1 = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk))
2013	if (fchd1.GT.rsmall) then
2014	fchd = xvpdT / (sqrt(fchd1) + tiny(fchd1))
2015	else
2016	fchd = 0.
2017	endif
2018	fjld = fchd * fald * xpar(ji,jj,jk) !! diatom J term
2019	fjlim_pd = fjld / xvpdT
2020
2021	!!----------------------------------------------------------------------
2022	!! Phytoplankton nutrient limitation
2023	!!----------------------------------------------------------------------
2024	!!
2025	!! non-diatoms (N, Fe)
2026	fnln = zdin / (zdin + xnln) !! non-diatom Qn term
2027	ffln = zfer / (zfer + xfln) !! non-diatom Qf term
2028	!!
2029	!! diatoms (N, Si, Fe)
2030	fnld = zdin / (zdin + xnld) !! diatom Qn term
2031	fsld = zsil / (zsil + xsld) !! diatom Qs term
2032	ffld = zfer / (zfer + xfld) !! diatom Qf term
2033
2034	!!----------------------------------------------------------------------
2035	!! Primary production (non-diatoms)
2036	!! (note: still needs multiplying by phytoplankton concentration)
2037	!!----------------------------------------------------------------------
2038	!!
2039	if (jliebig .eq. 0) then
2040	!! multiplicative nutrient limitation
2041	fpnlim = fnln * ffln
2042	elseif (jliebig .eq. 1) then
2043	!! Liebig Law (= most limiting) nutrient limitation
2044	fpnlim = min(fnln, ffln)
2045	endif
2046	fprn = fjln * fpnlim
2047
2048	!!----------------------------------------------------------------------
2049	!! Primary production (diatoms)
2050	!! (note: still needs multiplying by phytoplankton concentration)
2051	!!
2052	!! production here is split between nitrogen production and that of
2053	!! silicon; depending upon the "intracellular" ratio of Si:N, model
2054	!! diatoms will uptake nitrogen/silicon differentially; this borrows
2055	!! from the diatom model of Mongin et al. (2006)
2056	!!----------------------------------------------------------------------
2057	!!
2058	if (jliebig .eq. 0) then
2059	!! multiplicative nutrient limitation
2060	fpdlim = fnld * ffld
2061	elseif (jliebig .eq. 1) then
2062	!! Liebig Law (= most limiting) nutrient limitation
2063	fpdlim = min(fnld, ffld)
2064	endif
2065	!!
2066	if (zphd.GT.rsmall .AND. zpds.GT.rsmall) then
2067	!! "intracellular" elemental ratios
2068	! fsin = zpds / (zphd + tiny(zphd))
2069	! fnsi = zphd / (zpds + tiny(zpds))
2070	fsin = 0.0
2071	IF( zphd .GT. rsmall) fsin = zpds / zphd
2072	fnsi = 0.0
2073	IF( zpds .GT. rsmall) fnsi = zphd / zpds
2074	!! AXY (23/02/10): these next variables derive from Mongin et al. (2003)
2075	fsin1 = 3.0 * xsin0 !! = 0.6
2076	fnsi1 = 1.0 / fsin1 !! = 1.667
2077	fnsi2 = 1.0 / xsin0 !! = 5.0
2078	!!
2079	!! conditionalities based on ratios
2080	!! nitrogen (and iron and carbon)
2081	if (fsin.le.xsin0) then
2082	fprd = 0.0
2083	fsld2 = 0.0
2084	elseif (fsin.lt.fsin1) then
2085	fprd = xuif * ((fsin - xsin0) / (fsin + tiny(fsin))) * (fjld * fpdlim)
2086	fsld2 = xuif * ((fsin - xsin0) / (fsin + tiny(fsin)))
2087	elseif (fsin.ge.fsin1) then
2088	fprd = (fjld * fpdlim)
2089	fsld2 = 1.0
2090	endif
2091	!!
2092	!! silicon
2093	if (fsin.lt.fnsi1) then
2094	fprds = (fjld * fsld)
2095	elseif (fsin.lt.fnsi2) then
2096	fprds = xuif * ((fnsi - xnsi0) / (fnsi + tiny(fnsi))) * (fjld * fsld)
2097	else
2098	fprds = 0.0
2099	endif
2100	else
2101	fsin = 0.0
2102	fnsi = 0.0
2103	fprd = 0.0
2104	fsld2 = 0.0
2105	fprds = 0.0
2106	endif
2107
2108	!!----------------------------------------------------------------------
2109	!! Mixed layer primary production
2110	!! this block calculates the amount of primary production that occurs
2111	!! within the upper mixed layer; this allows the separate diagnosis
2112	!! of "sub-surface" primary production; it does assume that short-
2113	!! term variability in mixed layer depth doesn't mess with things
2114	!! though
2115	!!----------------------------------------------------------------------
2116	!!
2117	if (fdep1.le.hmld(ji,jj)) then
2118	!! this level is entirely in the mixed layer
2119	fq0 = 1.0
2120	elseif (fdep.ge.hmld(ji,jj)) then
2121	!! this level is entirely below the mixed layer
2122	fq0 = 0.0
2123	else
2124	!! this level straddles the mixed layer
2125	fq0 = (hmld(ji,jj) - fdep) / fthk
2126	endif
2127	!!
2128	fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn * zphn * fthk * fq0)
2129	fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd * zphd * fthk * fq0)
2130
2131	!!----------------------------------------------------------------------
2132	!! Vertical Integral --
2133	!!----------------------------------------------------------------------
2134	ftot_pn(ji,jj) = ftot_pn(ji,jj) + (zphn * fthk) !! vertical integral non-diatom phytoplankton
2135	ftot_pd(ji,jj) = ftot_pd(ji,jj) + (zphd * fthk) !! vertical integral diatom phytoplankton
2136	ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi * fthk) !! vertical integral microzooplankton
2137	ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme * fthk) !! vertical integral mesozooplankton
2138	ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet * fthk) !! vertical integral slow detritus, nitrogen
2139	ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc * fthk) !! vertical integral slow detritus, carbon
2140
2141	!!----------------------------------------------------------------------
2142	!! More chlorophyll calculations
2143	!!----------------------------------------------------------------------
2144	!!
2145	!! frn = (xthetam / fthetan) * (fprn / (fthetan * xpar(ji,jj,jk)))
2146	!! frd = (xthetam / fthetad) * (fprd / (fthetad * xpar(ji,jj,jk)))
2147	frn = (xthetam * fchn * fnln * ffln ) / (fthetan + tiny(fthetan))
2148	!! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid
2149	!! limitation, is used in the following line to regulate chlorophyll
2150	!! growth in a manner that is inconsistent with its use in the regulation
2151	!! of biomass growth; the Mongin term term used in growth is more complex
2152	!! than the simple multiplicative function used below
2153	!! frd = (xthetam * fchd * fnld * ffld * fsld) / (fthetad + tiny(fthetad))
2154	!! AXY (12/05/09): this replacement line uses the new variable, fsld2, to
2155	!! regulate chlorophyll growth
2156	frd = (xthetamd * fchd * fnld * ffld * fsld2) / (fthetad + tiny(fthetad))
2157
2158	!!======================================================================
2159	!! AXY (07/04/17): possible subroutine block; zooplankton grazing
2160	!!======================================================================
2161
2162	!!----------------------------------------------------------------------
2163	!! Zooplankton Grazing
2164	!! this code supplements the base grazing model with one that
2165	!! considers the C:N ratio of grazed food and balances this against
2166	!! the requirements of zooplankton growth; this model is derived
2167	!! from that of Anderson & Pondaven (2003)
2168	!!
2169	!! the current version of the code assumes a fixed C:N ratio for
2170	!! detritus (in contrast to Anderson & Pondaven, 2003), though the
2171	!! full equations are retained for future extension
2172	!!----------------------------------------------------------------------
2173	!!
2174	!!----------------------------------------------------------------------
2175	!! Microzooplankton first
2176	!!----------------------------------------------------------------------
2177	!!
2178	fmi1 = (xkmi * xkmi) + (xpmipn * zphn * zphn) + (xpmid * zdet * zdet)
2179	fmi = xgmi * zzmi / fmi1
2180	fgmipn = fmi * xpmipn * zphn * zphn !! grazing on non-diatoms
2181	fgmid = fmi * xpmid * zdet * zdet !! grazing on detrital nitrogen
2182	# if defined key_roam
2183	fgmidc = rsmall !acc
2184	IF ( zdet .GT. rsmall ) fgmidc = (zdtc / (zdet + tiny(zdet))) * fgmid !! grazing on detrital carbon
2185	# else
2186	!! AXY (26/11/08): implicit detrital carbon change
2187	fgmidc = xthetad * fgmid !! grazing on detrital carbon
2188	# endif
2189	!!
2190	!! which translates to these incoming N and C fluxes
2191	finmi = (1.0 - xphi) * (fgmipn + fgmid)
2192	ficmi = (1.0 - xphi) * ((xthetapn * fgmipn) + fgmidc)
2193	!!
2194	!! the ideal food C:N ratio for microzooplankton
2195	!! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
2196	fstarmi = (xbetan * xthetazmi) / (xbetac * xkc)
2197	!!
2198	!! process these to determine proportioning of grazed N and C
2199	!! (since there is no explicit consideration of respiration,
2200	!! only growth and excretion are calculated here)
2201	fmith = (ficmi / (finmi + tiny(finmi)))
2202	if (fmith.ge.fstarmi) then
2203	fmigrow = xbetan * finmi
2204	fmiexcr = 0.0
2205	else
2206	fmigrow = (xbetac * xkc * ficmi) / xthetazmi
2207	fmiexcr = ficmi * ((xbetan / (fmith + tiny(fmith))) - ((xbetac * xkc) / xthetazmi))
2208	endif
2209	# if defined key_roam
2210	fmiresp = (xbetac * ficmi) - (xthetazmi * fmigrow)
2211	# endif
2212
2213	!!----------------------------------------------------------------------
2214	!! Mesozooplankton second
2215	!!----------------------------------------------------------------------
2216	!!
2217	fme1 = (xkme * xkme) + (xpmepn * zphn * zphn) + (xpmepd * zphd * zphd) + &
2218	(xpmezmi * zzmi * zzmi) + (xpmed * zdet * zdet)
2219	fme = xgme * zzme / fme1
2220	fgmepn = fme * xpmepn * zphn * zphn !! grazing on non-diatoms
2221	fgmepd = fme * xpmepd * zphd * zphd !! grazing on diatoms
2222	fgmepds = fsin * fgmepd !! grazing on diatom silicon
2223	fgmezmi = fme * xpmezmi * zzmi * zzmi !! grazing on microzooplankton
2224	fgmed = fme * xpmed * zdet * zdet !! grazing on detrital nitrogen
2225	# if defined key_roam
2226	fgmedc = rsmall !acc
2227	IF ( zdet .GT. rsmall ) fgmedc = (zdtc / (zdet + tiny(zdet))) * fgmed !! grazing on detrital carbon
2228	# else
2229	!! AXY (26/11/08): implicit detrital carbon change
2230	fgmedc = xthetad * fgmed !! grazing on detrital carbon
2231	# endif
2232	!!
2233	!! which translates to these incoming N and C fluxes
2234	finme = (1.0 - xphi) * (fgmepn + fgmepd + fgmezmi + fgmed)
2235	ficme = (1.0 - xphi) * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + &
2236	(xthetazmi * fgmezmi) + fgmedc)
2237	!!
2238	!! the ideal food C:N ratio for mesozooplankton
2239	!! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
2240	fstarme = (xbetan * xthetazme) / (xbetac * xkc)
2241	!!
2242	!! process these to determine proportioning of grazed N and C
2243	!! (since there is no explicit consideration of respiration,
2244	!! only growth and excretion are calculated here)
2245	fmeth = (ficme / (finme + tiny(finme)))
2246	if (fmeth.ge.fstarme) then
2247	fmegrow = xbetan * finme
2248	fmeexcr = 0.0
2249	else
2250	fmegrow = (xbetac * xkc * ficme) / xthetazme
2251	fmeexcr = ficme * ((xbetan / (fmeth + tiny(fmeth))) - ((xbetac * xkc) / xthetazme))
2252	endif
2253	# if defined key_roam
2254	fmeresp = (xbetac * ficme) - (xthetazme * fmegrow)
2255	# endif
2256
2257	fzmi_i(ji,jj) = fzmi_i(ji,jj) + fthk * ( &
2258	fgmipn + fgmid )
2259	fzmi_o(ji,jj) = fzmi_o(ji,jj) + fthk * ( &
2260	fmigrow + (xphi * (fgmipn + fgmid)) + fmiexcr + ((1.0 - xbetan) * finmi) )
2261	fzme_i(ji,jj) = fzme_i(ji,jj) + fthk * ( &
2262	fgmepn + fgmepd + fgmezmi + fgmed )
2263	fzme_o(ji,jj) = fzme_o(ji,jj) + fthk * ( &
2264	fmegrow + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + fmeexcr + ((1.0 - xbetan) * finme) )
2265
2266	!!======================================================================
2267	!! AXY (07/04/17): possible subroutine block; miscellaneous plankton losses
2268	!!======================================================================
2269
2270	!!----------------------------------------------------------------------
2271	!! Plankton metabolic losses
2272	!! Linear loss processes assumed to be metabolic in origin
2273	!!----------------------------------------------------------------------
2274	!!
2275	fdpn2 = xmetapn * zphn
2276	fdpd2 = xmetapd * zphd
2277	fdpds2 = xmetapd * zpds
2278	fdzmi2 = xmetazmi * zzmi
2279	fdzme2 = xmetazme * zzme
2280
2281	!!----------------------------------------------------------------------
2282	!! Plankton mortality losses
2283	!! EKP (26/02/09): phytoplankton hyperbolic mortality term introduced
2284	!! to improve performance in gyres
2285	!!----------------------------------------------------------------------
2286	!!
2287	!! non-diatom phytoplankton
2288	if (jmpn.eq.1) fdpn = xmpn * zphn !! linear
2289	if (jmpn.eq.2) fdpn = xmpn * zphn * zphn !! quadratic
2290	if (jmpn.eq.3) fdpn = xmpn * zphn * & !! hyperbolic
2291	(zphn / (xkphn + zphn))
2292	if (jmpn.eq.4) fdpn = xmpn * zphn * & !! sigmoid
2293	((zphn * zphn) / (xkphn + (zphn * zphn)))
2294	!!
2295	!! diatom phytoplankton
2296	if (jmpd.eq.1) fdpd = xmpd * zphd !! linear
2297	if (jmpd.eq.2) fdpd = xmpd * zphd * zphd !! quadratic
2298	if (jmpd.eq.3) fdpd = xmpd * zphd * & !! hyperbolic
2299	(zphd / (xkphd + zphd))
2300	if (jmpd.eq.4) fdpd = xmpd * zphd * & !! sigmoid
2301	((zphd * zphd) / (xkphd + (zphd * zphd)))
2302	fdpds = fdpd * fsin
2303	!!
2304	!! microzooplankton
2305	if (jmzmi.eq.1) fdzmi = xmzmi * zzmi !! linear
2306	if (jmzmi.eq.2) fdzmi = xmzmi * zzmi * zzmi !! quadratic
2307	if (jmzmi.eq.3) fdzmi = xmzmi * zzmi * & !! hyperbolic
2308	(zzmi / (xkzmi + zzmi))
2309	if (jmzmi.eq.4) fdzmi = xmzmi * zzmi * & !! sigmoid
2310	((zzmi * zzmi) / (xkzmi + (zzmi * zzmi)))
2311	!!
2312	!! mesozooplankton
2313	if (jmzme.eq.1) fdzme = xmzme * zzme !! linear
2314	if (jmzme.eq.2) fdzme = xmzme * zzme * zzme !! quadratic
2315	if (jmzme.eq.3) fdzme = xmzme * zzme * & !! hyperbolic
2316	(zzme / (xkzme + zzme))
2317	if (jmzme.eq.4) fdzme = xmzme * zzme * & !! sigmoid
2318	((zzme * zzme) / (xkzme + (zzme * zzme)))
2319
2320	!!======================================================================
2321	!! AXY (07/04/17): possible subroutine block; detritus processes (fuse with later?)
2322	!!======================================================================
2323
2324	!!----------------------------------------------------------------------
2325	!! Detritus remineralisation
2326	!! Constant or temperature-dependent
2327	!!----------------------------------------------------------------------
2328	!!
2329	if (jmd.eq.1) then
2330	!! temperature-dependent
2331	fdd = xmd * fun_T * zdet
2332	# if defined key_roam
2333	fddc = xmdc * fun_T * zdtc
2334	# endif
2335	elseif (jmd.eq.2) then
2336	!! AXY (16/05/13): add in Q10-based parameterisation (def in nmlst)
2337	!! temperature-dependent
2338	fdd = xmd * fun_Q10 * zdet
2339	#if defined key_roam
2340	fddc = xmdc * fun_Q10 * zdtc
2341	#endif
2342	else
2343	!! temperature-independent
2344	fdd = xmd * zdet
2345	# if defined key_roam
2346	fddc = xmdc * zdtc
2347	# endif
2348	endif
2349	!!
2350	!! AXY (22/07/09): accelerate detrital remineralisation in the bottom box
2351	if ((jk.eq.jmbathy) .and. jsfd.eq.1) then
2352	fdd = 1.0 * zdet
2353	# if defined key_roam
2354	fddc = 1.0 * zdtc
2355	# endif
2356	endif
2357
2358	!!----------------------------------------------------------------------
2359	!! Detritus addition to benthos
2360	!! If activated, slow detritus in the bottom box will enter the
2361	!! benthic pool
2362	!!----------------------------------------------------------------------
2363	!!
2364	if ((jk.eq.jmbathy) .and. jorgben.eq.1) then
2365	!! this is the BOTTOM OCEAN BOX -> into the benthic pool!
2366	!!
2367	f_sbenin_n(ji,jj) = (zdet * vsed * 86400.)
2368	f_sbenin_fe(ji,jj) = (zdet * vsed * 86400. * xrfn)
2369	# if defined key_roam
2370	f_sbenin_c(ji,jj) = (zdtc * vsed * 86400.)
2371	# else
2372	f_sbenin_c(ji,jj) = (zdet * vsed * 86400. * xthetad)
2373	# endif
2374	endif
2375
2376	!!======================================================================
2377	!! AXY (07/04/17): possible subroutine block; iron chemistry and scavenging
2378	!!======================================================================
2379
2380	!!----------------------------------------------------------------------
2381	!! Iron chemistry and fractionation
2382	!! following the Parekh et al. (2004) scheme adopted by the Met.
2383	!! Office, Medusa models total iron but considers "free" and
2384	!! ligand-bound forms for the purposes of scavenging (only "free"
2385	!! iron can be scavenged
2386	!!----------------------------------------------------------------------
2387	!!
2388	!! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
2389	xFeT = zfer * 1.e3
2390	!!
2391	!! calculate fractionation (based on Diat-HadOCC; in turn based on Parekh et al., 2004)
2392	xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0
2393	xb2M4ac = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0)
2394	!!
2395	!! "free" ligand concentration
2396	xLgF = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL
2397	!!
2398	!! ligand-bound iron concentration
2399	xFeL = xLgT - xLgF
2400	!!
2401	!! "free" iron concentration (and convert to mmol Fe / m3)
2402	xFeF = (xFeT - xFeL) * 1.e-3
2403	xFree(ji,jj)= xFeF / (zfer + tiny(zfer))
2404	!!
2405	!! scavenging of iron
2406	!! AXY (05/04/17): formerly several schemes, now the only appropriate
2407	!! one, with all other options returning no scavenging (trc_nam_medusa
2408	!! reports on this)
2409	!!
2410	if (jiron.eq.1) then
2411	!!----------------------------------------------------------------------
2412	!! Scheme 1: Dutkiewicz et al. (2005)
2413	!! This scheme includes a single scavenging term based solely on a
2414	!! fixed rate and the availablility of "free" iron
2415	!!----------------------------------------------------------------------
2416	!!
2417	ffescav = xk_sc_Fe * xFeF ! = mmol/m3/d
2418	!!
2419	!!----------------------------------------------------------------------
2420	!!
2421	!! Mick's code contains a further (optional) implicit "scavenging" of
2422	!! iron that sets an upper bound on "free" iron concentration, and
2423	!! essentially caps the concentration of total iron as xFeL + "free"
2424	!! iron; since the former is constrained by a fixed total ligand
2425	!! concentration (= 1.0 umol/m3), and the latter isn't allowed above
2426	!! this upper bound, total iron is constrained to a maximum of ...
2427	!!
2428	!! xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3
2429	!!
2430	!! In Mick's code, the actual value of total iron is reset to this
2431	!! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't
2432	!! our favoured approach to tracer updating here (not least because
2433	!! of the leapfrog), so here the amount scavenged is augmented by an
2434	!! additional amount that serves to drag total iron back towards that
2435	!! expected from this limitation on iron concentration ...
2436	!!
2437	xmaxFeF = min((xFeF * 1.e3), 0.3) ! = umol/m3
2438	!!
2439	!! Here, the difference between current total Fe and (FeL + Fe') is
2440	!! calculated and added to the scavenging flux already calculated
2441	!! above ...
2442	!!
2443	fdeltaFe = (xFeT - (xFeL + xmaxFeF)) * 1.e-3 ! = mmol/m3
2444	!!
2445	!! This assumes that the "excess" iron is dissipated with a time-
2446	!! scale of 1 day; seems reasonable to me ... (famous last words)
2447	!!
2448	ffescav = ffescav + fdeltaFe ! = mmol/m3/d
2449	!!
2450	# if defined key_deep_fe_fix
2451	!! AXY (17/01/13)
2452	!! stop scavenging for iron concentrations below 0.5 umol / m3
2453	!! at depths greater than 1000 m; this aims to end MEDUSA's
2454	!! continual loss of iron at depth without impacting things
2455	!! at the surface too much; the justification for this is that
2456	!! it appears to be what Mick Follows et al. do in their work
2457	!! (as evidenced by the iron initial condition they supplied
2458	!! me with); to be honest, it looks like Follow et al. do this
2459	!! at shallower depths than 1000 m, but I'll stick with this
2460	!! for now; I suspect that this seemingly arbitrary approach
2461	!! effectively "parameterises" the particle-based scavenging
2462	!! rates that other models use (i.e. at depth there are no
2463	!! sinking particles, so scavenging stops); it might be fun
2464	!! justifying this in a paper though!
2465	!!
2466	if ((fdep.gt.1000.) .and. (xFeT.lt.0.5)) then
2467	ffescav = 0.
2468	endif
2469	# endif
2470	else
2471	!!----------------------------------------------------------------------
2472	!! No Scheme: you coward!
2473	!! This scheme puts its head in the sand and eskews any decision about
2474	!! how iron is removed from the ocean; prepare to get deluged in iron
2475	!! you fool!
2476	!!----------------------------------------------------------------------
2477	ffescav = 0.
2478	endif
2479
2480	!!----------------------------------------------------------------------
2481	!! Other iron cycle processes
2482	!!----------------------------------------------------------------------
2483	!!
2484	!! aeolian iron deposition
2485	if (jk.eq.1) then
2486	!! zirondep is in mmol-Fe / m2 / day
2487	!! ffetop is in mmol-dissolved-Fe / m3 / day
2488	ffetop = zirondep(ji,jj) * xfe_sol / fthk
2489	else
2490	ffetop = 0.0
2491	endif
2492	!!
2493	!! seafloor iron addition
2494	!! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down
2495	!! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then
2496	if ((jk.eq.jmbathy).AND.jk.le.i0500) then
2497	!! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a
2498	!! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value
2499	!! but apply it everywhere
2500	!! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 1-37)
2501	ffebot = (xfe_sed / fthk)
2502	else
2503	ffebot = 0.0
2504	endif
2505
2506	!!======================================================================
2507	!! AXY (07/04/17): possible subroutine block; miscellaneous processes (fuse?)
2508	!!======================================================================
2509
2510	!!----------------------------------------------------------------------
2511	!! Miscellaneous
2512	!!----------------------------------------------------------------------
2513	!!
2514	!! diatom frustule dissolution
2515	fsdiss = xsdiss * zpds
2516
2517	!!----------------------------------------------------------------------
2518	!! Slow detritus creation
2519	!!----------------------------------------------------------------------
2520	!! this variable integrates the creation of slow sinking detritus
2521	!! to allow the split between fast and slow detritus to be
2522	!! diagnosed
2523	fslown = fdpn + fdzmi + ((1.0 - xfdfrac1) * fdpd) + &
2524	((1.0 - xfdfrac2) * fdzme) + ((1.0 - xbetan) * (finmi + finme))
2525	!!
2526	!! this variable records the slow detrital sinking flux at this
2527	!! particular depth; it is used in the output of this flux at
2528	!! standard depths in the diagnostic outputs; needs to be
2529	!! adjusted from per second to per day because of parameter vsed
2530	fslownflux(ji,jj) = zdet * vsed * 86400.
2531	# if defined key_roam
2532	!!
2533	!! and the same for detrital carbon
2534	fslowc = (xthetapn * fdpn) + (xthetazmi * fdzmi) + &
2535	(xthetapd * (1.0 - xfdfrac1) * fdpd) + &
2536	(xthetazme * (1.0 - xfdfrac2) * fdzme) + &
2537	((1.0 - xbetac) * (ficmi + ficme))
2538	!!
2539	!! this variable records the slow detrital sinking flux at this
2540	!! particular depth; it is used in the output of this flux at
2541	!! standard depths in the diagnostic outputs; needs to be
2542	!! adjusted from per second to per day because of parameter vsed
2543	fslowcflux(ji,jj) = zdtc * vsed * 86400.
2544	# endif
2545
2546	!!----------------------------------------------------------------------
2547	!! Nutrient regeneration
2548	!! this variable integrates total nitrogen regeneration down the
2549	!! watercolumn; its value is stored and output as a 2D diagnostic;
2550	!! the corresponding dissolution flux of silicon (from sources
2551	!! other than fast detritus) is also integrated; note that,
2552	!! confusingly, the linear loss terms from plankton compartments
2553	!! are labelled as fdX2 when one might have expected fdX or fdX1
2554	!!----------------------------------------------------------------------
2555	!!
2556	!! nitrogen
2557	fregen = (( (xphi * (fgmipn + fgmid)) + & ! messy feeding
2558	(xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + & ! messy feeding
2559	fmiexcr + fmeexcr + fdd + & ! excretion + D remin.
2560	fdpn2 + fdpd2 + fdzmi2 + fdzme2) * fthk) ! linear mortality
2561	!!
2562	!! silicon
2563	fregensi = (( fsdiss + ((1.0 - xfdfrac1) * fdpds) + & ! dissolution + non-lin. mortality
2564	((1.0 - xfdfrac3) * fgmepds) + & ! egestion by zooplankton
2565	fdpds2) * fthk) ! linear mortality
2566	# if defined key_roam
2567	!!
2568	!! carbon
2569	fregenc = (( (xphi * ((xthetapn * fgmipn) + fgmidc)) + & ! messy feeding
2570	(xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & ! messy feeding
2571	(xthetazmi * fgmezmi) + fgmedc)) + & ! messy feeding
2572	fmiresp + fmeresp + fddc + & ! respiration + D remin.
2573	(xthetapn * fdpn2) + (xthetapd * fdpd2) + & ! linear mortality
2574	(xthetazmi * fdzmi2) + (xthetazme * fdzme2)) * fthk) ! linear mortality
2575	# endif
2576
2577
2578	!!======================================================================
2579	!! AXY (07/04/17): possible subroutine block; fast-sinking detritus
2580	!!======================================================================
2581
2582	!!----------------------------------------------------------------------
2583	!! Fast-sinking detritus terms
2584	!! "local" variables declared so that conservation can be checked;
2585	!! the calculated terms are added to the fast-sinking flux later on
2586	!! only after the flux entering this level has experienced some
2587	!! remineralisation
2588	!! note: these fluxes need to be scaled by the level thickness
2589	!!----------------------------------------------------------------------
2590	!!
2591	!! nitrogen: diatom and mesozooplankton mortality
2592	ftempn = b0 * ((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme))
2593	!!
2594	!! silicon: diatom mortality and grazed diatoms
2595	ftempsi = b0 * ((xfdfrac1 * fdpds) + (xfdfrac3 * fgmepds))
2596	!!
2597	!! iron: diatom and mesozooplankton mortality
2598	ftempfe = b0 * (((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) * xrfn)
2599	!!
2600	!! carbon: diatom and mesozooplankton mortality
2601	ftempc = b0 * ((xfdfrac1 * xthetapd * fdpd) + &
2602	(xfdfrac2 * xthetazme * fdzme))
2603	!!
2604	# if defined key_roam
2605	if (jrratio.eq.0) then
2606	!! CaCO3: latitudinally-based fraction of total primary production
2607	!! absolute latitude of current grid cell
2608	flat = abs(gphit(ji,jj))
2609	!! 0.10 at equator; 0.02 at pole
2610	fcaco3 = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0))
2611	elseif (jrratio.eq.1) then
2612	!! CaCO3: Ridgwell et al. (2007) submodel, version 1
2613	!! this uses SURFACE omega calcite to regulate rain ratio
2614	if (f_omcal(ji,jj).ge.1.0) then
2615	fq1 = (f_omcal(ji,jj) - 1.0)**0.81
2616	else
2617	fq1 = 0.
2618	endif
2619	fcaco3 = xridg_r0 * fq1
2620	elseif (jrratio.eq.2) then
2621	!! CaCO3: Ridgwell et al. (2007) submodel, version 2
2622	!! this uses FULL 3D omega calcite to regulate rain ratio
2623	if (f3_omcal(ji,jj,jk).ge.1.0) then
2624	fq1 = (f3_omcal(ji,jj,jk) - 1.0)**0.81
2625	else
2626	fq1 = 0.
2627	endif
2628	fcaco3 = xridg_r0 * fq1
2629	endif
2630	# else
2631	!! CaCO3: latitudinally-based fraction of total primary production
2632	!! absolute latitude of current grid cell
2633	flat = abs(gphit(ji,jj))
2634	!! 0.10 at equator; 0.02 at pole
2635	fcaco3 = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0))
2636	# endif
2637	!! AXY (09/03/09): convert CaCO3 production from function of
2638	!! primary production into a function of fast-sinking material;
2639	!! technically, this is what Dunne et al. (2007) do anyway; they
2640	!! convert total primary production estimated from surface
2641	!! chlorophyll to an export flux for which they apply conversion
2642	!! factors to estimate the various elemental fractions (Si, Ca)
2643	ftempca = ftempc * fcaco3
2644
2645	# if defined key_debug_medusa
2646	!! integrate total fast detritus production
2647	if (idf.eq.1) then
2648	fifd_n(ji,jj) = fifd_n(ji,jj) + (ftempn * fthk)
2649	fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi * fthk)
2650	fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe * fthk)
2651	# if defined key_roam
2652	fifd_c(ji,jj) = fifd_c(ji,jj) + (ftempc * fthk)
2653	# endif
2654	endif
2655	# endif
2656
2657	!!----------------------------------------------------------------------
2658	!! This version of MEDUSA offers a choice of three methods for
2659	!! handling the remineralisation of fast detritus. All three
2660	!! do so in broadly the same way:
2661	!!
2662	!! 1. Fast detritus is stored as a 2D array [ ffastX ]
2663	!! 2. Fast detritus is added level-by-level [ ftempX ]
2664	!! 3. Fast detritus is not remineralised in the top box [ freminX ]
2665	!! 4. Remaining fast detritus is remineralised in the bottom [ fsedX ]
2666	!! box
2667	!!
2668	!! The three remineralisation methods are:
2669	!!
2670	!! 1. Ballast model (i.e. that published in Yool et al., 2011)
2671	!! (1b. Ballast-sans-ballast model)
2672	!! 2. Martin et al. (1987)
2673	!! 3. Henson et al. (2011)
2674	!!
2675	!! The first of these couples C, N and Fe remineralisation to
2676	!! the remineralisation of particulate Si and CaCO3, but the
2677	!! latter two treat remineralisation of C, N, Fe, Si and CaCO3
2678	!! completely separately. At present a switch within the code
2679	!! regulates which submodel is used, but this should be moved
2680	!! to the namelist file.
2681	!!
2682	!! The ballast-sans-ballast submodel is an original development
2683	!! feature of MEDUSA in which the ballast submodel's general
2684	!! framework and parameterisation is used, but in which there
2685	!! is no protection of organic material afforded by ballasting
2686	!! minerals. While similar, it is not the same as the Martin
2687	!! et al. (1987) submodel.
2688	!!
2689	!! Since the three submodels behave the same in terms of
2690	!! accumulating sinking material and remineralising it all at
2691	!! the seafloor, these portions of the code below are common to
2692	!! all three.
2693	!!----------------------------------------------------------------------
2694
2695	if (jexport.eq.1) then
2696	!!======================================================================
2697	!! BALLAST SUBMODEL
2698	!!======================================================================
2699	!!
2700	!!----------------------------------------------------------------------
2701	!! Fast-sinking detritus fluxes, pt. 1: REMINERALISATION
2702	!! aside from explicitly modelled, slow-sinking detritus, the
2703	!! model includes an implicit representation of detrital
2704	!! particles that sink too quickly to be modelled with
2705	!! explicit state variables; this sinking flux is instead
2706	!! instantaneously remineralised down the water column using
2707	!! the version of Armstrong et al. (2002)'s ballast model
2708	!! used by Dunne et al. (2007); the version of this model
2709	!! here considers silicon and calcium carbonate ballast
2710	!! minerals; this section of the code redistributes the fast
2711	!! sinking material generated locally down the water column;
2712	!! this differs from Dunne et al. (2007) in that fast sinking
2713	!! material is distributed at every level below that it is
2714	!! generated, rather than at every level below some fixed
2715	!! depth; this scheme is also different in that sinking material
2716	!! generated in one level is aggregated with that generated by
2717	!! shallower levels; this should make the ballast model more
2718	!! self-consistent (famous last words)
2719	!!----------------------------------------------------------------------
2720	!!
2721	if (jk.eq.1) then
2722	!! this is the SURFACE OCEAN BOX (no remineralisation)
2723	!!
2724	freminc = 0.0
2725	freminn = 0.0
2726	freminfe = 0.0
2727	freminsi = 0.0
2728	freminca = 0.0
2729	elseif (jk.le.jmbathy) then
2730	!! this is an OCEAN BOX (remineralise some material)
2731	!!
2732	!! set up CCD depth to be used depending on user choice
2733	if (jocalccd.eq.0) then
2734	!! use default CCD field
2735	fccd_dep = ocal_ccd(ji,jj)
2736	elseif (jocalccd.eq.1) then
2737	!! use calculated CCD field
2738	fccd_dep = f2_ccd_cal(ji,jj)
2739	endif
2740	!!
2741	!! === organic carbon ===
2742	fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol)
2743	if (iball.eq.1) then
2744	fq1 = (fq0 * xmassc) !! how much it weighs (mass)
2745	fq2 = (ffastca(ji,jj) * xmassca) !! how much CaCO3 enters this box (mass)
2746	fq3 = (ffastsi(ji,jj) * xmasssi) !! how much opal enters this box (mass)
2747	fq4 = (fq2 * xprotca) + (fq3 * xprotsi) !! total protected organic C (mass)
2748	!! this next term is calculated for C but used for N and Fe as well
2749	!! it needs to be protected in case ALL C is protected
2750	if (fq4.lt.fq1) then
2751	fprotf = (fq4 / (fq1 + tiny(fq1))) !! protected fraction of total organic C (non-dim)
2752	else
2753	fprotf = 1.0 !! all organic C is protected (non-dim)
2754	endif
2755	fq5 = (1.0 - fprotf) !! unprotected fraction of total organic C (non-dim)
2756	fq6 = (fq0 * fq5) !! how much organic C is unprotected (mol)
2757	fq7 = (fq6 * exp(-(fthk / xfastc))) !! how much unprotected C leaves this box (mol)
2758	fq8 = (fq7 + (fq0 * fprotf)) !! how much total C leaves this box (mol)
2759	freminc = (fq0 - fq8) / fthk !! C remineralisation in this box (mol)
2760	ffastc(ji,jj) = fq8
2761	else
2762	fq1 = fq0 * exp(-(fthk / xfastc)) !! how much organic C leaves this box (mol)
2763	freminc = (fq0 - fq1) / fthk !! C remineralisation in this box (mol)
2764	ffastc(ji,jj) = fq1
2765	endif
2766	!!
2767	!! === organic nitrogen ===
2768	fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol)
2769	if (iball.eq.1) then
2770	fq5 = (1.0 - fprotf) !! unprotected fraction of total organic N (non-dim)
2771	fq6 = (fq0 * fq5) !! how much organic N is unprotected (mol)
2772	fq7 = (fq6 * exp(-(fthk / xfastc))) !! how much unprotected N leaves this box (mol)
2773	fq8 = (fq7 + (fq0 * fprotf)) !! how much total N leaves this box (mol)
2774	freminn = (fq0 - fq8) / fthk !! N remineralisation in this box (mol)
2775	ffastn(ji,jj) = fq8
2776	else
2777	fq1 = fq0 * exp(-(fthk / xfastc)) !! how much organic N leaves this box (mol)
2778	freminn = (fq0 - fq1) / fthk !! N remineralisation in this box (mol)
2779	ffastn(ji,jj) = fq1
2780	endif
2781	!!
2782	!! === organic iron ===
2783	fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol)
2784	if (iball.eq.1) then
2785	fq5 = (1.0 - fprotf) !! unprotected fraction of total organic Fe (non-dim)
2786	fq6 = (fq0 * fq5) !! how much organic Fe is unprotected (mol)
2787	fq7 = (fq6 * exp(-(fthk / xfastc))) !! how much unprotected Fe leaves this box (mol)
2788	fq8 = (fq7 + (fq0 * fprotf)) !! how much total Fe leaves this box (mol)
2789	freminfe = (fq0 - fq8) / fthk !! Fe remineralisation in this box (mol)
2790	ffastfe(ji,jj) = fq8
2791	else
2792	fq1 = fq0 * exp(-(fthk / xfastc)) !! how much total Fe leaves this box (mol)
2793	freminfe = (fq0 - fq1) / fthk !! Fe remineralisation in this box (mol)
2794	ffastfe(ji,jj) = fq1
2795	endif
2796	!!
2797	!! === biogenic silicon ===
2798	fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol)
2799	fq1 = fq0 * exp(-(fthk / xfastsi)) !! how much opal leaves this box (mol)
2800	freminsi = (fq0 - fq1) / fthk !! Si remineralisation in this box (mol)
2801	ffastsi(ji,jj) = fq1
2802	!!
2803	!! === biogenic calcium carbonate ===
2804	fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol)
2805	if (fdep.le.fccd_dep) then
2806	!! whole grid cell above CCD
2807	fq1 = fq0 !! above lysocline, no Ca dissolves (mol)
2808	freminca = 0.0 !! above lysocline, no Ca dissolves (mol)
2809	fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#)
2810	elseif (fdep.ge.fccd_dep) then
2811	!! whole grid cell below CCD
2812	fq1 = fq0 * exp(-(fthk / xfastca)) !! how much CaCO3 leaves this box (mol)
2813	freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol)
2814	else
2815	!! partial grid cell below CCD
2816	fq2 = fdep1 - fccd_dep !! amount of grid cell below CCD (m)
2817	fq1 = fq0 * exp(-(fq2 / xfastca)) !! how much CaCO3 leaves this box (mol)
2818	freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol)
2819	endif
2820	ffastca(ji,jj) = fq1
2821	else
2822	!! this is BELOW THE LAST OCEAN BOX (do nothing)
2823	freminc = 0.0
2824	freminn = 0.0
2825	freminfe = 0.0
2826	freminsi = 0.0
2827	freminca = 0.0
2828	endif
2829
2830	elseif (jexport.eq.2.or.jexport.eq.3) then
2831	if (jexport.eq.2) then
2832	!!======================================================================
2833	!! MARTIN ET AL. (1987) SUBMODEL
2834	!!======================================================================
2835	!!
2836	!!----------------------------------------------------------------------
2837	!! This submodel uses the classic Martin et al. (1987) curve
2838	!! to determine the attenuation of fast-sinking detritus down
2839	!! the water column. All three organic elements, C, N and Fe,
2840	!! are handled identically, and their quantities in sinking
2841	!! particles attenuate according to a power relationship
2842	!! governed by parameter "b". This is assigned a canonical
2843	!! value of -0.858. Biogenic opal and calcium carbonate are
2844	!! attentuated using the same function as in the ballast
2845	!! submodel
2846	!!----------------------------------------------------------------------
2847	!!
2848	fb_val = -0.858
2849	elseif (jexport.eq.3) then
2850	!!======================================================================
2851	!! HENSON ET AL. (2011) SUBMODEL
2852	!!======================================================================
2853	!!
2854	!!----------------------------------------------------------------------
2855	!! This submodel reconfigures the Martin et al. (1987) curve by
2856	!! allowing the "b" value to vary geographically. Its value is
2857	!! set, following Henson et al. (2011), as a function of local
2858	!! sea surface temperature:
2859	!! b = -1.06 + (0.024 * SST)
2860	!! This means that remineralisation length scales are longer in
2861	!! warm, tropical areas and shorter in cold, polar areas. This
2862	!! does seem back-to-front (i.e. one would expect GREATER
2863	!! remineralisation in warmer waters), but is an outcome of
2864	!! analysis of sediment trap data, and it may reflect details
2865	!! of ecosystem structure that pertain to particle production
2866	!! rather than simply Q10.
2867	!!----------------------------------------------------------------------
2868	!!
2869	fl_sst = tsn(ji,jj,1,jp_tem)
2870	fb_val = -1.06 + (0.024 * fl_sst)
2871	endif
2872	!!
2873	if (jk.eq.1) then
2874	!! this is the SURFACE OCEAN BOX (no remineralisation)
2875	!!
2876	freminc = 0.0
2877	freminn = 0.0
2878	freminfe = 0.0
2879	freminsi = 0.0
2880	freminca = 0.0
2881	elseif (jk.le.jmbathy) then
2882	!! this is an OCEAN BOX (remineralise some material)
2883	!!
2884	!! === organic carbon ===
2885	fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol)
2886	fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic C leaves this box (mol)
2887	freminc = (fq0 - fq1) / fthk !! C remineralisation in this box (mol)
2888	ffastc(ji,jj) = fq1
2889	!!
2890	!! === organic nitrogen ===
2891	fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol)
2892	fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic N leaves this box (mol)
2893	freminn = (fq0 - fq1) / fthk !! N remineralisation in this box (mol)
2894	ffastn(ji,jj) = fq1
2895	!!
2896	!! === organic iron ===
2897	fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol)
2898	fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic Fe leaves this box (mol)
2899	freminfe = (fq0 - fq1) / fthk !! Fe remineralisation in this box (mol)
2900	ffastfe(ji,jj) = fq1
2901	!!
2902	!! === biogenic silicon ===
2903	fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol)
2904	fq1 = fq0 * exp(-(fthk / xfastsi)) !! how much opal leaves this box (mol)
2905	freminsi = (fq0 - fq1) / fthk !! Si remineralisation in this box (mol)
2906	ffastsi(ji,jj) = fq1
2907	!!
2908	!! === biogenic calcium carbonate ===
2909	fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol)
2910	if (fdep.le.ocal_ccd(ji,jj)) then
2911	!! whole grid cell above CCD
2912	fq1 = fq0 !! above lysocline, no Ca dissolves (mol)
2913	freminca = 0.0 !! above lysocline, no Ca dissolves (mol)
2914	fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#)
2915	elseif (fdep.ge.ocal_ccd(ji,jj)) then
2916	!! whole grid cell below CCD
2917	fq1 = fq0 * exp(-(fthk / xfastca)) !! how much CaCO3 leaves this box (mol)
2918	freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol)
2919	else
2920	!! partial grid cell below CCD
2921	fq2 = fdep1 - ocal_ccd(ji,jj) !! amount of grid cell below CCD (m)
2922	fq1 = fq0 * exp(-(fq2 / xfastca)) !! how much CaCO3 leaves this box (mol)
2923	freminca = (fq0 - fq1) / fthk !! Ca remineralisation in this box (mol)
2924	endif
2925	ffastca(ji,jj) = fq1
2926	else
2927	!! this is BELOW THE LAST OCEAN BOX (do nothing)
2928	freminc = 0.0
2929	freminn = 0.0
2930	freminfe = 0.0
2931	freminsi = 0.0
2932	freminca = 0.0
2933	endif
2934
2935	endif
2936
2937	!!----------------------------------------------------------------------
2938	!! Fast-sinking detritus fluxes, pt. 2: UPDATE FAST FLUXES
2939	!! here locally calculated additions to the fast-sinking flux are added
2940	!! to the total fast-sinking flux; this is done here such that material
2941	!! produced in a particular layer is only remineralised below this
2942	!! layer
2943	!!----------------------------------------------------------------------
2944	!!
2945	!! add sinking material generated in this layer to running totals
2946	!!
2947	!! === organic carbon === (diatom and mesozooplankton mortality)
2948	ffastc(ji,jj) = ffastc(ji,jj) + (ftempc * fthk)
2949	!!
2950	!! === organic nitrogen === (diatom and mesozooplankton mortality)
2951	ffastn(ji,jj) = ffastn(ji,jj) + (ftempn * fthk)
2952	!!
2953	!! === organic iron === (diatom and mesozooplankton mortality)
2954	ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe * fthk)
2955	!!
2956	!! === biogenic silicon === (diatom mortality and grazed diatoms)
2957	ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi * fthk)
2958	!!
2959	!! === biogenic calcium carbonate === (latitudinally-based fraction of total primary production)
2960	ffastca(ji,jj) = ffastca(ji,jj) + (ftempca * fthk)
2961
2962	!!----------------------------------------------------------------------
2963	!! Fast-sinking detritus fluxes, pt. 3: SEAFLOOR
2964	!! remineralise all remaining fast-sinking detritus to dissolved
2965	!! nutrients; the sedimentation fluxes calculated here allow the
2966	!! separation of what's remineralised sinking through the final
2967	!! ocean box from that which is added to the final box by the
2968	!! remineralisation of material that reaches the seafloor (i.e.
2969	!! the model assumes that all material that hits the seafloor
2970	!! is remineralised and that none is permanently buried; hey,
2971	!! this is a giant GCM model that can't be run for long enough
2972	!! to deal with burial fluxes!)
2973	!!
2974	!! in a change to this process, in part so that MEDUSA behaves
2975	!! a little more like ERSEM et al., fast-sinking detritus (N, Fe
2976	!! and C) is converted to slow sinking detritus at the seafloor
2977	!! instead of being remineralised; the rationale is that in
2978	!! shallower shelf regions (... that are not fully mixed!) this
2979	!! allows the detrital material to return slowly to dissolved
2980	!! nutrient rather than instantaneously as now; the alternative
2981	!! would be to explicitly handle seafloor organic material - a
2982	!! headache I don't wish to experience at this point; note that
2983	!! fast-sinking Si and Ca detritus is just remineralised as
2984	!! per usual
2985	!!
2986	!! AXY (13/01/12)
2987	!! in a further change to this process, again so that MEDUSA is
2988	!! a little more like ERSEM et al., material that reaches the
2989	!! seafloor can now be added to sediment pools and stored for
2990	!! slow release; there are new 2D arrays for organic nitrogen,
2991	!! iron and carbon and inorganic silicon and carbon that allow
2992	!! fast and slow detritus that reaches the seafloor to be held
2993	!! and released back to the water column more slowly; these arrays
2994	!! are transferred via the tracer restart files between repeat
2995	!! submissions of the model
2996	!!----------------------------------------------------------------------
2997	!!
2998	ffast2slowc = 0.0
2999	ffast2slown = 0.0
3000	ffast2slowfe = 0.0
3001	!!
3002	if (jk.eq.jmbathy) then
3003	!! this is the BOTTOM OCEAN BOX (remineralise everything)
3004	!!
3005	!! AXY (17/01/12): tweaked to include benthos pools
3006	!!
3007	!! === organic carbon ===
3008	if (jfdfate.eq.0 .and. jorgben.eq.0) then
3009	freminc = freminc + (ffastc(ji,jj) / fthk) !! C remineralisation in this box (mol/m3)
3010	elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
3011	ffast2slowc = ffastc(ji,jj) / fthk !! fast C -> slow C (mol/m3)
3012	fslowc = fslowc + ffast2slowc
3013	elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
3014	f_fbenin_c(ji,jj) = ffastc(ji,jj) !! fast C -> benthic C (mol/m2)
3015	endif
3016	fsedc(ji,jj) = ffastc(ji,jj) !! record seafloor C (mol/m2)
3017	ffastc(ji,jj) = 0.0
3018	!!
3019	!! === organic nitrogen ===
3020	if (jfdfate.eq.0 .and. jorgben.eq.0) then
3021	freminn = freminn + (ffastn(ji,jj) / fthk) !! N remineralisation in this box (mol/m3)
3022	elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
3023	ffast2slown = ffastn(ji,jj) / fthk !! fast N -> slow N (mol/m3)
3024	fslown = fslown + ffast2slown
3025	elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
3026	f_fbenin_n(ji,jj) = ffastn(ji,jj) !! fast N -> benthic N (mol/m2)
3027	endif
3028	fsedn(ji,jj) = ffastn(ji,jj) !! record seafloor N (mol/m2)
3029	ffastn(ji,jj) = 0.0
3030	!!
3031	!! === organic iron ===
3032	if (jfdfate.eq.0 .and. jorgben.eq.0) then
3033	freminfe = freminfe + (ffastfe(ji,jj) / fthk) !! Fe remineralisation in this box (mol/m3)
3034	elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
3035	ffast2slowfe = ffastn(ji,jj) / fthk !! fast Fe -> slow Fe (mol/m3)
3036	elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
3037	f_fbenin_fe(ji,jj) = ffastfe(ji,jj) !! fast Fe -> benthic Fe (mol/m2)
3038	endif
3039	fsedfe(ji,jj) = ffastfe(ji,jj) !! record seafloor Fe (mol/m2)
3040	ffastfe(ji,jj) = 0.0
3041	!!
3042	!! === biogenic silicon ===
3043	if (jinorgben.eq.0) then
3044	freminsi = freminsi + (ffastsi(ji,jj) / fthk) !! Si remineralisation in this box (mol/m3)
3045	elseif (jinorgben.eq.1) then
3046	f_fbenin_si(ji,jj) = ffastsi(ji,jj) !! fast Si -> benthic Si (mol/m2)
3047	endif
3048	fsedsi(ji,jj) = ffastsi(ji,jj) !! record seafloor Si (mol/m2)
3049	ffastsi(ji,jj) = 0.0
3050	!!
3051	!! === biogenic calcium carbonate ===
3052	if (jinorgben.eq.0) then
3053	freminca = freminca + (ffastca(ji,jj) / fthk) !! Ca remineralisation in this box (mol/m3)
3054	elseif (jinorgben.eq.1) then
3055	f_fbenin_ca(ji,jj) = ffastca(ji,jj) !! fast Ca -> benthic Ca (mol/m2)
3056	endif
3057	fsedca(ji,jj) = ffastca(ji,jj) !! record seafloor Ca (mol/m2)
3058	ffastca(ji,jj) = 0.0
3059	endif
3060
3061	# if defined key_debug_medusa
3062	if (idf.eq.1) then
3063	!!----------------------------------------------------------------------
3064	!! Integrate total fast detritus remineralisation
3065	!!----------------------------------------------------------------------
3066	!!
3067	fofd_n(ji,jj) = fofd_n(ji,jj) + (freminn * fthk)
3068	fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi * fthk)
3069	fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe * fthk)
3070	# if defined key_roam
3071	fofd_c(ji,jj) = fofd_c(ji,jj) + (freminc * fthk)
3072	# endif
3073	endif
3074	# endif
3075
3076	!!----------------------------------------------------------------------
3077	!! Sort out remineralisation tally of fast-sinking detritus
3078	!!----------------------------------------------------------------------
3079	!!
3080	!! update fast-sinking regeneration arrays
3081	fregenfast(ji,jj) = fregenfast(ji,jj) + (freminn * fthk)
3082	fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi * fthk)
3083	# if defined key_roam
3084	fregenfastc(ji,jj) = fregenfastc(ji,jj) + (freminc * fthk)
3085	# endif
3086
3087	!!----------------------------------------------------------------------
3088	!! Benthic remineralisation fluxes
3089	!!----------------------------------------------------------------------
3090	!!
3091	if (jk.eq.jmbathy) then
3092	!!
3093	!! organic components
3094	if (jorgben.eq.1) then
3095	f_benout_n(ji,jj) = xsedn * zn_sed_n(ji,jj)
3096	f_benout_fe(ji,jj) = xsedfe * zn_sed_fe(ji,jj)
3097	f_benout_c(ji,jj) = xsedc * zn_sed_c(ji,jj)
3098	endif
3099	!!
3100	!! inorganic components
3101	if (jinorgben.eq.1) then
3102	f_benout_si(ji,jj) = xsedsi * zn_sed_si(ji,jj)
3103	f_benout_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
3104	!!
3105	!! account for CaCO3 that dissolves when it shouldn't
3106	if ( fdep .le. fccd_dep ) then
3107	f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
3108	endif
3109	endif
3110	endif
3111	CALL flush(numout)
3112
3113	!!======================================================================
3114	!! AXY (07/04/17): possible subroutine block; business and updating
3115	!!======================================================================
3116
3117	!!======================================================================
3118	!! LOCAL GRID CELL TRENDS
3119	!!======================================================================
3120	!!
3121	!!----------------------------------------------------------------------
3122	!! Determination of trends
3123	!!----------------------------------------------------------------------
3124	!!
3125	!!----------------------------------------------------------------------
3126	!! chlorophyll
3127	btra(jpchn) = b0 * ( &
3128	+ ((frn * fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2) * (fthetan / xxi) )
3129	btra(jpchd) = b0 * ( &
3130	+ ((frd * fprd * zphd) - fgmepd - fdpd - fdpd2) * (fthetad / xxi) )
3131	!!
3132	!!----------------------------------------------------------------------
3133	!! phytoplankton
3134	btra(jpphn) = b0 * ( &
3135	+ (fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2 )
3136	btra(jpphd) = b0 * ( &
3137	+ (fprd * zphd) - fgmepd - fdpd - fdpd2 )
3138	btra(jppds) = b0 * ( &
3139	+ (fprds * zpds) - fgmepds - fdpds - fsdiss - fdpds2 )
3140	!!
3141	!!----------------------------------------------------------------------
3142	!! zooplankton
3143	btra(jpzmi) = b0 * ( &
3144	+ fmigrow - fgmezmi - fdzmi - fdzmi2 )
3145	btra(jpzme) = b0 * ( &
3146	+ fmegrow - fdzme - fdzme2 )
3147	!!
3148	!!----------------------------------------------------------------------
3149	!! detritus
3150	btra(jpdet) = b0 * ( &
3151	+ fdpn + ((1.0 - xfdfrac1) * fdpd) & ! mort. losses
3152	+ fdzmi + ((1.0 - xfdfrac2) * fdzme) & ! mort. losses
3153	+ ((1.0 - xbetan) * (finmi + finme)) & ! assim. inefficiency
3154	- fgmid - fgmed - fdd & ! grazing and remin.
3155	+ ffast2slown ) ! seafloor fast->slow
3156	!!
3157	!!----------------------------------------------------------------------
3158	!! dissolved inorganic nitrogen nutrient
3159	fn_cons = 0.0 &
3160	- (fprn * zphn) - (fprd * zphd) ! primary production
3161	fn_prod = 0.0 &
3162	+ (xphi * (fgmipn + fgmid)) & ! messy feeding remin.
3163	+ (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) & ! messy feeding remin.
3164	+ fmiexcr + fmeexcr + fdd + freminn & ! excretion and remin.
3165	+ fdpn2 + fdpd2 + fdzmi2 + fdzme2 ! metab. losses
3166	!!
3167	!! riverine flux
3168	if ( jriver_n .gt. 0 ) then
3169	f_riv_loc_n = f_riv_n(ji,jj) * friver_dep(jk,jmbathy) / fthk
3170	fn_prod = fn_prod + f_riv_loc_n
3171	endif
3172	!!
3173	!! benthic remineralisation
3174	if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
3175	fn_prod = fn_prod + (f_benout_n(ji,jj) / fthk)
3176	endif
3177	!!
3178	btra(jpdin) = b0 * ( &
3179	fn_prod + fn_cons )
3180	!!
3181	fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved nitrogen
3182	fn_cons ) )
3183	fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fthk * ( & ! production of dissolved nitrogen
3184	fn_prod ) )
3185	!!
3186	!!----------------------------------------------------------------------
3187	!! dissolved silicic acid nutrient
3188	fs_cons = 0.0 &
3189	- (fprds * zpds) ! opal production
3190	fs_prod = 0.0 &
3191	+ fsdiss & ! opal dissolution
3192	+ ((1.0 - xfdfrac1) * fdpds) & ! mort. loss
3193	+ ((1.0 - xfdfrac3) * fgmepds) & ! egestion of grazed Si
3194	+ freminsi + fdpds2 ! fast diss. and metab. losses
3195	!!
3196	!! riverine flux
3197	if ( jriver_si .gt. 0 ) then
3198	f_riv_loc_si = f_riv_si(ji,jj) * friver_dep(jk,jmbathy) / fthk
3199	fs_prod = fs_prod + f_riv_loc_si
3200	endif
3201	!!
3202	!! benthic remineralisation
3203	if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
3204	fs_prod = fs_prod + (f_benout_si(ji,jj) / fthk)
3205	endif
3206	!!
3207	btra(jpsil) = b0 * ( &
3208	fs_prod + fs_cons )
3209	!!
3210	fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved silicon
3211	fs_cons ) )
3212	fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fthk * ( & ! production of dissolved silicon
3213	fs_prod ) )
3214	!!
3215	!!----------------------------------------------------------------------
3216	!! dissolved "iron" nutrient
3217	btra(jpfer) = b0 * ( &
3218	+ (xrfn * btra(jpdin)) + ffetop + ffebot - ffescav )
3219
3220	# if defined key_roam
3221	!!
3222	!!----------------------------------------------------------------------
3223	!! AXY (26/11/08): implicit detrital carbon change
3224	btra(jpdtc) = b0 * ( &
3225	+ (xthetapn * fdpn) + ((1.0 - xfdfrac1) * (xthetapd * fdpd)) & ! mort. losses
3226	+ (xthetazmi * fdzmi) + ((1.0 - xfdfrac2) * (xthetazme * fdzme)) & ! mort. losses
3227	+ ((1.0 - xbetac) * (ficmi + ficme)) & ! assim. inefficiency
3228	- fgmidc - fgmedc - fddc & ! grazing and remin.
3229	+ ffast2slowc ) ! seafloor fast->slow
3230	!!
3231	!!----------------------------------------------------------------------
3232	!! dissolved inorganic carbon
3233	fc_cons = 0.0 &
3234	- (xthetapn * fprn * zphn) - (xthetapd * fprd * zphd) ! primary production
3235	fc_prod = 0.0 &
3236	+ (xthetapn * xphi * fgmipn) + (xphi * fgmidc) & ! messy feeding remin
3237	+ (xthetapn * xphi * fgmepn) + (xthetapd * xphi * fgmepd) & ! messy feeding remin
3238	+ (xthetazmi * xphi * fgmezmi) + (xphi * fgmedc) & ! messy feeding remin
3239	+ fmiresp + fmeresp + fddc + freminc + (xthetapn * fdpn2) & ! resp., remin., losses
3240	+ (xthetapd * fdpd2) + (xthetazmi * fdzmi2) & ! losses
3241	+ (xthetazme * fdzme2) ! losses
3242	!!
3243	!! riverine flux
3244	if ( jriver_c .gt. 0 ) then
3245	f_riv_loc_c = f_riv_c(ji,jj) * friver_dep(jk,jmbathy) / fthk
3246	fc_prod = fc_prod + f_riv_loc_c
3247	endif
3248	!!
3249	!! benthic remineralisation
3250	if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
3251	fc_prod = fc_prod + (f_benout_c(ji,jj) / fthk)
3252	endif
3253	if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
3254	fc_prod = fc_prod + (f_benout_ca(ji,jj) / fthk)
3255	endif
3256	!!
3257	!! community respiration (does not include CaCO3 terms - obviously!)
3258	fcomm_resp(ji,jj) = fcomm_resp(ji,jj) + fc_prod
3259	!!
3260	!! CaCO3
3261	fc_prod = fc_prod - ftempca + freminca
3262	!!
3263	!! riverine flux
3264	if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then
3265	fc_prod = fc_prod + f_riv_c(ji,jj)
3266	endif
3267	!!
3268	btra(jpdic) = b0 * ( &
3269	fc_prod + fc_cons )
3270	!!
3271	fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved carbon
3272	fc_cons ) )
3273	fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fthk * ( & ! production of dissolved carbon
3274	fc_prod ) )
3275	!!
3276	!!----------------------------------------------------------------------
3277	!! alkalinity
3278	fa_prod = 0.0 &
3279	+ (2.0 * freminca) ! CaCO3 dissolution
3280	fa_cons = 0.0 &
3281	- (2.0 * ftempca) ! CaCO3 production
3282	!!
3283	!! riverine flux
3284	if ( jriver_alk .gt. 0 ) then
3285	f_riv_loc_alk = f_riv_alk(ji,jj) * friver_dep(jk,jmbathy) / fthk
3286	fa_prod = fa_prod + f_riv_loc_alk
3287	endif
3288	!!
3289	!! benthic remineralisation
3290	if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
3291	fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fthk)
3292	endif
3293	!!
3294	btra(jpalk) = b0 * ( &
3295	fa_prod + fa_cons )
3296	!!
3297	!!----------------------------------------------------------------------
3298	!! oxygen (has protection at low O2 concentrations; OCMIP-2 style)
3299	fo2_prod = 0.0 &
3300	+ (xthetanit * fprn * zphn) & ! Pn primary production, N
3301	+ (xthetanit * fprd * zphd) & ! Pd primary production, N
3302	+ (xthetarem * xthetapn * fprn * zphn) & ! Pn primary production, C
3303	+ (xthetarem * xthetapd * fprd * zphd) ! Pd primary production, C
3304	fo2_ncons = 0.0 &
3305	- (xthetanit * xphi * fgmipn) & ! Pn messy feeding remin., N
3306	- (xthetanit * xphi * fgmid) & ! D messy feeding remin., N
3307	- (xthetanit * xphi * fgmepn) & ! Pn messy feeding remin., N
3308	- (xthetanit * xphi * fgmepd) & ! Pd messy feeding remin., N
3309	- (xthetanit * xphi * fgmezmi) & ! Zi messy feeding remin., N
3310	- (xthetanit * xphi * fgmed) & ! D messy feeding remin., N
3311	- (xthetanit * fmiexcr) & ! microzoo excretion, N
3312	- (xthetanit * fmeexcr) & ! mesozoo excretion, N
3313	- (xthetanit * fdd) & ! slow detritus remin., N
3314	- (xthetanit * freminn) & ! fast detritus remin., N
3315	- (xthetanit * fdpn2) & ! Pn losses, N
3316	- (xthetanit * fdpd2) & ! Pd losses, N
3317	- (xthetanit * fdzmi2) & ! Zmi losses, N
3318	- (xthetanit * fdzme2) ! Zme losses, N
3319	!!
3320	!! benthic remineralisation
3321	if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
3322	fo2_ncons = fo2_ncons - (xthetanit * f_benout_n(ji,jj) / fthk)
3323	endif
3324	fo2_ccons = 0.0 &
3325	- (xthetarem * xthetapn * xphi * fgmipn) & ! Pn messy feeding remin., C
3326	- (xthetarem * xphi * fgmidc) & ! D messy feeding remin., C
3327	- (xthetarem * xthetapn * xphi * fgmepn) & ! Pn messy feeding remin., C
3328	- (xthetarem * xthetapd * xphi * fgmepd) & ! Pd messy feeding remin., C
3329	- (xthetarem * xthetazmi * xphi * fgmezmi) & ! Zi messy feeding remin., C
3330	- (xthetarem * xphi * fgmedc) & ! D messy feeding remin., C
3331	- (xthetarem * fmiresp) & ! microzoo respiration, C
3332	- (xthetarem * fmeresp) & ! mesozoo respiration, C
3333	- (xthetarem * fddc) & ! slow detritus remin., C
3334	- (xthetarem * freminc) & ! fast detritus remin., C
3335	- (xthetarem * xthetapn * fdpn2) & ! Pn losses, C
3336	- (xthetarem * xthetapd * fdpd2) & ! Pd losses, C
3337	- (xthetarem * xthetazmi * fdzmi2) & ! Zmi losses, C
3338	- (xthetarem * xthetazme * fdzme2) ! Zme losses, C
3339	!!
3340	!! benthic remineralisation
3341	if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
3342	fo2_ccons = fo2_ccons - (xthetarem * f_benout_c(ji,jj) / fthk)
3343	endif
3344	fo2_cons = fo2_ncons + fo2_ccons
3345	!!
3346	!! is this a suboxic zone?
3347	if (zoxy.lt.xo2min) then ! deficient O2; production fluxes only
3348	btra(jpoxy) = b0 * ( &
3349	fo2_prod )
3350	foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod )
3351	foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fthk * fo2_cons )
3352	else ! sufficient O2; production + consumption fluxes
3353	btra(jpoxy) = b0 * ( &
3354	fo2_prod + fo2_cons )
3355	foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod )
3356	foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fthk * fo2_cons )
3357	endif
3358	!!
3359	!! air-sea fluxes (if this is the surface box)
3360	if (jk.eq.1) then
3361	!!
3362	!! CO2 flux
3363	btra(jpdic) = btra(jpdic) + (b0 * f_co2flux)
3364	!!
3365	!! O2 flux (mol/m3/s -> mmol/m3/d)
3366	btra(jpoxy) = btra(jpoxy) + (b0 * f_o2flux)
3367	endif
3368	# endif
3369
3370	!!----------------------------------------------------------------------
3371	!! Integrate calculated fluxes for mass balance
3372	!!----------------------------------------------------------------------
3373	!!
3374	!! === nitrogen ===
3375	fflx_n(ji,jj) = fflx_n(ji,jj) + &
3376	fthk * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) )
3377	!! === silicon ===
3378	fflx_si(ji,jj) = fflx_si(ji,jj) + &
3379	fthk * ( btra(jppds) + btra(jpsil) )
3380	!! === iron ===
3381	fflx_fe(ji,jj) = fflx_fe(ji,jj) + &
3382	fthk * ( ( xrfn * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet)) ) + btra(jpfer) )
3383	# if defined key_roam
3384	!! === carbon ===
3385	fflx_c(ji,jj) = fflx_c(ji,jj) + &
3386	fthk * ( (xthetapn * btra(jpphn)) + (xthetapd * btra(jpphd)) + &
3387	(xthetazmi * btra(jpzmi)) + (xthetazme * btra(jpzme)) + btra(jpdtc) + btra(jpdic) )
3388	!! === alkalinity ===
3389	fflx_a(ji,jj) = fflx_a(ji,jj) + &
3390	fthk * ( btra(jpalk) )
3391	!! === oxygen ===
3392	fflx_o2(ji,jj) = fflx_o2(ji,jj) + &
3393	fthk * ( btra(jpoxy) )
3394	# endif
3395
3396	!!----------------------------------------------------------------------
3397	!! Apply calculated tracer fluxes
3398	!!----------------------------------------------------------------------
3399	!!
3400	!! units: [unit of tracer] per second (fluxes are calculated above per day)
3401	!!
3402	ibio_switch = 1
3403	# if defined key_gulf_finland
3404	!! AXY (17/05/13): fudge in a Gulf of Finland correction; uses longitude-
3405	!! latitude range to establish if this is a Gulf of Finland
3406	!! grid cell; if so, then BGC fluxes are ignored (though
3407	!! still calculated); for reference, this is meant to be a
3408	!! temporary fix to see if all of my problems can be done
3409	!! away with if I switch off BGC fluxes in the Gulf of
3410	!! Finland, which currently appears the source of trouble
3411	if ( flonx.gt.24.7 .and. flonx.lt.27.8 .and. &
3412	& flatx.gt.59.2 .and. flatx.lt.60.2 ) then
3413	ibio_switch = 0
3414	endif
3415	# endif
3416	if (ibio_switch.eq.1) then
3417	tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(jpchn) / 86400.)
3418	tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(jpchd) / 86400.)
3419	tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(jpphn) / 86400.)
3420	tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(jpphd) / 86400.)
3421	tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(jppds) / 86400.)
3422	tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(jpzmi) / 86400.)
3423	tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(jpzme) / 86400.)
3424	tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(jpdet) / 86400.)
3425	tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(jpdin) / 86400.)
3426	tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(jpsil) / 86400.)
3427	tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(jpfer) / 86400.)
3428	# if defined key_roam
3429	tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(jpdtc) / 86400.)
3430	tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(jpdic) / 86400.)
3431	tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(jpalk) / 86400.)
3432	tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(jpoxy) / 86400.)
3433	# endif
3434	endif
3435
3436	!! AXY (18/11/16): CMIP6 diagnostics
3437	IF( med_diag%FBDDTALK%dgsave ) THEN
3438	fbddtalk(ji,jj) = fbddtalk(ji,jj) + (btra(jpalk) * fthk)
3439	ENDIF
3440	IF( med_diag%FBDDTDIC%dgsave ) THEN
3441	fbddtdic(ji,jj) = fbddtdic(ji,jj) + (btra(jpdic) * fthk)
3442	ENDIF
3443	IF( med_diag%FBDDTDIFE%dgsave ) THEN
3444	fbddtdife(ji,jj) = fbddtdife(ji,jj) + (btra(jpfer) * fthk)
3445	ENDIF
3446	IF( med_diag%FBDDTDIN%dgsave ) THEN
3447	fbddtdin(ji,jj) = fbddtdin(ji,jj) + (btra(jpdin) * fthk)
3448	ENDIF
3449	IF( med_diag%FBDDTDISI%dgsave ) THEN
3450	fbddtdisi(ji,jj) = fbddtdisi(ji,jj) + (btra(jpsil) * fthk)
3451	ENDIF
3452	!!
3453	IF( med_diag%BDDTALK3%dgsave ) THEN
3454	bddtalk3(ji,jj,jk) = btra(jpalk)
3455	ENDIF
3456	IF( med_diag%BDDTDIC3%dgsave ) THEN
3457	bddtdic3(ji,jj,jk) = btra(jpdic)
3458	ENDIF
3459	IF( med_diag%BDDTDIFE3%dgsave ) THEN
3460	bddtdife3(ji,jj,jk) = btra(jpfer)
3461	ENDIF
3462	IF( med_diag%BDDTDIN3%dgsave ) THEN
3463	bddtdin3(ji,jj,jk) = btra(jpdin)
3464	ENDIF
3465	IF( med_diag%BDDTDISI3%dgsave ) THEN
3466	bddtdisi3(ji,jj,jk) = btra(jpsil)
3467	ENDIF
3468
3469	# if defined key_debug_medusa
3470	IF ( lwp ) write (numout,*) '------'
3471	IF ( lwp ) write (numout,*) 'trc_bio_medusa: end all calculations'
3472	IF ( lwp ) write (numout,*) 'trc_bio_medusa: now outputs'
3473	CALL flush(numout)
3474	# endif
3475
3476	# if defined key_axy_nancheck
3477	!!----------------------------------------------------------------------
3478	!! Check calculated tracer fluxes
3479	!!----------------------------------------------------------------------
3480	!!
3481	DO jn = 1,jptra
3482	fq0 = btra(jn)
3483	!! AXY (30/01/14): "isnan" problem on HECTOR
3484	!! if (fq0 /= fq0 ) then
3485	if ( ieee_is_nan( fq0 ) ) then
3486	!! there's a NaN here
3487	if (lwp) write(numout,*) 'NAN detected in btra(', ji, ',', &
3488	& jj, ',', jk, ',', jn, ') at time', kt
3489	CALL ctl_stop( 'trcbio_medusa, NAN in btra field' )
3490	endif
3491	ENDDO
3492	DO jn = 1,jptra
3493	fq0 = tra(ji,jj,jk,jn)
3494	!! AXY (30/01/14): "isnan" problem on HECTOR
3495	!! if (fq0 /= fq0 ) then
3496	if ( ieee_is_nan( fq0 ) ) then
3497	!! there's a NaN here
3498	if (lwp) write(numout,*) 'NAN detected in tra(', ji, ',', &
3499	& jj, ',', jk, ',', jn, ') at time', kt
3500	CALL ctl_stop( 'trcbio_medusa, NAN in tra field' )
3501	endif
3502	ENDDO
3503	CALL flush(numout)
3504	# endif
3505
3506	!!----------------------------------------------------------------------
3507	!! Check model conservation
3508	!! these terms merely sum up the tendency terms of the relevant
3509	!! state variables, which should sum to zero; the iron cycle is
3510	!! complicated by fluxes that add (aeolian deposition and seafloor
3511	!! remineralisation) and remove (scavenging) dissolved iron from
3512	!! the model (i.e. the sum of iron fluxes is unlikely to be zero)
3513	!!----------------------------------------------------------------------
3514	!!
3515	!! fnit0 = btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) ! + ftempn
3516	!! fsil0 = btra(jppds) + btra(jpsil) ! + ftempsi
3517	!! ffer0 = (xrfn * fnit0) + btra(jpfer)
3518	# if defined key_roam
3519	!! fcar0 = 0.
3520	!! falk0 = 0.
3521	!! foxy0 = 0.
3522	# endif
3523	!!
3524	!! if (kt/240*240.eq.kt) then
3525	!! if (ji.eq.2.and.jj.eq.2.and.jk.eq.1) then
3526	!! IF ( lwp ) write (,) '*****!MEDUSA Conservation!*****',kt
3527	# if defined key_roam
3528	!! IF ( lwp ) write (,) fnit0,fsil0,ffer0,fcar0,falk0,foxy0
3529	# else
3530	!! IF ( lwp ) write (,) fnit0,fsil0,ffer0
3531	# endif
3532	!! endif
3533	!! endif
3534
3535	IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN
3536	!!----------------------------------------------------------------------
3537	!! Add in XML diagnostics stuff
3538	!!----------------------------------------------------------------------
3539	!!
3540	!! ** 2D diagnostics
3541	# if defined key_debug_medusa
3542	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk loop'
3543	CALL flush(numout)
3544	# endif
3545	IF ( med_diag%PRN%dgsave ) THEN
3546	fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn * zphn * fthk)
3547	ENDIF
3548	IF ( med_diag%MPN%dgsave ) THEN
3549	fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn * fthk)
3550	ENDIF
3551	IF ( med_diag%PRD%dgsave ) THEN
3552	fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd * zphd * fthk)
3553	ENDIF
3554	IF( med_diag%MPD%dgsave ) THEN
3555	fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd * fthk)
3556	ENDIF
3557	! IF( med_diag%DSED%dgsave ) THEN
3558	! CALL iom_put( "DSED" , ftot_n )
3559	! ENDIF
3560	IF( med_diag%OPAL%dgsave ) THEN
3561	fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds * zpds * fthk)
3562	ENDIF
3563	IF( med_diag%OPALDISS%dgsave ) THEN
3564	fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss * fthk)
3565	ENDIF
3566	IF( med_diag%GMIPn%dgsave ) THEN
3567	fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn * fthk)
3568	ENDIF
3569	IF( med_diag%GMID%dgsave ) THEN
3570	fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid * fthk)
3571	ENDIF
3572	IF( med_diag%MZMI%dgsave ) THEN
3573	fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi * fthk)
3574	ENDIF
3575	IF( med_diag%GMEPN%dgsave ) THEN
3576	fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn * fthk)
3577	ENDIF
3578	IF( med_diag%GMEPD%dgsave ) THEN
3579	fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd * fthk)
3580	ENDIF
3581	IF( med_diag%GMEZMI%dgsave ) THEN
3582	fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi * fthk)
3583	ENDIF
3584	IF( med_diag%GMED%dgsave ) THEN
3585	fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed * fthk)
3586	ENDIF
3587	IF( med_diag%MZME%dgsave ) THEN
3588	fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme * fthk)
3589	ENDIF
3590	! IF( med_diag%DEXP%dgsave ) THEN
3591	! CALL iom_put( "DEXP" , ftot_n )
3592	! ENDIF
3593	IF( med_diag%DETN%dgsave ) THEN
3594	fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown * fthk)
3595	ENDIF
3596	IF( med_diag%MDET%dgsave ) THEN
3597	fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd * fthk)
3598	ENDIF
3599	IF( med_diag%AEOLIAN%dgsave ) THEN
3600	ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop * fthk)
3601	ENDIF
3602	IF( med_diag%BENTHIC%dgsave ) THEN
3603	ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot * fthk)
3604	ENDIF
3605	IF( med_diag%SCAVENGE%dgsave ) THEN
3606	ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav * fthk)
3607	ENDIF
3608	IF( med_diag%PN_JLIM%dgsave ) THEN
3609	! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln * zphn * fthk)
3610	fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn * zphn * fthk)
3611	ENDIF
3612	IF( med_diag%PN_NLIM%dgsave ) THEN
3613	fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln * zphn * fthk)
3614	ENDIF
3615	IF( med_diag%PN_FELIM%dgsave ) THEN
3616	ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln * zphn * fthk)
3617	ENDIF
3618	IF( med_diag%PD_JLIM%dgsave ) THEN
3619	! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld * zphd * fthk)
3620	fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd * zphd * fthk)
3621	ENDIF
3622	IF( med_diag%PD_NLIM%dgsave ) THEN
3623	fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld * zphd * fthk)
3624	ENDIF
3625	IF( med_diag%PD_FELIM%dgsave ) THEN
3626	ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld * zphd * fthk)
3627	ENDIF
3628	IF( med_diag%PD_SILIM%dgsave ) THEN
3629	fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2 * zphd * fthk)
3630	ENDIF
3631	IF( med_diag%PDSILIM2%dgsave ) THEN
3632	fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld * zphd * fthk)
3633	ENDIF
3634	!!
3635	IF( med_diag%TOTREG_N%dgsave ) THEN
3636	fregen2d(ji,jj) = fregen2d(ji,jj) + fregen
3637	ENDIF
3638	IF( med_diag%TOTRG_SI%dgsave ) THEN
3639	fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi
3640	ENDIF
3641	!!
3642	IF( med_diag%FASTN%dgsave ) THEN
3643	ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn * fthk)
3644	ENDIF
3645	IF( med_diag%FASTSI%dgsave ) THEN
3646	ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi * fthk)
3647	ENDIF
3648	IF( med_diag%FASTFE%dgsave ) THEN
3649	ftempfe2d(ji,jj) =ftempfe2d(ji,jj) + (ftempfe * fthk)
3650	ENDIF
3651	IF( med_diag%FASTC%dgsave ) THEN
3652	ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc * fthk)
3653	ENDIF
3654	IF( med_diag%FASTCA%dgsave ) THEN
3655	ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca * fthk)
3656	ENDIF
3657	!!
3658	IF( med_diag%REMINN%dgsave ) THEN
3659	freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn * fthk)
3660	ENDIF
3661	IF( med_diag%REMINSI%dgsave ) THEN
3662	freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi * fthk)
3663	ENDIF
3664	IF( med_diag%REMINFE%dgsave ) THEN
3665	freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe * fthk)
3666	ENDIF
3667	IF( med_diag%REMINC%dgsave ) THEN
3668	freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc * fthk)
3669	ENDIF
3670	IF( med_diag%REMINCA%dgsave ) THEN
3671	freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca * fthk)
3672	ENDIF
3673	!!
3674	# if defined key_roam
3675	!!
3676	!! AXY (09/11/16): CMIP6 diagnostics
3677	IF( med_diag%FD_NIT3%dgsave ) THEN
3678	fd_nit3(ji,jj,jk) = ffastn(ji,jj)
3679	ENDIF
3680	IF( med_diag%FD_SIL3%dgsave ) THEN
3681	fd_sil3(ji,jj,jk) = ffastsi(ji,jj)
3682	ENDIF
3683	IF( med_diag%FD_CAR3%dgsave ) THEN
3684	fd_car3(ji,jj,jk) = ffastc(ji,jj)
3685	ENDIF
3686	IF( med_diag%FD_CAL3%dgsave ) THEN
3687	fd_cal3(ji,jj,jk) = ffastca(ji,jj)
3688	ENDIF
3689	!!
3690	IF (jk.eq.i0100) THEN
3691	IF( med_diag%RR_0100%dgsave ) THEN
3692	ffastca2d(ji,jj) = &
3693	ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
3694	ENDIF
3695	ELSE IF (jk.eq.i0500) THEN
3696	IF( med_diag%RR_0500%dgsave ) THEN
3697	ffastca2d(ji,jj) = &
3698	ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
3699	ENDIF
3700	ELSE IF (jk.eq.i1000) THEN
3701	IF( med_diag%RR_1000%dgsave ) THEN
3702	ffastca2d(ji,jj) = &
3703	ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
3704	ENDIF
3705	ELSE IF (jk.eq.jmbathy) THEN
3706	IF( med_diag%IBEN_N%dgsave ) THEN
3707	iben_n2d(ji,jj) = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj)
3708	ENDIF
3709	IF( med_diag%IBEN_FE%dgsave ) THEN
3710	iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
3711	ENDIF
3712	IF( med_diag%IBEN_C%dgsave ) THEN
3713	iben_c2d(ji,jj) = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj)
3714	ENDIF
3715	IF( med_diag%IBEN_SI%dgsave ) THEN
3716	iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
3717	ENDIF
3718	IF( med_diag%IBEN_CA%dgsave ) THEN
3719	iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
3720	ENDIF
3721	IF( med_diag%OBEN_N%dgsave ) THEN
3722	oben_n2d(ji,jj) = f_benout_n(ji,jj)
3723	ENDIF
3724	IF( med_diag%OBEN_FE%dgsave ) THEN
3725	oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
3726	ENDIF
3727	IF( med_diag%OBEN_C%dgsave ) THEN
3728	oben_c2d(ji,jj) = f_benout_c(ji,jj)
3729	ENDIF
3730	IF( med_diag%OBEN_SI%dgsave ) THEN
3731	oben_si2d(ji,jj) = f_benout_si(ji,jj)
3732	ENDIF
3733	IF( med_diag%OBEN_CA%dgsave ) THEN
3734	oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
3735	ENDIF
3736	IF( med_diag%SFR_OCAL%dgsave ) THEN
3737	sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
3738	ENDIF
3739	IF( med_diag%SFR_OARG%dgsave ) THEN
3740	sfr_oarg2d(ji,jj) = f3_omarg(ji,jj,jk)
3741	ENDIF
3742	IF( med_diag%LYSO_CA%dgsave ) THEN
3743	lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
3744	ENDIF
3745	ENDIF
3746	!! end bathy-1 diags
3747	!!
3748	IF( med_diag%RIV_N%dgsave ) THEN
3749	rivn2d(ji,jj) = rivn2d(ji,jj) + (f_riv_loc_n * fthk)
3750	ENDIF
3751	IF( med_diag%RIV_SI%dgsave ) THEN
3752	rivsi2d(ji,jj) = rivsi2d(ji,jj) + (f_riv_loc_si * fthk)
3753	ENDIF
3754	IF( med_diag%RIV_C%dgsave ) THEN
3755	rivc2d(ji,jj) = rivc2d(ji,jj) + (f_riv_loc_c * fthk)
3756	ENDIF
3757	IF( med_diag%RIV_ALK%dgsave ) THEN
3758	rivalk2d(ji,jj) = rivalk2d(ji,jj) + (f_riv_loc_alk * fthk)
3759	ENDIF
3760	IF( med_diag%DETC%dgsave ) THEN
3761	fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc * fthk)
3762	ENDIF
3763	!!
3764	!!
3765	!!
3766	IF( med_diag%PN_LLOSS%dgsave ) THEN
3767	fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2 * fthk)
3768	ENDIF
3769	IF( med_diag%PD_LLOSS%dgsave ) THEN
3770	fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2 * fthk)
3771	ENDIF
3772	IF( med_diag%ZI_LLOSS%dgsave ) THEN
3773	fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2 * fthk)
3774	ENDIF
3775	IF( med_diag%ZE_LLOSS%dgsave ) THEN
3776	fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2 * fthk)
3777	ENDIF
3778	IF( med_diag%ZI_MES_N%dgsave ) THEN
3779	zimesn2d(ji,jj) = zimesn2d(ji,jj) + &
3780	(xphi * (fgmipn + fgmid) * fthk)
3781	ENDIF
3782	IF( med_diag%ZI_MES_D%dgsave ) THEN
3783	zimesd2d(ji,jj) = zimesd2d(ji,jj) + &
3784	((1. - xbetan) * finmi * fthk)
3785	ENDIF
3786	IF( med_diag%ZI_MES_C%dgsave ) THEN
3787	zimesc2d(ji,jj) = zimesc2d(ji,jj) + &
3788	(xphi * ((xthetapn * fgmipn) + fgmidc) * fthk)
3789	ENDIF
3790	IF( med_diag%ZI_MESDC%dgsave ) THEN
3791	zimesdc2d(ji,jj) = zimesdc2d(ji,jj) + &
3792	((1. - xbetac) * ficmi * fthk)
3793	ENDIF
3794	IF( med_diag%ZI_EXCR%dgsave ) THEN
3795	ziexcr2d(ji,jj) = ziexcr2d(ji,jj) + (fmiexcr * fthk)
3796	ENDIF
3797	IF( med_diag%ZI_RESP%dgsave ) THEN
3798	ziresp2d(ji,jj) = ziresp2d(ji,jj) + (fmiresp * fthk)
3799	ENDIF
3800	IF( med_diag%ZI_GROW%dgsave ) THEN
3801	zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow * fthk)
3802	ENDIF
3803	IF( med_diag%ZE_MES_N%dgsave ) THEN
3804	zemesn2d(ji,jj) = zemesn2d(ji,jj) + &
3805	(xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk)
3806	ENDIF
3807	IF( med_diag%ZE_MES_D%dgsave ) THEN
3808	zemesd2d(ji,jj) = zemesd2d(ji,jj) + &
3809	((1. - xbetan) * finme * fthk)
3810	ENDIF
3811	IF( med_diag%ZE_MES_C%dgsave ) THEN
3812	zemesc2d(ji,jj) = zemesc2d(ji,jj) + &
3813	(xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + &
3814	(xthetazmi * fgmezmi) + fgmedc) * fthk)
3815	ENDIF
3816	IF( med_diag%ZE_MESDC%dgsave ) THEN
3817	zemesdc2d(ji,jj) = zemesdc2d(ji,jj) + &
3818	((1. - xbetac) * ficme * fthk)
3819	ENDIF
3820	IF( med_diag%ZE_EXCR%dgsave ) THEN
3821	zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr * fthk)
3822	ENDIF
3823	IF( med_diag%ZE_RESP%dgsave ) THEN
3824	zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp * fthk)
3825	ENDIF
3826	IF( med_diag%ZE_GROW%dgsave ) THEN
3827	zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow * fthk)
3828	ENDIF
3829	IF( med_diag%MDETC%dgsave ) THEN
3830	mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc * fthk)
3831	ENDIF
3832	IF( med_diag%GMIDC%dgsave ) THEN
3833	gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc * fthk)
3834	ENDIF
3835	IF( med_diag%GMEDC%dgsave ) THEN
3836	gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc * fthk)
3837	ENDIF
3838	!!
3839	# endif
3840	!!
3841	!! ** 3D diagnostics
3842	IF( med_diag%TPP3%dgsave ) THEN
3843	tpp3d(ji,jj,jk) = (fprn * zphn) + (fprd * zphd)
3844	!CALL iom_put( "TPP3" , tpp3d )
3845	ENDIF
3846	IF( med_diag%TPPD3%dgsave ) THEN
3847	tppd3(ji,jj,jk) = (fprd * zphd)
3848	ENDIF
3849
3850	IF( med_diag%REMIN3N%dgsave ) THEN
3851	remin3dn(ji,jj,jk) = fregen + (freminn * fthk) !! remineralisation
3852	!CALL iom_put( "REMIN3N" , remin3dn )
3853	ENDIF
3854	!! IF( med_diag%PH3%dgsave ) THEN
3855	!! CALL iom_put( "PH3" , f3_pH )
3856	!! ENDIF
3857	!! IF( med_diag%OM_CAL3%dgsave ) THEN
3858	!! CALL iom_put( "OM_CAL3" , f3_omcal )
3859	!! ENDIF
3860	!!
3861	!! AXY (09/11/16): CMIP6 diagnostics
3862	IF ( med_diag%DCALC3%dgsave ) THEN
3863	dcalc3(ji,jj,jk) = freminca
3864	ENDIF
3865	IF ( med_diag%FEDISS3%dgsave ) THEN
3866	fediss3(ji,jj,jk) = ffetop
3867	ENDIF
3868	IF ( med_diag%FESCAV3%dgsave ) THEN
3869	fescav3(ji,jj,jk) = ffescav
3870	ENDIF
3871	IF ( med_diag%MIGRAZP3%dgsave ) THEN
3872	migrazp3(ji,jj,jk) = fgmipn * xthetapn
3873	ENDIF
3874	IF ( med_diag%MIGRAZD3%dgsave ) THEN
3875	migrazd3(ji,jj,jk) = fgmidc
3876	ENDIF
3877	IF ( med_diag%MEGRAZP3%dgsave ) THEN
3878	megrazp3(ji,jj,jk) = (fgmepn * xthetapn) + (fgmepd * xthetapd)
3879	ENDIF
3880	IF ( med_diag%MEGRAZD3%dgsave ) THEN
3881	megrazd3(ji,jj,jk) = fgmedc
3882	ENDIF
3883	IF ( med_diag%MEGRAZZ3%dgsave ) THEN
3884	megrazz3(ji,jj,jk) = (fgmezmi * xthetazmi)
3885	ENDIF
3886	IF ( med_diag%PBSI3%dgsave ) THEN
3887	pbsi3(ji,jj,jk) = (fprds * zpds)
3888	ENDIF
3889	IF ( med_diag%PCAL3%dgsave ) THEN
3890	pcal3(ji,jj,jk) = ftempca
3891	ENDIF
3892	IF ( med_diag%REMOC3%dgsave ) THEN
3893	remoc3(ji,jj,jk) = freminc
3894	ENDIF
3895	IF ( med_diag%PNLIMJ3%dgsave ) THEN
3896	! pnlimj3(ji,jj,jk) = fjln
3897	pnlimj3(ji,jj,jk) = fjlim_pn
3898	ENDIF
3899	IF ( med_diag%PNLIMN3%dgsave ) THEN
3900	pnlimn3(ji,jj,jk) = fnln
3901	ENDIF
3902	IF ( med_diag%PNLIMFE3%dgsave ) THEN
3903	pnlimfe3(ji,jj,jk) = ffln
3904	ENDIF
3905	IF ( med_diag%PDLIMJ3%dgsave ) THEN
3906	! pdlimj3(ji,jj,jk) = fjld
3907	pdlimj3(ji,jj,jk) = fjlim_pd
3908	ENDIF
3909	IF ( med_diag%PDLIMN3%dgsave ) THEN
3910	pdlimn3(ji,jj,jk) = fnld
3911	ENDIF
3912	IF ( med_diag%PDLIMFE3%dgsave ) THEN
3913	pdlimfe3(ji,jj,jk) = ffld
3914	ENDIF
3915	IF ( med_diag%PDLIMSI3%dgsave ) THEN
3916	pdlimsi3(ji,jj,jk) = fsld2
3917	ENDIF
3918	ENDIF ! end of ln_diatrc option
3919	!! CLOSE wet point IF..THEN loop
3920	endif
3921	!! CLOSE horizontal loops
3922	ENDDO
3923	ENDDO
3924	!!
3925	IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN
3926	!! first - 2D diag implemented
3927	!! on every K level
3928	!!-----------------------------------------
3929	!! --
3930	!!second - 2d specific k level diags
3931	!!
3932	!!-----------------------------------------
3933	IF (jk.eq.1) THEN
3934	# if defined key_debug_medusa
3935	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag jk = 1'
3936	CALL flush(numout)
3937	# endif
3938	IF( med_diag%MED_QSR%dgsave ) THEN
3939	CALL iom_put( "MED_QSR" , qsr ) !
3940	ENDIF
3941	IF( med_diag%MED_XPAR%dgsave ) THEN
3942	CALL iom_put( "MED_XPAR" , xpar(:,:,jk) ) !
3943	ENDIF
3944	IF( med_diag%OCAL_CCD%dgsave ) THEN
3945	CALL iom_put( "OCAL_CCD" , ocal_ccd ) !
3946	ENDIF
3947	IF( med_diag%FE_0000%dgsave ) THEN
3948	CALL iom_put( "FE_0000" , xFree ) !
3949	ENDIF
3950	IF( med_diag%MED_XZE%dgsave ) THEN
3951	CALL iom_put( "MED_XZE" , xze ) !
3952	ENDIF
3953	# if defined key_roam
3954	IF( med_diag%WIND%dgsave ) THEN
3955	CALL iom_put( "WIND" , wndm )
3956	ENDIF
3957	IF( med_diag%ATM_PCO2%dgsave ) THEN
3958	CALL iom_put( "ATM_PCO2" , f_pco2a2d )
3959	CALL wrk_dealloc( jpi, jpj, f_pco2a2d )
3960	ENDIF
3961	IF( med_diag%OCN_PH%dgsave ) THEN
3962	zw2d(:,:) = f3_pH(:,:,jk)
3963	CALL iom_put( "OCN_PH" , zw2d )
3964	ENDIF
3965	IF( med_diag%OCN_PCO2%dgsave ) THEN
3966	CALL iom_put( "OCN_PCO2" , f_pco2w2d )
3967	CALL wrk_dealloc( jpi, jpj, f_pco2w2d )
3968	ENDIF
3969	IF( med_diag%OCNH2CO3%dgsave ) THEN
3970	zw2d(:,:) = f3_h2co3(:,:,jk)
3971	CALL iom_put( "OCNH2CO3" , zw2d )
3972	ENDIF
3973	IF( med_diag%OCN_HCO3%dgsave ) THEN
3974	zw2d(:,:) = f3_hco3(:,:,jk)
3975	CALL iom_put( "OCN_HCO3" , zw2d )
3976	ENDIF
3977	IF( med_diag%OCN_CO3%dgsave ) THEN
3978	zw2d(:,:) = f3_co3(:,:,jk)
3979	CALL iom_put( "OCN_CO3" , zw2d )
3980	ENDIF
3981	IF( med_diag%CO2FLUX%dgsave ) THEN
3982	CALL iom_put( "CO2FLUX" , f_co2flux2d )
3983	CALL wrk_dealloc( jpi, jpj, f_co2flux2d )
3984	ENDIF
3985	!!
3986	!! AXY (10/11/16): repeat CO2 flux diagnostic in UKMO/CMIP6 units; this
3987	!! both outputs the CO2 flux in specified units and
3988	!! sends the resulting field to the coupler
3989	!! JPALM (17/11/16): put CO2 flux (fgco2) alloc/unalloc/pass to zn
3990	!! out of diag list request
3991	CALL lbc_lnk( fgco2(:,:),'T',1. )
3992	IF( med_diag%FGCO2%dgsave ) THEN
3993	CALL iom_put( "FGCO2" , fgco2 )
3994	ENDIF
3995	!! JPALM (17/11/16): should mv this fgco2 part
3996	!! out of lk_iomput loop
3997	zb_co2_flx = zn_co2_flx
3998	zn_co2_flx = fgco2
3999	IF ( lk_oasis ) THEN
4000	CO2Flux_out_cpl = zn_co2_flx
4001	ENDIF
4002	CALL wrk_dealloc( jpi, jpj, fgco2 )
4003	!! ---
4004	IF( med_diag%OM_CAL%dgsave ) THEN
4005	CALL iom_put( "OM_CAL" , f_omcal )
4006	ENDIF
4007	IF( med_diag%OM_ARG%dgsave ) THEN
4008	CALL iom_put( "OM_ARG" , f_omarg )
4009	ENDIF
4010	IF( med_diag%TCO2%dgsave ) THEN
4011	CALL iom_put( "TCO2" , f_TDIC2d )
4012	CALL wrk_dealloc( jpi, jpj, f_TDIC2d )
4013	ENDIF
4014	IF( med_diag%TALK%dgsave ) THEN
4015	CALL iom_put( "TALK" , f_TALK2d )
4016	CALL wrk_dealloc( jpi, jpj, f_TALK2d )
4017	ENDIF
4018	IF( med_diag%KW660%dgsave ) THEN
4019	CALL iom_put( "KW660" , f_kw6602d )
4020	CALL wrk_dealloc( jpi, jpj, f_kw6602d )
4021	ENDIF
4022	IF( med_diag%ATM_PP0%dgsave ) THEN
4023	CALL iom_put( "ATM_PP0" , f_pp02d )
4024	CALL wrk_dealloc( jpi, jpj, f_pp02d )
4025	ENDIF
4026	IF( med_diag%O2FLUX%dgsave ) THEN
4027	CALL iom_put( "O2FLUX" , f_o2flux2d )
4028	CALL wrk_dealloc( jpi, jpj, f_o2flux2d )
4029	ENDIF
4030	IF( med_diag%O2SAT%dgsave ) THEN
4031	CALL iom_put( "O2SAT" , f_o2sat2d )
4032	CALL wrk_dealloc( jpi, jpj, f_o2sat2d )
4033	ENDIF
4034	IF( med_diag%CAL_CCD%dgsave ) THEN
4035	CALL iom_put( "CAL_CCD" , f2_ccd_cal )
4036	ENDIF
4037	IF( med_diag%ARG_CCD%dgsave ) THEN
4038	CALL iom_put( "ARG_CCD" , f2_ccd_arg )
4039	ENDIF
4040	IF (jdms .eq. 1) THEN
4041	IF( med_diag%DMS_SURF%dgsave ) THEN
4042	CALL lbc_lnk(dms_surf2d(:,:),'T',1. )
4043	CALL iom_put( "DMS_SURF" , dms_surf2d )
4044	zb_dms_srf = zn_dms_srf
4045	zn_dms_srf = dms_surf2d
4046	IF ( lk_oasis ) THEN
4047	DMS_out_cpl = zn_dms_srf
4048	ENDIF
4049	CALL wrk_dealloc( jpi, jpj, dms_surf2d )
4050	ENDIF
4051	IF( med_diag%DMS_ANDR%dgsave ) THEN
4052	CALL iom_put( "DMS_ANDR" , dms_andr2d )
4053	CALL wrk_dealloc( jpi, jpj, dms_andr2d )
4054	ENDIF
4055	IF( med_diag%DMS_SIMO%dgsave ) THEN
4056	CALL iom_put( "DMS_SIMO" , dms_simo2d )
4057	CALL wrk_dealloc( jpi, jpj, dms_simo2d )
4058	ENDIF
4059	IF( med_diag%DMS_ARAN%dgsave ) THEN
4060	CALL iom_put( "DMS_ARAN" , dms_aran2d )
4061	CALL wrk_dealloc( jpi, jpj, dms_aran2d )
4062	ENDIF
4063	IF( med_diag%DMS_HALL%dgsave ) THEN
4064	CALL iom_put( "DMS_HALL" , dms_hall2d )
4065	CALL wrk_dealloc( jpi, jpj, dms_hall2d )
4066	ENDIF
4067	ENDIF
4068	!! AXY (24/11/16): extra MOCSY diagnostics
4069	IF( med_diag%ATM_XCO2%dgsave ) THEN
4070	CALL iom_put( "ATM_XCO2" , f_xco2a_2d )
4071	CALL wrk_dealloc( jpi, jpj, f_xco2a_2d )
4072	ENDIF
4073	IF( med_diag%OCN_FCO2%dgsave ) THEN
4074	CALL iom_put( "OCN_FCO2" , f_fco2w_2d )
4075	CALL wrk_dealloc( jpi, jpj, f_fco2w_2d )
4076	ENDIF
4077	IF( med_diag%ATM_FCO2%dgsave ) THEN
4078	CALL iom_put( "ATM_FCO2" , f_fco2a_2d )
4079	CALL wrk_dealloc( jpi, jpj, f_fco2a_2d )
4080	ENDIF
4081	IF( med_diag%OCN_RHOSW%dgsave ) THEN
4082	CALL iom_put( "OCN_RHOSW" , f_ocnrhosw_2d )
4083	CALL wrk_dealloc( jpi, jpj, f_ocnrhosw_2d )
4084	ENDIF
4085	IF( med_diag%OCN_SCHCO2%dgsave ) THEN
4086	CALL iom_put( "OCN_SCHCO2" , f_ocnschco2_2d )
4087	CALL wrk_dealloc( jpi, jpj, f_ocnschco2_2d )
4088	ENDIF
4089	IF( med_diag%OCN_KWCO2%dgsave ) THEN
4090	CALL iom_put( "OCN_KWCO2" , f_ocnkwco2_2d )
4091	CALL wrk_dealloc( jpi, jpj, f_ocnkwco2_2d )
4092	ENDIF
4093	IF( med_diag%OCN_K0%dgsave ) THEN
4094	CALL iom_put( "OCN_K0" , f_ocnk0_2d )
4095	CALL wrk_dealloc( jpi, jpj, f_ocnk0_2d )
4096	ENDIF
4097	IF( med_diag%CO2STARAIR%dgsave ) THEN
4098	CALL iom_put( "CO2STARAIR" , f_co2starair_2d )
4099	CALL wrk_dealloc( jpi, jpj, f_co2starair_2d )
4100	ENDIF
4101	IF( med_diag%OCN_DPCO2%dgsave ) THEN
4102	CALL iom_put( "OCN_DPCO2" , f_ocndpco2_2d )
4103	CALL wrk_dealloc( jpi, jpj, f_ocndpco2_2d )
4104	ENDIF
4105	# endif
4106	ELSE IF (jk.eq.i0100) THEN
4107	# if defined key_debug_medusa
4108	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag jk = 100'
4109	CALL flush(numout)
4110	# endif
4111	IF( med_diag%SDT__100%dgsave ) THEN
4112	zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
4113	CALL iom_put( "SDT__100" , zw2d )
4114	ENDIF
4115	IF( med_diag%REG__100%dgsave ) THEN
4116	CALL iom_put( "REG__100" , fregen2d )
4117	ENDIF
4118	IF( med_diag%FDT__100%dgsave ) THEN
4119	CALL iom_put( "FDT__100" , ffastn )
4120	ENDIF
4121	IF( med_diag%RG__100F%dgsave ) THEN
4122	CALL iom_put( "RG__100F" , fregenfast )
4123	ENDIF
4124	IF( med_diag%FDS__100%dgsave ) THEN
4125	CALL iom_put( "FDS__100" , ffastsi )
4126	ENDIF
4127	IF( med_diag%RGS_100F%dgsave ) THEN
4128	CALL iom_put( "RGS_100F" , fregenfastsi )
4129	ENDIF
4130	IF( med_diag%FE_0100%dgsave ) THEN
4131	CALL iom_put( "FE_0100" , xFree )
4132	ENDIF
4133	# if defined key_roam
4134	IF( med_diag%RR_0100%dgsave ) THEN
4135	CALL iom_put( "RR_0100" , ffastca2d )
4136	ENDIF
4137	IF( med_diag%SDC__100%dgsave ) THEN
4138	zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk)
4139	CALL iom_put( "SDC__100" , zw2d )
4140	ENDIF
4141	IF( med_diag%epC100%dgsave ) THEN
4142	zw2d(:,:) = (fslowcflux + ffastc) * tmask(:,:,jk)
4143	CALL iom_put( "epC100" , zw2d )
4144	ENDIF
4145	IF( med_diag%epCALC100%dgsave ) THEN
4146	CALL iom_put( "epCALC100" , ffastca )
4147	ENDIF
4148	IF( med_diag%epN100%dgsave ) THEN
4149	zw2d(:,:) = (fslownflux + ffastn) * tmask(:,:,jk)
4150	CALL iom_put( "epN100" , zw2d )
4151	ENDIF
4152	IF( med_diag%epSI100%dgsave ) THEN
4153	CALL iom_put( "epSI100" , ffastsi )
4154	ENDIF
4155	ELSE IF (jk.eq.i0150) THEN
4156	# if defined key_debug_medusa
4157	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag jk = 150'
4158	CALL flush(numout)
4159	# endif
4160	# endif
4161	ELSE IF (jk.eq.i0200) THEN
4162	# if defined key_debug_medusa
4163	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag jk = 200'
4164	CALL flush(numout)
4165	# endif
4166	IF( med_diag%SDT__200%dgsave ) THEN
4167	zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
4168	CALL iom_put( "SDT__200" , zw2d )
4169	ENDIF
4170	IF( med_diag%REG__200%dgsave ) THEN
4171	CALL iom_put( "REG__200" , fregen2d )
4172	ENDIF
4173	IF( med_diag%FDT__200%dgsave ) THEN
4174	CALL iom_put( "FDT__200" , ffastn )
4175	ENDIF
4176	IF( med_diag%RG__200F%dgsave ) THEN
4177	CALL iom_put( "RG__200F" , fregenfast )
4178	ENDIF
4179	IF( med_diag%FDS__200%dgsave ) THEN
4180	CALL iom_put( "FDS__200" , ffastsi )
4181	ENDIF
4182	IF( med_diag%RGS_200F%dgsave ) THEN
4183	CALL iom_put( "RGS_200F" , fregenfastsi )
4184	ENDIF
4185	IF( med_diag%FE_0200%dgsave ) THEN
4186	CALL iom_put( "FE_0200" , xFree )
4187	ENDIF
4188	# if defined key_roam
4189	IF( med_diag%SDC__200%dgsave ) THEN
4190	zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk)
4191	CALL iom_put( "SDC__200" , zw2d )
4192	ENDIF
4193	# endif
4194	ELSE IF (jk.eq.i0500) THEN
4195	# if defined key_debug_medusa
4196	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag jk = 500'
4197	CALL flush(numout)
4198	# endif
4199	IF( med_diag%SDT__500%dgsave ) THEN
4200	zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
4201	CALL iom_put( "SDT__500" , zw2d )
4202	ENDIF
4203	IF( med_diag%REG__500%dgsave ) THEN
4204	CALL iom_put( "REG__500" , fregen2d )
4205	ENDIF
4206	IF( med_diag%FDT__500%dgsave ) THEN
4207	CALL iom_put( "FDT__500" , ffastn )
4208	ENDIF
4209	IF( med_diag%RG__500F%dgsave ) THEN
4210	CALL iom_put( "RG__500F" , fregenfast )
4211	ENDIF
4212	IF( med_diag%FDS__500%dgsave ) THEN
4213	CALL iom_put( "FDS__500" , ffastsi )
4214	ENDIF
4215	IF( med_diag%RGS_500F%dgsave ) THEN
4216	CALL iom_put( "RGS_500F" , fregenfastsi )
4217	ENDIF
4218	IF( med_diag%FE_0500%dgsave ) THEN
4219	CALL iom_put( "FE_0500" , xFree )
4220	ENDIF
4221	# if defined key_roam
4222	IF( med_diag%RR_0500%dgsave ) THEN
4223	CALL iom_put( "RR_0500" , ffastca2d )
4224	ENDIF
4225	IF( med_diag%SDC__500%dgsave ) THEN
4226	zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk)
4227	CALL iom_put( "SDC__500" , zw2d )
4228	ENDIF
4229	# endif
4230	ELSE IF (jk.eq.i1000) THEN
4231	# if defined key_debug_medusa
4232	IF ( lwp ) write (numout,*) 'trc_bio_medusa: diag jk = 1000'
4233	CALL flush(numout)
4234	# endif
4235	IF( med_diag%SDT_1000%dgsave ) THEN
4236	zw2d(:,:) = fslownflux(:,:) * tmask(:,:,jk)
4237	CALL iom_put( "SDT_1000" , zw2d )
4238	ENDIF
4239	IF( med_diag%REG_1000%dgsave ) THEN
4240	CALL iom_put( "REG_1000" , fregen2d )
4241	ENDIF
4242	IF( med_diag%FDT_1000%dgsave ) THEN
4243	CALL iom_put( "FDT_1000" , ffastn )
4244	ENDIF
4245	IF( med_diag%RG_1000F%dgsave ) THEN
4246	CALL iom_put( "RG_1000F" , fregenfast )
4247	ENDIF
4248	IF( med_diag%FDS_1000%dgsave ) THEN
4249	CALL iom_put( "FDS_1000" , ffastsi )
4250	ENDIF
4251	IF( med_diag%RGS1000F%dgsave ) THEN
4252	CALL iom_put( "RGS1000F" , fregenfastsi )
4253	ENDIF
4254	IF( med_diag%FE_1000%dgsave ) THEN
4255	CALL iom_put( "FE_1000" , xFree )
4256	ENDIF
4257	# if defined key_roam
4258	IF( med_diag%RR_1000%dgsave ) THEN
4259	CALL iom_put( "RR_1000" , ffastca2d )
4260	CALL wrk_dealloc( jpi, jpj, ffastca2d )
4261	ENDIF
4262	IF( med_diag%SDC_1000%dgsave ) THEN
4263	zw2d(:,:) = fslowcflux(:,:) * tmask(:,:,jk)
4264	CALL iom_put( "SDC_1000" , zw2d )
4265	ENDIF
4266	# endif
4267	ENDIF
4268	!! to do on every k loop :
4269	IF( med_diag%DETFLUX3%dgsave ) THEN
4270	detflux3d(:,:,jk) = (fslownflux(:,:) + ffastn(:,:)) * tmask(:,:,jk) !! detrital flux
4271	!CALL iom_put( "DETFLUX3" , ftot_n )
4272	ENDIF
4273	# if defined key_roam
4274	IF( med_diag%EXPC3%dgsave ) THEN
4275	expc3(:,:,jk) = (fslowcflux(:,:) + ffastc(:,:)) * tmask(:,:,jk)
4276	ENDIF
4277	IF( med_diag%EXPN3%dgsave ) THEN
4278	expn3(:,:,jk) = (fslownflux(:,:) + ffastn(:,:)) * tmask(:,:,jk)
4279	ENDIF
4280	# endif
4281	ENDIF
4282	!! CLOSE vertical loop
4283	ENDDO
4284
4285	!!----------------------------------------------------------------------
4286	!! Process benthic in/out fluxes
4287	!! These can be handled outside of the 3D calculations since the
4288	!! benthic pools (and fluxes) are 2D in nature; this code is
4289	!! (shamelessly) borrowed from corresponding code in the LOBSTER
4290	!! model
4291	!!----------------------------------------------------------------------
4292	!!
4293	!! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt
4294	if (jorgben.eq.1) then
4295	za_sed_n(:,:) = zn_sed_n(:,:) + &
4296	& ( f_sbenin_n(:,:) + f_fbenin_n(:,:) - f_benout_n(:,:) ) * (rdt / 86400.)
4297	zn_sed_n(:,:) = za_sed_n(:,:)
4298	!!
4299	za_sed_fe(:,:) = zn_sed_fe(:,:) + &
4300	& ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * (rdt / 86400.)
4301	zn_sed_fe(:,:) = za_sed_fe(:,:)
4302	!!
4303	za_sed_c(:,:) = zn_sed_c(:,:) + &
4304	& ( f_sbenin_c(:,:) + f_fbenin_c(:,:) - f_benout_c(:,:) ) * (rdt / 86400.)
4305	zn_sed_c(:,:) = za_sed_c(:,:)
4306	endif
4307	if (jinorgben.eq.1) then
4308	za_sed_si(:,:) = zn_sed_si(:,:) + &
4309	& ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * (rdt / 86400.)
4310	zn_sed_si(:,:) = za_sed_si(:,:)
4311	!!
4312	za_sed_ca(:,:) = zn_sed_ca(:,:) + &
4313	& ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * (rdt / 86400.)
4314	zn_sed_ca(:,:) = za_sed_ca(:,:)
4315	endif
4316	!!
4317	if (ibenthic.eq.2) then
4318	!! The code below (in this if ... then ... endif loop) is
4319	!! effectively commented out because it does not work as
4320	!! anticipated; it can be deleted at a later date
4321	if (jorgben.eq.1) then
4322	za_sed_n(:,:) = ( f_sbenin_n(:,:) + f_fbenin_n(:,:) - f_benout_n(:,:) ) * rdt
4323	za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt
4324	za_sed_c(:,:) = ( f_sbenin_c(:,:) + f_fbenin_c(:,:) - f_benout_c(:,:) ) * rdt
4325	endif
4326	if (jinorgben.eq.1) then
4327	za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt
4328	za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt
4329	endif
4330	!!
4331	!! Leap-frog scheme - only in explicit case, otherwise the time stepping
4332	!! is already being done in trczdf
4333	!! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN
4334	!! zfact = 2. * rdttra(jk) * FLOAT( ndttrc )
4335	!! IF( neuler == 0 .AND. kt == nittrc000 ) zfact = rdttra(jk) * FLOAT(ndttrc)
4336	!! if (jorgben.eq.1) then
4337	!! za_sed_n(:,:) = zb_sed_n(:,:) + ( zfact * za_sed_n(:,:) )
4338	!! za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) )
4339	!! za_sed_c(:,:) = zb_sed_c(:,:) + ( zfact * za_sed_c(:,:) )
4340	!! endif
4341	!! if (jinorgben.eq.1) then
4342	!! za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) )
4343	!! za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) )
4344	!! endif
4345	!! ENDIF
4346	!!
4347	!! Time filter and swap of arrays
4348	IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN ! centred or tvd scheme
4349	IF( neuler == 0 .AND. kt == nittrc000 ) THEN
4350	if (jorgben.eq.1) then
4351	zb_sed_n(:,:) = zn_sed_n(:,:)
4352	zn_sed_n(:,:) = za_sed_n(:,:)
4353	za_sed_n(:,:) = 0.0
4354	!!
4355	zb_sed_fe(:,:) = zn_sed_fe(:,:)
4356	zn_sed_fe(:,:) = za_sed_fe(:,:)
4357	za_sed_fe(:,:) = 0.0
4358	!!
4359	zb_sed_c(:,:) = zn_sed_c(:,:)
4360	zn_sed_c(:,:) = za_sed_c(:,:)
4361	za_sed_c(:,:) = 0.0
4362	endif
4363	if (jinorgben.eq.1) then
4364	zb_sed_si(:,:) = zn_sed_si(:,:)
4365	zn_sed_si(:,:) = za_sed_si(:,:)
4366	za_sed_si(:,:) = 0.0
4367	!!
4368	zb_sed_ca(:,:) = zn_sed_ca(:,:)
4369	zn_sed_ca(:,:) = za_sed_ca(:,:)
4370	za_sed_ca(:,:) = 0.0
4371	endif
4372	ELSE
4373	if (jorgben.eq.1) then
4374	zb_sed_n(:,:) = (atfp * ( zb_sed_n(:,:) + za_sed_n(:,:) )) + (atfp1 * zn_sed_n(:,:) )
4375	zn_sed_n(:,:) = za_sed_n(:,:)
4376	za_sed_n(:,:) = 0.0
4377	!!
4378	zb_sed_fe(:,:) = (atfp * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:))
4379	zn_sed_fe(:,:) = za_sed_fe(:,:)
4380	za_sed_fe(:,:) = 0.0
4381	!!
4382	zb_sed_c(:,:) = (atfp * ( zb_sed_c(:,:) + za_sed_c(:,:) )) + (atfp1 * zn_sed_c(:,:) )
4383	zn_sed_c(:,:) = za_sed_c(:,:)
4384	za_sed_c(:,:) = 0.0
4385	endif
4386	if (jinorgben.eq.1) then
4387	zb_sed_si(:,:) = (atfp * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:))
4388	zn_sed_si(:,:) = za_sed_si(:,:)
4389	za_sed_si(:,:) = 0.0
4390	!!
4391	zb_sed_ca(:,:) = (atfp * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:))
4392	zn_sed_ca(:,:) = za_sed_ca(:,:)
4393	za_sed_ca(:,:) = 0.0
4394	endif
4395	ENDIF
4396	ELSE ! case of smolar scheme or muscl
4397	if (jorgben.eq.1) then
4398	zb_sed_n(:,:) = za_sed_n(:,:)
4399	zn_sed_n(:,:) = za_sed_n(:,:)
4400	za_sed_n(:,:) = 0.0
4401	!!
4402	zb_sed_fe(:,:) = za_sed_fe(:,:)
4403	zn_sed_fe(:,:) = za_sed_fe(:,:)
4404	za_sed_fe(:,:) = 0.0
4405	!!
4406	zb_sed_c(:,:) = za_sed_c(:,:)
4407	zn_sed_c(:,:) = za_sed_c(:,:)
4408	za_sed_c(:,:) = 0.0
4409	endif
4410	if (jinorgben.eq.1) then
4411	zb_sed_si(:,:) = za_sed_si(:,:)
4412	zn_sed_si(:,:) = za_sed_si(:,:)
4413	za_sed_si(:,:) = 0.0
4414	!!
4415	zb_sed_ca(:,:) = za_sed_ca(:,:)
4416	zn_sed_ca(:,:) = za_sed_ca(:,:)
4417	za_sed_ca(:,:) = 0.0
4418	endif
4419	ENDIF
4420	endif
4421
4422	!!-----------------------------------------------------------------------
4423	!! Add very last diag calculations
4424	!!-----------------------------------------------------------------------
4425	DO jj = 2,jpjm1
4426	DO ji = 2,jpim1
4427	!!
4428	IF( med_diag%PN_JLIM%dgsave ) THEN
4429	fjln2d(ji,jj) = fjln2d(ji,jj) / MAX(ftot_pn(ji,jj), rsmall)
4430	ENDIF
4431	IF( med_diag%PN_NLIM%dgsave ) THEN
4432	fnln2d(ji,jj) = fnln2d(ji,jj) / MAX(ftot_pn(ji,jj), rsmall)
4433	ENDIF
4434	IF( med_diag%PN_FELIM%dgsave ) THEN
4435	ffln2d(ji,jj) = ffln2d(ji,jj) / MAX(ftot_pn(ji,jj), rsmall)
4436	ENDIF
4437	IF( med_diag%PD_JLIM%dgsave ) THEN
4438	fjld2d(ji,jj) = fjld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
4439	ENDIF
4440	IF( med_diag%PD_NLIM%dgsave ) THEN
4441	fnld2d(ji,jj) = fnld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
4442	ENDIF
4443	IF( med_diag%PD_FELIM%dgsave ) THEN
4444	ffld2d(ji,jj) = ffld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
4445	ENDIF
4446	IF( med_diag%PD_SILIM%dgsave ) THEN
4447	fsld2d2(ji,jj) = fsld2d2(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
4448	ENDIF
4449	IF( med_diag%PDSILIM2%dgsave ) THEN
4450	fsld2d(ji,jj) = fsld2d(ji,jj) / MAX(ftot_pd(ji,jj), rsmall)
4451	ENDIF
4452	ENDDO
4453	ENDDO
4454
4455	!!----------------------------------------------------------------------
4456	!! Add in XML diagnostics stuff
4457	!!----------------------------------------------------------------------
4458	!!
4459	!! ** 2D diagnostics
4460	# if defined key_debug_medusa
4461	IF ( lwp ) write (numout,*) 'trc_bio_medusa: export all diag.'
4462	CALL flush(numout)
4463	# endif
4464	IF ( med_diag%INVTN%dgsave ) THEN
4465	CALL iom_put( "INVTN" , ftot_n )
4466	ENDIF
4467	IF ( med_diag%INVTSI%dgsave ) THEN
4468	CALL iom_put( "INVTSI" , ftot_si )
4469	ENDIF
4470	IF ( med_diag%INVTFE%dgsave ) THEN
4471	CALL iom_put( "INVTFE" , ftot_fe )
4472	ENDIF
4473	IF ( med_diag%ML_PRN%dgsave ) THEN
4474	CALL iom_put( "ML_PRN" , fprn_ml )
4475	ENDIF
4476	IF ( med_diag%ML_PRD%dgsave ) THEN
4477	CALL iom_put( "ML_PRD" , fprd_ml )
4478	ENDIF
4479	IF ( med_diag%OCAL_LVL%dgsave ) THEN
4480	CALL iom_put( "OCAL_LVL" , fccd )
4481	ENDIF
4482	IF ( med_diag%PN_JLIM%dgsave ) THEN
4483	CALL iom_put( "PN_JLIM" , fjln2d )
4484	CALL wrk_dealloc( jpi, jpj, fjln2d )
4485	ENDIF
4486	IF ( med_diag%PN_NLIM%dgsave ) THEN
4487	CALL iom_put( "PN_NLIM" , fnln2d )
4488	CALL wrk_dealloc( jpi, jpj, fnln2d )
4489	ENDIF
4490	IF ( med_diag%PN_FELIM%dgsave ) THEN
4491	CALL iom_put( "PN_FELIM" , ffln2d )
4492	CALL wrk_dealloc( jpi, jpj, ffln2d )
4493	ENDIF
4494	IF ( med_diag%PD_JLIM%dgsave ) THEN
4495	CALL iom_put( "PD_JLIM" , fjld2d )
4496	CALL wrk_dealloc( jpi, jpj, fjld2d )
4497	ENDIF
4498	IF ( med_diag%PD_NLIM%dgsave ) THEN
4499	CALL iom_put( "PD_NLIM" , fnld2d )
4500	CALL wrk_dealloc( jpi, jpj, fnld2d )
4501	ENDIF
4502	IF ( med_diag%PD_FELIM%dgsave ) THEN
4503	CALL iom_put( "PD_FELIM" , ffld2d )
4504	CALL wrk_dealloc( jpi, jpj, ffld2d )
4505	ENDIF
4506	IF ( med_diag%PD_SILIM%dgsave ) THEN
4507	CALL iom_put( "PD_SILIM" , fsld2d2 )
4508	CALL wrk_dealloc( jpi, jpj, fsld2d2 )
4509	ENDIF
4510	IF ( med_diag%PDSILIM2%dgsave ) THEN
4511	CALL iom_put( "PDSILIM2" , fsld2d )
4512	CALL wrk_dealloc( jpi, jpj, fsld2d )
4513	ENDIF
4514	IF ( med_diag%INTFLX_N%dgsave ) THEN
4515	CALL iom_put( "INTFLX_N" , fflx_n )
4516	ENDIF
4517	IF ( med_diag%INTFLX_SI%dgsave ) THEN
4518	CALL iom_put( "INTFLX_SI" , fflx_si )
4519	ENDIF
4520	IF ( med_diag%INTFLX_FE%dgsave ) THEN
4521	CALL iom_put( "INTFLX_FE" , fflx_fe )
4522	ENDIF
4523	IF ( med_diag%INT_PN%dgsave ) THEN
4524	CALL iom_put( "INT_PN" , ftot_pn )
4525	ENDIF
4526	IF ( med_diag%INT_PD%dgsave ) THEN
4527	CALL iom_put( "INT_PD" , ftot_pd )
4528	ENDIF
4529	IF ( med_diag%INT_ZMI%dgsave ) THEN
4530	CALL iom_put( "INT_ZMI" , ftot_zmi )
4531	ENDIF
4532	IF ( med_diag%INT_ZME%dgsave ) THEN
4533	CALL iom_put( "INT_ZME" , ftot_zme )
4534	ENDIF
4535	IF ( med_diag%INT_DET%dgsave ) THEN
4536	CALL iom_put( "INT_DET" , ftot_det )
4537	ENDIF
4538	IF ( med_diag%INT_DTC%dgsave ) THEN
4539	CALL iom_put( "INT_DTC" , ftot_dtc )
4540	ENDIF
4541	IF ( med_diag%BEN_N%dgsave ) THEN
4542	CALL iom_put( "BEN_N" , za_sed_n )
4543	ENDIF
4544	IF ( med_diag%BEN_FE%dgsave ) THEN
4545	CALL iom_put( "BEN_FE" , za_sed_fe )
4546	ENDIF
4547	IF ( med_diag%BEN_C%dgsave ) THEN
4548	CALL iom_put( "BEN_C" , za_sed_c )
4549	ENDIF
4550	IF ( med_diag%BEN_SI%dgsave ) THEN
4551	CALL iom_put( "BEN_SI" , za_sed_si )
4552	ENDIF
4553	IF ( med_diag%BEN_CA%dgsave ) THEN
4554	CALL iom_put( "BEN_CA" , za_sed_ca )
4555	ENDIF
4556	IF ( med_diag%RUNOFF%dgsave ) THEN
4557	CALL iom_put( "RUNOFF" , f_runoff )
4558	ENDIF
4559	# if defined key_roam
4560	IF ( med_diag%N_PROD%dgsave ) THEN
4561	CALL iom_put( "N_PROD" , fnit_prod )
4562	ENDIF
4563	IF ( med_diag%N_CONS%dgsave ) THEN
4564	CALL iom_put( "N_CONS" , fnit_cons )
4565	ENDIF
4566	IF ( med_diag%C_PROD%dgsave ) THEN
4567	CALL iom_put( "C_PROD" , fcar_prod )
4568	ENDIF
4569	IF ( med_diag%C_CONS%dgsave ) THEN
4570	CALL iom_put( "C_CONS" , fcar_cons )
4571	ENDIF
4572	IF ( med_diag%O2_PROD%dgsave ) THEN
4573	CALL iom_put( "O2_PROD" , foxy_prod )
4574	ENDIF
4575	IF ( med_diag%O2_CONS%dgsave ) THEN
4576	CALL iom_put( "O2_CONS" , foxy_cons )
4577	ENDIF
4578	IF ( med_diag%O2_ANOX%dgsave ) THEN
4579	CALL iom_put( "O2_ANOX" , foxy_anox )
4580	ENDIF
4581	IF ( med_diag%INVTC%dgsave ) THEN
4582	CALL iom_put( "INVTC" , ftot_c )
4583	ENDIF
4584	IF ( med_diag%INVTALK%dgsave ) THEN
4585	CALL iom_put( "INVTALK" , ftot_a )
4586	ENDIF
4587	IF ( med_diag%INVTO2%dgsave ) THEN
4588	CALL iom_put( "INVTO2" , ftot_o2 )
4589	ENDIF
4590	IF ( med_diag%COM_RESP%dgsave ) THEN
4591	CALL iom_put( "COM_RESP" , fcomm_resp )
4592	ENDIF
4593	# endif
4594	!!
4595	!! diagnostic filled in the i-j-k main loop
4596	!!--------------------------------------------
4597	IF ( med_diag%PRN%dgsave ) THEN
4598	CALL iom_put( "PRN" , fprn2d )
4599	CALL wrk_dealloc( jpi, jpj, fprn2d )
4600	ENDIF
4601	IF ( med_diag%MPN%dgsave ) THEN
4602	CALL iom_put( "MPN" ,fdpn2d )
4603	CALL wrk_dealloc( jpi, jpj, fdpn2d )
4604	ENDIF
4605	IF ( med_diag%PRD%dgsave ) THEN
4606	CALL iom_put( "PRD" ,fprd2d )
4607	CALL wrk_dealloc( jpi, jpj, fprd2d )
4608	ENDIF
4609	IF( med_diag%MPD%dgsave ) THEN
4610	CALL iom_put( "MPD" , fdpd2d )
4611	CALL wrk_dealloc( jpi, jpj, fdpd2d )
4612	ENDIF
4613	! IF( med_diag%DSED%dgsave ) THEN
4614	! CALL iom_put( "DSED" , ftot_n )
4615	! ENDIF
4616	IF( med_diag%OPAL%dgsave ) THEN
4617	CALL iom_put( "OPAL" , fprds2d )
4618	CALL wrk_dealloc( jpi, jpj, fprds2d )
4619	ENDIF
4620	IF( med_diag%OPALDISS%dgsave ) THEN
4621	CALL iom_put( "OPALDISS" , fsdiss2d )
4622	CALL wrk_dealloc( jpi, jpj, fsdiss2d )
4623	ENDIF
4624	IF( med_diag%GMIPn%dgsave ) THEN
4625	CALL iom_put( "GMIPn" , fgmipn2d )
4626	CALL wrk_dealloc( jpi, jpj, fgmipn2d )
4627	ENDIF
4628	IF( med_diag%GMID%dgsave ) THEN
4629	CALL iom_put( "GMID" , fgmid2d )
4630	CALL wrk_dealloc( jpi, jpj, fgmid2d )
4631	ENDIF
4632	IF( med_diag%MZMI%dgsave ) THEN
4633	CALL iom_put( "MZMI" , fdzmi2d )
4634	CALL wrk_dealloc( jpi, jpj, fdzmi2d )
4635	ENDIF
4636	IF( med_diag%GMEPN%dgsave ) THEN
4637	CALL iom_put( "GMEPN" , fgmepn2d )
4638	CALL wrk_dealloc( jpi, jpj, fgmepn2d )
4639	ENDIF
4640	IF( med_diag%GMEPD%dgsave ) THEN
4641	CALL iom_put( "GMEPD" , fgmepd2d )
4642	CALL wrk_dealloc( jpi, jpj, fgmepd2d )
4643	ENDIF
4644	IF( med_diag%GMEZMI%dgsave ) THEN
4645	CALL iom_put( "GMEZMI" , fgmezmi2d )
4646	CALL wrk_dealloc( jpi, jpj, fgmezmi2d )
4647	ENDIF
4648	IF( med_diag%GMED%dgsave ) THEN
4649	CALL iom_put( "GMED" , fgmed2d )
4650	CALL wrk_dealloc( jpi, jpj, fgmed2d )
4651	ENDIF
4652	IF( med_diag%MZME%dgsave ) THEN
4653	CALL iom_put( "MZME" , fdzme2d )
4654	CALL wrk_dealloc( jpi, jpj, fdzme2d )
4655	ENDIF
4656	! IF( med_diag%DEXP%dgsave ) THEN
4657	! CALL iom_put( "DEXP" , ftot_n )
4658	! ENDIF
4659	IF( med_diag%DETN%dgsave ) THEN
4660	CALL iom_put( "DETN" , fslown2d )
4661	CALL wrk_dealloc( jpi, jpj, fslown2d )
4662	ENDIF
4663	IF( med_diag%MDET%dgsave ) THEN
4664	CALL iom_put( "MDET" , fdd2d )
4665	CALL wrk_dealloc( jpi, jpj, fdd2d )
4666	ENDIF
4667	IF( med_diag%AEOLIAN%dgsave ) THEN
4668	CALL iom_put( "AEOLIAN" , ffetop2d )
4669	CALL wrk_dealloc( jpi, jpj, ffetop2d )
4670	ENDIF
4671	IF( med_diag%BENTHIC%dgsave ) THEN
4672	CALL iom_put( "BENTHIC" , ffebot2d )
4673	CALL wrk_dealloc( jpi, jpj, ffebot2d )
4674	ENDIF
4675	IF( med_diag%SCAVENGE%dgsave ) THEN
4676	CALL iom_put( "SCAVENGE" , ffescav2d )
4677	CALL wrk_dealloc( jpi, jpj, ffescav2d )
4678	ENDIF
4679	!!
4680	IF( med_diag%TOTREG_N%dgsave ) THEN
4681	CALL iom_put( "TOTREG_N" , fregen2d )
4682	CALL wrk_dealloc( jpi, jpj, fregen2d )
4683	ENDIF
4684	IF( med_diag%TOTRG_SI%dgsave ) THEN
4685	CALL iom_put( "TOTRG_SI" , fregensi2d )
4686	CALL wrk_dealloc( jpi, jpj, fregensi2d )
4687	ENDIF
4688	!!
4689	IF( med_diag%FASTN%dgsave ) THEN
4690	CALL iom_put( "FASTN" , ftempn2d )
4691	CALL wrk_dealloc( jpi, jpj, ftempn2d )
4692	ENDIF
4693	IF( med_diag%FASTSI%dgsave ) THEN
4694	CALL iom_put( "FASTSI" , ftempsi2d )
4695	CALL wrk_dealloc( jpi, jpj, ftempsi2d )
4696	ENDIF
4697	IF( med_diag%FASTFE%dgsave ) THEN
4698	CALL iom_put( "FASTFE" , ftempfe2d )
4699	CALL wrk_dealloc( jpi, jpj, ftempfe2d )
4700	ENDIF
4701	IF( med_diag%FASTC%dgsave ) THEN
4702	CALL iom_put( "FASTC" , ftempc2d )
4703	CALL wrk_dealloc( jpi, jpj, ftempc2d )
4704	ENDIF
4705	IF( med_diag%FASTCA%dgsave ) THEN
4706	CALL iom_put( "FASTCA" , ftempca2d )
4707	CALL wrk_dealloc( jpi, jpj, ftempca2d )
4708	ENDIF
4709	!!
4710	IF( med_diag%REMINN%dgsave ) THEN
4711	CALL iom_put( "REMINN" , freminn2d )
4712	CALL wrk_dealloc( jpi, jpj, freminn2d )
4713	ENDIF
4714	IF( med_diag%REMINSI%dgsave ) THEN
4715	CALL iom_put( "REMINSI" , freminsi2d )
4716	CALL wrk_dealloc( jpi, jpj, freminsi2d )
4717	ENDIF
4718	IF( med_diag%REMINFE%dgsave ) THEN
4719	CALL iom_put( "REMINFE" , freminfe2d )
4720	CALL wrk_dealloc( jpi, jpj, freminfe2d )
4721	ENDIF
4722	IF( med_diag%REMINC%dgsave ) THEN
4723	CALL iom_put( "REMINC" , freminc2d )
4724	CALL wrk_dealloc( jpi, jpj, freminc2d )
4725	ENDIF
4726	IF( med_diag%REMINCA%dgsave ) THEN
4727	CALL iom_put( "REMINCA" , freminca2d )
4728	CALL wrk_dealloc( jpi, jpj, freminca2d )
4729	ENDIF
4730	IF( med_diag%SEAFLRN%dgsave ) THEN
4731	CALL iom_put( "SEAFLRN" , fsedn )
4732	ENDIF
4733	IF( med_diag%SEAFLRSI%dgsave ) THEN
4734	CALL iom_put( "SEAFLRSI" , fsedsi )
4735	ENDIF
4736	IF( med_diag%SEAFLRFE%dgsave ) THEN
4737	CALL iom_put( "SEAFLRFE" , fsedfe )
4738	ENDIF
4739	IF( med_diag%SEAFLRC%dgsave ) THEN
4740	CALL iom_put( "SEAFLRC" , fsedc )
4741	ENDIF
4742	IF( med_diag%SEAFLRCA%dgsave ) THEN
4743	CALL iom_put( "SEAFLRCA" , fsedca )
4744	ENDIF
4745	!!
4746	# if defined key_roam
4747	!!
4748	IF( med_diag%RIV_N%dgsave ) THEN
4749	CALL iom_put( "RIV_N" , rivn2d )
4750	CALL wrk_dealloc( jpi, jpj, rivn2d )
4751	ENDIF
4752	IF( med_diag%RIV_SI%dgsave ) THEN
4753	CALL iom_put( "RIV_SI" , rivsi2d )
4754	CALL wrk_dealloc( jpi, jpj, rivsi2d )
4755	ENDIF
4756	IF( med_diag%RIV_C%dgsave ) THEN
4757	CALL iom_put( "RIV_C" , rivc2d )
4758	CALL wrk_dealloc( jpi, jpj, rivc2d )
4759	ENDIF
4760	IF( med_diag%RIV_ALK%dgsave ) THEN
4761	CALL iom_put( "RIV_ALK" , rivalk2d )
4762	CALL wrk_dealloc( jpi, jpj, rivalk2d )
4763	ENDIF
4764	IF( med_diag%DETC%dgsave ) THEN
4765	CALL iom_put( "DETC" , fslowc2d )
4766	CALL wrk_dealloc( jpi, jpj, fslowc2d )
4767	ENDIF
4768	!!
4769	IF( med_diag%PN_LLOSS%dgsave ) THEN
4770	CALL iom_put( "PN_LLOSS" , fdpn22d )
4771	CALL wrk_dealloc( jpi, jpj, fdpn22d )
4772	ENDIF
4773	IF( med_diag%PD_LLOSS%dgsave ) THEN
4774	CALL iom_put( "PD_LLOSS" , fdpd22d )
4775	CALL wrk_dealloc( jpi, jpj, fdpd22d )
4776	ENDIF
4777	IF( med_diag%ZI_LLOSS%dgsave ) THEN
4778	CALL iom_put( "ZI_LLOSS" , fdzmi22d )
4779	CALL wrk_dealloc( jpi, jpj, fdzmi22d )
4780	ENDIF
4781	IF( med_diag%ZE_LLOSS%dgsave ) THEN
4782	CALL iom_put( "ZE_LLOSS" , fdzme22d )
4783	CALL wrk_dealloc( jpi, jpj, fdzme22d )
4784	ENDIF
4785	IF( med_diag%ZI_MES_N%dgsave ) THEN
4786	CALL iom_put( "ZI_MES_N" , zimesn2d )
4787	CALL wrk_dealloc( jpi, jpj, zimesn2d )
4788	ENDIF
4789	IF( med_diag%ZI_MES_D%dgsave ) THEN
4790	CALL iom_put( "ZI_MES_D" , zimesd2d )
4791	CALL wrk_dealloc( jpi, jpj, zimesd2d )
4792	ENDIF
4793	IF( med_diag%ZI_MES_C%dgsave ) THEN
4794	CALL iom_put( "ZI_MES_C" , zimesc2d )
4795	CALL wrk_dealloc( jpi, jpj, zimesc2d )
4796	ENDIF
4797	IF( med_diag%ZI_MESDC%dgsave ) THEN
4798	CALL iom_put( "ZI_MESDC" ,zimesdc2d )
4799	CALL wrk_dealloc( jpi, jpj, zimesdc2d )
4800	ENDIF
4801	IF( med_diag%ZI_EXCR%dgsave ) THEN
4802	CALL iom_put( "ZI_EXCR" , ziexcr2d )
4803	CALL wrk_dealloc( jpi, jpj, ziexcr2d )
4804	ENDIF
4805	IF( med_diag%ZI_RESP%dgsave ) THEN
4806	CALL iom_put( "ZI_RESP" , ziresp2d )
4807	CALL wrk_dealloc( jpi, jpj, ziresp2d )
4808	ENDIF
4809	IF( med_diag%ZI_GROW%dgsave ) THEN
4810	CALL iom_put( "ZI_GROW" , zigrow2d )
4811	CALL wrk_dealloc( jpi, jpj, zigrow2d )
4812	ENDIF
4813	IF( med_diag%ZE_MES_N%dgsave ) THEN
4814	CALL iom_put( "ZE_MES_N" , zemesn2d )
4815	CALL wrk_dealloc( jpi, jpj, zemesn2d )
4816	ENDIF
4817	IF( med_diag%ZE_MES_D%dgsave ) THEN
4818	CALL iom_put( "ZE_MES_D" , zemesd2d )
4819	CALL wrk_dealloc( jpi, jpj, zemesd2d )
4820	ENDIF
4821	IF( med_diag%ZE_MES_C%dgsave ) THEN
4822	CALL iom_put( "ZE_MES_C" , zemesc2d )
4823	CALL wrk_dealloc( jpi, jpj, zemesc2d )
4824	ENDIF
4825	IF( med_diag%ZE_MESDC%dgsave ) THEN
4826	CALL iom_put( "ZE_MESDC" , zemesdc2d )
4827	CALL wrk_dealloc( jpi, jpj, zemesdc2d )
4828	ENDIF
4829	IF( med_diag%ZE_EXCR%dgsave ) THEN
4830	CALL iom_put( "ZE_EXCR" , zeexcr2d )
4831	CALL wrk_dealloc( jpi, jpj, zeexcr2d )
4832	ENDIF
4833	IF( med_diag%ZE_RESP%dgsave ) THEN
4834	CALL iom_put( "ZE_RESP" , zeresp2d )
4835	CALL wrk_dealloc( jpi, jpj, zeresp2d )
4836	ENDIF
4837	IF( med_diag%ZE_GROW%dgsave ) THEN
4838	CALL iom_put( "ZE_GROW" , zegrow2d )
4839	CALL wrk_dealloc( jpi, jpj, zegrow2d )
4840	ENDIF
4841	IF( med_diag%MDETC%dgsave ) THEN
4842	CALL iom_put( "MDETC" , mdetc2d )
4843	CALL wrk_dealloc( jpi, jpj, mdetc2d )
4844	ENDIF
4845	IF( med_diag%GMIDC%dgsave ) THEN
4846	CALL iom_put( "GMIDC" , gmidc2d )
4847	CALL wrk_dealloc( jpi, jpj, gmidc2d )
4848	ENDIF
4849	IF( med_diag%GMEDC%dgsave ) THEN
4850	CALL iom_put( "GMEDC" , gmedc2d )
4851	CALL wrk_dealloc( jpi, jpj, gmedc2d )
4852	ENDIF
4853	IF( med_diag%IBEN_N%dgsave ) THEN
4854	CALL iom_put( "IBEN_N" , iben_n2d )
4855	CALL wrk_dealloc( jpi, jpj, iben_n2d )
4856	ENDIF
4857	IF( med_diag%IBEN_FE%dgsave ) THEN
4858	CALL iom_put( "IBEN_FE" , iben_fe2d )
4859	CALL wrk_dealloc( jpi, jpj, iben_fe2d )
4860	ENDIF
4861	IF( med_diag%IBEN_C%dgsave ) THEN
4862	CALL iom_put( "IBEN_C" , iben_c2d )
4863	CALL wrk_dealloc( jpi, jpj, iben_c2d )
4864	ENDIF
4865	IF( med_diag%IBEN_SI%dgsave ) THEN
4866	CALL iom_put( "IBEN_SI" , iben_si2d )
4867	CALL wrk_dealloc( jpi, jpj, iben_si2d )
4868	ENDIF
4869	IF( med_diag%IBEN_CA%dgsave ) THEN
4870	CALL iom_put( "IBEN_CA" , iben_ca2d )
4871	CALL wrk_dealloc( jpi, jpj, iben_ca2d )
4872	ENDIF
4873	IF( med_diag%OBEN_N%dgsave ) THEN
4874	CALL iom_put( "OBEN_N" , oben_n2d )
4875	CALL wrk_dealloc( jpi, jpj, oben_n2d )
4876	ENDIF
4877	IF( med_diag%OBEN_FE%dgsave ) THEN
4878	CALL iom_put( "OBEN_FE" , oben_fe2d )
4879	CALL wrk_dealloc( jpi, jpj, oben_fe2d )
4880	ENDIF
4881	IF( med_diag%OBEN_C%dgsave ) THEN
4882	CALL iom_put( "OBEN_C" , oben_c2d )
4883	CALL wrk_dealloc( jpi, jpj, oben_c2d )
4884	ENDIF
4885	IF( med_diag%OBEN_SI%dgsave ) THEN
4886	CALL iom_put( "OBEN_SI" , oben_si2d )
4887	CALL wrk_dealloc( jpi, jpj, oben_si2d )
4888	ENDIF
4889	IF( med_diag%OBEN_CA%dgsave ) THEN
4890	CALL iom_put( "OBEN_CA" , oben_ca2d )
4891	CALL wrk_dealloc( jpi, jpj, oben_ca2d )
4892	ENDIF
4893	IF( med_diag%SFR_OCAL%dgsave ) THEN
4894	CALL iom_put( "SFR_OCAL" , sfr_ocal2d )
4895	CALL wrk_dealloc( jpi, jpj, sfr_ocal2d )
4896	ENDIF
4897	IF( med_diag%SFR_OARG%dgsave ) THEN
4898	CALL iom_put( "SFR_OARG" , sfr_oarg2d )
4899	CALL wrk_dealloc( jpi, jpj, sfr_oarg2d )
4900	ENDIF
4901	IF( med_diag%LYSO_CA%dgsave ) THEN
4902	CALL iom_put( "LYSO_CA" , lyso_ca2d )
4903	CALL wrk_dealloc( jpi, jpj, lyso_ca2d )
4904	ENDIF
4905	# endif
4906	!!
4907	!! ** 3D diagnostics
4908	IF( med_diag%TPP3%dgsave ) THEN
4909	CALL iom_put( "TPP3" , tpp3d )
4910	CALL wrk_dealloc( jpi, jpj, jpk, tpp3d )
4911	ENDIF
4912	IF( med_diag%DETFLUX3%dgsave ) THEN
4913	CALL iom_put( "DETFLUX3" , detflux3d )
4914	CALL wrk_dealloc( jpi, jpj, jpk, detflux3d )
4915	ENDIF
4916	IF( med_diag%REMIN3N%dgsave ) THEN
4917	CALL iom_put( "REMIN3N" , remin3dn )
4918	CALL wrk_dealloc( jpi, jpj, jpk, remin3dn )
4919	ENDIF
4920	# if defined key_roam
4921	IF( med_diag%PH3%dgsave ) THEN
4922	CALL iom_put( "PH3" , f3_pH )
4923	ENDIF
4924	IF( med_diag%OM_CAL3%dgsave ) THEN
4925	CALL iom_put( "OM_CAL3" , f3_omcal )
4926	ENDIF
4927	!!
4928	!! AXY (09/11/16): 2D CMIP6 diagnostics
4929	IF( med_diag%INTDISSIC%dgsave ) THEN
4930	CALL iom_put( "INTDISSIC" , intdissic )
4931	CALL wrk_dealloc( jpi, jpj, intdissic )
4932	ENDIF
4933	IF( med_diag%INTDISSIN%dgsave ) THEN
4934	CALL iom_put( "INTDISSIN" , intdissin )
4935	CALL wrk_dealloc( jpi, jpj, intdissin )
4936	ENDIF
4937	IF( med_diag%INTDISSISI%dgsave ) THEN
4938	CALL iom_put( "INTDISSISI" , intdissisi )
4939	CALL wrk_dealloc( jpi, jpj, intdissisi )
4940	ENDIF
4941	IF( med_diag%INTTALK%dgsave ) THEN
4942	CALL iom_put( "INTTALK" , inttalk )
4943	CALL wrk_dealloc( jpi, jpj, inttalk )
4944	ENDIF
4945	IF( med_diag%O2min%dgsave ) THEN
4946	CALL iom_put( "O2min" , o2min )
4947	CALL wrk_dealloc( jpi, jpj, o2min )
4948	ENDIF
4949	IF( med_diag%ZO2min%dgsave ) THEN
4950	CALL iom_put( "ZO2min" , zo2min )
4951	CALL wrk_dealloc( jpi, jpj, zo2min )
4952	ENDIF
4953	IF( med_diag%FBDDTALK%dgsave ) THEN
4954	CALL iom_put( "FBDDTALK" , fbddtalk )
4955	CALL wrk_dealloc( jpi, jpj, fbddtalk )
4956	ENDIF
4957	IF( med_diag%FBDDTDIC%dgsave ) THEN
4958	CALL iom_put( "FBDDTDIC" , fbddtdic )
4959	CALL wrk_dealloc( jpi, jpj, fbddtdic )
4960	ENDIF
4961	IF( med_diag%FBDDTDIFE%dgsave ) THEN
4962	CALL iom_put( "FBDDTDIFE" , fbddtdife )
4963	CALL wrk_dealloc( jpi, jpj, fbddtdife )
4964	ENDIF
4965	IF( med_diag%FBDDTDIN%dgsave ) THEN
4966	CALL iom_put( "FBDDTDIN" , fbddtdin )
4967	CALL wrk_dealloc( jpi, jpj, fbddtdin )
4968	ENDIF
4969	IF( med_diag%FBDDTDISI%dgsave ) THEN
4970	CALL iom_put( "FBDDTDISI" , fbddtdisi )
4971	CALL wrk_dealloc( jpi, jpj, fbddtdisi )
4972	ENDIF
4973	!!
4974	!! AXY (09/11/16): 3D CMIP6 diagnostics
4975	IF( med_diag%TPPD3%dgsave ) THEN
4976	CALL iom_put( "TPPD3" , tppd3 )
4977	CALL wrk_dealloc( jpi, jpj, jpk, tppd3 )
4978	ENDIF
4979	IF( med_diag%BDDTALK3%dgsave ) THEN
4980	CALL iom_put( "BDDTALK3" , bddtalk3 )
4981	CALL wrk_dealloc( jpi, jpj, jpk, bddtalk3 )
4982	ENDIF
4983	IF( med_diag%BDDTDIC3%dgsave ) THEN
4984	CALL iom_put( "BDDTDIC3" , bddtdic3 )
4985	CALL wrk_dealloc( jpi, jpj, jpk, bddtdic3 )
4986	ENDIF
4987	IF( med_diag%BDDTDIFE3%dgsave ) THEN
4988	CALL iom_put( "BDDTDIFE3" , bddtdife3 )
4989	CALL wrk_dealloc( jpi, jpj, jpk, bddtdife3 )
4990	ENDIF
4991	IF( med_diag%BDDTDIN3%dgsave ) THEN
4992	CALL iom_put( "BDDTDIN3" , bddtdin3 )
4993	CALL wrk_dealloc( jpi, jpj, jpk, bddtdin3 )
4994	ENDIF
4995	IF( med_diag%BDDTDISI3%dgsave ) THEN
4996	CALL iom_put( "BDDTDISI3" , bddtdisi3 )
4997	CALL wrk_dealloc( jpi, jpj, jpk, bddtdisi3 )
4998	ENDIF
4999	IF( med_diag%FD_NIT3%dgsave ) THEN
5000	CALL iom_put( "FD_NIT3" , fd_nit3 )
5001	CALL wrk_dealloc( jpi, jpj, jpk, fd_nit3 )
5002	ENDIF
5003	IF( med_diag%FD_SIL3%dgsave ) THEN
5004	CALL iom_put( "FD_SIL3" , fd_sil3 )
5005	CALL wrk_dealloc( jpi, jpj, jpk, fd_sil3 )
5006	ENDIF
5007	IF( med_diag%FD_CAL3%dgsave ) THEN
5008	CALL iom_put( "FD_CAL3" , fd_cal3 )
5009	CALL wrk_dealloc( jpi, jpj, jpk, fd_cal3 )
5010	ENDIF
5011	IF( med_diag%FD_CAR3%dgsave ) THEN
5012	CALL iom_put( "FD_CAR3" , fd_car3 )
5013	CALL wrk_dealloc( jpi, jpj, jpk, fd_car3 )
5014	ENDIF
5015	IF( med_diag%CO33%dgsave ) THEN
5016	CALL iom_put( "CO33" , f3_co3 )
5017	ENDIF
5018	IF( med_diag%CO3SATARAG3%dgsave ) THEN
5019	CALL iom_put( "CO3SATARAG3" , f3_omarg )
5020	ENDIF
5021	IF( med_diag%CO3SATCALC3%dgsave ) THEN
5022	CALL iom_put( "CO3SATCALC3" , f3_omcal )
5023	ENDIF
5024	IF( med_diag%EXPC3%dgsave ) THEN
5025	CALL iom_put( "EXPC3" , expc3 )
5026	CALL wrk_dealloc( jpi, jpj, jpk, expc3 )
5027	ENDIF
5028	IF( med_diag%EXPN3%dgsave ) THEN
5029	CALL iom_put( "EXPN3" , expn3 )
5030	CALL wrk_dealloc( jpi, jpj, jpk, expn3 )
5031	ENDIF
5032	IF( med_diag%DCALC3%dgsave ) THEN
5033	CALL iom_put( "DCALC3" , dcalc3 )
5034	CALL wrk_dealloc( jpi, jpj, jpk, dcalc3 )
5035	ENDIF
5036	IF( med_diag%FEDISS3%dgsave ) THEN
5037	CALL iom_put( "FEDISS3" , fediss3 )
5038	CALL wrk_dealloc( jpi, jpj, jpk, fediss3 )
5039	ENDIF
5040	IF( med_diag%FESCAV3%dgsave ) THEN
5041	CALL iom_put( "FESCAV3" , fescav3 )
5042	CALL wrk_dealloc( jpi, jpj, jpk, fescav3 )
5043	ENDIF
5044	IF( med_diag%MIGRAZP3%dgsave ) THEN
5045	CALL iom_put( "MIGRAZP3" , migrazp3 )
5046	CALL wrk_dealloc( jpi, jpj, jpk, migrazp3 )
5047	ENDIF
5048	IF( med_diag%MIGRAZD3%dgsave ) THEN
5049	CALL iom_put( "MIGRAZD3" , migrazd3 )
5050	CALL wrk_dealloc( jpi, jpj, jpk, migrazd3 )
5051	ENDIF
5052	IF( med_diag%MEGRAZP3%dgsave ) THEN
5053	CALL iom_put( "MEGRAZP3" , megrazp3 )
5054	CALL wrk_dealloc( jpi, jpj, jpk, megrazp3 )
5055	ENDIF
5056	IF( med_diag%MEGRAZD3%dgsave ) THEN
5057	CALL iom_put( "MEGRAZD3" , megrazd3 )
5058	CALL wrk_dealloc( jpi, jpj, jpk, megrazd3 )
5059	ENDIF
5060	IF( med_diag%MEGRAZZ3%dgsave ) THEN
5061	CALL iom_put( "MEGRAZZ3" , megrazz3 )
5062	CALL wrk_dealloc( jpi, jpj, jpk, megrazz3 )
5063	ENDIF
5064	IF( med_diag%O2SAT3%dgsave ) THEN
5065	CALL iom_put( "O2SAT3" , o2sat3 )
5066	CALL wrk_dealloc( jpi, jpj, jpk, o2sat3 )
5067	ENDIF
5068	IF( med_diag%PBSI3%dgsave ) THEN
5069	CALL iom_put( "PBSI3" , pbsi3 )
5070	CALL wrk_dealloc( jpi, jpj, jpk, pbsi3 )
5071	ENDIF
5072	IF( med_diag%PCAL3%dgsave ) THEN
5073	CALL iom_put( "PCAL3" , pcal3 )
5074	CALL wrk_dealloc( jpi, jpj, jpk, pcal3 )
5075	ENDIF
5076	IF( med_diag%REMOC3%dgsave ) THEN
5077	CALL iom_put( "REMOC3" , remoc3 )
5078	CALL wrk_dealloc( jpi, jpj, jpk, remoc3 )
5079	ENDIF
5080	IF( med_diag%PNLIMJ3%dgsave ) THEN
5081	CALL iom_put( "PNLIMJ3" , pnlimj3 )
5082	CALL wrk_dealloc( jpi, jpj, jpk, pnlimj3 )
5083	ENDIF
5084	IF( med_diag%PNLIMN3%dgsave ) THEN
5085	CALL iom_put( "PNLIMN3" , pnlimn3 )
5086	CALL wrk_dealloc( jpi, jpj, jpk, pnlimn3 )
5087	ENDIF
5088	IF( med_diag%PNLIMFE3%dgsave ) THEN
5089	CALL iom_put( "PNLIMFE3" , pnlimfe3 )
5090	CALL wrk_dealloc( jpi, jpj, jpk, pnlimfe3 )
5091	ENDIF
5092	IF( med_diag%PDLIMJ3%dgsave ) THEN
5093	CALL iom_put( "PDLIMJ3" , pdlimj3 )
5094	CALL wrk_dealloc( jpi, jpj, jpk, pdlimj3 )
5095	ENDIF
5096	IF( med_diag%PDLIMN3%dgsave ) THEN
5097	CALL iom_put( "PDLIMN3" , pdlimn3 )
5098	CALL wrk_dealloc( jpi, jpj, jpk, pdlimn3 )
5099	ENDIF
5100	IF( med_diag%PDLIMFE3%dgsave ) THEN
5101	CALL iom_put( "PDLIMFE3" , pdlimfe3 )
5102	CALL wrk_dealloc( jpi, jpj, jpk, pdlimfe3 )
5103	ENDIF
5104	IF( med_diag%PDLIMSI3%dgsave ) THEN
5105	CALL iom_put( "PDLIMSI3" , pdlimsi3 )
5106	CALL wrk_dealloc( jpi, jpj, jpk, pdlimsi3 )
5107	ENDIF
5108	# endif
5109
5110	CALL wrk_dealloc( jpi, jpj, zw2d )
5111
5112	# if defined key_debug_medusa
5113	IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
5114	CALL flush(numout)
5115	# endif
5116
5117	END SUBROUTINE trc_bio_medusa
5118
5119	#else
5120	!!======================================================================
5121	!! Dummy module : No MEDUSA bio-model
5122	!!======================================================================
5123	CONTAINS
5124	SUBROUTINE trc_bio_medusa( kt ) ! Empty routine
5125	INTEGER, INTENT( in ) :: kt
5126	WRITE(,) 'trc_bio_medusa: You should not have seen this print! error?', kt
5127	END SUBROUTINE trc_bio_medusa
5128	#endif
5129
5130	!!======================================================================
5131	END MODULE trcbio_medusa

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

Download in other formats: