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

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

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

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

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