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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

trcbio_medusa.F90 in branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 @ 8495

Last change on this file since 8495 was 8495, checked in by dford, 7 years ago
Merge in changes from dev_r5518_GO6_package_asm_surf_bgc, and adapt to the updated MEDUSA structure.
File size: 31.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: