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

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

JPALM —14-11-2015 — add MEDUSA CMIP6 diags-update

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