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

Last change on this file since 7703 was 7703, checked in by marc, 4 years ago

Marc 20/2/17. Fix two bugs.

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