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

source: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 9073

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

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