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

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

Marc 6/4/17. Fix bug where flatx is not set in call to mocsy_interface

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