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

source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 10047

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

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