source: branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8643

Last change on this file since 8643 was 8643, checked in by jpalmier, 4 years ago

JPALM — 19-10-17 — Add ctl_warn for master node output

File size: 36.5 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   !!----------------------------------------------------------------------
26   !!
27#if defined key_roam
28   !!----------------------------------------------------------------------
29   !! Updates for the ROAM project include:
30   !!   - addition of DIC, alkalinity, detrital carbon and oxygen tracers
31   !!   - addition of air-sea fluxes of CO2 and oxygen
32   !!   - periodic (monthly) calculation of full 3D carbonate chemistry
33   !!   - detrital C:N ratio now free to evolve dynamically
34   !!   - benthic storage pools
35   !!
36   !! Opportunity also taken to add functionality:
37   !!   - switch for Liebig Law (= most-limiting) nutrient uptake
38   !!   - switch for accelerated seafloor detritus remineralisation
39   !!   - switch for fast -> slow detritus transfer at seafloor
40   !!   - switch for ballast vs. Martin vs. Henson fast detritus remin.
41   !!   - per GMD referee remarks, xfdfrac3 introduced for grazed PDS
42   !!----------------------------------------------------------------------
43#endif
44   !!
45#if defined key_mocsy
46   !!----------------------------------------------------------------------
47   !! Updates with the addition of MOCSY include:
48   !!   - option to use PML or MOCSY carbonate chemistry (the latter is
49   !!     preferred)
50   !!   - central calculation of gas transfer velocity, f_kw660; previously
51   !!     this was done separately for CO2 and O2 with predictable results
52   !!   - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux
53   !!     calculations and to those for O2 air-sea flux
54   !!   - extra diagnostics included for MOCSY
55   !!----------------------------------------------------------------------
56#endif
57   !!
58#if defined key_medusa
59   !!----------------------------------------------------------------------
60   !!                                        MEDUSA bio-model
61   !!----------------------------------------------------------------------
62   !!   trc_bio_medusa        : 
63   !!----------------------------------------------------------------------
64      !! AXY (30/01/14): necessary to find NaNs on HECTOR
65      USE, INTRINSIC :: ieee_arithmetic 
66
67      USE bio_medusa_mod,             ONLY: b0, fdep1,                      & 
68                                            ibenthic, idf, idfval,          &
69# if defined key_roam
70                                            f_xco2a,                        &
71                                            zalk, zdic, zoxy, zsal, ztmp,   &
72# endif
73# if defined key_mocsy
74                                            zpho,                           &
75# endif
76                                            zchd, zchn, zdet, zdin, zdtc,   &
77                                            zfer, zpds, zphd, zphn, zsil,   &
78                                            zzme, zzmi
79      USE dom_oce,                    ONLY: e3t_0, e3t_n,                   &
80                                            gdept_0, gdept_n,               &
81                                            gdepw_0, gdepw_n,               &
82                                            nday_year, nsec_day, nyear,     &
83                                            rdt, tmask, mig, mjg
84      USE in_out_manager,             ONLY: lwp, numout, nn_date0
85# if defined key_iomput
86      USE iom,                        ONLY: lk_iomput
87# endif
88      USE lbclnk,                     ONLY: lbc_lnk
89      USE lib_mpp,                    ONLY: ctl_stop, ctl_warn
90      USE oce,                        ONLY: tsb, tsn
91      USE par_kind,                   ONLY: wp
92      USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     &
93                                            jpdic, jpdin, jpdtc, jpfer,     &
94                                            jpoxy, jppds, jpphd, jpphn,     &
95                                            jpsil, jpzme, jpzmi
96      USE par_oce,                    ONLY: jp_sal, jp_tem, jpi, jpim1,     &
97                                            jpj, jpjm1, jpk
98      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
99      USE sbc_oce,                    ONLY: lk_oasis
100      USE sms_medusa,                 ONLY: hist_pco2
101      USE trc,                        ONLY: ln_rsttr, nittrc000, trn
102      USE bio_medusa_init_mod,        ONLY: bio_medusa_init
103      USE carb_chem_mod,              ONLY: carb_chem
104      USE air_sea_mod,                ONLY: air_sea
105      USE plankton_mod,               ONLY: plankton
106      USE iron_chem_scav_mod,         ONLY: iron_chem_scav
107      USE detritus_mod,               ONLY: detritus
108      USE bio_medusa_update_mod,      ONLY: bio_medusa_update
109      USE bio_medusa_diag_mod,        ONLY: bio_medusa_diag
110      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
111      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
112
113      IMPLICIT NONE
114      PRIVATE
115     
116      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
117
118   !!* Substitution
119#  include "domzgr_substitute.h90"
120   !!----------------------------------------------------------------------
121   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
122   !! $Id$
123   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
124   !!----------------------------------------------------------------------
125
126CONTAINS
127
128   SUBROUTINE trc_bio_medusa( kt )
129      !!------------------------------------------------------------------
130      !!                     ***  ROUTINE trc_bio  ***
131      !!
132      !! ** Purpose : compute the now trend due to biogeochemical processes
133      !!              and add it to the general trend of passive tracers
134      !!              equations
135      !!
136      !! ** Method  : each now biological flux is calculated in function of
137      !!              now concentrations of tracers.
138      !!              depending on the tracer, these fluxes are sources or
139      !!              sinks.
140      !!              The total of the sources and sinks for each tracer
141      !!              is added to the general trend.
142      !!       
143      !!                      tra = tra + zf...tra - zftra...
144      !!                                     |         |
145      !!                                     |         |
146      !!                                  source      sink
147      !!       
148      !!              IF 'key_trc_diabio' defined , the biogeochemical trends
149      !!              for passive tracers are saved for futher diagnostics.
150      !!------------------------------------------------------------------
151      !!
152      !!
153      !!------------------------------------------------------------------
154      !! Variable conventions
155      !!------------------------------------------------------------------
156      !!
157      !! names: z*** - state variable
158      !!        f*** - function (or temporary variable used in part of
159      !!               a function)
160      !!        x*** - parameter
161      !!        b*** - right-hand part (sources and sinks)
162      !!        i*** - integer variable (usually used in yes/no flags)
163      !!
164      !! time (integer timestep)
165      INTEGER, INTENT( in ) ::    kt
166      !!
167      !! spatial array indices
168      INTEGER  ::    ji,jj,jk,jn
169      !!
170      INTEGER  ::    iball
171# if defined key_roam
172      !!
173      INTEGER  ::    iyr1, iyr2
174      !!
175# endif
176      !!
177      !! temporary variables
178      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
179      !!
180      !! T and S check temporary variable
181      REAL(wp) ::    sumtsn, tsnavg
182      INTEGER  ::    summask
183      !!
184      !!------------------------------------------------------------------
185
186# if defined key_debug_medusa
187      IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined'
188      CALL flush(numout)
189# endif 
190
191      !! AXY (20/11/14): alter this to report on first MEDUSA call
192      !! IF( kt == nit000 ) THEN
193      IF( kt == nittrc000 ) THEN
194         IF(lwp) WRITE(numout,*)
195         IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model'
196         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
197    IF(lwp) WRITE(numout,*) ' kt =',kt
198      ENDIF
199
200      !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
201      ibenthic = 1
202
203      !!------------------------------------------------------------------
204      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency
205      !! terms of all biological equations to 0.
206      !!------------------------------------------------------------------
207      !!
208      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit
209      !!                 the bill for now; another item on the things-to-sort-
210      !!     out-in-the-future list ...
211# if defined key_kill_medusa
212      b0 = 0.
213# else
214      b0 = 1.
215# endif
216      !!------------------------------------------------------------------
217      !! fast detritus ballast scheme (0 = no; 1 = yes)
218      !! alternative to ballast scheme is same scheme but with no ballast
219      !! protection (not dissimilar to Martin et al., 1987)
220      !!------------------------------------------------------------------
221      !!
222      iball = 1
223
224      !!------------------------------------------------------------------
225      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
226      !! these should *only* be used in 1D since they give comprehensive
227      !! output for ecological functions in the model; primarily used in
228      !! debugging
229      !!------------------------------------------------------------------
230      !!
231      idf    = 0
232      !!
233      !! timer mechanism
234      if (kt/120*120.eq.kt) then
235         idfval = 1
236      else
237         idfval = 0
238      endif
239
240      !!------------------------------------------------------------------
241      !! Initialise arrays to zero and set up arrays for diagnostics
242      !!------------------------------------------------------------------
243      CALL bio_medusa_init( kt )
244
245# if defined key_axy_nancheck
246       DO jn = jp_msa0,jp_msa1
247          !! fq0 = MINVAL(trn(:,:,:,jn))
248          !! fq1 = MAXVAL(trn(:,:,:,jn))
249          fq2 = SUM(trn(:,:,:,jn))
250          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK',     &
251          !!                kt, jn, fq0, fq1, fq2
252          !! AXY (30/01/14): much to our surprise, the next line doesn't
253          !!                 work on HECTOR and has been replaced here with
254          !!                 a specialist routine
255          !! if (fq2 /= fq2 ) then
256          if ( ieee_is_nan( fq2 ) ) then
257             !! there's a NaN here
258             if (lwp) write(numout,*) 'NAN detected in field', jn,           &
259                                      'at time', kt, 'at position:'
260             DO jk = 1,jpk
261                DO jj = 1,jpj
262                   DO ji = 1,jpi
263                      !! AXY (30/01/14): "isnan" problem on HECTOR
264                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then
265                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
266                         if (lwp) write (numout,'(a,1pe12.2,4i6)')           &
267                            'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn
268                      endif
269                   enddo
270                enddo
271             enddo
272             CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
273          endif
274       ENDDO
275       CALL flush(numout)
276# endif
277
278# if defined key_debug_medusa
279      IF (lwp) write (numout,*)                                              &
280                     'trc_bio_medusa: variables initialised and checked'
281      CALL flush(numout)
282# endif 
283
284# if defined key_roam
285      !!------------------------------------------------------------------
286      !! calculate atmospheric pCO2
287      !!------------------------------------------------------------------
288      !!
289      !! what's atmospheric pCO2 doing? (data start in 1859)
290      iyr1  = nyear - 1859 + 1
291      iyr2  = iyr1 + 1
292      if (iyr1 .le. 1) then
293         !! before 1860
294         f_xco2a(:,:) = hist_pco2(1)
295      elseif (iyr2 .ge. 242) then
296         !! after 2099
297         f_xco2a(:,:) = hist_pco2(242)
298      else
299         !! just right
300         fq0 = hist_pco2(iyr1)
301         fq1 = hist_pco2(iyr2)
302         fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0)
303         !! AXY (14/06/12): tweaked to make more sense (and be correct)
304#  if defined key_bs_axy_yrlen
305         !! bugfix: for 360d year with HadGEM2-ES forcing
306         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 
307#  else
308         !! original use of 365 days (not accounting for leap year or
309         !! 360d year)
310         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0
311#  endif
312         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
313         f_xco2a(:,:) = fq4
314      endif
315#  if defined key_axy_pi_co2
316      !! OCMIP pre-industrial pCO2
317      !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2
318      f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2
319#  endif
320      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
321      !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day)
322      !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)
323      !! AXY (29/01/14): remove surplus diagnostics
324      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0
325      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1
326      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
327      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
328      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
329# endif
330
331# if defined key_debug_medusa
332      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry'
333      IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt
334      IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000
335      CALL flush(numout)
336# endif 
337
338# if defined key_roam
339      !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every
340      !!                 month (this is hardwired as 960 timesteps but should
341      !!                 be calculated and done properly
342      !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN
343      !! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN
344      !!=============================
345      !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call :
346      !!          we don't want to call on the first time-step of all run
347      !!          submission, but only on the very first time-step, and
348      !!          then every month. So we call on nittrc000 if not
349      !!          restarted run, else if one month after last call.
350      !!          Assume one month is 30d --> 3600*24*30 : 2592000s
351      !!          try to call carb-chem at 1st month's tm-stp :
352      !!          x * 30d + 1*rdt(i.e: mod = rdt)   
353      !!          ++ need to pass carb-chem output var through restarts
354      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
355           ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN
356         !!---------------------------------------------------------------
357         !! Calculate the carbonate chemistry for the whole ocean on the first
358         !! simulation timestep and every month subsequently; the resulting 3D
359         !! field of omega calcite is used to determine the depth of the CCD
360         !!---------------------------------------------------------------
361         CALL carb_chem( kt )
362
363      ENDIF
364# endif
365
366# if defined key_debug_medusa
367      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations'
368      CALL flush(numout)
369# endif 
370
371      !!------------------------------------------------------------------
372      !! MEDUSA has unified equation through the water column
373      !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)
374      !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1         
375      !!------------------------------------------------------------------
376      !!
377      !! NOTE: the ordering of the loops below differs from that of some other
378      !! models; looping over the vertical dimension is the outermost loop and
379      !! this complicates some calculations (e.g. storage of vertical fluxes
380      !! that can otherwise be done via a singular variable require 2D fields
381      !! here); however, these issues are relatively easily resolved, but the
382      !! loops CANNOT be reordered without potentially causing code efficiency
383      !! problems (e.g. array indexing means that reordering the loops would
384      !! require skipping between widely-spaced memory location; potentially
385      !! outside those immediately cached)
386      !!
387      !! OPEN vertical loop
388      DO jk = 1,jpk
389         !! OPEN horizontal loops
390         DO jj = 2,jpjm1
391            DO ji = 2,jpim1
392               !! OPEN wet point IF..THEN loop
393               if (tmask(ji,jj,jk) == 1) then               
394                  !!======================================================
395                  !! SETUP LOCAL GRID CELL
396                  !!======================================================
397                  !!
398                  !!------------------------------------------------------
399                  !! Some notes on grid vertical structure
400                  !! - fsdepw(ji,jj,jk) is the depth of the upper surface of
401                  !!   level jk
402                  !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of
403                  !!   level jk
404                  !! - fse3t(ji,jj,jk)  is the thickness of level jk
405                  !!------------------------------------------------------
406                  !!
407                  !! AXY (01/03/10): set up level depth (bottom of level)
408                  fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
409                  !!
410                  !! set up model tracers
411                  !! negative values of state variables are not allowed to
412                  !! contribute to the calculated fluxes
413                  !! non-diatom chlorophyll
414                  zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
415                  !! diatom chlorophyll
416                  zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
417                  !! non-diatoms
418                  zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
419                  !! diatoms
420                  zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
421                  !! diatom silicon
422                  zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
423                  !! AXY (28/01/10): probably need to take account of
424                  !! chl/biomass connection
425                  if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
426                  if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
427                  if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
428                  if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
429             !! AXY (23/01/14): duh - why did I forget diatom silicon?
430             if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
431             if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
432               ENDIF
433            ENDDO
434         ENDDO
435
436         DO jj = 2,jpjm1
437            DO ji = 2,jpim1
438               if (tmask(ji,jj,jk) == 1) then
439                  !! microzooplankton
440                  zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
441                  !! mesozooplankton
442                  zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
443                  !! detrital nitrogen
444                  zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
445                  !! dissolved inorganic nitrogen
446                  zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
447                  !! dissolved silicic acid
448                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
449                  !! dissolved "iron"
450                  zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
451               ENDIF
452            ENDDO
453         ENDDO
454
455# if defined key_roam
456         DO jj = 2,jpjm1
457            DO ji = 2,jpim1
458               if (tmask(ji,jj,jk) == 1) then
459                  !! detrital carbon
460                  zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
461                  !! dissolved inorganic carbon
462                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
463                  !! alkalinity
464                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
465                  !! oxygen
466                  zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
467#  if defined key_axy_carbchem && defined key_mocsy
468                  !! phosphate via DIN and Redfield
469                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
470#  endif
471                  !!
472                  !! also need physical parameters for gas exchange
473                  !! calculations
474                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
475                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
476                  !!
477             !! AXY (28/02/14): check input fields
478                  if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
479                     Call ctl_warn(' trc_bio_medusa: T WARNING 3D : ',       &
480                                   tsn(ji,jj,jk,jp_tem), 'at I, J, K, kt :', & 
481                                   mig(ji), mjg(jj), jk, kt )
482# if defined key_debug_medusa
483                     !! temperature
484                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
485                      IF(lwp) WRITE(numout,*)                                 &
486                     ' trc_bio_medusa: Abnormal T suroundings Temp'
487                      IF(lwp) WRITE(numout,*)                                 &
488                     tsn(ji-1,jj+1,jk,jp_tem), tsn(ji,jj+1,jk,jp_tem), tsn(ji+1,jj+1,jk,jp_tem)
489                      IF(lwp) WRITE(numout,*)                                 &
490                     tsn(ji-1,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), tsn(ji+1,jj,jk,jp_tem)
491                      IF(lwp) WRITE(numout,*)                                 &
492                     tsn(ji-1,jj-1,jk,jp_tem), tsn(ji,jj-1,jk,jp_tem), tsn(ji+1,jj-1,jk,jp_tem)
493                      IF(lwp) WRITE(numout,*)                                 &
494                     ' trc_bio_medusa: Abnormal T suroundings Sal'
495                      IF(lwp) WRITE(numout,*)                                 &
496                     tsn(ji-1,jj+1,jk,jp_sal), tsn(ji,jj+1,jk,jp_sal), tsn(ji+1,jj+1,jk,jp_sal)
497                      IF(lwp) WRITE(numout,*)                                 &
498                     tsn(ji-1,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), tsn(ji+1,jj,jk,jp_sal)
499                      IF(lwp) WRITE(numout,*)                                 &
500                     tsn(ji-1,jj-1,jk,jp_sal), tsn(ji,jj-1,jk,jp_sal), tsn(ji+1,jj-1,jk,jp_sal)
501                     !!
502                      IF(lwp) WRITE(numout,*)                                 &
503                     ' trc_bio_medusa: Abnormal T suroundings DIC'
504                      IF(lwp) WRITE(numout,*)                                 &
505                     trn(ji-1,jj+1,jk,jpdic), trn(ji,jj+1,jk,jpdic), trn(ji+1,jj+1,jk,jpdic)
506                      IF(lwp) WRITE(numout,*)                                 &
507                     trn(ji-1,jj,jk,jpdic), trn(ji,jj,jk,jpdic), trn(ji+1,jj,jk,jpdic)
508                      IF(lwp) WRITE(numout,*)                                 &
509                     trn(ji-1,jj-1,jk,jpdic), trn(ji,jj-1,jk,jpdic), trn(ji+1,jj-1,jk,jpdic)
510                      IF(lwp) WRITE(numout,*)                                 &
511                     ' trc_bio_medusa: Abnormal T suroundings Alk'
512                      IF(lwp) WRITE(numout,*)                                 &
513                     trn(ji-1,jj+1,jk,jpalk), trn(ji,jj+1,jk,jpalk), trn(ji+1,jj+1,jk,jpalk)
514                      IF(lwp) WRITE(numout,*)                                 &
515                     trn(ji-1,jj,jk,jpalk), trn(ji,jj,jk,jpalk), trn(ji+1,jj,jk,jpalk)
516                      IF(lwp) WRITE(numout,*)                                 &
517                     trn(ji-1,jj-1,jk,jpalk), trn(ji,jj-1,jk,jpalk), trn(ji+1,jj-1,jk,jpalk)
518                      IF(lwp) WRITE(numout,*)                                 &
519                     ' trc_bio_medusa: Abnormal T suroundings tmask'
520                      IF(lwp) WRITE(numout,*)                                 &
521                     tmask(ji-1,jj+1,jk), tmask(ji,jj+1,jk), tmask(ji+1,jj+1,jk)
522                      IF(lwp) WRITE(numout,*)                                 &
523                     tmask(ji-1,jj,jk), tmask(ji,jj,jk), tmask(ji+1,jj,jk)
524                      IF(lwp) WRITE(numout,*)                                 &
525                     tmask(ji-1,jj-1,jk), tmask(ji,jj-1,jk), tmask(ji+1,jj-1,jk)
526# endif
527                     !! Correct out of range values
528                     sumtsn = ( tmask(ji-1,jj+1,jk) * tsn(ji-1,jj+1,jk,jp_tem) ) + &
529                              ( tmask(ji  ,jj+1,jk) * tsn(ji  ,jj+1,jk,jp_tem) ) + &
530                              ( tmask(ji+1,jj+1,jk) * tsn(ji+1,jj+1,jk,jp_tem) ) + &
531                              ( tmask(ji-1,jj  ,jk) * tsn(ji-1,jj  ,jk,jp_tem) ) + &
532                              ( tmask(ji+1,jj  ,jk) * tsn(ji+1,jj  ,jk,jp_tem) ) + &
533                              ( tmask(ji-1,jj-1,jk) * tsn(ji-1,jj-1,jk,jp_tem) ) + &
534                              ( tmask(ji  ,jj-1,jk) * tsn(ji  ,jj-1,jk,jp_tem) ) + &
535                              ( tmask(ji+1,jj-1,jk) * tsn(ji+1,jj-1,jk,jp_tem) )
536                     summask = tmask(ji-1,jj+1,jk) + tmask(ji  ,jj+1,jk)   +  &
537                               tmask(ji+1,jj+1,jk) + tmask(ji-1,jj  ,jk)   +  &
538                               tmask(ji+1,jj  ,jk) + tmask(ji-1,jj-1,jk)   +  &
539                               tmask(ji  ,jj-1,jk) + tmask(ji+1,jj-1,jk)
540                     tsnavg = ( sumtsn / summask )
541                     !!
542                     IF ( ( summask .EQ. 0.0 ) .OR. (tsnavg .LT. -3.0 ) .OR.   &
543                          ( tsnavg .GT. 40.0 ) ) THEN   
544                        IF (ztmp(ji,jj) .LT. -3.0 ) THEN
545                           CALL ctl_warn(                                    &
546                           ' trc_bio_medusa: T SWITCHING 3D, ',              &
547                           tsn(ji,jj,jk,jp_tem), ' -> -3.0 ' )
548                           ztmp(ji,jj) = -3.0
549                        ENDIF
550                        IF (ztmp(ji,jj) .GT. 40.0 ) THEN
551                           CALL ctl_warn(                                    &
552                           ' trc_bio_medusa: T SWITCHING 3D, ',              &
553                           tsn(ji,jj,jk,jp_tem), ' -> 40.0 ' )
554                           ztmp(ji,jj) = 40.0
555                        ENDIF
556                     ELSE
557                        CALL ctl_warn(                                       &
558                        ' trc_bio_medusa: T SWITCHING 3D, ',                 &
559                        tsn(ji,jj,jk,jp_tem), ' -> surounding avg : ', tsnavg )
560                        ztmp(ji,jj) = tsnavg
561                     ENDIF
562                  endif
563                  !! end T chack
564                  if (zsal(ji,jj) .lt. 1.0 .or. zsal(ji,jj) .gt. 47.0 ) then
565                        CALL ctl_warn(                                       &
566                        ' trc_bio_medusa: S WARNING 2D, ',                   &
567                        tsn(ji,jj,jk,jp_sal),                                &
568                        ' at I, J, K, kt : ', mig(ji), mjg(jj), jk, kt)
569                     !! Correct out of range values
570                     IF (zsal(ji,jj) .LT. 1.0 ) THEN
571                        CALL ctl_warn(                                       &
572                        ' trc_bio_medusa: S SWITCHING 3D, ',                 &
573                        tsn(ji,jj,jk,jp_sal), ' -> 1.0 ')
574                        zsal(ji,jj) = 1.0
575                     ENDIF
576                     IF (zsal(ji,jj) .GT. 47.0 ) THEN
577                        CALL ctl_warn(                                       &
578                        ' trc_bio_medusa: T SWITCHING 3D, ',                 &
579                        tsn(ji,jj,jk,jp_sal), ' -> 47.0 ')
580                        zsal(ji,jj) = 47.0
581                     ENDIF
582                  endif
583                  !! end S check
584               ENDIF
585            ENDDO
586         ENDDO
587# else
588         DO jj = 2,jpjm1
589            DO ji = 2,jpim1
590               if (tmask(ji,jj,jk) == 1) then
591                  !! implicit detrital carbon
592                  zdtc(ji,jj) = zdet(ji,jj) * xthetad
593               ENDIF
594            ENDDO
595         ENDDO
596# endif
597# if defined key_debug_medusa
598         DO jj = 2,jpjm1
599            DO ji = 2,jpim1
600               if (tmask(ji,jj,jk) == 1) then
601                  if (idf.eq.1) then
602                     !! AXY (15/01/10)
603                     if (trn(ji,jj,jk,jpdin).lt.0.) then
604                        IF (lwp) write (numout,*)                            &
605                           '------------------------------'
606                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',    &
607                           trn(ji,jj,jk,jpdin)
608                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',    &
609                           ji, jj, jk, kt
610                     endif
611                     if (trn(ji,jj,jk,jpsil).lt.0.) then
612                        IF (lwp) write (numout,*)                            &
613                           '------------------------------'
614                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',    &
615                           trn(ji,jj,jk,jpsil)
616                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',    &
617                           ji, jj, jk, kt
618                     endif
619#  if defined key_roam
620                     if (trn(ji,jj,jk,jpdic).lt.0.) then
621                        IF (lwp) write (numout,*)                            &
622                           '------------------------------'
623                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',    &
624                           trn(ji,jj,jk,jpdic)
625                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',    &
626                           ji, jj, jk, kt
627                     endif
628                     if (trn(ji,jj,jk,jpalk).lt.0.) then
629                        IF (lwp) write (numout,*)                            &
630                           '------------------------------'
631                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',    &
632                           trn(ji,jj,jk,jpalk)
633                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',    &
634                           ji, jj, jk, kt
635                     endif
636                     if (trn(ji,jj,jk,jpoxy).lt.0.) then
637                        IF (lwp) write (numout,*)                            &
638                           '------------------------------'
639                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',    &
640                           trn(ji,jj,jk,jpoxy)
641                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',    &
642                           ji, jj, jk, kt
643                     endif
644#  endif
645                  endif
646               ENDIF
647            ENDDO
648         ENDDO
649# endif
650# if defined key_debug_medusa
651! I'M NOT SURE THIS IS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
652!         if (idf.eq.1.AND.idfval.eq.1) then
653!            DO jj = 2,jpjm1
654!               DO ji = 2,jpim1
655!                  if (tmask(ji,jj,jk) == 1) then
656!                     !! report state variable values
657!                     IF (lwp) write (numout,*)                               &
658!                        '------------------------------'
659!                     IF (lwp) write (numout,*) 'fthk(',jk,') = ',            &
660!                        fse3t(ji,jj,jk)
661!                     IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
662!                     IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
663!                     IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
664!                     IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
665!                     IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
666!                     IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
667!                     IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
668!                     IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
669!                     IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
670#  if defined key_roam
671!                     IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
672!                     IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
673!                     IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
674!                     IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)
675#  endif
676!                  ENDIF
677!               ENDDO
678!            ENDDO
679!         endif
680# endif
681
682# if defined key_debug_medusa
683! I'M NOT SURE THIS IS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17
684!         if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then
685!            DO jj = 2,jpjm1
686!               DO ji = 2,jpim1
687!                  if (tmask(ji,jj,jk) == 1) then
688!                     IF (lwp) write (numout,*)                               &
689!                       '------------------------------'
690!                     IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj)
691!                  ENDIF
692!               ENDDO
693!            ENDDO
694!         endif
695# endif
696
697         !!---------------------------------------------------------------
698         !! Calculate air-sea gas exchange and river inputs
699         !!---------------------------------------------------------------
700         IF ( jk == 1 ) THEN
701            CALL air_sea( kt )
702         ENDIF
703
704         !!---------------------------------------------------------------
705         !! Phytoplankton growth, zooplankton grazing and miscellaneous
706         !! plankton losses.
707         !!---------------------------------------------------------------
708         CALL plankton( jk )
709
710         !!---------------------------------------------------------------
711         !! Iron chemistry and scavenging
712         !!---------------------------------------------------------------
713         CALL iron_chem_scav( jk )
714
715         !!---------------------------------------------------------------
716         !! Detritus processes
717         !!---------------------------------------------------------------
718         CALL detritus( jk, iball )
719
720         !!---------------------------------------------------------------
721         !! Updating tracers
722         !!---------------------------------------------------------------
723         CALL bio_medusa_update( kt, jk )
724
725         !!---------------------------------------------------------------
726         !! Diagnostics
727         !!---------------------------------------------------------------
728         CALL bio_medusa_diag( jk )
729
730         !!-------------------------------------------------------
731         !! 2d specific k level diags
732         !!-------------------------------------------------------
733         IF( lk_iomput ) THEN
734            CALL bio_medusa_diag_slice( jk )
735         ENDIF
736
737      !! CLOSE vertical loop
738      ENDDO
739
740      !!------------------------------------------------------------------
741      !! Final calculations for diagnostics
742      !!------------------------------------------------------------------
743      CALL bio_medusa_fin( kt )
744
745# if defined key_debug_medusa
746       IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
747       CALL flush(numout)
748# endif
749
750   END SUBROUTINE trc_bio_medusa
751
752#else
753   !!=====================================================================
754   !!  Dummy module :                                   No MEDUSA bio-model
755   !!=====================================================================
756CONTAINS
757   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
758      INTEGER, INTENT( in ) ::   kt
759      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
760   END SUBROUTINE trc_bio_medusa
761#endif 
762
763   !!=====================================================================
764END MODULE  trcbio_medusa
Note: See TracBrowser for help on using the repository browser.