source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8495

Last change on this file since 8495 was 8495, checked in by dford, 3 years ago

Merge in changes from dev_r5518_GO6_package_asm_surf_bgc, and adapt to the updated MEDUSA structure.

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