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.
bio_medusa_diag.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/bio_medusa_diag.F90 @ 8023

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

Tidying up of headers for several files

File size: 83.6 KB
Line 
1MODULE bio_medusa_diag_mod
2   !!======================================================================
3   !!                         ***  MODULE bio_medusa_diag_mod  ***
4   !! Calculates diagnostics
5   !!======================================================================
6   !! History :
7   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
8   !!----------------------------------------------------------------------
9#if defined key_medusa
10   !!----------------------------------------------------------------------
11   !!                                                   MEDUSA bio-model
12   !!----------------------------------------------------------------------
13
14   IMPLICIT NONE
15   PRIVATE
16     
17   PUBLIC   bio_medusa_diag        ! Called in trcbio_medusa.F90
18
19   !!----------------------------------------------------------------------
20   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
21   !! $Id$
22   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
23   !!----------------------------------------------------------------------
24
25CONTAINS
26
27   SUBROUTINE bio_medusa_diag( kt, jk )
28      !!-------------------------------------------------------------------
29      !!                     ***  ROUTINE bio_medusa_diag  ***
30      !! This called from TRC_BIO_MEDUSA and calculates diagnostics
31      !!-------------------------------------------------------------------
32      USE bio_medusa_mod
33      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,       &
34                                   mbathy, tmask
35      USE in_out_manager,    ONLY: lwp, numout
36# if defined key_iomput
37      USE iom,               ONLY: lk_iomput
38# endif
39      USE par_kind,          ONLY: wp
40      USE par_oce,           ONLY: jpim1, jpjm1
41      USE phycst,            ONLY: rsmall
42      USE sbc_oce,           ONLY: qsr, wndm
43      USE sms_medusa,        ONLY: f2_ccd_arg, f2_ccd_cal,               &
44                                   f3_omarg, f3_omcal, f3_pH,            &
45                                   i0100, i0150, i0200, i0500, i1000,    &
46                                   jdms, ocal_ccd,                       &
47                                   xbetac, xbetan, xpar, xphi, xrfn,     &
48                                   xthetapd, xthetapn, xthetazme, xthetazmi, xze
49      USE trc,               ONLY: ln_diatrc, med_diag, trc2d, trc3d 
50# if defined key_roam
51      USE trcoxy_medusa,     ONLY: oxy_sato
52# endif
53
54   !!* Substitution
55#  include "domzgr_substitute.h90"
56
57      !! time (integer timestep)
58      INTEGER, INTENT( in ) :: kt
59      !! level
60      INTEGER, INTENT( in ) :: jk
61
62      !! Loop avariables
63      INTEGER :: ji, jj, jn
64
65# if defined key_trc_diabio
66      !!==========================================================
67      !! LOCAL GRID CELL DIAGNOSTICS
68      !!==========================================================
69      !!
70      !!----------------------------------------------------------
71      !! Full diagnostics key_trc_diabio:
72      !! LOBSTER and PISCES support full diagnistics option
73      !! key_trc_diabio which gives an option of FULL output of
74      !! biological sourses and sinks. I cannot see any reason
75      !! for doing this. If needed, it can be done as shown
76      !! below.
77      !!----------------------------------------------------------
78      !!
79      IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio'
80# endif
81
82      !! The section below, down to calculation of zo2min, was moved
83      !! from before the call to AIR_SEA in trcbio_medusa.F90 - marc 9/5/17
84      IF( lk_iomput ) THEN
85         DO jj = 2,jpjm1
86            DO ji = 2,jpim1
87               if (tmask(ji,jj,1) == 1) then
88                  !! sum tracers for inventory checks
89                  IF ( med_diag%INVTN%dgsave )   THEN
90                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        &
91                        (fse3t(ji,jj,jk) * (zphn(ji,jj) + zphd(ji,jj) +      &
92                                            zzmi(ji,jj) + zzme(ji,jj) +      &
93                                            zdet(ji,jj) + zdin(ji,jj)))
94                  ENDIF
95                  IF ( med_diag%INVTSI%dgsave )  THEN
96                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       & 
97                       (fse3t(ji,jj,jk) * (zpds(ji,jj) + zsil(ji,jj)))
98                  ENDIF
99                  IF ( med_diag%INVTFE%dgsave )  THEN
100                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       & 
101                        (fse3t(ji,jj,jk) * (xrfn *                           &
102                                            (zphn(ji,jj) + zphd(ji,jj) +     &
103                                             zzmi(ji,jj) + zzme(ji,jj) +     &
104                                             zdet(ji,jj)) +                  &
105                                            zfer(ji,jj)))
106                  ENDIF
107               ENDIF
108            ENDDO
109         ENDDO
110
111# if defined key_roam
112         DO jj = 2,jpjm1
113            DO ji = 2,jpim1
114               if (tmask(ji,jj,1) == 1) then
115                  IF ( med_diag%INVTC%dgsave )  THEN
116                     ftot_c(ji,jj)  = ftot_c(ji,jj) +                        & 
117                        (fse3t(ji,jj,jk) * ((xthetapn * zphn(ji,jj)) +       &
118                                            (xthetapd * zphd(ji,jj)) +       &
119                                            (xthetazmi * zzmi(ji,jj)) +      &
120                                            (xthetazme * zzme(ji,jj)) +      &
121                                            zdtc(ji,jj) + zdic(ji,jj)))
122                  ENDIF
123                  IF ( med_diag%INVTALK%dgsave ) THEN
124                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     &
125                                                       zalk(ji,jj))
126                  ENDIF
127                  IF ( med_diag%INVTO2%dgsave )  THEN
128                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    &
129                                                        zoxy(ji,jj))
130                  ENDIF
131               ENDIF
132            ENDDO
133         ENDDO
134
135         DO jj = 2,jpjm1
136            DO ji = 2,jpim1
137               if (tmask(ji,jj,1) == 1) then
138                  IF ( med_diag%INVTC%dgsave )  THEN
139                     !!
140                     !! AXY (10/11/16): CMIP6 diagnostics
141                     IF ( med_diag%INTDISSIC%dgsave ) THEN
142                        intdissic(ji,jj) = intdissic(ji,jj) +                &
143                                           (fse3t(ji,jj,jk) * zdic(ji,jj))
144                     ENDIF
145                     IF ( med_diag%INTDISSIN%dgsave ) THEN
146                        intdissin(ji,jj) = intdissin(ji,jj) +                &
147                                           (fse3t(ji,jj,jk) * zdin(ji,jj))
148                     ENDIF
149                     IF ( med_diag%INTDISSISI%dgsave ) THEN
150                        intdissisi(ji,jj) = intdissisi(ji,jj) +              &
151                                            (fse3t(ji,jj,jk) * zsil(ji,jj))
152                     ENDIF
153                     IF ( med_diag%INTTALK%dgsave ) THEN
154                        inttalk(ji,jj) = inttalk(ji,jj) +                    &
155                                         (fse3t(ji,jj,jk) * zalk(ji,jj))
156                     ENDIF
157                  ENDIF
158               ENDIF
159            ENDDO
160         ENDDO
161
162         DO jj = 2,jpjm1
163            DO ji = 2,jpim1
164               if (tmask(ji,jj,1) == 1) then
165                  IF ( med_diag%O2min%dgsave ) THEN
166                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then
167                        o2min(ji,jj)  = zoxy(ji,jj)
168                        IF ( med_diag%ZO2min%dgsave ) THEN
169                           !! layer midpoint
170                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               &
171                                            fdep1(ji,jj)) / 2.0
172                        ENDIF
173                     endif
174                  ENDIF
175               ENDIF
176            ENDDO
177         ENDDO
178# endif
179      ENDIF
180
181# if defined key_roam
182      !! This section is moved from just below CALL to AIR_SEA in
183      !! trcbio_medusa.F90 - marc 9/5/17
184      DO jj = 2,jpjm1
185         DO ji = 2,jpim1
186            !! OPEN wet point IF..THEN loop
187            if (tmask(ji,jj,jk) == 1) then
188
189               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
190               IF ( med_diag%O2SAT3%dgsave ) THEN
191! Remove f_o2sat3 - marc 9/5/17
192!                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 )
193!                  o2sat3(ji, jj, jk) = f_o2sat3
194                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj),                   &
195                                 o2sat3(ji,jj,jk) )
196               ENDIF
197            ENDIF
198         ENDDO
199      ENDDO
200# endif
201
202      IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN
203
204         DO jj = 2,jpjm1
205            DO ji = 2,jpim1
206               !! OPEN wet point IF..THEN loop
207               IF (tmask(ji,jj,jk) == 1) THEN
208                  !!-------------------------------------------------------
209                  !! Add in XML diagnostics stuff
210                  !!-------------------------------------------------------
211                  !!
212                  !! ** 2D diagnostics
213#   if defined key_debug_medusa
214                  IF (lwp) write (numout,*)                                  &
215                     'trc_bio_medusa: diag in ij-jj-jk loop'
216                  CALL flush(numout)
217#   endif
218                  IF ( med_diag%PRN%dgsave ) THEN
219                      fprn2d(ji,jj) = fprn2d(ji,jj) +                        &
220                                      (fprn(ji,jj)  * zphn(ji,jj) *          &
221                                       fse3t(ji,jj,jk)) 
222                  ENDIF
223                  IF ( med_diag%MPN%dgsave ) THEN
224                      fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn(ji,jj) *         &
225                                                       fse3t(ji,jj,jk))
226                  ENDIF
227                  IF ( med_diag%PRD%dgsave ) THEN
228                      fprd2d(ji,jj) = fprd2d(ji,jj) +                        &
229                                      (fprd(ji,jj)  * zphd(ji,jj) *          &
230                                       fse3t(ji,jj,jk))
231                  ENDIF
232                  IF( med_diag%MPD%dgsave ) THEN
233                      fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd(ji,jj) *         &
234                                                       fse3t(ji,jj,jk)) 
235                  ENDIF
236                  !  IF( med_diag%DSED%dgsave ) THEN
237                  !      CALL iom_put( "DSED"  , ftot_n )
238                  !  ENDIF
239                  IF( med_diag%OPAL%dgsave ) THEN
240                      fprds2d(ji,jj) = fprds2d(ji,jj) +                      &
241                                       (fprds(ji,jj) * zpds(ji,jj) *         &
242                                        fse3t(ji,jj,jk)) 
243                  ENDIF
244               ENDIF
245            ENDDO
246         ENDDO
247
248         DO jj = 2,jpjm1
249            DO ji = 2,jpim1
250               IF (tmask(ji,jj,jk) == 1) THEN
251                  IF( med_diag%OPALDISS%dgsave ) THEN
252                      fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss(ji,jj) *   &
253                                                           fse3t(ji,jj,jk)) 
254                  ENDIF
255                  IF( med_diag%GMIPn%dgsave ) THEN
256                      fgmipn2d(ji,jj) = fgmipn2d(ji,jj) +                    &
257                                        (fgmipn(ji,jj)  * fse3t(ji,jj,jk)) 
258                  ENDIF
259                  IF( med_diag%GMID%dgsave ) THEN
260                      fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid(ji,jj) *      &
261                                                         fse3t(ji,jj,jk)) 
262                  ENDIF
263                  IF( med_diag%MZMI%dgsave ) THEN
264                      fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi(ji,jj) *      &
265                                                         fse3t(ji,jj,jk)) 
266                  ENDIF
267               ENDIF
268            ENDDO
269         ENDDO
270
271         DO jj = 2,jpjm1
272            DO ji = 2,jpim1
273               IF (tmask(ji,jj,jk) == 1) THEN
274                  IF( med_diag%GMEPN%dgsave ) THEN
275                      fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn(ji,jj) *   &
276                                                           fse3t(ji,jj,jk))
277                  ENDIF
278                  IF( med_diag%GMEPD%dgsave ) THEN
279                      fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd(ji,jj) *   &
280                                                           fse3t(ji,jj,jk)) 
281                  ENDIF
282                  IF( med_diag%GMEZMI%dgsave ) THEN
283                      fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) +                  &
284                                         (fgmezmi(ji,jj) * fse3t(ji,jj,jk)) 
285                  ENDIF
286                  IF( med_diag%GMED%dgsave ) THEN
287                      fgmed2d(ji,jj) = fgmed2d(ji,jj) +                      &
288                                       (fgmed(ji,jj) * fse3t(ji,jj,jk)) 
289                  ENDIF
290                  IF( med_diag%MZME%dgsave ) THEN
291                      fdzme2d(ji,jj) = fdzme2d(ji,jj) +                      &
292                                       (fdzme(ji,jj) * fse3t(ji,jj,jk)) 
293                  ENDIF
294                  !  IF( med_diag%DEXP%dgsave ) THEN
295                  !      CALL iom_put( "DEXP"  , ftot_n )
296                  !  ENDIF
297               ENDIF
298            ENDDO
299         ENDDO
300
301         DO jj = 2,jpjm1
302            DO ji = 2,jpim1
303               IF (tmask(ji,jj,jk) == 1) THEN
304                  IF( med_diag%DETN%dgsave ) THEN
305                      fslown2d(ji,jj) = fslown2d(ji,jj) +                    &
306                                        (fslown(ji,jj) * fse3t(ji,jj,jk)) 
307                  ENDIF
308                  IF( med_diag%MDET%dgsave ) THEN
309                      fdd2d(ji,jj) = fdd2d(ji,jj) +                          &
310                                     (fdd(ji,jj) * fse3t(ji,jj,jk)) 
311                  ENDIF
312               ENDIF
313            ENDDO
314         ENDDO
315
316         DO jj = 2,jpjm1
317            DO ji = 2,jpim1
318               IF (tmask(ji,jj,jk) == 1) THEN
319                  IF( med_diag%AEOLIAN%dgsave ) THEN
320                      ffetop2d(ji,jj) = ffetop2d(ji,jj) +                    &
321                                        (ffetop(ji,jj) * fse3t(ji,jj,jk)) 
322                  ENDIF
323                  IF( med_diag%BENTHIC%dgsave ) THEN
324                      ffebot2d(ji,jj) = ffebot2d(ji,jj) +                    &
325                                        (ffebot(ji,jj) * fse3t(ji,jj,jk)) 
326                  ENDIF
327                  IF( med_diag%SCAVENGE%dgsave ) THEN
328                      ffescav2d(ji,jj) = ffescav2d(ji,jj) +                  &
329                                         (ffescav(ji,jj) * fse3t(ji,jj,jk)) 
330                  ENDIF
331               ENDIF
332            ENDDO
333         ENDDO
334
335         DO jj = 2,jpjm1
336            DO ji = 2,jpim1
337               IF (tmask(ji,jj,jk) == 1) THEN
338                  IF( med_diag%PN_JLIM%dgsave ) THEN
339                      ! fjln2d(ji,jj) = fjln2d(ji,jj) +                       &
340                      !                 (fjln(ji,jj)  * zphn(ji,jj) *         &
341                      !                  fse3t(ji,jj,jk))
342                      fjln2d(ji,jj) = fjln2d(ji,jj) +                        &
343                                      (fjlim_pn(ji,jj) * zphn(ji,jj) *       &
344                                       fse3t(ji,jj,jk)) 
345                  ENDIF
346                  IF( med_diag%PN_NLIM%dgsave ) THEN
347                      fnln2d(ji,jj) = fnln2d(ji,jj) +                        &
348                                      (fnln(ji,jj) * zphn(ji,jj) *           &
349                                       fse3t(ji,jj,jk)) 
350                  ENDIF
351                  IF( med_diag%PN_FELIM%dgsave ) THEN
352                      ffln2d(ji,jj) = ffln2d(ji,jj) +                        &
353                                      (ffln2(ji,jj) * zphn(ji,jj) *          &
354                                       fse3t(ji,jj,jk)) 
355                  ENDIF
356               ENDIF
357            ENDDO
358         ENDDO
359
360         DO jj = 2,jpjm1
361            DO ji = 2,jpim1
362               IF (tmask(ji,jj,jk) == 1) THEN
363                  IF( med_diag%PD_JLIM%dgsave ) THEN
364                      ! fjld2d(ji,jj) = fjld2d(ji,jj) +                       &
365                      !                 (fjld(ji,jj)  * zphd(ji,jj) *         &
366                      !                  fse3t(ji,jj,jk))
367                      fjld2d(ji,jj) = fjld2d(ji,jj) +                        &
368                                      (fjlim_pd(ji,jj) * zphd(ji,jj) *       &
369                                       fse3t(ji,jj,jk)) 
370                  ENDIF
371                  IF( med_diag%PD_NLIM%dgsave ) THEN
372                      fnld2d(ji,jj) = fnld2d(ji,jj) +                        &
373                                      (fnld(ji,jj) * zphd(ji,jj) *           &
374                                       fse3t(ji,jj,jk)) 
375                  ENDIF
376                  IF( med_diag%PD_FELIM%dgsave ) THEN
377                      ffld2d(ji,jj) = ffld2d(ji,jj) +                        &
378                                      (ffld(ji,jj) * zphd(ji,jj) *           &
379                                       fse3t(ji,jj,jk)) 
380                  ENDIF
381                  IF( med_diag%PD_SILIM%dgsave ) THEN
382                      fsld2d2(ji,jj) = fsld2d2(ji,jj) +                      &
383                                       (fsld2(ji,jj) * zphd(ji,jj) *         &
384                                        fse3t(ji,jj,jk)) 
385                  ENDIF
386                  IF( med_diag%PDSILIM2%dgsave ) THEN
387                      fsld2d(ji,jj) = fsld2d(ji,jj) +                        &
388                                      (fsld(ji,jj) * zphd(ji,jj) *           &
389                                       fse3t(ji,jj,jk))
390                  ENDIF
391               ENDIF
392            ENDDO
393         ENDDO
394
395         DO jj = 2,jpjm1
396            DO ji = 2,jpim1
397               IF (tmask(ji,jj,jk) == 1) THEN
398                  !!
399                  IF( med_diag%TOTREG_N%dgsave ) THEN
400                      fregen2d(ji,jj) = fregen2d(ji,jj) + fregen(ji,jj)
401                  ENDIF
402                  IF( med_diag%TOTRG_SI%dgsave ) THEN
403                      fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi(ji,jj)
404                  ENDIF
405               ENDIF
406            ENDDO
407         ENDDO
408
409         DO jj = 2,jpjm1
410            DO ji = 2,jpim1
411               IF (tmask(ji,jj,jk) == 1) THEN
412                  !!
413                  IF( med_diag%FASTN%dgsave ) THEN
414                      ftempn2d(ji,jj) = ftempn2d(ji,jj) +                    &
415                                        (ftempn(ji,jj)  * fse3t(ji,jj,jk))
416                  ENDIF
417                  IF( med_diag%FASTSI%dgsave ) THEN
418                      ftempsi2d(ji,jj) = ftempsi2d(ji,jj) +                  &
419                                         (ftempsi(ji,jj) * fse3t(ji,jj,jk))
420                  ENDIF
421                  IF( med_diag%FASTFE%dgsave ) THEN
422                      ftempfe2d(ji,jj) = ftempfe2d(ji,jj) +                  &
423                                         (ftempfe(ji,jj) * fse3t(ji,jj,jk)) 
424                  ENDIF
425                  IF( med_diag%FASTC%dgsave ) THEN
426                      ftempc2d(ji,jj) = ftempc2d(ji,jj) +                    &
427                                        (ftempc(ji,jj) * fse3t(ji,jj,jk))
428                  ENDIF
429                  IF( med_diag%FASTCA%dgsave ) THEN
430                      ftempca2d(ji,jj) = ftempca2d(ji,jj) +                  &
431                                         (ftempca(ji,jj) * fse3t(ji,jj,jk))
432                  ENDIF
433               ENDIF
434            ENDDO
435         ENDDO
436
437         DO jj = 2,jpjm1
438            DO ji = 2,jpim1
439               IF (tmask(ji,jj,jk) == 1) THEN
440                  !!
441                  IF( med_diag%REMINN%dgsave ) THEN
442                      freminn2d(ji,jj) = freminn2d(ji,jj) +                  &
443                                         (freminn(ji,jj)  * fse3t(ji,jj,jk))
444                  ENDIF
445                  IF( med_diag%REMINSI%dgsave ) THEN
446                      freminsi2d(ji,jj) = freminsi2d(ji,jj) +                &
447                                          (freminsi(ji,jj) * fse3t(ji,jj,jk))
448                  ENDIF
449                  IF( med_diag%REMINFE%dgsave ) THEN
450                      freminfe2d(ji,jj) = freminfe2d(ji,jj) +                &
451                                          (freminfe(ji,jj) * fse3t(ji,jj,jk)) 
452                  ENDIF
453                  IF( med_diag%REMINC%dgsave ) THEN
454                      freminc2d(ji,jj) = freminc2d(ji,jj) +                  &
455                                         (freminc(ji,jj)  * fse3t(ji,jj,jk)) 
456                  ENDIF
457                  IF( med_diag%REMINCA%dgsave ) THEN
458                      freminca2d(ji,jj) = freminca2d(ji,jj) +                &
459                                          (freminca(ji,jj) * fse3t(ji,jj,jk)) 
460                  ENDIF
461                  !!
462               ENDIF
463            ENDDO
464         ENDDO
465
466# if defined key_roam
467         DO jj = 2,jpjm1
468            DO ji = 2,jpim1
469               IF (tmask(ji,jj,jk) == 1) THEN
470                  !!
471                  !! AXY (09/11/16): CMIP6 diagnostics
472                  IF( med_diag%FD_NIT3%dgsave ) THEN
473                     fd_nit3(ji,jj,jk) = ffastn(ji,jj)
474                  ENDIF
475                  IF( med_diag%FD_SIL3%dgsave ) THEN
476                     fd_sil3(ji,jj,jk) = ffastsi(ji,jj)
477                  ENDIF
478                  IF( med_diag%FD_CAR3%dgsave ) THEN
479                     fd_car3(ji,jj,jk) = ffastc(ji,jj)
480                  ENDIF
481                  IF( med_diag%FD_CAL3%dgsave ) THEN
482                     fd_cal3(ji,jj,jk) = ffastca(ji,jj)
483                  ENDIF
484               ENDIF
485            ENDDO
486         ENDDO
487
488         IF (jk.eq.i0100) THEN
489            DO jj = 2,jpjm1
490               DO ji = 2,jpim1
491                  IF (tmask(ji,jj,jk) == 1) THEN
492                     IF( med_diag%RR_0100%dgsave ) THEN
493                        ffastca2d(ji,jj) =                                   &
494                           ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
495                     ENDIF
496                  ENDIF
497               ENDDO
498            ENDDO
499         ELSE IF (jk.eq.i0500) THEN
500            DO jj = 2,jpjm1
501               DO ji = 2,jpim1
502                  IF (tmask(ji,jj,jk) == 1) THEN
503                     IF( med_diag%RR_0500%dgsave ) THEN
504                        ffastca2d(ji,jj) =                                   &
505                           ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
506                     ENDIF
507                  ENDIF
508               ENDDO
509            ENDDO
510         ELSE IF (jk.eq.i1000) THEN
511            DO jj = 2,jpjm1
512               DO ji = 2,jpim1
513                  IF (tmask(ji,jj,jk) == 1) THEN
514                     IF( med_diag%RR_1000%dgsave ) THEN
515                        ffastca2d(ji,jj) =                                   &
516                           ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
517                     ENDIF
518                  ENDIF
519               ENDDO
520            ENDDO
521         ELSE IF (jk.eq.mbathy(ji,jj)) THEN
522            DO jj = 2,jpjm1
523               DO ji = 2,jpim1
524                  IF (tmask(ji,jj,jk) == 1) THEN
525                     IF( med_diag%IBEN_N%dgsave ) THEN
526                        iben_n2d(ji,jj) = f_sbenin_n(ji,jj) +                &
527                                          f_fbenin_n(ji,jj)
528                     ENDIF
529                     IF( med_diag%IBEN_FE%dgsave ) THEN
530                        iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) +              &
531                                           f_fbenin_fe(ji,jj)
532                     ENDIF
533                     IF( med_diag%IBEN_C%dgsave ) THEN
534                        iben_c2d(ji,jj) = f_sbenin_c(ji,jj) +                &
535                                          f_fbenin_c(ji,jj)
536                     ENDIF
537                     IF( med_diag%IBEN_SI%dgsave ) THEN
538                        iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
539                     ENDIF
540                     IF( med_diag%IBEN_CA%dgsave ) THEN
541                        iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
542                     ENDIF
543                     IF( med_diag%OBEN_N%dgsave ) THEN
544                        oben_n2d(ji,jj) = f_benout_n(ji,jj)
545                     ENDIF
546                     IF( med_diag%OBEN_FE%dgsave ) THEN
547                        oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
548                     ENDIF
549                     IF( med_diag%OBEN_C%dgsave ) THEN
550                        oben_c2d(ji,jj) = f_benout_c(ji,jj)
551                     ENDIF
552                     IF( med_diag%OBEN_SI%dgsave ) THEN
553                        oben_si2d(ji,jj) = f_benout_si(ji,jj)
554                     ENDIF
555                     IF( med_diag%OBEN_CA%dgsave ) THEN
556                        oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
557                     ENDIF
558                     IF( med_diag%SFR_OCAL%dgsave ) THEN
559                        sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
560                     ENDIF
561                     IF( med_diag%SFR_OARG%dgsave ) THEN
562                        sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
563                     ENDIF
564                     IF( med_diag%LYSO_CA%dgsave ) THEN
565                        lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
566                     ENDIF
567                  ENDIF
568               ENDDO
569            ENDDO
570         ENDIF
571         !! end bathy-1 diags
572
573         DO jj = 2,jpjm1
574            DO ji = 2,jpim1
575               IF (tmask(ji,jj,jk) == 1) THEN
576                  !!
577                  IF( med_diag%RIV_N%dgsave ) THEN
578                     rivn2d(ji,jj) = rivn2d(ji,jj) +                         &
579                                     (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
580                  ENDIF
581                  IF( med_diag%RIV_SI%dgsave ) THEN
582                     rivsi2d(ji,jj) = rivsi2d(ji,jj) +                       &
583                                      (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
584                  ENDIF
585                  IF( med_diag%RIV_C%dgsave ) THEN
586                     rivc2d(ji,jj) = rivc2d(ji,jj) +                         &
587                                     (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
588                  ENDIF
589                  IF( med_diag%RIV_ALK%dgsave ) THEN
590                     rivalk2d(ji,jj) = rivalk2d(ji,jj) +                     &
591                                       (f_riv_loc_alk(ji,jj) *               &
592                                        fse3t(ji,jj,jk))
593                  ENDIF
594                  IF( med_diag%DETC%dgsave ) THEN
595                     fslowc2d(ji,jj) = fslowc2d(ji,jj) +                     &
596                                       (fslowc(ji,jj)  * fse3t(ji,jj,jk))   
597                  ENDIF
598               ENDIF
599            ENDDO
600         ENDDO
601
602         DO jj = 2,jpjm1
603            DO ji = 2,jpim1
604               IF (tmask(ji,jj,jk) == 1) THEN
605                  !!
606                  IF( med_diag%PN_LLOSS%dgsave ) THEN
607                     fdpn22d(ji,jj) = fdpn22d(ji,jj) +                       &
608                                      (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
609                  ENDIF
610                  IF( med_diag%PD_LLOSS%dgsave ) THEN
611                     fdpd22d(ji,jj) = fdpd22d(ji,jj) +                       &
612                                      (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
613                  ENDIF
614               ENDIF
615            ENDDO
616         ENDDO
617
618         DO jj = 2,jpjm1
619            DO ji = 2,jpim1
620               IF (tmask(ji,jj,jk) == 1) THEN
621                  IF( med_diag%ZI_LLOSS%dgsave ) THEN
622                     fdzmi22d(ji,jj) = fdzmi22d(ji,jj) +                     &
623                                       (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
624                  ENDIF
625                  IF( med_diag%ZE_LLOSS%dgsave ) THEN
626                     fdzme22d(ji,jj) = fdzme22d(ji,jj) +                     &
627                                       (fdzme2(ji,jj) * fse3t(ji,jj,jk))
628                  ENDIF
629               ENDIF
630            ENDDO
631         ENDDO
632
633         DO jj = 2,jpjm1
634            DO ji = 2,jpim1
635               IF (tmask(ji,jj,jk) == 1) THEN
636                  IF( med_diag%ZI_MES_N%dgsave ) THEN
637                     zimesn2d(ji,jj) = zimesn2d(ji,jj) +                     &
638                                       (xphi * (fgmipn(ji,jj) +              &
639                                                fgmid(ji,jj)) *              &
640                                        fse3t(ji,jj,jk))
641                  ENDIF
642                  IF( med_diag%ZI_MES_D%dgsave ) THEN
643                     zimesd2d(ji,jj) = zimesd2d(ji,jj) +                     & 
644                                       ((1. - xbetan) * finmi(ji,jj) *       &
645                                        fse3t(ji,jj,jk))
646                  ENDIF
647                  IF( med_diag%ZI_MES_C%dgsave ) THEN
648                     zimesc2d(ji,jj) = zimesc2d(ji,jj) +                     &
649                                       (xphi * ((xthetapn * fgmipn(ji,jj)) + &
650                                                fgmidc(ji,jj)) *             &
651                                                fse3t(ji,jj,jk))
652                  ENDIF
653                  IF( med_diag%ZI_MESDC%dgsave ) THEN
654                     zimesdc2d(ji,jj) = zimesdc2d(ji,jj) +                   &
655                                        ((1. - xbetac) * ficmi(ji,jj) *      &
656                                         fse3t(ji,jj,jk))
657                  ENDIF
658                  IF( med_diag%ZI_EXCR%dgsave ) THEN
659                     ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +                     &
660                                       (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
661                  ENDIF
662                  IF( med_diag%ZI_RESP%dgsave ) THEN
663                     ziresp2d(ji,jj) = ziresp2d(ji,jj) +                     &
664                                       (fmiresp(ji,jj) * fse3t(ji,jj,jk))
665                  ENDIF
666                  IF( med_diag%ZI_GROW%dgsave ) THEN
667                     zigrow2d(ji,jj) = zigrow2d(ji,jj) +                     &
668                                       (fmigrow(ji,jj) * fse3t(ji,jj,jk))
669                  ENDIF
670               ENDIF
671            ENDDO
672         ENDDO
673
674         DO jj = 2,jpjm1
675            DO ji = 2,jpim1
676               IF (tmask(ji,jj,jk) == 1) THEN
677                  IF( med_diag%ZE_MES_N%dgsave ) THEN
678                     zemesn2d(ji,jj) = zemesn2d(ji,jj) +                     &
679                                       (xphi *                               &
680                                        (fgmepn(ji,jj) + fgmepd(ji,jj) +     &
681                                         fgmezmi(ji,jj) + fgmed(ji,jj)) *    &
682                                        fse3t(ji,jj,jk))
683                  ENDIF
684                  IF( med_diag%ZE_MES_D%dgsave ) THEN
685                     zemesd2d(ji,jj) = zemesd2d(ji,jj) +                     &
686                                       ((1. - xbetan) * finme(ji,jj) *       &
687                                        fse3t(ji,jj,jk))
688                  ENDIF
689                  IF( med_diag%ZE_MES_C%dgsave ) THEN
690                     zemesc2d(ji,jj) = zemesc2d(ji,jj) +                     & 
691                                       (xphi *                               &
692                                        ((xthetapn * fgmepn(ji,jj)) +        &
693                                         (xthetapd * fgmepd(ji,jj)) +        &
694                                         (xthetazmi * fgmezmi(ji,jj)) +      &
695                                         fgmedc(ji,jj)) * fse3t(ji,jj,jk))
696                  ENDIF
697                  IF( med_diag%ZE_MESDC%dgsave ) THEN
698                     zemesdc2d(ji,jj) = zemesdc2d(ji,jj) +                   &
699                                        ((1. - xbetac) * ficme(ji,jj) *      &
700                                         fse3t(ji,jj,jk))
701                  ENDIF
702                  IF( med_diag%ZE_EXCR%dgsave ) THEN
703                     zeexcr2d(ji,jj) = zeexcr2d(ji,jj) +                     &
704                                       (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
705                  ENDIF
706                  IF( med_diag%ZE_RESP%dgsave ) THEN
707                     zeresp2d(ji,jj) = zeresp2d(ji,jj) +                     &
708                                       (fmeresp(ji,jj) * fse3t(ji,jj,jk))
709                  ENDIF
710                  IF( med_diag%ZE_GROW%dgsave ) THEN
711                     zegrow2d(ji,jj) = zegrow2d(ji,jj) +                     &
712                                       (fmegrow(ji,jj) * fse3t(ji,jj,jk))
713                  ENDIF
714               ENDIF
715            ENDDO
716         ENDDO
717
718         DO jj = 2,jpjm1
719            DO ji = 2,jpim1
720               IF (tmask(ji,jj,jk) == 1) THEN
721                  IF( med_diag%MDETC%dgsave ) THEN
722                     mdetc2d(ji,jj) = mdetc2d(ji,jj) +                       &
723                                      (fddc(ji,jj) * fse3t(ji,jj,jk))
724                  ENDIF
725                  IF( med_diag%GMIDC%dgsave ) THEN
726                     gmidc2d(ji,jj) = gmidc2d(ji,jj) +                       &
727                                      (fgmidc(ji,jj) * fse3t(ji,jj,jk))
728                  ENDIF
729                  IF( med_diag%GMEDC%dgsave ) THEN
730                     gmedc2d(ji,jj) = gmedc2d(ji,jj) +                       &
731                                      (fgmedc(ji,jj) * fse3t(ji,jj,jk))
732                  ENDIF
733               ENDIF
734            ENDDO
735         ENDDO
736# endif                   
737
738         DO jj = 2,jpjm1
739            DO ji = 2,jpim1
740               IF (tmask(ji,jj,jk) == 1) THEN
741                  !!
742                  !! ** 3D diagnostics
743                  IF( med_diag%TPP3%dgsave ) THEN
744                     tpp3d(ji,jj,jk) = (fprn(ji,jj) * zphn(ji,jj)) +         &
745                                       (fprd(ji,jj) * zphd(ji,jj))
746                     !CALL iom_put( "TPP3"  , tpp3d )
747                  ENDIF
748                  IF( med_diag%TPPD3%dgsave ) THEN
749                     tppd3(ji,jj,jk) = (fprd(ji,jj) * zphd(ji,jj))
750                  ENDIF
751               ENDIF
752            ENDDO
753         ENDDO
754
755         DO jj = 2,jpjm1
756            DO ji = 2,jpim1
757               IF (tmask(ji,jj,jk) == 1) THEN
758                 
759                  IF( med_diag%REMIN3N%dgsave ) THEN
760                     !! remineralisation
761                     remin3dn(ji,jj,jk) = fregen(ji,jj) +                    &
762                                          (freminn(ji,jj) * fse3t(ji,jj,jk))
763                     !CALL iom_put( "REMIN3N"  , remin3dn )
764                  ENDIF
765                  !! IF( med_diag%PH3%dgsave ) THEN
766                  !!     CALL iom_put( "PH3"  , f3_pH )
767                  !! ENDIF
768                  !! IF( med_diag%OM_CAL3%dgsave ) THEN
769                  !!     CALL iom_put( "OM_CAL3"  , f3_omcal )
770                  !! ENDIF
771        !!
772        !! AXY (09/11/16): CMIP6 diagnostics
773        IF ( med_diag%DCALC3%dgsave   ) THEN
774                     dcalc3(ji,jj,jk) = freminca(ji,jj)
775                  ENDIF
776               ENDIF
777            ENDDO
778         ENDDO
779
780         DO jj = 2,jpjm1
781            DO ji = 2,jpim1
782               IF (tmask(ji,jj,jk) == 1) THEN
783        IF ( med_diag%FEDISS3%dgsave  ) THEN
784                     fediss3(ji,jj,jk) = ffetop(ji,jj)
785                  ENDIF
786        IF ( med_diag%FESCAV3%dgsave  ) THEN
787                     fescav3(ji,jj,jk) = ffescav(ji,jj)
788                  ENDIF
789               ENDIF
790            ENDDO
791         ENDDO
792
793         DO jj = 2,jpjm1
794            DO ji = 2,jpim1
795               IF (tmask(ji,jj,jk) == 1) THEN
796        IF ( med_diag%MIGRAZP3%dgsave ) THEN
797                     migrazp3(ji,jj,jk) = fgmipn(ji,jj) * xthetapn
798                  ENDIF
799        IF ( med_diag%MIGRAZD3%dgsave ) THEN
800                     migrazd3(ji,jj,jk) = fgmidc(ji,jj)
801                  ENDIF
802        IF ( med_diag%MEGRAZP3%dgsave ) THEN
803                     megrazp3(ji,jj,jk) = (fgmepn(ji,jj) * xthetapn) +       &
804                                          (fgmepd(ji,jj) * xthetapd)
805                  ENDIF
806        IF ( med_diag%MEGRAZD3%dgsave ) THEN
807                     megrazd3(ji,jj,jk) = fgmedc(ji,jj)
808                  ENDIF
809        IF ( med_diag%MEGRAZZ3%dgsave ) THEN
810                     megrazz3(ji,jj,jk) = (fgmezmi(ji,jj) * xthetazmi)
811                  ENDIF
812               ENDIF
813            ENDDO
814         ENDDO
815
816         DO jj = 2,jpjm1
817            DO ji = 2,jpim1
818               IF (tmask(ji,jj,jk) == 1) THEN
819        IF ( med_diag%PBSI3%dgsave    ) THEN
820                     pbsi3(ji,jj,jk)    = (fprds(ji,jj) * zpds(ji,jj))
821                  ENDIF
822        IF ( med_diag%PCAL3%dgsave    ) THEN
823                     pcal3(ji,jj,jk)    = ftempca(ji,jj)
824                  ENDIF
825        IF ( med_diag%REMOC3%dgsave   ) THEN
826                     remoc3(ji,jj,jk)   = freminc(ji,jj)
827                  ENDIF
828               ENDIF
829            ENDDO
830         ENDDO
831
832         DO jj = 2,jpjm1
833            DO ji = 2,jpim1
834               IF (tmask(ji,jj,jk) == 1) THEN
835        IF ( med_diag%PNLIMJ3%dgsave  ) THEN
836                     ! pnlimj3(ji,jj,jk)  = fjln(ji,jj)
837                     pnlimj3(ji,jj,jk)  = fjlim_pn(ji,jj)
838                  ENDIF
839        IF ( med_diag%PNLIMN3%dgsave  ) THEN
840                     pnlimn3(ji,jj,jk)  = fnln(ji,jj)
841                  ENDIF
842        IF ( med_diag%PNLIMFE3%dgsave ) THEN
843                     pnlimfe3(ji,jj,jk) = ffln2(ji,jj)
844                  ENDIF
845        IF ( med_diag%PDLIMJ3%dgsave  ) THEN
846                     ! pdlimj3(ji,jj,jk)  = fjld(ji,jj)
847                     pdlimj3(ji,jj,jk)  = fjlim_pd(ji,jj)
848                  ENDIF
849        IF ( med_diag%PDLIMN3%dgsave  ) THEN
850                     pdlimn3(ji,jj,jk)  = fnld(ji,jj)
851                  ENDIF
852        IF ( med_diag%PDLIMFE3%dgsave ) THEN
853                     pdlimfe3(ji,jj,jk) = ffld(ji,jj)
854                  ENDIF
855        IF ( med_diag%PDLIMSI3%dgsave ) THEN
856                     pdlimsi3(ji,jj,jk) = fsld2(ji,jj)
857                  ENDIF
858               ENDIF
859            ENDDO
860         ENDDO
861
862      ELSE IF( ln_diatrc ) THEN
863
864         !!
865         !! ** Without using iom_use
866#   if defined key_debug_medusa
867         IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ij-jj-jk ln_diatrc'
868         CALL flush(numout)
869#   endif
870         DO jj = 2,jpjm1
871            DO ji = 2,jpim1
872               IF (tmask(ji,jj,jk) == 1) then
873                  !!-------------------------------------------------------
874                  !! Prepare 2D diagnostics
875                  !!-------------------------------------------------------
876                  !!
877                  !! if ((kt / 240*240).eq.kt) then
878                  !!    IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt
879                  !! endif     
880                  !! nitrogen inventory
881                  trc2d(ji,jj,1)  =  ftot_n(ji,jj)
882                  !! silicon  inventory
883                  trc2d(ji,jj,2)  =  ftot_si(ji,jj)
884                  !! iron     inventory
885                  trc2d(ji,jj,3)  =  ftot_fe(ji,jj)
886               ENDIF
887            ENDDO
888         ENDDO
889
890         DO jj = 2,jpjm1
891            DO ji = 2,jpim1
892               IF (tmask(ji,jj,jk) == 1) THEN
893                  !! non-diatom production
894                  trc2d(ji,jj,4)  = trc2d(ji,jj,4)  +                        &
895                                    (fprn(ji,jj)  * zphn(ji,jj) *            &
896                                     fse3t(ji,jj,jk))
897                  !! non-diatom non-grazing losses
898                  trc2d(ji,jj,5)  = trc2d(ji,jj,5)  +                        &
899                                    (fdpn(ji,jj) * fse3t(ji,jj,jk))
900                  !! diatom production
901                  trc2d(ji,jj,6)  = trc2d(ji,jj,6)  +                        &
902                                    (fprd(ji,jj) * zphd(ji,jj) *             &
903                                     fse3t(ji,jj,jk))
904                  !! diatom non-grazing losses
905                  !! diagnostic field  8 is (ostensibly) supplied by trcsed.F
906                  trc2d(ji,jj,7)  = trc2d(ji,jj,7)  +                        &
907                                    (fdpd(ji,jj) * fse3t(ji,jj,jk))
908                  !! diatom silicon production
909                  trc2d(ji,jj,9)  = trc2d(ji,jj,9)  +                        &
910                                    (fprds(ji,jj) * zpds(ji,jj) *            &
911                                     fse3t(ji,jj,jk))
912                  !! diatom silicon dissolution
913                  trc2d(ji,jj,10) = trc2d(ji,jj,10) +                        &
914                                    (fsdiss(ji,jj)  * fse3t(ji,jj,jk))
915               ENDIF
916            ENDDO
917         ENDDO
918
919         DO jj = 2,jpjm1
920            DO ji = 2,jpim1
921               IF (tmask(ji,jj,jk) == 1) THEN
922                  !! microzoo grazing on non-diatoms
923                  trc2d(ji,jj,11) = trc2d(ji,jj,11) +                        &
924                                    (fgmipn(ji,jj)  * fse3t(ji,jj,jk))
925                  !! microzoo grazing on detrital nitrogen
926                  trc2d(ji,jj,12) = trc2d(ji,jj,12) +                        &
927                                    (fgmid(ji,jj) * fse3t(ji,jj,jk))
928                  !! microzoo non-grazing losses
929                  trc2d(ji,jj,13) = trc2d(ji,jj,13) +                        &
930                                    (fdzmi(ji,jj) * fse3t(ji,jj,jk))
931               ENDIF
932            ENDDO
933         ENDDO
934
935         DO jj = 2,jpjm1
936            DO ji = 2,jpim1
937               IF (tmask(ji,jj,jk) == 1) THEN
938                  !! mesozoo  grazing on non-diatoms
939                  trc2d(ji,jj,14) = trc2d(ji,jj,14) +                        &
940                                    (fgmepn(ji,jj)  * fse3t(ji,jj,jk))
941                  !! mesozoo  grazing on diatoms
942                  trc2d(ji,jj,15) = trc2d(ji,jj,15) +                        &
943                                    (fgmepd(ji,jj)  * fse3t(ji,jj,jk))
944                  !! mesozoo  grazing on microzoo
945                  trc2d(ji,jj,16) = trc2d(ji,jj,16) +                        &
946                                    (fgmezmi(ji,jj) * fse3t(ji,jj,jk))
947                  !! mesozoo  grazing on detrital nitrogen
948                  trc2d(ji,jj,17) = trc2d(ji,jj,17) +                        &
949                                    (fgmed(ji,jj) * fse3t(ji,jj,jk))
950                  !! mesozoo  non-grazing losses
951                  trc2d(ji,jj,18) = trc2d(ji,jj,18) +                        &
952                                    (fdzme(ji,jj)   * fse3t(ji,jj,jk))
953               ENDIF
954            ENDDO
955         ENDDO
956
957         DO jj = 2,jpjm1
958            DO ji = 2,jpim1
959               IF (tmask(ji,jj,jk) == 1) THEN
960                  !! diagnostic field 19 is (ostensibly) supplied by trcexp.F
961                  !! slow sinking detritus N production
962                  trc2d(ji,jj,20) = trc2d(ji,jj,20) +                        &
963                                    (fslown(ji,jj) * fse3t(ji,jj,jk))
964                  !! detrital remineralisation
965                  trc2d(ji,jj,21) = trc2d(ji,jj,21) +                        &
966                                    (fdd(ji,jj) * fse3t(ji,jj,jk))
967                  !! aeolian  iron addition
968                  trc2d(ji,jj,22) = trc2d(ji,jj,22) +                        &
969                                    (ffetop(ji,jj) * fse3t(ji,jj,jk))
970                  !! seafloor iron addition
971                  trc2d(ji,jj,23) = trc2d(ji,jj,23) +                        &
972                                    (ffebot(ji,jj) * fse3t(ji,jj,jk))
973                  !! "free" iron scavenging
974                  trc2d(ji,jj,24) = trc2d(ji,jj,24) +                        &
975                                    (ffescav(ji,jj) * fse3t(ji,jj,jk))
976               ENDIF
977            ENDDO
978         ENDDO
979
980         DO jj = 2,jpjm1
981            DO ji = 2,jpim1
982               IF (tmask(ji,jj,jk) == 1) THEN
983                  !! non-diatom J  limitation term
984                  trc2d(ji,jj,25) = trc2d(ji,jj,25) +                        &
985                                    (fjlim_pn(ji,jj) * zphn(ji,jj) *         &
986                                     fse3t(ji,jj,jk))
987                  !! non-diatom N  limitation term
988                  trc2d(ji,jj,26) = trc2d(ji,jj,26) +                        &
989                                    (fnln(ji,jj) * zphn(ji,jj) *             &
990                                     fse3t(ji,jj,jk))
991                  !! non-diatom Fe limitation term
992                  trc2d(ji,jj,27) = trc2d(ji,jj,27) +                        &
993                                    (ffln2(ji,jj) * zphn(ji,jj) *            &
994                                     fse3t(ji,jj,jk))
995                  !! diatom     J  limitation term
996                  trc2d(ji,jj,28) = trc2d(ji,jj,28) +                        &
997                                    (fjlim_pd(ji,jj) * zphd(ji,jj) *         &
998                                     fse3t(ji,jj,jk))
999                  !! diatom     N  limitation term
1000                  trc2d(ji,jj,29) = trc2d(ji,jj,29) +                        &
1001                                    (fnld(ji,jj) * zphd(ji,jj) *             &
1002                                     fse3t(ji,jj,jk))
1003                  !! diatom     Fe limitation term
1004                  trc2d(ji,jj,30) = trc2d(ji,jj,30) +                        &
1005                                    (ffld(ji,jj) * zphd(ji,jj) *             &
1006                                     fse3t(ji,jj,jk))
1007                  !! diatom     Si limitation term
1008                  trc2d(ji,jj,31) = trc2d(ji,jj,31) +                        &
1009                                    (fsld2(ji,jj) * zphd(ji,jj) *            &
1010                                     fse3t(ji,jj,jk))
1011                  !! diatom     Si uptake limitation term
1012                  trc2d(ji,jj,32) = trc2d(ji,jj,32) +                        &
1013                                    (fsld(ji,jj) * zphd(ji,jj) *             &
1014                                     fse3t(ji,jj,jk))
1015               ENDIF
1016            ENDDO
1017         ENDDO
1018
1019         IF (jk.eq.i0100) THEN
1020            DO jj = 2,jpjm1
1021               DO ji = 2,jpim1
1022                  IF (tmask(ji,jj,jk) == 1) THEN
1023                     !! slow detritus flux at  100 m
1024                     trc2d(ji,jj,33) = fslownflux(ji,jj)
1025                  ENDIF
1026               ENDDO
1027            ENDDO
1028         ENDIF
1029
1030         IF (jk.eq.i0200) THEN
1031            DO jj = 2,jpjm1
1032               DO ji = 2,jpim1
1033                  IF (tmask(ji,jj,jk) == 1) THEN
1034                     !! slow detritus flux at  200 m
1035                     trc2d(ji,jj,34) = fslownflux(ji,jj)
1036                  ENDIF
1037               ENDDO
1038            ENDDO
1039         ENDIF
1040
1041         IF (jk.eq.i0500) THEN
1042            DO jj = 2,jpjm1
1043               DO ji = 2,jpim1
1044                  IF (tmask(ji,jj,jk) == 1) THEN
1045                     !! slow detritus flux at  500 m
1046                     trc2d(ji,jj,35) = fslownflux(ji,jj)
1047                  ENDIF
1048               ENDDO
1049            ENDDO
1050         ENDIF
1051
1052         IF (jk.eq.i1000) THEN
1053            DO jj = 2,jpjm1
1054               DO ji = 2,jpim1
1055                  IF (tmask(ji,jj,jk) == 1) THEN
1056                     !! slow detritus flux at 1000 m
1057                     trc2d(ji,jj,36) = fslownflux(ji,jj)
1058                  ENDIF
1059               ENDDO
1060            ENDDO
1061         ENDIF
1062
1063         DO jj = 2,jpjm1
1064            DO ji = 2,jpim1
1065               IF (tmask(ji,jj,jk) == 1) THEN
1066                  !! non-fast N  full column regeneration
1067                  trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen(ji,jj)
1068                  !! non-fast Si full column regeneration
1069                  trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi(ji,jj)
1070                  !! non-fast N  regeneration to  100 m
1071               ENDIF
1072            ENDDO
1073         ENDDO
1074
1075         IF (jk.eq.i0100) THEN
1076            DO jj = 2,jpjm1
1077               DO ji = 2,jpim1
1078                  IF (tmask(ji,jj,jk) == 1) THEN
1079                     trc2d(ji,jj,39) = trc2d(ji,jj,37)
1080                  ENDIF
1081               ENDDO
1082            ENDDO
1083         ENDIF
1084
1085         IF (jk.eq.i0200) THEN
1086            DO jj = 2,jpjm1
1087               DO ji = 2,jpim1
1088                  IF (tmask(ji,jj,jk) == 1) THEN
1089                     !! non-fast N  regeneration to  200 m
1090                     trc2d(ji,jj,40) = trc2d(ji,jj,37)
1091                  ENDIF
1092               ENDDO
1093            ENDDO
1094         ENDIF
1095
1096         IF (jk.eq.i0500) THEN
1097            DO jj = 2,jpjm1
1098               DO ji = 2,jpim1
1099                  IF (tmask(ji,jj,jk) == 1) THEN
1100                     !! non-fast N  regeneration to  500 m
1101                     trc2d(ji,jj,41) = trc2d(ji,jj,37)
1102                  ENDIF
1103               ENDDO
1104            ENDDO
1105         ENDIF
1106
1107         IF (jk.eq.i1000) THEN
1108            DO jj = 2,jpjm1
1109               DO ji = 2,jpim1
1110                  IF (tmask(ji,jj,jk) == 1) THEN
1111                     !! non-fast N  regeneration to 1000 m
1112                     trc2d(ji,jj,42) = trc2d(ji,jj,37)
1113                  ENDIF
1114               ENDDO
1115            ENDDO
1116         ENDIF
1117
1118         DO jj = 2,jpjm1
1119            DO ji = 2,jpim1
1120               IF (tmask(ji,jj,jk) == 1) THEN
1121                  !! fast sinking detritus N production
1122                  trc2d(ji,jj,43) = trc2d(ji,jj,43) +                        &
1123                                    (ftempn(ji,jj) * fse3t(ji,jj,jk))
1124                  !! fast sinking detritus Si production
1125                  trc2d(ji,jj,44) = trc2d(ji,jj,44) +                        &
1126                                    (ftempsi(ji,jj) * fse3t(ji,jj,jk))
1127                  !! fast sinking detritus Fe production
1128                  trc2d(ji,jj,45) = trc2d(ji,jj,45) +                        &
1129                                    (ftempfe(ji,jj) * fse3t(ji,jj,jk))
1130                  !! fast sinking detritus C production
1131                  trc2d(ji,jj,46) = trc2d(ji,jj,46) +                        &
1132                                    (ftempc(ji,jj)  * fse3t(ji,jj,jk))
1133                  !! fast sinking detritus CaCO3 production
1134                  trc2d(ji,jj,47) = trc2d(ji,jj,47) +                        &
1135                                    (ftempca(ji,jj) * fse3t(ji,jj,jk))
1136               ENDIF
1137            ENDDO
1138         ENDDO
1139
1140         IF (jk.eq.i0100) THEN
1141            DO jj = 2,jpjm1
1142               DO ji = 2,jpim1
1143                  IF (tmask(ji,jj,jk) == 1) THEN
1144                     !! fast detritus N  flux at  100 m
1145                     trc2d(ji,jj,48) = ffastn(ji,jj)
1146                  ENDIF
1147               ENDDO
1148            ENDDO
1149         ENDIF
1150
1151         IF (jk.eq.i0200) THEN
1152            DO jj = 2,jpjm1
1153               DO ji = 2,jpim1
1154                  IF (tmask(ji,jj,jk) == 1) THEN
1155                     !! fast detritus N  flux at  200 m
1156                     trc2d(ji,jj,49) = ffastn(ji,jj)
1157                  ENDIF
1158               ENDDO
1159            ENDDO
1160         ENDIF
1161
1162         IF (jk.eq.i0500) THEN
1163            DO jj = 2,jpjm1
1164               DO ji = 2,jpim1
1165                  IF (tmask(ji,jj,jk) == 1) THEN
1166                     !! fast detritus N  flux at  500 m
1167                     trc2d(ji,jj,50) = ffastn(ji,jj)
1168                  ENDIF
1169               ENDDO
1170            ENDDO
1171         ENDIF
1172
1173         IF (jk.eq.i1000) THEN
1174            DO jj = 2,jpjm1
1175               DO ji = 2,jpim1
1176                  IF (tmask(ji,jj,jk) == 1) THEN
1177                     !! fast detritus N  flux at 1000 m
1178                     trc2d(ji,jj,51) = ffastn(ji,jj)
1179                  ENDIF
1180               ENDDO
1181            ENDDO
1182         ENDIF
1183
1184         IF (jk.eq.i0100) THEN
1185            DO jj = 2,jpjm1
1186               DO ji = 2,jpim1
1187                  IF (tmask(ji,jj,jk) == 1) THEN
1188                     !! N  regeneration to  100 m
1189                     trc2d(ji,jj,52) = fregenfast(ji,jj)
1190                  ENDIF
1191               ENDDO
1192            ENDDO
1193         ENDIF
1194
1195         IF (jk.eq.i0200) THEN
1196            DO jj = 2,jpjm1
1197               DO ji = 2,jpim1
1198                  IF (tmask(ji,jj,jk) == 1) THEN
1199                     !! N  regeneration to  200 m
1200                     trc2d(ji,jj,53) = fregenfast(ji,jj)
1201                  ENDIF
1202               ENDDO
1203            ENDDO
1204         ENDIF
1205
1206         IF (jk.eq.i0500) THEN
1207            DO jj = 2,jpjm1
1208               DO ji = 2,jpim1
1209                  IF (tmask(ji,jj,jk) == 1) THEN
1210                     !! N  regeneration to  500 m
1211                     trc2d(ji,jj,54) = fregenfast(ji,jj)
1212                  ENDIF
1213               ENDDO
1214            ENDDO
1215         ENDIF
1216
1217         IF (jk.eq.i1000) THEN
1218            DO jj = 2,jpjm1
1219               DO ji = 2,jpim1
1220                  IF (tmask(ji,jj,jk) == 1) THEN
1221                     !! N  regeneration to 1000 m
1222                     trc2d(ji,jj,55) = fregenfast(ji,jj)
1223                  ENDIF
1224               ENDDO
1225            ENDDO
1226         ENDIF
1227
1228         IF (jk.eq.i0100) THEN
1229            DO jj = 2,jpjm1
1230               DO ji = 2,jpim1
1231                  IF (tmask(ji,jj,jk) == 1) THEN
1232                     !! fast detritus Si flux at  100 m
1233                     trc2d(ji,jj,56) = ffastsi(ji,jj)
1234                  ENDIF
1235               ENDDO
1236            ENDDO
1237         ENDIF
1238
1239         IF (jk.eq.i0200) THEN
1240            DO jj = 2,jpjm1
1241               DO ji = 2,jpim1
1242                  IF (tmask(ji,jj,jk) == 1) THEN
1243                     !! fast detritus Si flux at  200 m
1244                     trc2d(ji,jj,57) = ffastsi(ji,jj)
1245                  ENDIF
1246               ENDDO
1247            ENDDO
1248         ENDIF
1249
1250         IF (jk.eq.i0500) THEN
1251            DO jj = 2,jpjm1
1252               DO ji = 2,jpim1
1253                  IF (tmask(ji,jj,jk) == 1) THEN
1254                     !! fast detritus Si flux at  500 m
1255                     trc2d(ji,jj,58) = ffastsi(ji,jj)
1256                  ENDIF
1257               ENDDO
1258            ENDDO
1259         ENDIF
1260
1261         IF (jk.eq.i1000) THEN
1262            DO jj = 2,jpjm1
1263               DO ji = 2,jpim1
1264                  IF (tmask(ji,jj,jk) == 1) THEN
1265                     !! fast detritus Si flux at 1000 m
1266                     trc2d(ji,jj,59) = ffastsi(ji,jj)
1267                  ENDIF
1268               ENDDO
1269            ENDDO
1270         ENDIF
1271
1272         IF (jk.eq.i0100) THEN
1273            DO jj = 2,jpjm1
1274               DO ji = 2,jpim1
1275                  IF (tmask(ji,jj,jk) == 1) THEN
1276                     !! Si regeneration to  100 m
1277                     trc2d(ji,jj,60) = fregenfastsi(ji,jj)
1278                  ENDIF
1279               ENDDO
1280            ENDDO
1281         ENDIF
1282
1283         IF (jk.eq.i0200) THEN
1284            DO jj = 2,jpjm1
1285               DO ji = 2,jpim1
1286                  IF (tmask(ji,jj,jk) == 1) THEN
1287                     !! Si regeneration to  200 m
1288                     trc2d(ji,jj,61) = fregenfastsi(ji,jj)
1289                  ENDIF
1290               ENDDO
1291            ENDDO
1292         ENDIF
1293
1294         IF (jk.eq.i0500) THEN
1295            DO jj = 2,jpjm1
1296               DO ji = 2,jpim1
1297                  IF (tmask(ji,jj,jk) == 1) THEN
1298                     !! Si regeneration to  500 m
1299                     trc2d(ji,jj,62) = fregenfastsi(ji,jj)
1300                  ENDIF
1301               ENDDO
1302            ENDDO
1303         ENDIF
1304
1305         IF (jk.eq.i1000) THEN
1306            DO jj = 2,jpjm1
1307               DO ji = 2,jpim1
1308                  IF (tmask(ji,jj,jk) == 1) THEN
1309                     !! Si regeneration to 1000 m
1310                     trc2d(ji,jj,63) = fregenfastsi(ji,jj)
1311                  ENDIF
1312               ENDDO
1313            ENDDO
1314         ENDIF
1315
1316         DO jj = 2,jpjm1
1317            DO ji = 2,jpim1
1318               IF (tmask(ji,jj,jk) == 1) THEN
1319                  !! sum of fast-sinking N  fluxes
1320                  trc2d(ji,jj,64) = trc2d(ji,jj,64) +                        &
1321                                    (freminn(ji,jj) * fse3t(ji,jj,jk))
1322                  !! sum of fast-sinking Si fluxes
1323                  trc2d(ji,jj,65) = trc2d(ji,jj,65) +                        &
1324                                    (freminsi(ji,jj) * fse3t(ji,jj,jk))
1325                  !! sum of fast-sinking Fe fluxes
1326                  trc2d(ji,jj,66) = trc2d(ji,jj,66) +                        &
1327                                    (freminfe(ji,jj) * fse3t(ji,jj,jk))
1328                  !! sum of fast-sinking C  fluxes
1329                  trc2d(ji,jj,67) = trc2d(ji,jj,67) +                        &
1330                                    (freminc(ji,jj) * fse3t(ji,jj,jk))
1331                  !! sum of fast-sinking Ca fluxes
1332                  trc2d(ji,jj,68) = trc2d(ji,jj,68) +                        &
1333                                    (freminca(ji,jj) * fse3t(ji,jj,jk))
1334               ENDIF
1335            ENDDO
1336         ENDDO
1337
1338
1339         if (jk.eq.mbathy(ji,jj)) then
1340            DO jj = 2,jpjm1
1341               DO ji = 2,jpim1
1342                  IF (tmask(ji,jj,jk) == 1) THEN
1343                     !! N  sedimentation flux
1344                     trc2d(ji,jj,69) = fsedn(ji,jj)
1345                     !! Si sedimentation flux
1346                     trc2d(ji,jj,70) = fsedsi(ji,jj)
1347                     !! Fe sedimentation flux
1348                     trc2d(ji,jj,71) = fsedfe(ji,jj)
1349                     !! C  sedimentation flux
1350                     trc2d(ji,jj,72) = fsedc(ji,jj)
1351                     !! Ca sedimentation flux
1352                     trc2d(ji,jj,73) = fsedca(ji,jj)
1353                  ENDIF
1354               ENDDO
1355            ENDDO
1356         endif
1357
1358         if (jk.eq.1) then
1359            DO jj = 2,jpjm1
1360               DO ji = 2,jpim1
1361                  IF (tmask(ji,jj,jk) == 1) THEN
1362                     trc2d(ji,jj,74) = qsr(ji,jj)
1363                     trc2d(ji,jj,75) = xpar(ji,jj,jk)
1364                     !! trc2d(ji,jj,75) = real(iters(ji,jj))
1365                  ENDIF
1366               ENDDO
1367            ENDDO
1368         endif
1369
1370         DO jj = 2,jpjm1
1371            DO ji = 2,jpim1
1372               IF (tmask(ji,jj,jk) == 1) THEN
1373                  !! diagnostic fields 76 to 80 calculated below
1374                  !! mixed layer non-diatom production
1375                  trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj)
1376                  !! mixed layer     diatom production
1377                  trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj)
1378               ENDIF
1379            ENDDO
1380         ENDDO
1381
1382# if defined key_gulf_finland
1383         if (jk.eq.1) then
1384            DO jj = 2,jpjm1
1385               DO ji = 2,jpim1
1386                  IF (tmask(ji,jj,jk) == 1) THEN
1387                     !! Gulf of Finland check
1388                     trc2d(ji,jj,83) = real(ibio_switch)
1389                  ENDIF
1390               ENDDO
1391            ENDDO
1392         endif
1393# else
1394         DO jj = 2,jpjm1
1395            DO ji = 2,jpim1
1396               IF (tmask(ji,jj,jk) == 1) THEN
1397                  !! calcite CCD depth
1398                  trc2d(ji,jj,83) = ocal_ccd(ji,jj)
1399               ENDIF
1400            ENDDO
1401         ENDDO
1402# endif
1403         DO jj = 2,jpjm1
1404            DO ji = 2,jpim1
1405               IF (tmask(ji,jj,jk) == 1) THEN
1406                  !! last model level above calcite CCD depth
1407                  trc2d(ji,jj,84) = fccd(ji,jj)
1408               ENDIF
1409            ENDDO
1410         ENDDO
1411
1412         IF (jk.eq.1) THEN
1413            DO jj = 2,jpjm1
1414               DO ji = 2,jpim1
1415                  IF (tmask(ji,jj,jk) == 1) THEN
1416                     !! surface "free" iron
1417                     trc2d(ji,jj,85) = xFree(ji,jj)
1418                  ENDIF
1419               ENDDO
1420            ENDDO
1421         ENDIF
1422
1423! I'm keeping this the same as before, but it looks like it should
1424! be i0100 and not i0200 - marc 8/5/17
1425         IF (jk.eq.i0200) THEN
1426            DO jj = 2,jpjm1
1427               DO ji = 2,jpim1
1428                  IF (tmask(ji,jj,jk) == 1) THEN
1429                     !! "free" iron at  100 m
1430                     trc2d(ji,jj,86) = xFree(ji,jj)
1431                  ENDIF
1432               ENDDO
1433            ENDDO
1434         ENDIF
1435
1436
1437         IF (jk.eq.i0200) THEN
1438            DO jj = 2,jpjm1
1439               DO ji = 2,jpim1
1440                  IF (tmask(ji,jj,jk) == 1) THEN
1441                     !! "free" iron at  200 m
1442                     trc2d(ji,jj,87) = xFree(ji,jj)
1443                  ENDIF
1444               ENDDO
1445            ENDDO
1446         ENDIF
1447
1448
1449         IF (jk.eq.i0500) THEN
1450            DO jj = 2,jpjm1
1451               DO ji = 2,jpim1
1452                  IF (tmask(ji,jj,jk) == 1) THEN
1453                     !! "free" iron at  500 m
1454                     trc2d(ji,jj,88) = xFree(ji,jj)
1455                  ENDIF
1456               ENDDO
1457            ENDDO
1458         ENDIF
1459
1460
1461         IF (jk.eq.i1000) THEN
1462            DO jj = 2,jpjm1
1463               DO ji = 2,jpim1
1464                  IF (tmask(ji,jj,jk) == 1) THEN
1465                     !! "free" iron at 1000 m
1466                     trc2d(ji,jj,89) = xFree(ji,jj)
1467                  ENDIF
1468               ENDDO
1469            ENDDO
1470         ENDIF
1471
1472
1473         IF (jk.eq.1) THEN
1474            DO jj = 2,jpjm1
1475               DO ji = 2,jpim1
1476                  IF (tmask(ji,jj,jk) == 1) THEN
1477                     !! AXY (27/06/12): extract "euphotic depth"
1478                     trc2d(ji,jj,90) = xze(ji,jj)
1479                  ENDIF
1480               ENDDO
1481            ENDDO
1482         ENDIF
1483
1484# if defined key_roam
1485         if (jk .eq. 1) then
1486            DO jj = 2,jpjm1
1487               DO ji = 2,jpim1
1488                  IF (tmask(ji,jj,jk) == 1) THEN
1489                     !! ROAM provisionally has access to a further 20 2D
1490                     !! diagnostics
1491                     !! surface wind
1492                     trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + wndm(ji,jj)
1493                     !! atmospheric pCO2
1494                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm(ji,jj)
1495                     !! ocean pH
1496                     trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph(ji,jj)
1497                     !! ocean pCO2
1498                     trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w(ji,jj)
1499                     !! ocean H2CO3 conc.
1500                     trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3(ji,jj)
1501                     !! ocean HCO3 conc.
1502                     trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3(ji,jj)
1503                     !! ocean CO3 conc.
1504                     trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3(ji,jj)
1505                     !! air-sea CO2 flux
1506                     trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux(ji,jj)
1507                  ENDIF
1508               ENDDO
1509           ENDDO
1510
1511            DO jj = 2,jpjm1
1512               DO ji = 2,jpim1
1513                  IF (tmask(ji,jj,jk) == 1) THEN
1514                     !! ocean omega calcite
1515                     trc2d(ji,jj,99)  = trc2d(ji,jj,99)  + f_omcal(ji,jj)
1516                     !! ocean omega aragonite
1517                     trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj)
1518                     !! ocean TDIC
1519                     trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC(ji,jj)
1520                     !! ocean TALK
1521                     trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK(ji,jj)
1522                     !! surface kw660
1523                     trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660(ji,jj)
1524                     !! surface pressure
1525                     trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0(ji,jj)
1526                     !! air-sea O2 flux
1527                     trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux(ji,jj)
1528                     !! ocean O2 saturation
1529                     trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat(ji,jj)
1530                     !! depth calcite CCD
1531                     trc2d(ji,jj,107) = f2_ccd_cal(ji,jj)
1532                     !! depth aragonite CCD
1533                     trc2d(ji,jj,108) = f2_ccd_arg(ji,jj)
1534                  ENDIF
1535               ENDDO
1536            ENDDO
1537         endif
1538
1539         if (jk .eq. mbathy(ji,jj)) then
1540            DO jj = 2,jpjm1
1541               DO ji = 2,jpim1
1542                  IF (tmask(ji,jj,jk) == 1) THEN
1543                     !! seafloor omega calcite
1544                     trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)
1545                     !! seafloor omega aragonite
1546                     trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)
1547                  ENDIF
1548               ENDDO
1549            ENDDO
1550         endif
1551
1552         if (jk.eq.i0100) then
1553            DO jj = 2,jpjm1
1554               DO ji = 2,jpim1
1555                  IF (tmask(ji,jj,jk) == 1) THEN
1556                     !! diagnostic fields 111 to 117 calculated below
1557                     !! rain ratio at  100 m
1558                     trc2d(ji,jj,118) =                                      &
1559                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
1560                  ENDIF
1561               ENDDO
1562            ENDDO
1563         endif
1564
1565         if (jk.eq.i0500) then
1566            DO jj = 2,jpjm1
1567               DO ji = 2,jpim1
1568                  IF (tmask(ji,jj,jk) == 1) THEN
1569                     !! rain ratio at  500 m
1570                     trc2d(ji,jj,119) =                                      &
1571                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
1572                  ENDIF
1573               ENDDO
1574            ENDDO
1575         endif
1576
1577         if (jk.eq.i1000) then
1578            DO jj = 2,jpjm1
1579               DO ji = 2,jpim1
1580                  IF (tmask(ji,jj,jk) == 1) THEN
1581                     !! rain ratio at 1000 m
1582                     trc2d(ji,jj,120) =                                      &
1583                                   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
1584                  ENDIF
1585               ENDDO
1586            ENDDO
1587         endif
1588
1589         if (jk.eq.mbathy(ji,jj)) then
1590            DO jj = 2,jpjm1
1591               DO ji = 2,jpim1
1592                  IF (tmask(ji,jj,jk) == 1) THEN
1593                     !! AXY (18/01/12): benthic flux diagnostics
1594                     trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
1595                     trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
1596                     trc2d(ji,jj,123) = f_sbenin_c(ji,jj)  + f_fbenin_c(ji,jj)
1597                     trc2d(ji,jj,124) = f_fbenin_si(ji,jj)
1598                     trc2d(ji,jj,125) = f_fbenin_ca(ji,jj)
1599                     trc2d(ji,jj,126) = f_benout_n(ji,jj)
1600                     trc2d(ji,jj,127) = f_benout_fe(ji,jj)
1601                     trc2d(ji,jj,128) = f_benout_c(ji,jj)
1602                     trc2d(ji,jj,129) = f_benout_si(ji,jj)
1603                     trc2d(ji,jj,130) = f_benout_ca(ji,jj)
1604                  ENDIF
1605               ENDDO
1606            ENDDO
1607         endif
1608
1609         DO jj = 2,jpjm1
1610            DO ji = 2,jpim1
1611               IF (tmask(ji,jj,jk) == 1) THEN
1612                  !! diagnostics fields 131 to 135 calculated below
1613                  trc2d(ji,jj,136) = f_runoff(ji,jj)
1614                  !! AXY (19/07/12): amended to allow for riverine
1615                  !! nutrient addition below surface
1616                  trc2d(ji,jj,137) = trc2d(ji,jj,137) +                      &
1617                                     (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
1618                  trc2d(ji,jj,138) = trc2d(ji,jj,138) +                      &
1619                                     (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
1620                  trc2d(ji,jj,139) = trc2d(ji,jj,139) +                      &
1621                                     (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
1622                  trc2d(ji,jj,140) = trc2d(ji,jj,140) +                      &
1623                                     (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
1624                  !! slow sinking detritus C production
1625                  trc2d(ji,jj,141) = trc2d(ji,jj,141) +                      &
1626                                     (fslowc(ji,jj)  * fse3t(ji,jj,jk))
1627               ENDIF
1628            ENDDO
1629         ENDDO
1630
1631         if (jk.eq.i0100) then
1632            DO jj = 2,jpjm1
1633               DO ji = 2,jpim1
1634                  IF (tmask(ji,jj,jk) == 1) THEN
1635                     !! slow detritus flux at  100 m
1636                     trc2d(ji,jj,142) = fslowcflux(ji,jj)
1637                  ENDIF
1638               ENDDO
1639            ENDDO
1640         endif
1641
1642         if (jk.eq.i0200) then
1643            DO jj = 2,jpjm1
1644               DO ji = 2,jpim1
1645                  IF (tmask(ji,jj,jk) == 1) THEN
1646                     !! slow detritus flux at  200 m
1647                     trc2d(ji,jj,143) = fslowcflux(ji,jj)
1648                  ENDIF
1649               ENDDO
1650            ENDDO
1651         endif
1652
1653
1654         if (jk.eq.i0500) then
1655            DO jj = 2,jpjm1
1656               DO ji = 2,jpim1
1657                  IF (tmask(ji,jj,jk) == 1) THEN
1658                     !! slow detritus flux at  500 m
1659                     trc2d(ji,jj,144) = fslowcflux(ji,jj)
1660                  ENDIF
1661               ENDDO
1662            ENDDO
1663         endif
1664
1665
1666         if (jk.eq.i1000) then
1667            DO jj = 2,jpjm1
1668               DO ji = 2,jpim1
1669                  IF (tmask(ji,jj,jk) == 1) THEN
1670                     !! slow detritus flux at 1000 m
1671                     trc2d(ji,jj,145) = fslowcflux(ji,jj)
1672                  ENDIF
1673               ENDDO
1674            ENDDO
1675         endif
1676
1677         DO jj = 2,jpjm1
1678            DO ji = 2,jpim1
1679               IF (tmask(ji,jj,jk) == 1) THEN
1680                  !! carbon     inventory
1681                  trc2d(ji,jj,146)  = trc2d(ji,jj,146)  + ftot_c(ji,jj)
1682                  !! alkalinity inventory
1683                  trc2d(ji,jj,147)  = trc2d(ji,jj,147)  + ftot_a(ji,jj)
1684                  !! oxygen     inventory
1685                  trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)
1686               ENDIF
1687            ENDDO
1688         ENDDO
1689
1690         if (jk.eq.mbathy(ji,jj)) then
1691            DO jj = 2,jpjm1
1692               DO ji = 2,jpim1
1693                  IF (tmask(ji,jj,jk) == 1) THEN
1694                     trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
1695                  ENDIF
1696               ENDDO
1697            ENDDO
1698         endif
1699
1700         DO jj = 2,jpjm1
1701            DO ji = 2,jpim1
1702               IF (tmask(ji,jj,jk) == 1) THEN
1703                  !! community respiration
1704                  trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)
1705               ENDIF
1706            ENDDO
1707         ENDDO
1708
1709         DO jj = 2,jpjm1
1710            DO ji = 2,jpim1
1711               IF (tmask(ji,jj,jk) == 1) THEN
1712        !!
1713        !! AXY (14/02/14): a Valentines Day gift to BASIN - a
1714                  !!                 shedload of new diagnostics that
1715                  !!                 they'll most likely never need!
1716                  !!                 (actually, as with all such gifts,
1717                  !!                 I'm giving them some things I'd like
1718                  !!                 myself!)
1719                  !!
1720                  !! ------------------------------------------------------
1721                  !! linear losses
1722                  !! non-diatom
1723                  trc2d(ji,jj,151) = trc2d(ji,jj,151) +                      &
1724                                     (fdpn2(ji,jj) * fse3t(ji,jj,jk))
1725                  !! diatom
1726                  trc2d(ji,jj,152) = trc2d(ji,jj,152) +                      &
1727                                     (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
1728                  !! microzooplankton
1729                  trc2d(ji,jj,153) = trc2d(ji,jj,153) +                      &
1730                                     (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
1731                  !! mesozooplankton
1732                  trc2d(ji,jj,154) = trc2d(ji,jj,154) +                      &
1733                                     (fdzme2(ji,jj) * fse3t(ji,jj,jk))
1734               ENDIF
1735            ENDDO
1736         ENDDO
1737
1738         DO jj = 2,jpjm1
1739            DO ji = 2,jpim1
1740               IF (tmask(ji,jj,jk) == 1) THEN
1741                  !! ------------------------------------------------------
1742                  !! microzooplankton grazing
1743                  !! microzooplankton messy -> N
1744                  trc2d(ji,jj,155) = trc2d(ji,jj,155) +                      &
1745                                     (xphi * (fgmipn(ji,jj) +                &
1746                                              fgmid(ji,jj)) * fse3t(ji,jj,jk))
1747                  !! microzooplankton messy -> D
1748                  trc2d(ji,jj,156) = trc2d(ji,jj,156) +                      &
1749                                     ((1. - xbetan) * finmi(ji,jj) *         &
1750                                      fse3t(ji,jj,jk))
1751                  !! microzooplankton messy -> DIC
1752                  trc2d(ji,jj,157) = trc2d(ji,jj,157) +                      &
1753                                     (xphi * ((xthetapn * fgmipn(ji,jj)) +   &
1754                                              fgmidc(ji,jj)) *               &
1755                                      fse3t(ji,jj,jk))
1756                  !! microzooplankton messy -> Dc
1757                  trc2d(ji,jj,158) = trc2d(ji,jj,158) +                      &
1758                                     ((1. - xbetac) * ficmi(ji,jj) *         &
1759                                      fse3t(ji,jj,jk))
1760                  !! microzooplankton excretion
1761                  trc2d(ji,jj,159) = trc2d(ji,jj,159) +                      &
1762                                     (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
1763                  !! microzooplankton respiration
1764                  trc2d(ji,jj,160) = trc2d(ji,jj,160) +                      &
1765                                     (fmiresp(ji,jj) * fse3t(ji,jj,jk))
1766                  !! microzooplankton growth
1767                  trc2d(ji,jj,161) = trc2d(ji,jj,161) +                      &
1768                                     (fmigrow(ji,jj) * fse3t(ji,jj,jk))
1769               ENDIF
1770            ENDDO
1771         ENDDO
1772
1773         DO jj = 2,jpjm1
1774            DO ji = 2,jpim1
1775               IF (tmask(ji,jj,jk) == 1) THEN
1776                  !! ------------------------------------------------------
1777                  !! mesozooplankton grazing
1778                  !! mesozooplankton messy -> N
1779                  trc2d(ji,jj,162) = trc2d(ji,jj,162) +                      &
1780                                     (xphi *                                 &
1781                                      (fgmepn(ji,jj) + fgmepd(ji,jj) +       &
1782                                       fgmezmi(ji,jj) + fgmed(ji,jj)) *      &
1783                                      fse3t(ji,jj,jk))
1784                  !! mesozooplankton messy -> D
1785                  trc2d(ji,jj,163) = trc2d(ji,jj,163) +                      &
1786                                     ((1. - xbetan) * finme(ji,jj) *         &
1787                                      fse3t(ji,jj,jk))
1788                  !! mesozooplankton messy -> DIC
1789                  trc2d(ji,jj,164) = trc2d(ji,jj,164) +                      &
1790                                     (xphi *                                 &
1791                                      ((xthetapn * fgmepn(ji,jj)) +          &
1792                                       (xthetapd * fgmepd(ji,jj)) +          &
1793                                       (xthetazmi * fgmezmi(ji,jj)) +        &
1794                                      fgmedc(ji,jj)) * fse3t(ji,jj,jk))
1795                  !! mesozooplankton messy -> Dc
1796                  trc2d(ji,jj,165) = trc2d(ji,jj,165) +                      &
1797                                     ((1. - xbetac) * ficme(ji,jj) *         &
1798                                      fse3t(ji,jj,jk))
1799                  !! mesozooplankton excretion
1800                  trc2d(ji,jj,166) = trc2d(ji,jj,166) +                      &
1801                                     (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
1802                  !! mesozooplankton respiration
1803                  trc2d(ji,jj,167) = trc2d(ji,jj,167) +                      &
1804                                     (fmeresp(ji,jj) * fse3t(ji,jj,jk))
1805                  !! mesozooplankton growth
1806                  trc2d(ji,jj,168) = trc2d(ji,jj,168) +                      &
1807                                     (fmegrow(ji,jj) * fse3t(ji,jj,jk))
1808               ENDIF
1809            ENDDO
1810         ENDDO
1811
1812         DO jj = 2,jpjm1
1813            DO ji = 2,jpim1
1814               IF (tmask(ji,jj,jk) == 1) THEN
1815                  !! ------------------------------------------------------
1816                  !! miscellaneous
1817                  !! detrital C remineralisation
1818                  trc2d(ji,jj,169) = trc2d(ji,jj,169) +                      &
1819                                     (fddc(ji,jj) * fse3t(ji,jj,jk))
1820                  !! microzoo grazing on detrital carbon
1821                  trc2d(ji,jj,170) = trc2d(ji,jj,170) +                      &
1822                                     (fgmidc(ji,jj)  * fse3t(ji,jj,jk))
1823                  !! mesozoo  grazing on detrital carbon
1824                  trc2d(ji,jj,171) = trc2d(ji,jj,171) +                      &
1825                                     (fgmedc(ji,jj)  * fse3t(ji,jj,jk))
1826                  !!
1827               ENDIF
1828            ENDDO
1829         ENDDO
1830
1831         !! ------------------------------------------------------
1832    !!
1833         !! AXY (23/10/14): extract primary production related
1834         !!                 surface fields to deal with diel
1835         !!                 cycle issues; hijacking BASIN 150m
1836         !!                 diagnostics to do so (see commented
1837         !!                 out diagnostics below this section)
1838         !!
1839         !! extract relevant BASIN fields at 150m
1840         if (jk .eq. i0150) then
1841            DO jj = 2,jpjm1
1842               DO ji = 2,jpim1
1843                  IF (tmask(ji,jj,jk) == 1) THEN
1844                     !! Pn PP
1845                     trc2d(ji,jj,172) = trc2d(ji,jj,4)
1846                     !! Pn linear loss
1847                     trc2d(ji,jj,173) = trc2d(ji,jj,151)
1848                     !! Pn non-linear loss
1849                     trc2d(ji,jj,174) = trc2d(ji,jj,5)
1850                     !! Pn grazing to Zmi
1851                     trc2d(ji,jj,175) = trc2d(ji,jj,11)
1852                     !! Pn grazing to Zme
1853                     trc2d(ji,jj,176) = trc2d(ji,jj,14)
1854                     !! Pd PP
1855                     trc2d(ji,jj,177) = trc2d(ji,jj,6)
1856                     !! Pd linear loss
1857                     trc2d(ji,jj,178) = trc2d(ji,jj,152)
1858                     !! Pd non-linear loss
1859                     trc2d(ji,jj,179) = trc2d(ji,jj,7)
1860                     !! Pd grazing to Zme
1861                     trc2d(ji,jj,180) = trc2d(ji,jj,15)
1862                     !! Zmi grazing on D
1863                     trc2d(ji,jj,181) = trc2d(ji,jj,12)
1864                     !! Zmi grazing on Dc
1865                     trc2d(ji,jj,182) = trc2d(ji,jj,170)
1866                     !! Zmi messy feeding loss to N
1867                     trc2d(ji,jj,183) = trc2d(ji,jj,155)
1868                     !! Zmi messy feeding loss to D
1869                     trc2d(ji,jj,184) = trc2d(ji,jj,156)
1870                     !! Zmi messy feeding loss to DIC
1871                     trc2d(ji,jj,185) = trc2d(ji,jj,157)
1872                     !! Zmi messy feeding loss to Dc
1873                     trc2d(ji,jj,186) = trc2d(ji,jj,158)
1874                     !! Zmi excretion
1875                     trc2d(ji,jj,187) = trc2d(ji,jj,159)
1876                     !! Zmi respiration
1877                     trc2d(ji,jj,188) = trc2d(ji,jj,160)
1878                     !! Zmi growth
1879                     trc2d(ji,jj,189) = trc2d(ji,jj,161)
1880                     !! Zmi linear loss
1881                     trc2d(ji,jj,190) = trc2d(ji,jj,153)
1882                     !! Zmi non-linear loss
1883                     trc2d(ji,jj,191) = trc2d(ji,jj,13)
1884                     !! Zmi grazing to Zme
1885                     trc2d(ji,jj,192) = trc2d(ji,jj,16)
1886                     !! Zme grazing on D
1887                     trc2d(ji,jj,193) = trc2d(ji,jj,17)
1888                     !! Zme grazing on Dc
1889                     trc2d(ji,jj,194) = trc2d(ji,jj,171)
1890                     !! Zme messy feeding loss to N
1891                     trc2d(ji,jj,195) = trc2d(ji,jj,162)
1892                     !! Zme messy feeding loss to D
1893                     trc2d(ji,jj,196) = trc2d(ji,jj,163)
1894                     !! Zme messy feeding loss to DIC
1895                     trc2d(ji,jj,197) = trc2d(ji,jj,164)
1896                     !! Zme messy feeding loss to Dc
1897                     trc2d(ji,jj,198) = trc2d(ji,jj,165)
1898                     !! Zme excretion
1899                     trc2d(ji,jj,199) = trc2d(ji,jj,166)
1900                     !! Zme respiration
1901                     trc2d(ji,jj,200) = trc2d(ji,jj,167)
1902                     !! Zme growth
1903                     trc2d(ji,jj,201) = trc2d(ji,jj,168)
1904                     !! Zme linear loss
1905                     trc2d(ji,jj,202) = trc2d(ji,jj,154)
1906                     !! Zme non-linear loss
1907                     trc2d(ji,jj,203) = trc2d(ji,jj,18)
1908                     !! Slow detritus production, N
1909                     trc2d(ji,jj,204) = trc2d(ji,jj,20)
1910                     !! Slow detritus remineralisation, N
1911                     trc2d(ji,jj,205) = trc2d(ji,jj,21)
1912                     !! Slow detritus production, C
1913                     trc2d(ji,jj,206) = trc2d(ji,jj,141)
1914                     !! Slow detritus remineralisation, C
1915                     trc2d(ji,jj,207) = trc2d(ji,jj,169)
1916                     !! Fast detritus production, N
1917                     trc2d(ji,jj,208) = trc2d(ji,jj,43)
1918                     !! Fast detritus remineralisation, N
1919                     trc2d(ji,jj,209) = trc2d(ji,jj,21)
1920                     !! Fast detritus production, C
1921                     trc2d(ji,jj,210) = trc2d(ji,jj,64)
1922                     !! Fast detritus remineralisation, C
1923                     trc2d(ji,jj,211) = trc2d(ji,jj,67)
1924                     !! Community respiration
1925                     trc2d(ji,jj,212) = trc2d(ji,jj,150)
1926                     !! Slow detritus N flux at 150 m
1927                     trc2d(ji,jj,213) = fslownflux(ji,jj)
1928                     !! Slow detritus C flux at 150 m
1929                     trc2d(ji,jj,214) = fslowcflux(ji,jj)
1930                     !! Fast detritus N flux at 150 m
1931                     trc2d(ji,jj,215) = ffastn(ji,jj)
1932                     !! Fast detritus C flux at 150 m
1933                     trc2d(ji,jj,216) = ffastc(ji,jj)
1934                  ENDIF
1935               ENDDO
1936            ENDDO
1937         endif
1938
1939         !!
1940         !! Jpalm (11-08-2014)
1941         !! Add UKESM1 diagnoatics
1942         !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1943         if ((jk .eq. 1) .and.( jdms.eq.1)) then
1944            DO jj = 2,jpjm1
1945               DO ji = 2,jpim1
1946                  IF (tmask(ji,jj,jk) == 1) THEN
1947                     !! DMS surface concentration
1948                     trc2d(ji,jj,221) = dms_surf(ji,jj)
1949                     !! AXY (13/03/15): add in other DMS estimates
1950                     !! DMS surface concentration
1951                     trc2d(ji,jj,222) = dms_andr(ji,jj)
1952                     !! DMS surface concentration
1953                     trc2d(ji,jj,223) = dms_simo(ji,jj)
1954                     !! DMS surface concentration
1955                     trc2d(ji,jj,224) = dms_aran(ji,jj)
1956                     !! DMS surface concentration
1957                     trc2d(ji,jj,225) = dms_hall(ji,jj)
1958                  ENDIF
1959               ENDDO
1960            ENDDO
1961         endif
1962# endif
1963
1964         DO jj = 2,jpjm1
1965            DO ji = 2,jpim1
1966               IF (tmask(ji,jj,jk) == 1) THEN
1967                  !! other possible future diagnostics include:
1968                  !!   - integrated tracer values (esp. biological)
1969                  !!   - mixed layer tracer values
1970                  !!   - sub-surface chlorophyll maxima (plus depth)
1971                  !!   - different mixed layer depth criteria (T, sigma,
1972                  !!     var. sigma)
1973                  !!-------------------------------------------------------
1974                  !! Prepare 3D diagnostics
1975                  !!-------------------------------------------------------
1976                  !!
1977                  !! primary production 
1978                  trc3d(ji,jj,jk,1)  = ((fprn(ji,jj) + fprd(ji,jj)) *        &
1979                                        zphn(ji,jj))
1980                  !! detrital flux
1981                  trc3d(ji,jj,jk,2)  = fslownflux(ji,jj) + ffastn(ji,jj)
1982                  !! remineralisation
1983                  trc3d(ji,jj,jk,3)  = fregen(ji,jj) +                       &
1984                                       (freminn(ji,jj) * fse3t(ji,jj,jk))
1985               ENDIF
1986            ENDDO
1987         ENDDO
1988# if defined key_roam
1989         DO jj = 2,jpjm1
1990            DO ji = 2,jpim1
1991               IF (tmask(ji,jj,jk) == 1) THEN
1992                  !! pH
1993                  trc3d(ji,jj,jk,4)  = f3_pH(ji,jj,jk)
1994                  !! omega calcite
1995                  trc3d(ji,jj,jk,5)  = f3_omcal(ji,jj,jk)
1996               ENDIF
1997            ENDDO
1998         ENDDO
1999# else
2000         DO jj = 2,jpjm1
2001            DO ji = 2,jpim1
2002               IF (tmask(ji,jj,jk) == 1) THEN
2003                  !! fast Si flux
2004                  trc3d(ji,jj,jk,4)  = ffastsi(ji,jj)
2005               ENDIF
2006            ENDDO
2007         ENDDO
2008# endif
2009
2010      ENDIF   ! end of ln_diatrc option
2011
2012   END SUBROUTINE bio_medusa_diag
2013
2014#else
2015   !!======================================================================
2016   !!  Dummy module :                                   No MEDUSA bio-model
2017   !!======================================================================
2018CONTAINS
2019   SUBROUTINE bio_medusa_diag( )                    ! Empty routine
2020      WRITE(*,*) 'bio_medusa_diag: You should not have seen this print! error?'
2021   END SUBROUTINE bio_medusa_diag
2022#endif 
2023
2024   !!======================================================================
2025END MODULE bio_medusa_diag_mod
Note: See TracBrowser for help on using the repository browser.