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

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

JPALM — 17-11-2016 — CMIP6 diags for MEDUSA bugfix

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