New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcbio_medusa.F90 in branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

Last change on this file was 11738, checked in by marc, 5 years ago

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