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/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

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

Last change on this file since 8642 was 8642, checked in by jpalmier, 7 years ago

jpalm --19-10-17-- Carb failure due to 1-cell exceptionnal and ephemeral T increase - update T passed to MEDUSA and add kill switch

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