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

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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8442

Last change on this file since 8442 was 8442, checked in by frrh, 7 years ago

Commit changes relating to Met Office GMED ticket 340 for the
tidying of MEDUSA related code and debugging statements in the TOP code.

Only code introduced at revision 8434 of branch
http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
is included here, all previous revisions of that branch having been dealt with
under GMED ticket 339.

File size: 30.3 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_rsttr, nittrc000, trn
102      USE bio_medusa_init_mod,        ONLY: bio_medusa_init
103      USE carb_chem_mod,              ONLY: carb_chem
104      USE air_sea_mod,                ONLY: air_sea
105      USE plankton_mod,               ONLY: plankton
106      USE iron_chem_scav_mod,         ONLY: iron_chem_scav
107      USE detritus_mod,               ONLY: detritus
108      USE bio_medusa_update_mod,      ONLY: bio_medusa_update
109      USE bio_medusa_diag_mod,        ONLY: bio_medusa_diag
110      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
111      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
112
113      IMPLICIT NONE
114      PRIVATE
115     
116      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
117
118   !!* Substitution
119#  include "domzgr_substitute.h90"
120   !!----------------------------------------------------------------------
121   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
122   !! $Id$
123   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
124   !!----------------------------------------------------------------------
125
126CONTAINS
127
128   SUBROUTINE trc_bio_medusa( kt )
129      !!------------------------------------------------------------------
130      !!                     ***  ROUTINE trc_bio  ***
131      !!
132      !! ** Purpose : compute the now trend due to biogeochemical processes
133      !!              and add it to the general trend of passive tracers
134      !!              equations
135      !!
136      !! ** Method  : each now biological flux is calculated in function of
137      !!              now concentrations of tracers.
138      !!              depending on the tracer, these fluxes are sources or
139      !!              sinks.
140      !!              The total of the sources and sinks for each tracer
141      !!              is added to the general trend.
142      !!       
143      !!                      tra = tra + zf...tra - zftra...
144      !!                                     |         |
145      !!                                     |         |
146      !!                                  source      sink
147      !!       
148      !!              IF 'key_trc_diabio' defined , the biogeochemical trends
149      !!              for passive tracers are saved for futher diagnostics.
150      !!------------------------------------------------------------------
151      !!
152      !!
153      !!------------------------------------------------------------------
154      !! Variable conventions
155      !!------------------------------------------------------------------
156      !!
157      !! names: z*** - state variable
158      !!        f*** - function (or temporary variable used in part of
159      !!               a function)
160      !!        x*** - parameter
161      !!        b*** - right-hand part (sources and sinks)
162      !!        i*** - integer variable (usually used in yes/no flags)
163      !!
164      !! time (integer timestep)
165      INTEGER, INTENT( in ) ::    kt
166      !!
167      !! spatial array indices
168      INTEGER  ::    ji,jj,jk,jn
169      !!
170      INTEGER  ::    iball
171# if defined key_roam
172      !!
173      INTEGER  ::    iyr1, iyr2
174      !!
175# endif
176      !!
177      !! temporary variables
178      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
179      !!
180      !!------------------------------------------------------------------
181
182# if defined key_debug_medusa
183      IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined'
184      CALL flush(numout)
185# endif 
186
187      !! AXY (20/11/14): alter this to report on first MEDUSA call
188      !! IF( kt == nit000 ) THEN
189      IF( kt == nittrc000 ) THEN
190         IF(lwp) WRITE(numout,*)
191         IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA bio-model'
192         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
193    IF(lwp) WRITE(numout,*) ' kt =',kt
194      ENDIF
195
196      !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
197      ibenthic = 1
198
199      !!------------------------------------------------------------------
200      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency
201      !! terms of all biological equations to 0.
202      !!------------------------------------------------------------------
203      !!
204      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit
205      !!                 the bill for now; another item on the things-to-sort-
206      !!     out-in-the-future list ...
207# if defined key_kill_medusa
208      b0 = 0.
209# else
210      b0 = 1.
211# endif
212      !!------------------------------------------------------------------
213      !! fast detritus ballast scheme (0 = no; 1 = yes)
214      !! alternative to ballast scheme is same scheme but with no ballast
215      !! protection (not dissimilar to Martin et al., 1987)
216      !!------------------------------------------------------------------
217      !!
218      iball = 1
219
220      !!------------------------------------------------------------------
221      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output
222      !! these should *only* be used in 1D since they give comprehensive
223      !! output for ecological functions in the model; primarily used in
224      !! debugging
225      !!------------------------------------------------------------------
226      !!
227      idf    = 0
228      !!
229      !! timer mechanism
230      if (kt/120*120.eq.kt) then
231         idfval = 1
232      else
233         idfval = 0
234      endif
235
236      !!------------------------------------------------------------------
237      !! Initialise arrays to zero and set up arrays for diagnostics
238      !!------------------------------------------------------------------
239      CALL bio_medusa_init( kt )
240
241# if defined key_axy_nancheck
242       DO jn = jp_msa0,jp_msa1
243          !! fq0 = MINVAL(trn(:,:,:,jn))
244          !! fq1 = MAXVAL(trn(:,:,:,jn))
245          fq2 = SUM(trn(:,:,:,jn))
246          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK',     &
247          !!                kt, jn, fq0, fq1, fq2
248          !! AXY (30/01/14): much to our surprise, the next line doesn't
249          !!                 work on HECTOR and has been replaced here with
250          !!                 a specialist routine
251          !! if (fq2 /= fq2 ) then
252          if ( ieee_is_nan( fq2 ) ) then
253             !! there's a NaN here
254             if (lwp) write(numout,*) 'NAN detected in field', jn,           &
255                                      'at time', kt, 'at position:'
256             DO jk = 1,jpk
257                DO jj = 1,jpj
258                   DO ji = 1,jpi
259                      !! AXY (30/01/14): "isnan" problem on HECTOR
260                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then
261                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then
262                         if (lwp) write (numout,'(a,1pe12.2,4i6)')           &
263                            'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn
264                      endif
265                   enddo
266                enddo
267             enddo
268             CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
269          endif
270       ENDDO
271       CALL flush(numout)
272# endif
273
274# if defined key_debug_medusa
275      IF (lwp) write (numout,*)                                              &
276                     'trc_bio_medusa: variables initialised and checked'
277      CALL flush(numout)
278# endif 
279
280# if defined key_roam
281      !!------------------------------------------------------------------
282      !! calculate atmospheric pCO2
283      !!------------------------------------------------------------------
284      !!
285      !! what's atmospheric pCO2 doing? (data start in 1859)
286      iyr1  = nyear - 1859 + 1
287      iyr2  = iyr1 + 1
288      if (iyr1 .le. 1) then
289         !! before 1860
290         f_xco2a(:,:) = hist_pco2(1)
291      elseif (iyr2 .ge. 242) then
292         !! after 2099
293         f_xco2a(:,:) = hist_pco2(242)
294      else
295         !! just right
296         fq0 = hist_pco2(iyr1)
297         fq1 = hist_pco2(iyr2)
298         fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0)
299         !! AXY (14/06/12): tweaked to make more sense (and be correct)
300#  if defined key_bs_axy_yrlen
301         !! bugfix: for 360d year with HadGEM2-ES forcing
302         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 
303#  else
304         !! original use of 365 days (not accounting for leap year or
305         !! 360d year)
306         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0
307#  endif
308         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
309         f_xco2a(:,:) = fq4
310      endif
311#  if defined key_axy_pi_co2
312      !! OCMIP pre-industrial pCO2
313      !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2
314      f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2
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           ((86400*mod(nn_date0,100) + 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,jk) == 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,jk) == 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,jk) == 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,jk) == 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 IS 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,jk) == 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 IS 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,jk) == 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( jk )
636
637         !!-------------------------------------------------------
638         !! 2d specific k level diags
639         !!-------------------------------------------------------
640         IF( lk_iomput ) 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_debug_medusa
653       IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
654       CALL flush(numout)
655# endif
656
657   END SUBROUTINE trc_bio_medusa
658
659#else
660   !!=====================================================================
661   !!  Dummy module :                                   No MEDUSA bio-model
662   !!=====================================================================
663CONTAINS
664   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
665      INTEGER, INTENT( in ) ::   kt
666      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
667   END SUBROUTINE trc_bio_medusa
668#endif 
669
670   !!=====================================================================
671END MODULE  trcbio_medusa
Note: See TracBrowser for help on using the repository browser.