New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcbio_medusa.F90 in branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8023

Last change on this file since 8023 was 8023, checked in by marc, 7 years ago

Tidying up of headers for several files

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