source: branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8649

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

JPALM — 20-10-17 — ctl_warn text fix

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