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

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

JPALM — 14-12-2016 — MEDUSA Diags update ; remove extra lbc-lnk

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