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

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

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

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

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