Context Navigation

trcbio_medusa.F90 @ 8649

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

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

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

Context Navigation

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

Download in other formats: