source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 11990

Last change on this file since 11990 was 11990, checked in by dford, 8 months ago

Get the nitrogen balancing working with 3D chlorophyll increments.

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