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_fix_cpp_keys/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_fix_cpp_keys/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 10004

Last change on this file since 10004 was 10004, checked in by frrh, 4 years ago

Various fixes

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