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

Last change on this file since 8182 was 8182, checked in by jpalmier, 5 years ago

JPALM -- update MEDUSA xCO2 to CMIP6 value

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