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 @ 10045

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

Andrew's changes to add the OMIP double_DIC (activated with key_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,                           &
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, e3t_n,                   &
86                                            gdept_0, gdept_n,               &
87                                            gdepw_0, gdepw_n,               &
88                                            nday_year, nsec_day,            &
89                                            nyear, nyear_len, ndastp,       &
90                                            nsec_month,                     &
91                                            rdt, tmask, mig, mjg, nimpp,    &
92                                            njmpp 
93      USE in_out_manager,             ONLY: lwp, numout, nn_date0
94# if defined key_iomput
95      USE iom,                        ONLY: lk_iomput
96# endif
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_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            ENDDO
723         ENDDO
724# endif
725# if defined key_debug_medusa
726! I'M NOT SURE THIS IS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
727!         if (idf.eq.1.AND.idfval.eq.1) then
728!            DO jj = 2,jpjm1
729!               DO ji = 2,jpim1
730!                  if (tmask(ji,jj,jk) == 1) then
731!                     !! report state variable values
732!                     IF (lwp) write (numout,*)                               &
733!                        '------------------------------'
734!                     IF (lwp) write (numout,*) 'fthk(',jk,') = ',            &
735!                        fse3t(ji,jj,jk)
736!                     IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
737!                     IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
738!                     IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
739!                     IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
740!                     IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
741!                     IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
742!                     IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
743!                     IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
744!                     IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
745#  if defined key_roam
746!                     IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
747!                     IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
748!                     IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
749!                     IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)
750#  endif
751!                  ENDIF
752!               ENDDO
753!            ENDDO
754!         endif
755# endif
756
757# if defined key_debug_medusa
758! I'M NOT SURE THIS IS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
759!         if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
760!            DO jj = 2,jpjm1
761!               DO ji = 2,jpim1
762!                  if (tmask(ji,jj,jk) == 1) then
763!                     IF (lwp) write (numout,*)                               &
764!                       '------------------------------'
765!                     IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
766!                  ENDIF
767!               ENDDO
768!            ENDDO
769!         endif
770# endif
771
772         !!---------------------------------------------------------------
773         !! Calculate air-sea gas exchange and river inputs
774         !!---------------------------------------------------------------
775         IF ( jk == 1 ) THEN
776            CALL air_sea( kt )
777         ENDIF
778
779         !!---------------------------------------------------------------
780         !! Phytoplankton growth, zooplankton grazing and miscellaneous
781         !! plankton losses.
782         !!---------------------------------------------------------------
783         CALL plankton( jk )
784
785         !!---------------------------------------------------------------
786         !! Iron chemistry and scavenging
787         !!---------------------------------------------------------------
788         CALL iron_chem_scav( jk )
789
790         !!---------------------------------------------------------------
791         !! Detritus processes
792         !!---------------------------------------------------------------
793         CALL detritus( jk, iball )
794
795         !!---------------------------------------------------------------
796         !! Updating tracers
797         !!---------------------------------------------------------------
798         CALL bio_medusa_update( kt, jk )
799
800         !!---------------------------------------------------------------
801         !! Diagnostics
802         !!---------------------------------------------------------------
803         CALL bio_medusa_diag( jk )
804
805         !!-------------------------------------------------------
806         !! 2d specific k level diags
807         !!-------------------------------------------------------
808         IF( lk_iomput ) THEN
809            CALL bio_medusa_diag_slice( jk )
810         ENDIF
811
812      !! CLOSE vertical loop
813      ENDDO
814
815      !!------------------------------------------------------------------
816      !! Final calculations for diagnostics
817      !!------------------------------------------------------------------
818      CALL bio_medusa_fin( kt )
819
820# if defined key_debug_medusa
821       IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
822       CALL flush(numout)
823# endif
824
825   END SUBROUTINE trc_bio_medusa
826
827
828
829   SUBROUTINE trc_bio_exceptionnal_fix(tiny_var, tiny_mask, var_nm, mini, maxi, varout)
830      !! JPALM (27/10/17): This function is called only when abnormal values that
831      !! could break the model's carbonate system are fed to MEDUSA
832      !!
833      !! The function calculates an average tracer value based on the 3x3 cell
834      !! neighbourhood around the abnormal cell, and reports this back
835      !!
836      !! Tracer array values are not modified, but MEDUSA uses "corrected" values
837      !! in its calculations
838      !!
839      !! temporary variables
840      REAL(wp), INTENT( in ), DIMENSION(3,3) :: tiny_var, tiny_mask
841      CHARACTER(25), INTENT( in )            :: var_nm
842      REAL(wp), INTENT( in )                 :: mini, maxi
843      REAL(wp), INTENT( out )                :: varout
844      REAL(wp)      :: sumtsn, tsnavg
845      INTEGER       :: summask
846      CHARACTER(25) :: charout1, charout2
847      CHARACTER(60) :: charout3, charout4
848      INTEGER       :: ii,ij
849   
850      !! point to the center of the 3*3 zoom-grid, to check around
851      ii = 2
852      ij = 2
853      !! Print surounding values to check if isolated Crazy value or
854      !! General error
855      IF(lwp) THEN
856          WRITE(numout,*)                                 &
857            '----------------------------------------------------------------------'
858          WRITE(numout,*)                                 &
859            'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm)
860          WRITE(numout,9100)                              &
861            3, tiny_var(ii-1,ij+1), tiny_var(ii  ,ij+1), tiny_var(ii+1,ij+1)
862          WRITE(numout,9100)                              &
863            2, tiny_var(ii-1,ij  ), tiny_var(ii  ,ij  ), tiny_var(ii+1,ij  )
864          WRITE(numout,9100)                              &
865            1, tiny_var(ii-1,ij-1), tiny_var(ii  ,ij-1), tiny_var(ii+1,ij-1)
866          WRITE(numout,*)                                 &
867            'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask'
868          WRITE(numout,9100)                              &
869            3, tiny_mask(ii-1,ij+1), tiny_mask(ii  ,ij+1), tiny_mask(ii+1,ij+1)
870          WRITE(numout,9100)                              &
871            2, tiny_mask(ii-1,ij  ), tiny_mask(ii  ,ij  ), tiny_mask(ii+1,ij  )
872          WRITE(numout,9100)                              &
873            1, tiny_mask(ii-1,ij-1), tiny_mask(ii  ,ij-1), tiny_mask(ii+1,ij-1)
874      ENDIF
875      !! Correct out of range values
876      sumtsn = ( tiny_mask(ii-1,ij+1) * tiny_var(ii-1,ij+1) ) +  &
877               ( tiny_mask(ii  ,ij+1) * tiny_var(ii  ,ij+1) ) +  &
878               ( tiny_mask(ii+1,ij+1) * tiny_var(ii+1,ij+1) ) +  &
879               ( tiny_mask(ii-1,ij  ) * tiny_var(ii-1,ij  ) ) +  &
880               ( tiny_mask(ii+1,ij  ) * tiny_var(ii+1,ij  ) ) +  &
881               ( tiny_mask(ii-1,ij-1) * tiny_var(ii-1,ij-1) ) +  &
882               ( tiny_mask(ii  ,ij-1) * tiny_var(ii  ,ij-1) ) +  &
883               ( tiny_mask(ii+1,ij-1) * tiny_var(ii+1,ij-1) )
884      !!
885      summask = tiny_mask(ii-1,ij+1) + tiny_mask(ii  ,ij+1)   +  &
886                tiny_mask(ii+1,ij+1) + tiny_mask(ii-1,ij  )   +  &
887                tiny_mask(ii+1,ij  ) + tiny_mask(ii-1,ij-1)   +  &
888                tiny_mask(ii  ,ij-1) + tiny_mask(ii+1,ij-1)
889      !!
890      IF ( summask .GT. 0 ) THEN
891         tsnavg = ( sumtsn / summask )
892         varout = MAX( MIN( maxi, tsnavg), mini )
893      ELSE   
894         IF (ztmp(ii,ij) .LT. mini )  varout = mini
895         IF (ztmp(ii,ij) .GT. maxi )  varout = maxi
896      ENDIF
897      !!
898      IF (lwp) THEN
899          WRITE(charout1,9200) tiny_var(ii,ij)
900          WRITE(charout2,9200) varout
901          WRITE(charout3,*) ' ', charout1, ' -> ', charout2
902          WRITE(charout4,*) ' Tracer: ', trim(var_nm)
903      !!
904          WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **'
905          WRITE(numout,*) charout4 
906          WRITE(numout,*) charout3
907          WRITE(numout,*) '----------------------------------------------------------------------'
908          WRITE(numout,*) ' '
909      ENDIF
910
9119100  FORMAT('Row:', i1, '  ', e16.6, ' ', e16.6, ' ', e16.6)
9129200  FORMAT(e16.6)
913
914   END SUBROUTINE trc_bio_exceptionnal_fix 
915
916   SUBROUTINE trc_bio_check(kt, jk)
917      !!-----------------------------------
918      !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure
919      !!                     problem. The model is now able to correct a local
920      !!                     crazy value. but does it silently.
921      !!                     We need to spread the word to the master processor. we
922      !!                     don't want the model  to correct values without telling us
923      !!                     This module will tell at least when crazy DIC or
924      !!                     ALK appears.
925      !!-------------------------------------
926      REAL(wp)              :: zmax, zmin    ! temporary scalars
927      INTEGER               :: ji,jj         ! dummy loop
928      INTEGER               :: ii,ij         ! temporary scalars
929      INTEGER, DIMENSION(2) :: ilocs         !
930      INTEGER, INTENT( in ) :: kt, jk
931      !!
932      !!==========================
933      !! DIC Check
934      zmax =  -5.0  ! arbitrary  low maximum value
935      zmin =  4.0E4  ! arbitrary high minimum value
936      DO jj = 2, jpjm1
937         DO ji = 2,jpim1
938            IF( tmask(ji,jj,1) == 1) THEN
939               zmax = MAX(zmax,zdic(ji,jj))     ! find local maximum
940               zmin = MIN(zmin,zdic(ji,jj))     ! find local minimum
941            ENDIF
942         END DO
943      END DO
944      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain
945      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain
946      !
947      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem
948         IF (lk_mpp) THEN
949            CALL mpp_maxloc ( zdic(:,:),tmask(:,:,1), zmax, ii,ij )
950         ELSE
951            ilocs = MAXLOC( zdic(:,:), mask = tmask(:,:,1) == 1. )
952            ii = ilocs(1) + nimpp - 1
953            ij = ilocs(2) + njmpp - 1
954         ENDIF
955         !
956         IF(lwp) THEN
957            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
958            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration > 4000 '
959            WRITE(numout,9600) kt, zmax, ii, ij, jk
960            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
961         ENDIF
962      ENDIF
963      !
964      IF( zmin .LE. 0.0) THEN  ! we've got a problem
965         IF (lk_mpp) THEN
966            CALL mpp_minloc ( zdic(:,:),tmask(:,:,1), zmin, ii,ij )
967         ELSE
968            ilocs = MINLOC( zdic(:,:), mask = tmask(:,:,1) == 1. )
969            ii = ilocs(1) + nimpp - 1
970            ij = ilocs(2) + njmpp - 1
971         ENDIF
972         !
973         IF(lwp) THEN
974            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
975            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration <= 0 '
976            WRITE(numout,9700) kt, zmin, ii, ij, jk
977            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
978         ENDIF
979      ENDIF
980      !!
981      !!==========================
982      !! ALKALINITY Check
983      zmax =  -5.0  ! arbitrary  low maximum value
984      zmin =  4.0E4  ! arbitrary high minimum value
985      DO jj = 2, jpjm1
986         DO ji = 2,jpim1
987            IF( tmask(ji,jj,1) == 1) THEN
988               zmax = MAX(zmax,zalk(ji,jj))     ! find local maximum
989               zmin = MIN(zmin,zalk(ji,jj))     ! find local minimum
990            ENDIF
991         END DO
992      END DO
993      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain
994      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain
995      !
996      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem
997         IF (lk_mpp) THEN
998            CALL mpp_maxloc ( zalk(:,:),tmask(:,:,1), zmax, ii,ij )
999         ELSE
1000            ilocs = MAXLOC( zalk(:,:), mask = tmask(:,:,1) == 1. )
1001            ii = ilocs(1) + nimpp - 1
1002            ij = ilocs(2) + njmpp - 1
1003         ENDIF
1004         !
1005         IF(lwp) THEN
1006            WRITE(numout,*) 'trc_bio:tracer anomaly: *****     WARNING    *****'
1007            WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration > 4000 '
1008            WRITE(numout,9800) kt, zmax, ii, ij, jk
1009            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
1010         ENDIF
1011      ENDIF
1012      !
1013      IF( zmin .LE. 0.0) THEN  ! we've got a problem
1014         IF (lk_mpp) THEN
1015            CALL mpp_minloc ( zalk(:,:),tmask(:,:,1), zmin, ii,ij )
1016         ELSE
1017            ilocs = MINLOC( zalk(:,:), mask = tmask(:,:,1) == 1. )
1018            ii = ilocs(1) + nimpp - 1
1019            ij = ilocs(2) + njmpp - 1
1020         ENDIF
1021         !
1022         IF(lwp) THEN
1023            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
1024            WRITE(numout,*) 'trc_bio:tracer anomaly:  ALK concentration <= 0 '
1025            WRITE(numout,9900) kt, zmin, ii, ij, jk
1026            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
1027         ENDIF
1028      ENDIF
1029
1030
10319600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5)
10329700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5)
10339800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5)
10349900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5)
1035
1036   END SUBROUTINE trc_bio_check
1037
1038
1039#else
1040   !!=====================================================================
1041   !!  Dummy module :                                   No MEDUSA bio-model
1042   !!=====================================================================
1043CONTAINS
1044   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
1045      INTEGER, INTENT( in ) ::   kt
1046      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
1047   END SUBROUTINE trc_bio_medusa
1048#endif 
1049
1050   !!=====================================================================
1051END MODULE  trcbio_medusa
Note: See TracBrowser for help on using the repository browser.