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

source: branches/UKMO/dev_r5518_GO6_fix_key_comp/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 9991

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

Fixes to allow MEDUSA to compile with C1D without
the need for multiple (apparently) unrelated CPP keys
merely to satisfy spurious code interdependencies.

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