source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 5707

Last change on this file since 5707 was 5707, checked in by acc, 5 years ago

JPALM —25-08-2015 — add MEDUSA in the branch. MEDUSA version already up-to-date with this trunk revision

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