source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8457

Last change on this file since 8457 was 8457, checked in by dford, 3 years ago

Option to read in a value for atmospheric xCO2.

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