source: branches/NERC/dev_r5518_GO6_split_trcbiomedusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8395

Last change on this file since 8395 was 8395, checked in by jpalmier, 3 years ago

JPALM — GMED #339 - split trcbio_medusa only

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