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

Rev	Line
[5726]	1	MODULE trcbio_medusa
	2	!!======================================================================
	3	!! * MODULE trcbio *
	4	!! TOP : MEDUSA
	5	!!======================================================================
	6	!! History :
	7	!! - ! 1999-07 (M. Levy) original code
[8441]	8	!! - ! 2000-12 (E. Kestenare) assign parameters to name
	9	!! individual tracers
[5726]	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
[5841]	19	!! - ! 2015-06 (A. Yool) Update to include MOCSY
	20	!! - ! 2015-07 (A. Yool) Update for rolling averages
[8441]	21	!! - ! 2015-10 (J. Palm) Update for diag outputs through
	22	!! iom_use
[7224]	23	!! - ! 2016-11 (A. Yool) Updated diags for CMIP6
[8131]	24	!! - ! 2017-05 (A. Yool) Added extra DMS calculation
[5726]	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
[8074]	31	!! - addition of air-sea fluxes of CO2 and oxygen
[5726]	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
[8074]	42	!!----------------------------------------------------------------------
	43	#endif
[5726]	44	!!
[8074]	45	#if defined key_mocsy
	46	!!----------------------------------------------------------------------
[5841]	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	!!
[5726]	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
[8441]	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
[6744]	98	!! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
[8441]	99	USE sbc_oce, ONLY: lk_oasis
[8495]	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
[8441]	105	USE sms_medusa, ONLY: hist_pco2
[8495]	106	# endif
[8442]	107	USE trc, ONLY: ln_rsttr, nittrc000, trn
[8441]	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
[6744]	118
[5726]	119	IMPLICIT NONE
	120	PRIVATE
	121
[8441]	122	PUBLIC trc_bio_medusa ! called in trcsms_medusa.F90
[5726]	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
[8074]	134	SUBROUTINE trc_bio_medusa( kt )
[8441]	135	!!------------------------------------------------------------------
[8074]	136	!! * ROUTINE trc_bio *
	137	!!
[8441]	138	!! ** Purpose : compute the now trend due to biogeochemical processes
	139	!! and add it to the general trend of passive tracers
	140	!! equations
[8074]	141	!!
[8441]	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
[8074]	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.
[8441]	156	!!------------------------------------------------------------------
[8074]	157	!!
	158	!!
[8441]	159	!!------------------------------------------------------------------
[8074]	160	!! Variable conventions
[8441]	161	!!------------------------------------------------------------------
[8074]	162	!!
	163	!! names: z*** - state variable
[8441]	164	!! f*** - function (or temporary variable used in part of
	165	!! a function)
[8074]	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
[5726]	177	# if defined key_roam
[8074]	178	!!
	179	INTEGER :: iyr1, iyr2
	180	!!
[5726]	181	# endif
[8074]	182	!!
[8441]	183	!! temporary variables
	184	REAL(wp) :: fq0,fq1,fq2,fq3,fq4
[8074]	185	!!
[8441]	186	!!------------------------------------------------------------------
[5726]	187
[5841]	188	# if defined key_debug_medusa
[8074]	189	IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined'
	190	CALL flush(numout)
[5841]	191	# endif
	192
[8074]	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
[5726]	201
[8074]	202	!! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes
	203	ibenthic = 1
[5726]	204
[8441]	205	!!------------------------------------------------------------------
[8074]	206	!! b0 is present for debugging purposes; using b0 = 0 sets the tendency
	207	!! terms of all biological equations to 0.
[8441]	208	!!------------------------------------------------------------------
[8074]	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 ...
[5726]	213	# if defined key_kill_medusa
[8074]	214	b0 = 0.
[5726]	215	# else
[8074]	216	b0 = 1.
[5726]	217	# endif
[8441]	218	!!------------------------------------------------------------------
[8074]	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)
[8441]	222	!!------------------------------------------------------------------
[8074]	223	!!
	224	iball = 1
[5726]	225
[8441]	226	!!------------------------------------------------------------------
[8074]	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
[8441]	231	!!------------------------------------------------------------------
[8074]	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
[5726]	241
[8441]	242	!!------------------------------------------------------------------
	243	!! Initialise arrays to zero and set up arrays for diagnostics
	244	!!------------------------------------------------------------------
	245	CALL bio_medusa_init( kt )
[5726]	246
	247	# if defined key_axy_nancheck
[8441]	248	DO jn = jp_msa0,jp_msa1
[8074]	249	!! fq0 = MINVAL(trn(:,:,:,jn))
	250	!! fq1 = MAXVAL(trn(:,:,:,jn))
	251	fq2 = SUM(trn(:,:,:,jn))
[8441]	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
[8074]	257	!! if (fq2 /= fq2 ) then
	258	if ( ieee_is_nan( fq2 ) ) then
	259	!! there's a NaN here
[8441]	260	if (lwp) write(numout,*) 'NAN detected in field', jn, &
	261	'at time', kt, 'at position:'
[8074]	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
[8441]	268	if (lwp) write (numout,'(a,1pe12.2,4i6)') &
	269	'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn
[8074]	270	endif
	271	enddo
[7224]	272	enddo
	273	enddo
[8074]	274	CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' )
	275	endif
	276	ENDDO
	277	CALL flush(numout)
[5726]	278	# endif
	279
[5841]	280	# if defined key_debug_medusa
[8441]	281	IF (lwp) write (numout,*) &
	282	'trc_bio_medusa: variables initialised and checked'
[8074]	283	CALL flush(numout)
[5841]	284	# endif
	285
[5726]	286	# if defined key_roam
[8441]	287	!!------------------------------------------------------------------
[8074]	288	!! calculate atmospheric pCO2
[8441]	289	!!------------------------------------------------------------------
[8074]	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
[8441]	296	f_xco2a(:,:) = hist_pco2(1)
[8074]	297	elseif (iyr2 .ge. 242) then
	298	!! after 2099
[8441]	299	f_xco2a(:,:) = hist_pco2(242)
[8074]	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
[8441]	307	!! bugfix: for 360d year with HadGEM2-ES forcing
	308	fq3 = (real(nday_year) - 1.0 + fq2) / 360.0
[5726]	309	# else
[8441]	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
[5726]	313	# endif
[8074]	314	fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
[8441]	315	f_xco2a(:,:) = fq4
[8074]	316	endif
[8495]	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
[8074]	325	# if defined key_axy_pi_co2
[8441]	326	!! OCMIP pre-industrial pCO2
	327	!! f_xco2a(:,:) = 284.725 !! CMIP5 pre-industrial pCO2
	328	f_xco2a = 284.317 !! CMIP6 pre-industrial pCO2
[8074]	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
[8441]	338	IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a(1,1)
[5726]	339	# endif
	340
[5841]	341	# if defined key_debug_medusa
[8074]	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)
[5841]	346	# endif
	347
[5726]	348	# if defined key_roam
[8074]	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 :
[8441]	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)
[8074]	363	!! ++ need to pass carb-chem output var through restarts
[8495]	364	#if defined key_foam_medusa
	365	!! DAF (Aug 2017): For FOAM we want to run daily
[8224]	366	If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR. &
[8495]	367	(mod(kt*rdt,86400.) == rdt) ) THEN
	368	#else
	369	If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR. &
[8224]	370	((86400mod(nn_date0,100) + mod(ktrdt,2592000.)) == rdt) ) THEN
[8495]	371	#endif
	372	!!----------------------------------------------------------------------
[8074]	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
[8441]	376	!!---------------------------------------------------------------
	377	CALL carb_chem( kt )
	378
[8074]	379	ENDIF
[5726]	380	# endif
	381
[5841]	382	# if defined key_debug_medusa
[8074]	383	IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations'
	384	CALL flush(numout)
[5841]	385	# endif
	386
[8441]	387	!!------------------------------------------------------------------
[8074]	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
[8441]	391	!!------------------------------------------------------------------
[8074]	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
[8441]	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	!!======================================================
[8074]	413	!!
[8441]	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	!!------------------------------------------------------
[8074]	422	!!
[8441]	423	!! AXY (01/03/10): set up level depth (bottom of level)
	424	fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
[8074]	425	!!
[8441]	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
[7224]	451
[8441]	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))
[8074]	467	ENDIF
[8441]	468	ENDDO
	469	ENDDO
[7224]	470
[8441]	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))
[8074]	479	!! alkalinity
[8441]	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
[8074]	486	# endif
	487	!!
[8441]	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)
[8074]	492	!!
[8441]	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)
[8074]	504	endif
[8441]	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
[8074]	510	endif
[8441]	511	ENDIF
	512	ENDDO
	513	ENDDO
[5726]	514	# else
[8441]	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
[5726]	523	# endif
	524	# if defined key_debug_medusa
[8441]	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
[8074]	537	endif
[8441]	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
[8074]	545	endif
[8441]	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
[8074]	554	endif
[8441]	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
[8074]	562	endif
[8441]	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
[8074]	570	endif
[5726]	571	# endif
[8074]	572	endif
[8441]	573	ENDIF
	574	ENDDO
	575	ENDDO
[5726]	576	# endif
[8074]	577	# if defined key_debug_medusa
[8441]	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)
[8074]	597	# if defined key_roam
[8441]	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)
[8074]	602	# endif
[8441]	603	! ENDIF
	604	! ENDDO
	605	! ENDDO
	606	! endif
[8074]	607	# endif
	608
[8441]	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
[5726]	622	# endif
	623
[8441]	624	!!---------------------------------------------------------------
	625	!! Calculate air-sea gas exchange and river inputs
	626	!!---------------------------------------------------------------
	627	IF ( jk == 1 ) THEN
	628	CALL air_sea( kt )
	629	ENDIF
[5726]	630
[8441]	631	!!---------------------------------------------------------------
	632	!! Phytoplankton growth, zooplankton grazing and miscellaneous
	633	!! plankton losses.
	634	!!---------------------------------------------------------------
	635	CALL plankton( jk )
[7224]	636
[8441]	637	!!---------------------------------------------------------------
	638	!! Iron chemistry and scavenging
	639	!!---------------------------------------------------------------
	640	CALL iron_chem_scav( jk )
[6022]	641
[8441]	642	!!---------------------------------------------------------------
	643	!! Detritus processes
	644	!!---------------------------------------------------------------
	645	CALL detritus( jk, iball )
[5726]	646
[8441]	647	!!---------------------------------------------------------------
	648	!! Updating tracers
	649	!!---------------------------------------------------------------
	650	CALL bio_medusa_update( kt, jk )
[5726]	651
[8441]	652	!!---------------------------------------------------------------
	653	!! Diagnostics
	654	!!---------------------------------------------------------------
	655	CALL bio_medusa_diag( jk )
[8074]	656
[8441]	657	!!-------------------------------------------------------
	658	!! 2d specific k level diags
	659	!!-------------------------------------------------------
[8442]	660	IF( lk_iomput ) THEN
[8441]	661	CALL bio_medusa_diag_slice( jk )
	662	ENDIF
[5726]	663
[8074]	664	!! CLOSE vertical loop
	665	ENDDO
	666
[8441]	667	!!------------------------------------------------------------------
	668	!! Final calculations for diagnostics
	669	!!------------------------------------------------------------------
	670	CALL bio_medusa_fin( kt )
[8074]	671
[5726]	672	# if defined key_debug_medusa
[8074]	673	IF(lwp) WRITE(numout,*) ' MEDUSA exiting trc_bio_medusa at kt =', kt
	674	CALL flush(numout)
[5726]	675	# endif
	676
[8074]	677	END SUBROUTINE trc_bio_medusa
[5726]	678
	679	#else
[8441]	680	!!=====================================================================
[5726]	681	!! Dummy module : No MEDUSA bio-model
[8441]	682	!!=====================================================================
[5726]	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
[8441]	690	!!=====================================================================
[5726]	691	END MODULE trcbio_medusa

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

Download in other formats: