source: branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8224

Last change on this file since 8224 was 8224, checked in by frrh, 4 years ago

Merge Marc's fix for MEDUSA NRUN reproducibility since this is the cleanest
way to include this fix in testing and will be required as a permanent fix
anyway.

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