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

Last change on this file since 6466 was 6466, checked in by jpalmier, 6 years ago

JPALM — 11-04-2016 — add dust deposition input through namelist

— relax time to IDTRA surface flux

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