source: branches/NERC/dev_r5518_GO6_CarbFix_DMSMin_3DMOCSY/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 9335

Last change on this file since 9335 was 9335, checked in by jpalmier, 2 years ago

Import 3D carb-chem calling conditions for tests in the coupled model

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