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

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

Implement initial version of surface chlorophyll assimilation for MEDUSA.

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