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

Last change on this file since 9309 was 9309, checked in by jpalmier, 2 years ago

JPALM — make clearer in the code and output which atm xCO2 source is used

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