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

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

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 7927

Last change on this file since 7927 was 7927, checked in by marc, 7 years ago

Moving the diagnostics at end of the big DO loop in trcbio_medusa.F90 into bio_medusa_diag_slice.F90

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