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

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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 10020

Last change on this file since 10020 was 10020, checked in by marc, 6 years ago

GMED ticket 406. CPP key fixes.
File size: 46.2 KB
Line 
1MODULE trcbio_medusa
2   !!======================================================================
3   !!                         ***  MODULE trcbio  ***
4   !! TOP :   MEDUSA
5   !!======================================================================
6   !! History :
7   !!  -   !  1999-07  (M. Levy)              original code
8   !!  -   !  2000-12  (E. Kestenare)         assign parameters to name
9   !!                                         individual tracers
10   !!  -   !  2001-03  (M. Levy)              LNO3 + dia2d
11   !! 2.0  !  2007-12  (C. Deltel, G. Madec)  F90
12   !!  -   !  2008-08  (K. Popova)            adaptation for MEDUSA
13   !!  -   !  2008-11  (A. Yool)              continuing adaptation for MEDUSA
14   !!  -   !  2010-03  (A. Yool)              updated for branch inclusion
15   !!  -   !  2011-08  (A. Yool)              updated for ROAM (see below)
16   !!  -   !  2013-03  (A. Yool)              updated for iMARNET
17   !!  -   !  2013-05  (A. Yool)              updated for v3.5
18   !!  -   !  2014-08  (A. Yool, J. Palm)     Add DMS module for UKESM1 model
19   !!  -   !  2015-06  (A. Yool)              Update to include MOCSY
20   !!  -   !  2015-07  (A. Yool)              Update for rolling averages
21   !!  -   !  2015-10  (J. Palm)              Update for diag outputs through
22   !!                                         iom_use 
23   !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6
24   !!  -   !  2017-05  (A. Yool)              Added extra DMS calculation
25   !!  -   !  2017-11  (J. Palm, A. Yool)     Diagnose tracer excursions
26   !!----------------------------------------------------------------------
27   !!
28#if defined key_roam
29   !!----------------------------------------------------------------------
30   !! Updates for the ROAM project include:
31   !!   - addition of DIC, alkalinity, detrital carbon and oxygen tracers
32   !!   - addition of air-sea fluxes of CO2 and oxygen
33   !!   - periodic (monthly) calculation of full 3D carbonate chemistry
34   !!   - detrital C:N ratio now free to evolve dynamically
35   !!   - benthic storage pools
36   !!
37   !! Opportunity also taken to add functionality:
38   !!   - switch for Liebig Law (= most-limiting) nutrient uptake
39   !!   - switch for accelerated seafloor detritus remineralisation
40   !!   - switch for fast -> slow detritus transfer at seafloor
41   !!   - switch for ballast vs. Martin vs. Henson fast detritus remin.
42   !!   - per GMD referee remarks, xfdfrac3 introduced for grazed PDS
43   !!----------------------------------------------------------------------
44#endif
45   !!
46#if defined key_mocsy
47   !!----------------------------------------------------------------------
48   !! Updates with the addition of MOCSY include:
49   !!   - option to use PML or MOCSY carbonate chemistry (the latter is
50   !!     preferred)
51   !!   - central calculation of gas transfer velocity, f_kw660; previously
52   !!     this was done separately for CO2 and O2 with predictable results
53   !!   - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux
54   !!     calculations and to those for O2 air-sea flux
55   !!   - extra diagnostics included for MOCSY
56   !!----------------------------------------------------------------------
57#endif
58   !!
59#if defined key_medusa
60   !!----------------------------------------------------------------------
61   !!                                        MEDUSA bio-model
62   !!----------------------------------------------------------------------
63   !!   trc_bio_medusa        : 
64   !!----------------------------------------------------------------------
65      !! AXY (30/01/14): necessary to find NaNs on HECTOR
66      USE, INTRINSIC :: ieee_arithmetic 
67
68      USE bio_medusa_mod,             ONLY: b0, fdep1,                      & 
69                                            ibenthic, idf, idfval,          &
70# if defined key_roam
71                                            f_xco2a,                        &
72                                            zalk, zdic, zoxy, zsal, ztmp,   &
73# endif
74# if defined key_mocsy
75                                            zpho,                           &
76# endif
77                                            zchd, zchn, zdet, zdin, zdtc,   &
78                                            zfer, zpds, zphd, zphn, zsil,   &
79                                            zzme, zzmi
80      USE dom_oce,                    ONLY: e3t_0, gdept_0,  gdepw_0,       &
81                                            nday_year, nsec_day,            & 
82                                            nyear, nyear_len, ndastp,       &
83                                            nsec_month,                     &
84                                            rdt, tmask, mig, mjg, nimpp,    &
85                                            njmpp 
86#if defined key_vvl
87      USE dom_oce,                    ONLY: e3t_n, gdept_n,  gdepw_n
88#endif
89
90      USE in_out_manager,             ONLY: lwp, numout, nn_date0
91      USE iom,                        ONLY: lk_iomput
92      USE lbclnk,                     ONLY: lbc_lnk
93      USE lib_mpp,                    ONLY: mpp_max, mpp_maxloc,            & 
94                                            mpp_min, mpp_minloc,            &
95                                            ctl_stop, ctl_warn, lk_mpp 
96      USE oce,                        ONLY: tsb, tsn, PCO2a_in_cpl
97      USE par_kind,                   ONLY: wp
98      USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     &
99                                            jpdic, jpdin, jpdtc, jpfer,     &
100                                            jpoxy, jppds, jpphd, jpphn,     &
101                                            jpsil, jpzme, jpzmi
102      USE par_oce,                    ONLY: jp_sal, jp_tem, jpi, jpim1,     &
103                                            jpj, jpjm1, jpk
104      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
105      USE sbc_oce,                    ONLY: lk_oasis
106      USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend, &
107                                            lk_pi_co2
108      USE trc,                        ONLY: ln_rsttr, nittrc000, trn
109      USE bio_medusa_init_mod,        ONLY: bio_medusa_init
110      USE carb_chem_mod,              ONLY: carb_chem
111      USE air_sea_mod,                ONLY: air_sea
112      USE plankton_mod,               ONLY: plankton
113      USE iron_chem_scav_mod,         ONLY: iron_chem_scav
114      USE detritus_mod,               ONLY: detritus
115      USE bio_medusa_update_mod,      ONLY: bio_medusa_update
116      USE bio_medusa_diag_mod,        ONLY: bio_medusa_diag
117      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
118      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
119      USE trcstat,                    ONLY: trc_rst_dia_stat
120
121      IMPLICIT NONE
122      PRIVATE
123     
124      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
125      PUBLIC   trc_bio_exceptionnal_fix ! here
126
127   !!* Substitution
128#  include "domzgr_substitute.h90"
129   !!----------------------------------------------------------------------
130   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
131   !! $Id$
132   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
133   !!----------------------------------------------------------------------
134
135CONTAINS
136
137   SUBROUTINE trc_bio_medusa( kt )
138      !!------------------------------------------------------------------
139      !!                     ***  ROUTINE trc_bio  ***
140      !!
141      !! ** Purpose : compute the now trend due to biogeochemical processes
142      !!              and add it to the general trend of passive tracers
143      !!              equations
144      !!
145      !! ** Method  : each now biological flux is calculated in function of
146      !!              now concentrations of tracers.
147      !!              depending on the tracer, these fluxes are sources or
148      !!              sinks.
149      !!              The total of the sources and sinks for each tracer
150      !!              is added to the general trend.
151      !!       
152      !!                      tra = tra + zf...tra - zftra...
153      !!                                     |         |
154      !!                                     |         |
155      !!                                  source      sink
156      !!       
157      !!              IF 'key_trc_diabio' defined , the biogeochemical trends
158      !!              for passive tracers are saved for futher diagnostics.
159      !!------------------------------------------------------------------
160      !!
161      !!
162      !!------------------------------------------------------------------
163      !! Variable conventions
164      !!------------------------------------------------------------------
165      !!
166      !! names: z*** - state variable
167      !!        f*** - function (or temporary variable used in part of
168      !!               a function)
169      !!        x*** - parameter
170      !!        b*** - right-hand part (sources and sinks)
171      !!        i*** - integer variable (usually used in yes/no flags)
172      !!
173      !! time (integer timestep)
174      INTEGER, INTENT( in ) ::    kt
175      !!
176      !! spatial array indices
177      INTEGER  ::    ji,jj,jk,jn
178      !!
179      INTEGER  ::    iball
180# if defined key_roam
181      !!
182      INTEGER  ::    iyr1, iyr2
183      !!
184# endif
185      !!
186      !! temporary variables
187      REAL(wp) ::    fq3,fq4
188      REAL(wp) ::    this_y, this_d, this_s, fyear
189      !!
190      !! T and S check temporary variable
191      REAL(wp) ::    sumtsn, tsnavg
192      INTEGER  ::    summask
193      CHARACTER(40) :: charout, charout2, charout3, charout4, charout5
194      !!
195      !!------------------------------------------------------------------
196
197# if defined key_debug_medusa
198      IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined'
199      CALL flush(numout)
200# endif 
201
202      !! AXY (20/11/14): alter this to report on first MEDUSA call
203      !! IF( kt == nit000 ) THEN
204      IF( kt == nittrc000 ) THEN
205         IF(lwp) WRITE(numout,*)
206         IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model'
207         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
208    IF(lwp) WRITE(numout,*) ' kt =',kt
209      ENDIF
210
211      !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
212      ibenthic = 1
213
214      !!------------------------------------------------------------------
215      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency
216      !! terms of all biological equations to 0.
217      !!------------------------------------------------------------------
218      !!
219      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit
220      !!                 the bill for now; another item on the things-to-sort-
221      !!     out-in-the-future list ...
222# if defined key_kill_medusa
223      b0 = 0.
224# else
225      b0 = 1.
226# endif
227      !!------------------------------------------------------------------
228      !! fast detritus ballast scheme (0 = no; 1 = yes)
229      !! alternative to ballast scheme is same scheme but with no ballast
230      !! protection (not dissimilar to Martin et al., 1987)
231      !!------------------------------------------------------------------
232      !!
233      iball = 1
234
235      !!------------------------------------------------------------------
236      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
237      !! these should *only* be used in 1D since they give comprehensive
238      !! output for ecological functions in the model; primarily used in
239      !! debugging
240      !!------------------------------------------------------------------
241      !!
242      idf    = 0
243      !!
244      !! timer mechanism
245      if (kt/120*120.eq.kt) then
246         idfval = 1
247      else
248         idfval = 0
249      endif
250
251      !!------------------------------------------------------------------
252      !! Initialise arrays to zero and set up arrays for diagnostics
253      !!------------------------------------------------------------------
254      CALL bio_medusa_init( kt )
255
256# if defined key_axy_nancheck
257       DO jn = jp_msa0,jp_msa1
258          !! fq0 = MINVAL(trn(:,:,:,jn))
259          !! fq1 = MAXVAL(trn(:,:,:,jn))
260          fq2 = SUM(trn(:,:,:,jn))
261          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK',     &
262          !!                kt, jn, fq0, fq1, fq2
263          !! AXY (30/01/14): much to our surprise, the next line doesn't
264          !!                 work on HECTOR and has been replaced here with
265          !!                 a specialist routine
266          !! if (fq2 /= fq2 ) then
267          if ( ieee_is_nan( fq2 ) ) then
268             !! there's a NaN here
269             if (lwp) write(numout,*) 'NAN detected in field', jn,           &
270                                      'at time', kt, 'at position:'
271             DO jk = 1,jpk
272                DO jj = 1,jpj
273                   DO ji = 1,jpi
274                      !! AXY (30/01/14): "isnan" problem on HECTOR
275                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then
276                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
277                         if (lwp) write (numout,'(a,1pe12.2,4i6)')           &
278                            'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn
279                      endif
280                   enddo
281                enddo
282             enddo
283             CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
284          endif
285       ENDDO
286       CALL flush(numout)
287# endif
288
289# if defined key_debug_medusa
290      IF (lwp) write (numout,*)                                              &
291                     'trc_bio_medusa: variables initialised and checked'
292      CALL flush(numout)
293# endif 
294
295# if defined key_roam
296      !!------------------------------------------------------------------
297      !! calculate atmospheric pCO2
298      !!------------------------------------------------------------------
299      !!
300      IF (lk_oasis) THEN 
301         !! xCO2 from coupled
302         IF ( ( kt == nittrc000 ) .AND. lwp )     &
303             WRITE(numout,*) '** MEDUSA Atm xCO2 given by the UM **'
304         f_xco2a(:,:) = PCO2a_in_cpl(:,:)
305         !! Check the xCO2 from the UM is OK
306         !! piece of code moved from air-sea.F90
307         !!---
308         DO jj = 2,jpjm1
309            DO ji = 2,jpim1
310               !! OPEN wet point IF..THEN loop
311               IF (tmask(ji,jj,1) == 1) then
312                  !!! Jpalm test on atm xCO2
313                  IF ( (f_xco2a(ji,jj) .GT. 10000.0 ).OR.     &
314                       (f_xco2a(ji,jj) .LE. 0.0 ) ) THEN
315                     IF(lwp) THEN
316                        WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),       &
317                                        ' -- ji =', mig(ji),' jj = ', mjg(jj)
318                     ENDIF
319                     CALL ctl_stop( 'MEDUSA - trc_bio :',     &
320                                    'unrealistic coupled atm xCO2 ' )
321                  ENDIF
322               ENDIF
323            ENDDO
324         ENDDO
325         !!---
326      ELSEIF (lk_pi_co2) THEN
327         !! OCMIP pre-industrial xCO2
328         IF ( ( kt == nittrc000 ) .AND. lwp )     &
329             WRITE(numout,*) '** MEDUSA Atm xCO2 fixed to pre-industrial value **'
330         !! f_xco2a(:,:) = 284.725       !! CMIP5 pre-industrial pCO2
331         f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2
332      ELSE     
333         !! xCO2 from file
334         !! AXY - JPALM new interpolation scheme usinf nyear_len
335         this_y = real(nyear)
336         this_d = real(nday_year)
337         this_s = real(nsec_day)
338         !!
339         fyear = this_y + ((this_d - 1) + (this_s / (60. * 60. * 24.))) / real(nyear_len(1))
340         !!
341         IF ( ( kt == nittrc000 ) .AND. lwp ) THEN
342            WRITE(numout,*) '** MEDUSA Atm xCO2 from file **'
343            WRITE(numout,*) ' MEDUSA year      =', this_y
344            WRITE(numout,*) ' Year length      =', real(nyear_len(1))
345            WRITE(numout,*) ' MEDUSA nday_year =', this_d
346            WRITE(numout,*) ' MEDUSA nsec_day  =', this_s
347         ENDIF
348         !!
349         !! different case test
350         IF (fyear .LE. co2_yinit) THEN
351            !! before first record -- pre-industrial value
352            f_xco2a(:,:) = hist_pco2(1) 
353         ELSEIF (fyear .GE. co2_yend) THEN
354            !! after last record - continue to use the last value
355            f_xco2a(:,:) = hist_pco2(int(co2_yend - co2_yinit + 1.) )
356         ELSE
357            !! just right
358            iyr1 = int(fyear - co2_yinit) + 1
359            iyr2 = iyr1 + 1 
360            fq3 = fyear - real(iyr1) - co2_yinit + 1.
361            fq4 = ((1 - fq3) * hist_pco2(iyr1)) + (fq3 * hist_pco2(iyr2))
362            f_xco2a(:,:) = fq4
363            !!
364            IF ( ( kt == nittrc000 ) .AND. lwp ) THEN
365               WRITE(numout,*) ' MEDUSA year1     =', iyr1
366               WRITE(numout,*) ' MEDUSA year2     =', iyr2
367               WRITE(numout,*) ' xCO2 year1       =', hist_pco2(iyr1)
368               WRITE(numout,*) ' xCO2 year2       =', hist_pco2(iyr2)
369               WRITE(numout,*) ' Year2 weight     =', fq3
370            ENDIF
371         ENDIF
372      ENDIF
373 
374      !! Writing xCO2 in output on start and on the 1st tsp of each  month
375      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
376           ( nsec_month .LE. INT(rdt) ) )  THEN
377           IF ( lwp )  WRITE(numout,*) ' *** Atm xCO2 *** -- kt:', kt,       &
378                                       ';  current date:', ndastp
379          call trc_rst_dia_stat(f_xco2a(:,:), 'atm xCO2')
380      ENDIF
381# endif
382
383# if defined key_debug_medusa
384      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry'
385      IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt
386      IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000
387      CALL flush(numout)
388# endif 
389
390# if defined key_roam
391      !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every
392      !!                 month (this is hardwired as 960 timesteps but should
393      !!                 be calculated and done properly
394      !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN
395      !! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN
396      !!=============================
397      !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call :
398      !!          we don't want to call on the first time-step of all run
399      !!          submission, but only on the very first time-step, and
400      !!          then every month. So we call on nittrc000 if not
401      !!          restarted run, else if one month after last call.
402      !!          Assume one month is 30d --> 3600*24*30 : 2592000s
403      !!          try to call carb-chem at 1st month's tm-stp :
404      !!          x * 30d + 1*rdt(i.e: mod = rdt)   
405      !!          ++ need to pass carb-chem output var through restarts
406      !!If ( (kt == nitt8rc000 .AND. .NOT.ln_rsttr) .OR.          &
407      !!     ( (mod(kt*rdt,2592000.)) == rdt) THEN
408      !!=============================
409      !! (Jpalm -- updated for restartability issues)
410      !! We want this to be start of month or if starting afresh from 
411      !! climatology - marc 20/6/17
412      !!If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                         &
413      !!     ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN
414      !!=============================
415      !! Jpalm -- 15-02-2018 -- need to change 3D carb-chem call freq again.
416      !! previous call did not work, probably the (86400*mod(nn_date0,100) part
417      !! should not be in...
418      !! now use the NEMO calendar tool : nsec_month to be sure to call
419      !! at the beginning of a new month .
420      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
421           ( nsec_month .LE. INT(rdt) ) )  THEN
422           IF ( lwp )  WRITE(numout,*)                                       &         
423                              ' *** 3D carb chem call *** -- kt:', kt,       &
424                              ';  current date:', ndastp 
425         !!---------------------------------------------------------------
426         !! Calculate the carbonate chemistry for the whole ocean on the first
427         !! simulation timestep and every month subsequently; the resulting 3D
428         !! field of omega calcite is used to determine the depth of the CCD
429         !!---------------------------------------------------------------
430         CALL carb_chem( kt )
431
432      ENDIF
433# endif
434
435# if defined key_debug_medusa
436      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations'
437      CALL flush(numout)
438# endif 
439
440      !!------------------------------------------------------------------
441      !! MEDUSA has unified equation through the water column
442      !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)
443      !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1         
444      !!------------------------------------------------------------------
445      !!
446      !! NOTE: the ordering of the loops below differs from that of some other
447      !! models; looping over the vertical dimension is the outermost loop and
448      !! this complicates some calculations (e.g. storage of vertical fluxes
449      !! that can otherwise be done via a singular variable require 2D fields
450      !! here); however, these issues are relatively easily resolved, but the
451      !! loops CANNOT be reordered without potentially causing code efficiency
452      !! problems (e.g. array indexing means that reordering the loops would
453      !! require skipping between widely-spaced memory location; potentially
454      !! outside those immediately cached)
455      !!
456      !! OPEN vertical loop
457      DO jk = 1,jpk
458         !! OPEN horizontal loops
459         DO jj = 2,jpjm1
460            DO ji = 2,jpim1
461               !! OPEN wet point IF..THEN loop
462               if (tmask(ji,jj,jk) == 1) then               
463                  !!======================================================
464                  !! SETUP LOCAL GRID CELL
465                  !!======================================================
466                  !!
467                  !!------------------------------------------------------
468                  !! Some notes on grid vertical structure
469                  !! - fsdepw(ji,jj,jk) is the depth of the upper surface of
470                  !!   level jk
471                  !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of
472                  !!   level jk
473                  !! - fse3t(ji,jj,jk)  is the thickness of level jk
474                  !!------------------------------------------------------
475                  !!
476                  !! AXY (01/03/10): set up level depth (bottom of level)
477                  fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
478                  !!
479                  !! set up model tracers
480                  !! negative values of state variables are not allowed to
481                  !! contribute to the calculated fluxes
482                  !! non-diatom chlorophyll
483                  zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
484                  !! diatom chlorophyll
485                  zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
486                  !! non-diatoms
487                  zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
488                  !! diatoms
489                  zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
490                  !! diatom silicon
491                  zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
492                  !! AXY (28/01/10): probably need to take account of
493                  !! chl/biomass connection
494                  if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
495                  if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
496                  if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
497                  if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
498             !! AXY (23/01/14): duh - why did I forget diatom silicon?
499             if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
500             if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
501               ENDIF
502            ENDDO
503         ENDDO
504
505         DO jj = 2,jpjm1
506            DO ji = 2,jpim1
507               if (tmask(ji,jj,jk) == 1) then
508                  !! microzooplankton
509                  zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
510                  !! mesozooplankton
511                  zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
512                  !! detrital nitrogen
513                  zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
514                  !! dissolved inorganic nitrogen
515                  zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
516                  !! dissolved silicic acid
517                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
518                  !! dissolved "iron"
519                  zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
520               ENDIF
521            ENDDO
522         ENDDO
523
524# if defined key_roam
525         !! extra MEDUSA-2 tracers
526         DO jj = 2,jpjm1
527            DO ji = 2,jpim1
528               if (tmask(ji,jj,jk) == 1) then
529                  !! detrital carbon
530                  zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
531                  !! dissolved inorganic carbon
532                  zdic(ji,jj) = trn(ji,jj,jk,jpdic)
533                  !! alkalinity
534                  zalk(ji,jj) = trn(ji,jj,jk,jpalk)
535                  !! oxygen
536                  zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
537#  if defined key_axy_carbchem && defined key_mocsy
538                  !! phosphate via DIN and Redfield
539                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
540#  endif
541                  !!
542                  !! also need physical parameters for gas exchange
543                  !! calculations
544                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
545                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
546               ENDIF
547            ENDDO
548         ENDDO
549# else
550         !! diagnostic MEDUSA-1 detrital carbon tracer
551         DO jj = 2,jpjm1
552            DO ji = 2,jpim1
553               IF (tmask(ji,jj,jk) == 1) THEN
554                  !! implicit detrital carbon
555                  zdtc(ji,jj) = zdet(ji,jj) * xthetad
556               ENDIF
557            ENDDO
558         ENDDO
559# endif
560
561# if defined key_roam
562         !! ---------------------------------------------
563         !! JPALM -- 14-12-2017 -- Here, before any exeptionnal crazy value is
564         !!              removed, we want to tell to the Master Processor that this
565         !!              Exceptionnal value did exist
566         !!
567         Call trc_bio_check(kt, jk)
568
569         !!================================================================
570    !! AXY (03/11/17): check input fields
571         !! tracer values that exceed thresholds can cause carbonate system
572         !! failures when passed to MOCSY; temporary temperature excursions
573         !! in recent UKESM0.8 runs trigger such failures but are too short
574         !! to have physical consequences; this section checks for such
575         !! values and amends them using neighbouring values
576         !!
577         !! the check on temperature here is also carried out at the end of
578         !! each model time step and anomalies are reported in the master
579         !! ocean.output file; the error reporting below is strictly local
580         !! to the relevant ocean.output_XXXX file so will not be visible
581         !! unless all processors are reporting output
582         !!================================================================
583         !!
584         DO jj = 2,jpjm1
585            DO ji = 2,jpim1
586               if (tmask(ji,jj,jk) == 1) then
587                  !! the thresholds for the four tracers are ...
588                  IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT. 40.0  ) .OR.   &
589                       (zsal(ji,jj) .LE.  0.0) .OR. (zsal(ji,jj) .GT. 50.0  ) .OR.   &
590                       (zdic(ji,jj) .LE.  0.0) .OR. (zdic(ji,jj) .GT. 4.0E3 ) .OR.   &
591                       (zalk(ji,jj) .LE.  0.0) .OR. (zalk(ji,jj) .GT. 4.0E3 ) ) THEN
592                     !!
593                     !! all tracer values are reported in the event of any excursion
594                     IF (lwp) THEN
595                         WRITE(charout,*)  ' Tmp = ', ztmp(ji,jj)
596                         WRITE(charout2,*) ' Sal = ', zsal(ji,jj)
597                         WRITE(charout3,*) ' DIC = ', zdic(ji,jj)
598                         WRITE(charout4,*) ' Alk = ', zalk(ji,jj)
599                         WRITE(charout5,*) mig(ji), mjg(jj), jk, kt 
600                         CALL ctl_warn( 'trc_bio_medusa: carbonate chemistry WARNING:',  &
601                            TRIM(charout),TRIM(charout2),TRIM(charout3),TRIM(charout4),  & 
602                            'at i, j, k, kt:', TRIM(charout5) )
603                     ENDIF
604                     !!
605                     !! Detect, report and correct tracer excursions
606                     IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT.  40.0) )             &
607                        CALL trc_bio_exceptionnal_fix(                                         &
608                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_tem), tmask(ji-1:ji+1,jj-1:jj+1,jk),    &
609                        'Tmp', -3.0, 40.0, ztmp(ji,jj) )
610                     !!
611                     IF ( (zsal(ji,jj) .LE. 0.0) .OR. (zsal(ji,jj) .GT.  50.0)  )             &
612                        CALL trc_bio_exceptionnal_fix(                                         &
613                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_sal), tmask(ji-1:ji+1,jj-1:jj+1,jk),    &
614                        'Sal', 1.0, 50.0, zsal(ji,jj) )
615                     !!
616                     IF ( (zdic(ji,jj) .LE. 0.0) .OR. (zdic(ji,jj) .GT. 4.0E3)  )             &
617                        CALL trc_bio_exceptionnal_fix(                                         &
618                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpdic), tmask(ji-1:ji+1,jj-1:jj+1,jk),     &
619                        'DIC', 100.0, 4.0E3, zdic(ji,jj) )
620                     !!
621                     IF ( (zalk(ji,jj) .LE. 0.0) .OR. (zalk(ji,jj) .GT. 4.0E3)  )             &
622                        CALL trc_bio_exceptionnal_fix(                                         &
623                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpalk), tmask(ji-1:ji+1,jj-1:jj+1,jk),     &
624                        'Alk', 100.0, 4.0E3, zalk(ji,jj) )
625                  ENDIF
626               ENDIF
627            ENDDO
628         ENDDO
629# endif
630
631# if defined key_debug_medusa
632         DO jj = 2,jpjm1
633            DO ji = 2,jpim1
634               if (tmask(ji,jj,jk) == 1) then
635                  if (idf.eq.1) then
636                     !! AXY (15/01/10)
637                     if (trn(ji,jj,jk,jpdin).lt.0.) then
638                        IF (lwp) write (numout,*)                            &
639                           '------------------------------'
640                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',    &
641                           trn(ji,jj,jk,jpdin)
642                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',    &
643                           ji, jj, jk, kt
644                     endif
645                     if (trn(ji,jj,jk,jpsil).lt.0.) then
646                        IF (lwp) write (numout,*)                            &
647                           '------------------------------'
648                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',    &
649                           trn(ji,jj,jk,jpsil)
650                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',    &
651                           ji, jj, jk, kt
652                     endif
653#  if defined key_roam
654                     if (trn(ji,jj,jk,jpdic).lt.0.) then
655                        IF (lwp) write (numout,*)                            &
656                           '------------------------------'
657                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',    &
658                           trn(ji,jj,jk,jpdic)
659                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',    &
660                           ji, jj, jk, kt
661                     endif
662                     if (trn(ji,jj,jk,jpalk).lt.0.) then
663                        IF (lwp) write (numout,*)                            &
664                           '------------------------------'
665                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',    &
666                           trn(ji,jj,jk,jpalk)
667                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',    &
668                           ji, jj, jk, kt
669                     endif
670                     if (trn(ji,jj,jk,jpoxy).lt.0.) then
671                        IF (lwp) write (numout,*)                            &
672                           '------------------------------'
673                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',    &
674                           trn(ji,jj,jk,jpoxy)
675                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',    &
676                           ji, jj, jk, kt
677                     endif
678#  endif
679                  endif
680               ENDIF
681            ENDDO
682         ENDDO
683# endif
684# if defined key_debug_medusa
685! I'M NOT SURE THIS IS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
686!         if (idf.eq.1.AND.idfval.eq.1) then
687!            DO jj = 2,jpjm1
688!               DO ji = 2,jpim1
689!                  if (tmask(ji,jj,jk) == 1) then
690!                     !! report state variable values
691!                     IF (lwp) write (numout,*)                               &
692!                        '------------------------------'
693!                     IF (lwp) write (numout,*) 'fthk(',jk,') = ',            &
694!                        fse3t(ji,jj,jk)
695!                     IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
696!                     IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
697!                     IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
698!                     IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
699!                     IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
700!                     IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
701!                     IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
702!                     IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
703!                     IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
704#  if defined key_roam
705!                     IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
706!                     IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
707!                     IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
708!                     IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)
709#  endif
710!                  ENDIF
711!               ENDDO
712!            ENDDO
713!         endif
714# endif
715
716# if defined key_debug_medusa
717! I'M NOT SURE THIS IS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
718!         if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
719!            DO jj = 2,jpjm1
720!               DO ji = 2,jpim1
721!                  if (tmask(ji,jj,jk) == 1) then
722!                     IF (lwp) write (numout,*)                               &
723!                       '------------------------------'
724!                     IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
725!                  ENDIF
726!               ENDDO
727!            ENDDO
728!         endif
729# endif
730
731         !!---------------------------------------------------------------
732         !! Calculate air-sea gas exchange and river inputs
733         !!---------------------------------------------------------------
734         IF ( jk == 1 ) THEN
735            CALL air_sea( kt )
736         ENDIF
737
738         !!---------------------------------------------------------------
739         !! Phytoplankton growth, zooplankton grazing and miscellaneous
740         !! plankton losses.
741         !!---------------------------------------------------------------
742         CALL plankton( jk )
743
744         !!---------------------------------------------------------------
745         !! Iron chemistry and scavenging
746         !!---------------------------------------------------------------
747         CALL iron_chem_scav( jk )
748
749         !!---------------------------------------------------------------
750         !! Detritus processes
751         !!---------------------------------------------------------------
752         CALL detritus( jk, iball )
753
754         !!---------------------------------------------------------------
755         !! Updating tracers
756         !!---------------------------------------------------------------
757         CALL bio_medusa_update( kt, jk )
758
759         !!---------------------------------------------------------------
760         !! Diagnostics
761         !!---------------------------------------------------------------
762         CALL bio_medusa_diag( jk )
763
764         !!-------------------------------------------------------
765         !! 2d specific k level diags
766         !!-------------------------------------------------------
767         IF( lk_iomput ) THEN
768            CALL bio_medusa_diag_slice( jk )
769         ENDIF
770
771      !! CLOSE vertical loop
772      ENDDO
773
774      !!------------------------------------------------------------------
775      !! Final calculations for diagnostics
776      !!------------------------------------------------------------------
777      CALL bio_medusa_fin( kt )
778
779# if defined key_debug_medusa
780       IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
781       CALL flush(numout)
782# endif
783
784   END SUBROUTINE trc_bio_medusa
785
786
787
788   SUBROUTINE trc_bio_exceptionnal_fix(tiny_var, tiny_mask, var_nm, mini, maxi, varout)
789      !! JPALM (27/10/17): This function is called only when abnormal values that
790      !! could break the model's carbonate system are fed to MEDUSA
791      !!
792      !! The function calculates an average tracer value based on the 3x3 cell
793      !! neighbourhood around the abnormal cell, and reports this back
794      !!
795      !! Tracer array values are not modified, but MEDUSA uses "corrected" values
796      !! in its calculations
797      !!
798      !! temporary variables
799      REAL(wp), INTENT( in ), DIMENSION(3,3) :: tiny_var, tiny_mask
800      CHARACTER(3), INTENT( in )            :: var_nm
801      REAL(wp), INTENT( in )                 :: mini, maxi
802      REAL(wp), INTENT( out )                :: varout
803      REAL(wp)      :: sumtsn, tsnavg
804      INTEGER       :: summask
805      CHARACTER(25) :: charout1, charout2
806      CHARACTER(60) :: charout3, charout4
807      INTEGER       :: ii,ij
808   
809      !! point to the center of the 3*3 zoom-grid, to check around
810      ii = 2
811      ij = 2
812      !! Print surounding values to check if isolated Crazy value or
813      !! General error
814      IF(lwp) THEN
815          WRITE(numout,*)                                 &
816            '----------------------------------------------------------------------'
817          WRITE(numout,*)                                 &
818            'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm)
819          WRITE(numout,9100)                              &
820            3, tiny_var(ii-1,ij+1), tiny_var(ii  ,ij+1), tiny_var(ii+1,ij+1)
821          WRITE(numout,9100)                              &
822            2, tiny_var(ii-1,ij  ), tiny_var(ii  ,ij  ), tiny_var(ii+1,ij  )
823          WRITE(numout,9100)                              &
824            1, tiny_var(ii-1,ij-1), tiny_var(ii  ,ij-1), tiny_var(ii+1,ij-1)
825          WRITE(numout,*)                                 &
826            'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask'
827          WRITE(numout,9100)                              &
828            3, tiny_mask(ii-1,ij+1), tiny_mask(ii  ,ij+1), tiny_mask(ii+1,ij+1)
829          WRITE(numout,9100)                              &
830            2, tiny_mask(ii-1,ij  ), tiny_mask(ii  ,ij  ), tiny_mask(ii+1,ij  )
831          WRITE(numout,9100)                              &
832            1, tiny_mask(ii-1,ij-1), tiny_mask(ii  ,ij-1), tiny_mask(ii+1,ij-1)
833      ENDIF
834      !! Correct out of range values
835      sumtsn = ( tiny_mask(ii-1,ij+1) * tiny_var(ii-1,ij+1) ) +  &
836               ( tiny_mask(ii  ,ij+1) * tiny_var(ii  ,ij+1) ) +  &
837               ( tiny_mask(ii+1,ij+1) * tiny_var(ii+1,ij+1) ) +  &
838               ( tiny_mask(ii-1,ij  ) * tiny_var(ii-1,ij  ) ) +  &
839               ( tiny_mask(ii+1,ij  ) * tiny_var(ii+1,ij  ) ) +  &
840               ( tiny_mask(ii-1,ij-1) * tiny_var(ii-1,ij-1) ) +  &
841               ( tiny_mask(ii  ,ij-1) * tiny_var(ii  ,ij-1) ) +  &
842               ( tiny_mask(ii+1,ij-1) * tiny_var(ii+1,ij-1) )
843      !!
844      summask = tiny_mask(ii-1,ij+1) + tiny_mask(ii  ,ij+1)   +  &
845                tiny_mask(ii+1,ij+1) + tiny_mask(ii-1,ij  )   +  &
846                tiny_mask(ii+1,ij  ) + tiny_mask(ii-1,ij-1)   +  &
847                tiny_mask(ii  ,ij-1) + tiny_mask(ii+1,ij-1)
848      !!
849      IF ( summask .GT. 0 ) THEN
850         tsnavg = ( sumtsn / summask )
851         varout = MAX( MIN( maxi, tsnavg), mini )
852      ELSE   
853         IF (ztmp(ii,ij) .LT. mini )  varout = mini
854         IF (ztmp(ii,ij) .GT. maxi )  varout = maxi
855      ENDIF
856      !!
857      IF (lwp) THEN
858          WRITE(charout1,9200) tiny_var(ii,ij)
859          WRITE(charout2,9200) varout
860          WRITE(charout3,*) ' ', charout1, ' -> ', charout2
861          WRITE(charout4,*) ' Tracer: ', trim(var_nm)
862      !!
863          WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **'
864          WRITE(numout,*) charout4 
865          WRITE(numout,*) charout3
866          WRITE(numout,*) '----------------------------------------------------------------------'
867          WRITE(numout,*) ' '
868      ENDIF
869
8709100  FORMAT('Row:', i1, '  ', e16.6, ' ', e16.6, ' ', e16.6)
8719200  FORMAT(e16.6)
872
873   END SUBROUTINE trc_bio_exceptionnal_fix 
874
875   SUBROUTINE trc_bio_check(kt, jk)
876      !!-----------------------------------
877      !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure
878      !!                     problem. The model is now able to correct a local
879      !!                     crazy value. but does it silently.
880      !!                     We need to spread the word to the master processor. we
881      !!                     don't want the model  to correct values without telling us
882      !!                     This module will tell at least when crazy DIC or
883      !!                     ALK appears.
884      !!-------------------------------------
885      REAL(wp)              :: zmax, zmin    ! temporary scalars
886      INTEGER               :: ji,jj         ! dummy loop
887      INTEGER               :: ii,ij         ! temporary scalars
888      INTEGER, DIMENSION(2) :: ilocs         !
889      INTEGER, INTENT( in ) :: kt, jk
890      !!
891      !!==========================
892      !! DIC Check
893      zmax =  -5.0  ! arbitrary  low maximum value
894      zmin =  4.0E4  ! arbitrary high minimum value
895      DO jj = 2, jpjm1
896         DO ji = 2,jpim1
897            IF( tmask(ji,jj,1) == 1) THEN
898               zmax = MAX(zmax,zdic(ji,jj))     ! find local maximum
899               zmin = MIN(zmin,zdic(ji,jj))     ! find local minimum
900            ENDIF
901         END DO
902      END DO
903      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain
904      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain
905      !
906      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem
907         IF (lk_mpp) THEN
908            CALL mpp_maxloc ( zdic(:,:),tmask(:,:,1), zmax, ii,ij )
909         ELSE
910            ilocs = MAXLOC( zdic(:,:), mask = tmask(:,:,1) == 1. )
911            ii = ilocs(1) + nimpp - 1
912            ij = ilocs(2) + njmpp - 1
913         ENDIF
914         !
915         IF(lwp) THEN
916            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
917            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration > 4000 '
918            WRITE(numout,9600) kt, zmax, ii, ij, jk
919            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
920         ENDIF
921      ENDIF
922      !
923      IF( zmin .LE. 0.0) THEN  ! we've got a problem
924         IF (lk_mpp) THEN
925            CALL mpp_minloc ( zdic(:,:),tmask(:,:,1), zmin, ii,ij )
926         ELSE
927            ilocs = MINLOC( zdic(:,:), mask = tmask(:,:,1) == 1. )
928            ii = ilocs(1) + nimpp - 1
929            ij = ilocs(2) + njmpp - 1
930         ENDIF
931         !
932         IF(lwp) THEN
933            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
934            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration <= 0 '
935            WRITE(numout,9700) kt, zmin, ii, ij, jk
936            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
937         ENDIF
938      ENDIF
939      !!
940      !!==========================
941      !! ALKALINITY Check
942      zmax =  -5.0  ! arbitrary  low maximum value
943      zmin =  4.0E4  ! arbitrary high minimum value
944      DO jj = 2, jpjm1
945         DO ji = 2,jpim1
946            IF( tmask(ji,jj,1) == 1) THEN
947               zmax = MAX(zmax,zalk(ji,jj))     ! find local maximum
948               zmin = MIN(zmin,zalk(ji,jj))     ! find local minimum
949            ENDIF
950         END DO
951      END DO
952      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain
953      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain
954      !
955      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem
956         IF (lk_mpp) THEN
957            CALL mpp_maxloc ( zalk(:,:),tmask(:,:,1), zmax, ii,ij )
958         ELSE
959            ilocs = MAXLOC( zalk(:,:), mask = tmask(:,:,1) == 1. )
960            ii = ilocs(1) + nimpp - 1
961            ij = ilocs(2) + njmpp - 1
962         ENDIF
963         !
964         IF(lwp) THEN
965            WRITE(numout,*) 'trc_bio:tracer anomaly: *****     WARNING    *****'
966            WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration > 4000 '
967            WRITE(numout,9800) kt, zmax, ii, ij, jk
968            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
969         ENDIF
970      ENDIF
971      !
972      IF( zmin .LE. 0.0) THEN  ! we've got a problem
973         IF (lk_mpp) THEN
974            CALL mpp_minloc ( zalk(:,:),tmask(:,:,1), zmin, ii,ij )
975         ELSE
976            ilocs = MINLOC( zalk(:,:), mask = tmask(:,:,1) == 1. )
977            ii = ilocs(1) + nimpp - 1
978            ij = ilocs(2) + njmpp - 1
979         ENDIF
980         !
981         IF(lwp) THEN
982            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
983            WRITE(numout,*) 'trc_bio:tracer anomaly:  ALK concentration <= 0 '
984            WRITE(numout,9900) kt, zmin, ii, ij, jk
985            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
986         ENDIF
987      ENDIF
988
989
9909600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5)
9919700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5)
9929800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5)
9939900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5)
994
995   END SUBROUTINE trc_bio_check
996
997
998#else
999   !!=====================================================================
1000   !!  Dummy module :                                   No MEDUSA bio-model
1001   !!=====================================================================
1002CONTAINS
1003   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
1004      IMPLICIT NONE
1005      INTEGER, INTENT( in ) ::   kt
1006      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
1007   END SUBROUTINE trc_bio_medusa
1008#endif 
1009
1010   !!=====================================================================
1011END MODULE  trcbio_medusa
Note: See TracBrowser for help on using the repository browser.