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/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 9292

Last change on this file since 9292 was 9292, checked in by dford, 6 years ago

Merge in changes from dev_r5518_GO6_package_asm_surf_bgc_v2.

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