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

source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 9327

Last change on this file since 9327 was 9327, checked in by jpalmier, 6 years ago

Jpalm -- GMED 380 -- the test used to call the 3D carb chem every month was not working.

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