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

Last change on this file since 9285 was 9285, checked in by jpalmier, 3 years ago

JPALM —26-01-2018 — change MEDUSA atm co2 input method. now reads a file a variable length

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