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

Last change on this file since 7894 was 7894, checked in by jpalmier, 4 years ago

JPALM — 11-04-2017 — MEDUSA spring tidy-up refreshning session

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