New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcbio_medusa.F90 in branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

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

Last change on this file since 5841 was 5841, checked in by jpalmier, 8 years ago

JPALM --30-10-2015-- Add MOCSY and DMS to MEDUSA-NEMO3.6

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