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

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

source: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 9247

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

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