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

Last change on this file was 10149, checked in by frrh, 6 years ago

Met Office GMED ticket 379: Merged David Ford's MEDUSA assimilation changes
using command:

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