Context Navigation

initracer.F @ 3757

Last change on this file since 3757 was 3726, checked in by emillour, 3 months ago
Mars PCM: Turn "callkeys.h" into module "callkeys_mod.F90" EM
File size: 38.3 KB

Rev	Line
[1224]	1	SUBROUTINE initracer(ngrid,nq,qsurf)
[38]	2
[1036]	3	use tracer_mod
[2321]	4	use comcstfi_h, only: pi
[2409]	5	use dust_param_mod, only: doubleq, submicron, dustbin
[3726]	6	use callkeys_mod, only: water, microphys, scavenging, rdstorm,
	7	& topflows, photochem, callthermos, hdo,
	8	& callnlte, nltemodel, dustinjection,
	9	& co2clouds, co2useh2o, meteo_flux
[38]	10	IMPLICIT NONE
	11	c=======================================================================
	12	c subject:
	13	c --------
	14	c Initialization related to tracer
	15	c (transported dust, water, chemical species, ice...)
	16	c
	17	c Name of the tracer
	18	c
	19	c Test of dimension :
[1036]	20	c Initialize tracer related data in tracer_mod, using tracer names provided
	21	c by the dynamics in "infotrac"
[38]	22	c
	23	c
	24	c author: F.Forget
	25	c ------
	26	c Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003)
	27	c Ehouarn Millour (oct. 2008) identify tracers by their names
	28	c=======================================================================
	29
	30
[1036]	31	integer,intent(in) :: ngrid ! number of atmospheric columns
	32	integer,intent(in) :: nq ! number of tracers
	33	real,intent(out) :: qsurf(ngrid,nq) ! tracer on surface (e.g. kg.m-2)
	34
[283]	35	integer iq,ig,count
[38]	36	real r0_lift , reff_lift, nueff_lift
[1974]	37	real r0_storm,reff_storm
[38]	38	c Ratio of small over large dust particles (used when both
	39	c doubleq and the submicron mode are active); In Montmessin
	40	c et al. (2002), a value of 25 has been deduced;
	41	real, parameter :: popratio = 25.
[1974]	42	character(len=30) :: txt ! to store some text
[3027]	43	character(len=30), dimension(nq) :: idtracers
[38]	44
	45	c-----------------------------------------------------------------------
[1036]	46	c radius(nq) ! aerosol particle radius (m)
	47	c rho_q(nq) ! tracer densities (kg.m-3)
	48	c alpha_lift(nq) ! saltation vertical flux/horiz flux ratio (m-1)
	49	c alpha_devil(nq) ! lifting coeeficient by dust devil
[38]	50	c rho_dust ! Mars dust density
	51	c rho_ice ! Water ice density
	52	c doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio
	53	c varian ! Characteristic variance of log-normal distribution
	54	c-----------------------------------------------------------------------
	55
[2823]	56	! Sanity check: check that we are running with at least 1 tracer:
	57	if (nq<1) then
	58	write(,) "initracer error, nq=",nq," must be >=1!"
	59	call abort_physic("initracer","nq<1",1)
	60	endif
[38]	61	c------------------------------------------------------------
	62	c NAME and molar mass of the tracer
	63	c------------------------------------------------------------
	64
	65	! Identify tracers by their names: (and set corresponding values of mmol)
	66	! 0. initialize tracer indexes to zero:
[1036]	67	igcm_dustbin(1:nq)=0
[1617]	68	igcm_co2_ice=0
	69	igcm_ccnco2_mass=0
	70	igcm_ccnco2_number=0
[2589]	71	igcm_ccnco2_meteor_mass=0
	72	igcm_ccnco2_meteor_number=0
[2562]	73	igcm_ccnco2_h2o_mass_ice=0
	74	igcm_ccnco2_h2o_mass_ccn=0
	75	igcm_ccnco2_h2o_number=0
[38]	76	igcm_dust_mass=0
	77	igcm_dust_number=0
[358]	78	igcm_ccn_mass=0
	79	igcm_ccn_number=0
[38]	80	igcm_dust_submicron=0
	81	igcm_h2o_vap=0
	82	igcm_h2o_ice=0
[2312]	83	igcm_hdo_vap=0
	84	igcm_hdo_ice=0
[1974]	85	igcm_stormdust_mass=0
	86	igcm_stormdust_number=0
[2199]	87	igcm_topdust_mass=0
	88	igcm_topdust_number=0
[38]	89	igcm_co2=0
	90	igcm_co=0
	91	igcm_o=0
	92	igcm_o1d=0
	93	igcm_o2=0
	94	igcm_o3=0
	95	igcm_h=0
[2461]	96	igcm_d=0
	97	igcm_hd=0
[38]	98	igcm_h2=0
[2461]	99	igcm_od=0
	100	igcm_do2=0
	101	igcm_hdo2=0
[38]	102	igcm_oh=0
	103	igcm_ho2=0
	104	igcm_h2o2=0
[459]	105	igcm_ch4=0
[38]	106	igcm_n2=0
	107	igcm_ar=0
	108	igcm_ar_n2=0
[171]	109	igcm_n=0
	110	igcm_no=0
	111	igcm_no2=0
	112	igcm_n2d=0
[1660]	113	igcm_he=0
[171]	114	igcm_co2plus=0
	115	igcm_oplus=0
	116	igcm_o2plus=0
	117	igcm_coplus=0
	118	igcm_cplus=0
	119	igcm_nplus=0
	120	igcm_noplus=0
	121	igcm_n2plus=0
	122	igcm_hplus=0
[635]	123	igcm_hco2plus=0
[2284]	124	igcm_hcoplus=0
[2302]	125	igcm_h2oplus=0
[2321]	126	igcm_h3oplus=0
	127	igcm_ohplus=0
[171]	128	igcm_elec=0
[2616]	129
[38]	130	! 1. find dust tracers
	131	count=0
	132	if (dustbin.gt.0) then
[1036]	133	do iq=1,nq
[38]	134	txt=" "
	135	write(txt,'(a4,i2.2)')'dust',count+1
	136	if (noms(iq).eq.txt) then
	137	count=count+1
[3027]	138	idtracers(count) = txt
[38]	139	igcm_dustbin(count)=iq
	140	mmol(iq)=100.
	141	endif
[1036]	142	enddo !do iq=1,nq
[38]	143	endif ! of if (dustbin.gt.0)
	144	if (doubleq) then
[1036]	145	do iq=1,nq
[38]	146	if (noms(iq).eq."dust_mass") then
	147	igcm_dust_mass=iq
	148	count=count+1
[3027]	149	idtracers(count) = "dust_mass"
[38]	150	endif
	151	if (noms(iq).eq."dust_number") then
	152	igcm_dust_number=iq
	153	count=count+1
[3027]	154	idtracers(count) = "dust_number"
[38]	155	endif
	156	enddo
	157	endif ! of if (doubleq)
[740]	158	if (microphys) then
[1036]	159	do iq=1,nq
[358]	160	if (noms(iq).eq."ccn_mass") then
	161	igcm_ccn_mass=iq
	162	count=count+1
[3027]	163	idtracers(count) = "ccn_mass"
[358]	164	endif
	165	if (noms(iq).eq."ccn_number") then
	166	igcm_ccn_number=iq
	167	count=count+1
[3027]	168	idtracers(count) = "ccn_number"
[358]	169	endif
	170	enddo
[740]	171	endif ! of if (microphys)
[38]	172	if (submicron) then
[1036]	173	do iq=1,nq
[38]	174	if (noms(iq).eq."dust_submicron") then
	175	igcm_dust_submicron=iq
	176	mmol(iq)=100.
	177	count=count+1
[3027]	178	idtracers(count) = "dust_submicron"
[38]	179	endif
	180	enddo
	181	endif ! of if (submicron)
[1974]	182	if (rdstorm) then
	183	do iq=1,nq
	184	if (noms(iq).eq."stormdust_mass") then
	185	igcm_stormdust_mass=iq
	186	count=count+1
[3027]	187	idtracers(count) = "stormdust_mass"
[1974]	188	endif
	189	if (noms(iq).eq."stormdust_number") then
	190	igcm_stormdust_number=iq
	191	count=count+1
[3027]	192	idtracers(count) = "stormdust_number"
[1974]	193	endif
	194	enddo
[2199]	195	endif ! of if (rdstorm)
[2628]	196	if (topflows) then
[2199]	197	do iq=1,nq
	198	if (noms(iq).eq."topdust_mass") then
	199	igcm_topdust_mass=iq
	200	count=count+1
[3027]	201	idtracers(count) = "topdust_mass"
[2199]	202	endif
	203	if (noms(iq).eq."topdust_number") then
	204	igcm_topdust_number=iq
	205	count=count+1
[3027]	206	idtracers(count) = "topdust_number"
[2199]	207	endif
	208	enddo
[2628]	209	endif ! of if (topflows)
[38]	210	! 2. find chemistry and water tracers
[1036]	211	do iq=1,nq
[38]	212	if (noms(iq).eq."co2") then
	213	igcm_co2=iq
	214	mmol(igcm_co2)=44.
	215	count=count+1
[3027]	216	idtracers(count) = "co2"
[38]	217	endif
	218	if (noms(iq).eq."co") then
	219	igcm_co=iq
	220	mmol(igcm_co)=28.
	221	count=count+1
[3027]	222	idtracers(count) = "co"
[38]	223	endif
	224	if (noms(iq).eq."o") then
	225	igcm_o=iq
	226	mmol(igcm_o)=16.
	227	count=count+1
[3027]	228	idtracers(count) = "o"
[38]	229	endif
	230	if (noms(iq).eq."o1d") then
	231	igcm_o1d=iq
	232	mmol(igcm_o1d)=16.
	233	count=count+1
[3027]	234	idtracers(count) = "o1d"
[38]	235	endif
	236	if (noms(iq).eq."o2") then
	237	igcm_o2=iq
	238	mmol(igcm_o2)=32.
	239	count=count+1
[3027]	240	idtracers(count) = "o2"
[38]	241	endif
	242	if (noms(iq).eq."o3") then
	243	igcm_o3=iq
	244	mmol(igcm_o3)=48.
	245	count=count+1
[3027]	246	idtracers(count) = "o3"
[38]	247	endif
	248	if (noms(iq).eq."h") then
	249	igcm_h=iq
	250	mmol(igcm_h)=1.
	251	count=count+1
[3027]	252	idtracers(count) = "h"
[38]	253	endif
	254	if (noms(iq).eq."h2") then
	255	igcm_h2=iq
	256	mmol(igcm_h2)=2.
	257	count=count+1
[3027]	258	idtracers(count) = "h2"
[38]	259	endif
	260	if (noms(iq).eq."oh") then
	261	igcm_oh=iq
	262	mmol(igcm_oh)=17.
	263	count=count+1
[3027]	264	idtracers(count) = "oh"
[38]	265	endif
	266	if (noms(iq).eq."ho2") then
	267	igcm_ho2=iq
	268	mmol(igcm_ho2)=33.
	269	count=count+1
[3027]	270	idtracers(count) = "ho2"
[38]	271	endif
	272	if (noms(iq).eq."h2o2") then
	273	igcm_h2o2=iq
	274	mmol(igcm_h2o2)=34.
	275	count=count+1
[3027]	276	idtracers(count) = "h2o2"
[38]	277	endif
	278	if (noms(iq).eq."n2") then
	279	igcm_n2=iq
	280	mmol(igcm_n2)=28.
	281	count=count+1
[3027]	282	idtracers(count) = "n2"
[38]	283	endif
[324]	284	if (noms(iq).eq."ch4") then
	285	igcm_ch4=iq
	286	mmol(igcm_ch4)=16.
	287	count=count+1
[3027]	288	idtracers(count) = "ch4"
[324]	289	endif
[38]	290	if (noms(iq).eq."ar") then
	291	igcm_ar=iq
	292	mmol(igcm_ar)=40.
	293	count=count+1
[3027]	294	idtracers(count) = "ar"
[38]	295	endif
[324]	296	if (noms(iq).eq."n") then
[171]	297	igcm_n=iq
	298	mmol(igcm_n)=14.
	299	count=count+1
[3027]	300	idtracers(count) = "n"
[171]	301	endif
	302	if (noms(iq).eq."no") then
	303	igcm_no=iq
	304	mmol(igcm_no)=30.
	305	count=count+1
[3027]	306	idtracers(count) = "no"
[171]	307	endif
	308	if (noms(iq).eq."no2") then
	309	igcm_no2=iq
	310	mmol(igcm_no2)=46.
	311	count=count+1
[3027]	312	idtracers(count) = "no2"
[171]	313	endif
	314	if (noms(iq).eq."n2d") then
	315	igcm_n2d=iq
	316	mmol(igcm_n2d)=28.
	317	count=count+1
[3027]	318	idtracers(count) = "n2d"
[171]	319	endif
[1660]	320	if (noms(iq).eq."he") then
	321	igcm_he=iq
	322	mmol(igcm_he)=4.
	323	count=count+1
[3027]	324	idtracers(count) = "he"
[1660]	325	endif
[171]	326	if (noms(iq).eq."co2plus") then
	327	igcm_co2plus=iq
	328	mmol(igcm_co2plus)=44.
	329	count=count+1
[3027]	330	idtracers(count) = "co2plus"
[171]	331	endif
	332	if (noms(iq).eq."oplus") then
	333	igcm_oplus=iq
	334	mmol(igcm_oplus)=16.
	335	count=count+1
[3027]	336	idtracers(count) = "oplus"
[171]	337	endif
	338	if (noms(iq).eq."o2plus") then
	339	igcm_o2plus=iq
	340	mmol(igcm_o2plus)=32.
	341	count=count+1
[3027]	342	idtracers(count) = "o2plus"
[171]	343	endif
	344	if (noms(iq).eq."coplus") then
	345	igcm_coplus=iq
	346	mmol(igcm_coplus)=28.
	347	count=count+1
[3027]	348	idtracers(count) = "coplus"
[171]	349	endif
	350	if (noms(iq).eq."cplus") then
	351	igcm_cplus=iq
	352	mmol(igcm_cplus)=12.
	353	count=count+1
[3027]	354	idtracers(count) = "cplus"
[171]	355	endif
	356	if (noms(iq).eq."nplus") then
	357	igcm_nplus=iq
	358	mmol(igcm_nplus)=14.
	359	count=count+1
[3027]	360	idtracers(count) = "nplus"
[171]	361	endif
	362	if (noms(iq).eq."noplus") then
	363	igcm_noplus=iq
	364	mmol(igcm_noplus)=30.
	365	count=count+1
[3027]	366	idtracers(count) = "noplus"
[171]	367	endif
[324]	368	if (noms(iq).eq."n2plus") then
[171]	369	igcm_n2plus=iq
	370	mmol(igcm_n2plus)=28.
	371	count=count+1
[3027]	372	idtracers(count) = "n2plus"
[171]	373	endif
	374	if (noms(iq).eq."hplus") then
	375	igcm_hplus=iq
	376	mmol(igcm_hplus)=1.
	377	count=count+1
[3027]	378	idtracers(count) = "hplus"
[171]	379	endif
[635]	380	if (noms(iq).eq."hco2plus") then
	381	igcm_hco2plus=iq
	382	mmol(igcm_hco2plus)=45.
	383	count=count+1
[3027]	384	idtracers(count) = "hco2plus"
[635]	385	endif
[2284]	386	if (noms(iq).eq."hcoplus") then
	387	igcm_hcoplus=iq
	388	mmol(igcm_hcoplus)=29.
	389	count=count+1
[3027]	390	idtracers(count) = "hcoplus"
[2284]	391	endif
[2302]	392	if (noms(iq).eq."h2oplus") then
	393	igcm_h2oplus=iq
	394	mmol(igcm_h2oplus)=18.
	395	count=count+1
[3027]	396	idtracers(count) = "h2oplus"
[2302]	397	endif
[2321]	398	if (noms(iq).eq."h3oplus") then
	399	igcm_h3oplus=iq
	400	mmol(igcm_h3oplus)=19.
	401	count=count+1
[3027]	402	idtracers(count) = "h3oplus"
[2321]	403	endif
	404	if (noms(iq).eq."ohplus") then
	405	igcm_ohplus=iq
	406	mmol(igcm_ohplus)=17.
	407	count=count+1
[3027]	408	idtracers(count) = "ohplus"
[2321]	409	endif
[171]	410	if (noms(iq).eq."elec") then
	411	igcm_elec=iq
	412	mmol(igcm_elec)=1./1822.89
	413	count=count+1
[3027]	414	idtracers(count) = "elec"
[171]	415	endif
[38]	416	if (noms(iq).eq."h2o_vap") then
	417	igcm_h2o_vap=iq
	418	mmol(igcm_h2o_vap)=18.
	419	count=count+1
[3027]	420	idtracers(count) = "h2o_vap"
[38]	421	endif
[2312]	422	if (noms(iq).eq."hdo_vap") then
	423	igcm_hdo_vap=iq
[2462]	424	mmol(igcm_hdo_vap)=19.
[2312]	425	count=count+1
[3027]	426	idtracers(count) = "hdo_vap"
[2312]	427	endif
[2461]	428	if (noms(iq).eq."od") then
	429	igcm_od=iq
	430	mmol(igcm_od)=18.
	431	count=count+1
[3027]	432	idtracers(count) = "od"
[2461]	433	endif
	434	if (noms(iq).eq."d") then
	435	igcm_d=iq
	436	mmol(igcm_d)=2.
	437	count=count+1
[3027]	438	idtracers(count) = "d"
[2461]	439	endif
	440	if (noms(iq).eq."hd") then
	441	igcm_hd=iq
	442	mmol(igcm_hd)=3.
	443	count=count+1
[3027]	444	idtracers(count) = "hd"
[2461]	445	endif
	446	if (noms(iq).eq."do2") then
	447	igcm_do2=iq
	448	mmol(igcm_do2)=34.
	449	count=count+1
[3027]	450	idtracers(count) = "do2"
[2461]	451	endif
	452	if (noms(iq).eq."hdo2") then
	453	igcm_hdo2=iq
	454	mmol(igcm_hdo2)=35.
	455	count=count+1
[3027]	456	idtracers(count) = "hdo2"
[2461]	457	endif
[1617]	458	if (noms(iq).eq."co2_ice") then
	459	igcm_co2_ice=iq
	460	mmol(igcm_co2_ice)=44.
	461	count=count+1
[3027]	462	idtracers(count) = "co2_ice"
[1617]	463	endif
[38]	464	if (noms(iq).eq."h2o_ice") then
	465	igcm_h2o_ice=iq
	466	mmol(igcm_h2o_ice)=18.
	467	count=count+1
[3027]	468	idtracers(count) = "h2o_ice"
[38]	469	endif
[2312]	470	if (noms(iq).eq."hdo_ice") then
	471	igcm_hdo_ice=iq
[2463]	472	mmol(igcm_hdo_ice)=19.
[2312]	473	count=count+1
[3027]	474	idtracers(count) = "hdo_ice"
[2312]	475	endif
[38]	476	! Other stuff: e.g. for simulations using co2 + neutral gaz
	477	if (noms(iq).eq."Ar_N2") then
	478	igcm_ar_n2=iq
	479	mmol(igcm_ar_n2)=30.
	480	count=count+1
[3027]	481	idtracers(count) = "Ar_N2"
[38]	482	endif
[1720]	483	if (co2clouds) then
[1617]	484	if (noms(iq).eq."ccnco2_mass") then
	485	igcm_ccnco2_mass=iq
	486	count=count+1
[3027]	487	idtracers(count) = "ccnco2_mass"
[1617]	488	endif
	489	if (noms(iq).eq."ccnco2_number") then
	490	igcm_ccnco2_number=iq
	491	count=count+1
[3027]	492	idtracers(count) = "ccnco2_number"
[1617]	493	endif
[2589]	494	if (meteo_flux) then
	495	if (noms(iq).eq."ccnco2_meteor_mass") then
	496	igcm_ccnco2_meteor_mass=iq
	497	count=count+1
[3027]	498	idtracers(count) = "ccnco2_meteor_mass"
[2589]	499	endif
	500	if (noms(iq).eq."ccnco2_meteor_number") then
	501	igcm_ccnco2_meteor_number=iq
	502	count=count+1
[3027]	503	idtracers(count) = "ccnco2_meteor_number"
[2589]	504	endif
	505	end if
[2562]	506	if (co2useh2o) then
	507	if (noms(iq).eq."ccnco2_h2o_number") then
	508	igcm_ccnco2_h2o_number=iq
	509	count=count+1
[3027]	510	idtracers(count) = "ccnco2_h2o_number"
[2562]	511	endif
	512	if (noms(iq).eq."ccnco2_h2o_mass_ice") then
	513	igcm_ccnco2_h2o_mass_ice=iq
	514	count=count+1
[3027]	515	idtracers(count) = "ccnco2_h2o_mass_ice"
[2562]	516	endif
	517	if (noms(iq).eq."ccnco2_h2o_mass_ccn") then
	518	igcm_ccnco2_h2o_mass_ccn=iq
	519	count=count+1
[3027]	520	idtracers(count) = "ccnco2_h2o_mass_ccn"
[2562]	521	endif
	522	end if
[1617]	523	endif
	524	enddo ! of do iq=1,nq
	525
[38]	526	! check that we identified all tracers:
[3027]	527	if (count /= nq) then
	528	write(,) "initracer: found only",count,"tracers but",
	529	&nq,"are expected according to traceur.def!"
	530	write(,) "initracer: the only identified tracers are:"
	531	do iq = 1,count
	532	write(,) ' ',iq,trim(idtracers(iq))
	533	enddo
[2302]	534	call abort_physic("initracer","tracer mismatch",1)
[38]	535	else
	536	write(,) "initracer: found all expected tracers, namely:"
[3027]	537	do iq = 1,nq
	538	write(,) ' ',iq,' ',trim(noms(iq))
[38]	539	enddo
	540	endif
	541
	542	! if water cycle but iceparty=.false., there will nevertheless be
	543	! water ice at the surface (iceparty is not used anymore, but this
	544	! part is still relevant, as we want to stay compatible with the
	545	! older versions).
	546	if (water.and.(igcm_h2o_ice.eq.0)) then
	547	igcm_h2o_ice=igcm_h2o_vap ! so that qsurf(i_h2o_ice) is identified
	548	! even though there is no q(i_h2o_ice)
	549	else
	550	! surface ice qsurf(i_h2o_ice) was loaded twice by phyetat0,
	551	! as qsurf(i_h2o_vap) & as qsurf(i_h2o_ice), so to be clean:
	552	if (igcm_h2o_vap.ne.0) then
[1036]	553	qsurf(1:ngrid,igcm_h2o_vap)=0
[38]	554	endif
	555	endif
	556
[2312]	557	! Additional test required for HDO
	558	! We need to compute some things for H2O before HDO
	559	if (hdo) then
	560	if (igcm_h2o_vap.gt.igcm_hdo_vap) then
[2398]	561	call abort_physic("initracer",
	562	& "Tracer H2O must be initialized before HDO",1)
[2324]	563	else if ((nqfils(igcm_h2o_ice).lt.1)
	564	& .or. (nqfils(igcm_h2o_vap).lt.1)) then
	565	write(,) "HDO must be transported as a son",
	566	& " of H2O: complete traceur.def"
[2398]	567	call abort_physic("initracer","adapt your tracer.def",1)
[2324]	568	else if ((igcm_hdo_ice.lt.nq-2)
	569	& .or. (igcm_hdo_vap.lt.nq-2)) then
	570	write(,) "The isotopes (HDO) must be placed at",
	571	& " the end of the list in traceur.def"
[2398]	572	call abort_physic("initracer","adapt your tracer.def",1)
[2312]	573	endif
	574	endif
	575
[38]	576	c------------------------------------------------------------
[1036]	577	c Initialize tracers .... (in tracer_mod)
[38]	578	c------------------------------------------------------------
[1005]	579	! start by setting everything to (default) zero
[1036]	580	rho_q(1:nq)=0 ! tracer density (kg.m-3)
	581	radius(1:nq)=0. ! tracer particle radius (m)
	582	alpha_lift(1:nq) =0. ! tracer saltation vertical flux/horiz flux ratio (m-1)
	583	alpha_devil(1:nq)=0. ! tracer lifting coefficient by dust devils
[38]	584
[1005]	585
	586	! some reference values
[38]	587	rho_dust=2500. ! Mars dust density (kg.m-3)
	588	rho_ice=920. ! Water ice density (kg.m-3)
[1685]	589	rho_ice_co2=1650.
[38]	590
	591	if (doubleq) then
	592	c "doubleq" technique
	593	c -------------------
	594	c (transport of mass and number mixing ratio)
	595	c iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
	596
[2312]	597	if( (nq.lt.2).or.(water.and.(nq.lt.4))
	598	* .or.(hdo.and.(nq.lt.6) )) then
[1036]	599	write(,)'initracer: nq is too low : nq=', nq
[38]	600	write(,)'water= ',water,' doubleq= ',doubleq
	601	end if
	602
	603	nueff_lift = 0.5
	604	varian=sqrt(log(1.+nueff_lift))
	605
	606	rho_q(igcm_dust_mass)=rho_dust
	607	rho_q(igcm_dust_number)=rho_dust
	608
	609	c Intermediate calcul for computing geometric mean radius r0
	610	c as a function of mass and number mixing ratio Q and N
	611	c (r0 = (r3n_q * Q/ N)^(1/3))
	612	r3n_q = exp(-4.5varian2)(3./4.)/(pi*rho_dust)
	613
	614	c Intermediate calcul for computing effective radius reff
	615	c from geometric mean radius r0
	616	c (reff = ref_r0 * r0)
	617	ref_r0 = exp(2.5varian*2)
	618
	619	c lifted dust :
	620	c '''''''''''
	621	reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m)
	622	alpha_devil(igcm_dust_mass)=9.e-9 ! dust devil lift mass coeff
	623	c alpha_lift(igcm_dust_mass)=3.0e-15 ! Lifted mass coeff
[1455]	624
	625	!! default lifting settings
	626	!! -- GCM: alpha_lift not zero because large-scale lifting by default
	627	!! -- MESOSCALE: alpha_lift zero because no lifting at all in mesoscale by default
	628	#ifdef MESOSCALE
	629	alpha_lift(igcm_dust_mass)=0.0
	630	#else
[38]	631	alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
[1974]	632	IF (dustinjection.ge.1) THEN
	633	reff_lift = 3.0e-6 ! Effective radius of lifted dust (m)
	634	alpha_lift(igcm_dust_mass)=(4/3.)reff_liftrho_dust
	635	& /2.4
	636	ENDIF
[1455]	637	#endif
[38]	638
	639	r0_lift = reff_lift/ref_r0
	640	alpha_devil(igcm_dust_number)=r3n_q*
	641	& alpha_devil(igcm_dust_mass)/r0_lift**3
	642	alpha_lift(igcm_dust_number)=r3n_q*
	643	& alpha_lift(igcm_dust_mass)/r0_lift**3
	644
	645	radius(igcm_dust_mass) = reff_lift
	646	radius(igcm_dust_number) = reff_lift
	647
	648	write(,) "initracer: doubleq_param reff_lift:", reff_lift
	649	write(,) "initracer: doubleq_param nueff_lift:", nueff_lift
	650	write(,) "initracer: doubleq_param alpha_lift:",
	651	& alpha_lift(igcm_dust_mass)
[1974]	652	!c ----------------------------------------------------------------------
	653	!c rocket dust storm scheme
	654	!c lifting tracer stormdust using same distribution than
	655	!c normal dust
	656	if (rdstorm) then
	657	reff_storm=3.e-6 ! reff_lift !3.e-6
	658	r0_storm=reff_storm/ref_r0
	659	rho_q(igcm_stormdust_mass)=rho_dust
	660	rho_q(igcm_stormdust_number)=rho_dust
	661
	662	alpha_devil(igcm_stormdust_mass)=9.e-9 ! dust devil lift mass coeff
	663	alpha_lift(igcm_stormdust_mass)=4./3./2.4reff_stormrho_dust
	664
	665	write(,) 'alpha_lift(rds):',alpha_lift(igcm_stormdust_mass)
	666
	667	alpha_devil(igcm_stormdust_number)=r3n_q*
	668	& alpha_devil(igcm_stormdust_mass)/r0_storm**3
	669	alpha_lift(igcm_stormdust_number)=r3n_q*
	670	& alpha_lift(igcm_stormdust_mass)/r0_storm**3
	671
	672	radius(igcm_stormdust_mass) = reff_storm
	673	radius(igcm_stormdust_number) = reff_storm
	674	end if !(rdstorm)
	675	!c ----------------------------------------------------------------------
[2628]	676	!c mountain top dust flows scheme
[2199]	677	!c you need a radius value for topdust to active its sedimentation
	678	!c we take the same value as for the normal dust
[2628]	679	if (topflows) then
[2199]	680	rho_q(igcm_topdust_mass)=rho_dust
	681	rho_q(igcm_topdust_number)=rho_dust
	682	radius(igcm_topdust_mass) = 3.e-6
	683	radius(igcm_topdust_number) = 3.e-6
[2628]	684	end if !(topflows)
[2199]	685	!c ----------------------------------------------------------------------
[1974]	686
[38]	687	else
	688
[648]	689	! initialize varian, which may be used (e.g. by surfacearea)
	690	! even with conrath dust
	691	nueff_lift = 0.5
	692	varian=sqrt(log(1.+nueff_lift))
	693
[38]	694	if (dustbin.gt.1) then
	695	print*,'initracer: STOP!',
	696	$ ' properties of dust need to be set in initracer !!!'
[2302]	697	call abort_physic("initracer","dustbin properties issue",1)
[38]	698
	699	else if (dustbin.eq.1) then
	700
	701	c This will be used for 1 dust particle size:
	702	c ------------------------------------------
	703	radius(igcm_dustbin(1))=3.e-6
	704	alpha_lift(igcm_dustbin(1))=0.0e-6
	705	alpha_devil(igcm_dustbin(1))=7.65e-9
	706	rho_q(igcm_dustbin(1))=rho_dust
	707
	708	endif
	709	end if ! (doubleq)
	710
[358]	711
	712	c Scavenging of dust particles by H2O clouds:
	713	c ------------------------------------------
	714	c Initialize the two tracers used for the CCNs
	715	if (water.AND.doubleq.AND.scavenging) then
	716	radius(igcm_ccn_mass) = radius(igcm_dust_mass)
	717	alpha_lift(igcm_ccn_mass) = 1e-30
	718	alpha_devil(igcm_ccn_mass) = 1e-30
	719	rho_q(igcm_ccn_mass) = rho_dust
	720
	721	radius(igcm_ccn_number) = radius(igcm_ccn_mass)
	722	alpha_lift(igcm_ccn_number) = alpha_lift(igcm_ccn_mass)
	723	alpha_devil(igcm_ccn_number) = alpha_devil(igcm_ccn_mass)
	724	rho_q(igcm_ccn_number) = rho_q(igcm_ccn_mass)
	725	endif ! of if (water.AND.doubleq.AND.scavenging)
	726
[38]	727	c Submicron dust mode:
	728	c --------------------
	729
	730	if (submicron) then
	731	radius(igcm_dust_submicron)=0.1e-6
	732	rho_q(igcm_dust_submicron)=rho_dust
	733	if (doubleq) then
	734	c If doubleq is also active, we use the population ratio:
	735	alpha_lift(igcm_dust_submicron) =
	736	& alpha_lift(igcm_dust_number)popratio
	737	& rho_q(igcm_dust_submicron)4./3.pi*
	738	& radius(igcm_dust_submicron)**3.
	739	alpha_devil(igcm_dust_submicron)=1.e-30
	740	else
	741	alpha_lift(igcm_dust_submicron)=1e-6
	742	alpha_devil(igcm_dust_submicron)=1.e-30
	743	endif ! (doubleq)
	744	end if ! (submicron)
	745
	746	c Initialization for water vapor
	747	c ------------------------------
	748	if(water) then
	749	radius(igcm_h2o_vap)=0.
	750	alpha_lift(igcm_h2o_vap) =0.
	751	alpha_devil(igcm_h2o_vap)=0.
[1036]	752	if(water.and.(nq.ge.2)) then
[38]	753	radius(igcm_h2o_ice)=3.e-6
	754	rho_q(igcm_h2o_ice)=rho_ice
	755	alpha_lift(igcm_h2o_ice) =0.
	756	alpha_devil(igcm_h2o_ice)=0.
[1036]	757	elseif(water.and.(nq.lt.2)) then
	758	write(,) 'nq is too low : nq=', nq
[38]	759	write(,) 'water= ',water
	760	endif
	761
	762	end if ! (water)
[2312]	763
	764	c Initialization for hdo vapor
	765	c ------------------------------
	766	if (hdo) then
	767	radius(igcm_hdo_vap)=0.
	768	alpha_lift(igcm_hdo_vap) =0.
	769	alpha_devil(igcm_hdo_vap)=0.
	770	if(water.and.(nq.ge.2)) then
	771	radius(igcm_hdo_ice)=3.e-6
	772	rho_q(igcm_hdo_ice)=rho_ice
	773	alpha_lift(igcm_hdo_ice) =0.
	774	alpha_devil(igcm_hdo_ice)=0.
	775	elseif(hdo.and.(nq.lt.6)) then
	776	write(,) 'nq is too low : nq=', nq
	777	write(,) 'hdo= ',hdo
	778	endif
	779
	780	end if ! (hdo)
	781
[1617]	782
[38]	783
	784	c Output for records:
	785	c ~~~~~~~~~~~~~~~~~~
	786	write(,)
	787	Write(,) '******** initracer : dust transport parameters :'
	788	write(,) 'alpha_lift = ', alpha_lift
	789	write(,) 'alpha_devil = ', alpha_devil
	790	write(,) 'radius = ', radius
	791	if(doubleq) then
	792	write(,) 'reff_lift (um) = ', reff_lift
	793	write(,) 'size distribution variance = ', varian
	794	write(,) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
	795	end if
	796
	797	!
	798	! some extra (possibly redundant) sanity checks for tracers:
	799	! ---------------------------------------------------------
	800
	801	if (doubleq) then
	802	! verify that we indeed have dust_mass and dust_number tracers
	803	if (igcm_dust_mass.eq.0) then
	804	write(,) "initracer: error !!"
	805	write(,) " cannot use doubleq option without ",
	806	& "a dust_mass tracer !"
[2302]	807	call abort_physic("initracer","doubleq issue",1)
[38]	808	endif
	809	if (igcm_dust_number.eq.0) then
	810	write(,) "initracer: error !!"
	811	write(,) " cannot use doubleq option without ",
	812	& "a dust_number tracer !"
[2302]	813	call abort_physic("initracer","doubleq issue",1)
[38]	814	endif
	815	endif
	816
	817	if ((.not.doubleq).and.(dustbin.gt.0)) then
	818	! verify that we indeed have 'dustbin' dust tracers
	819	count=0
	820	do iq=1,dustbin
	821	if (igcm_dustbin(iq).ne.0) then
	822	count=count+1
	823	endif
	824	enddo
	825	if (count.ne.dustbin) then
	826	write(,) "initracer: error !!"
[2302]	827	write(,) " dustbin is set to ",dustbin,
[38]	828	& " but we only have the following dust tracers:"
	829	do iq=1,count
	830	write(,)" ",trim(noms(igcm_dustbin(iq)))
	831	enddo
[2302]	832	call abort_physic("initracer","dustbin issue",1)
[38]	833	endif
	834	endif
	835
	836	if (water) then
	837	! verify that we indeed have h2o_vap and h2o_ice tracers
	838	if (igcm_h2o_vap.eq.0) then
	839	write(,) "initracer: error !!"
	840	write(,) " cannot use water option without ",
	841	& "an h2o_vap tracer !"
[2302]	842	call abort_physic("initracer","water cycle issue",1)
[38]	843	endif
	844	if (igcm_h2o_ice.eq.0) then
	845	write(,) "initracer: error !!"
	846	write(,) " cannot use water option without ",
	847	& "an h2o_ice tracer !"
[2302]	848	call abort_physic("initracer","water cycle issue",1)
[38]	849	endif
	850	endif
	851
[2312]	852	if (hdo) then
	853	! verify that we indeed have hdo_vap and hdo_ice tracers
	854	if (igcm_hdo_vap.eq.0) then
	855	write(,) "initracer: error !!"
	856	write(,) " cannot use hdo option without ",
	857	& "an hdo_vap tracer !"
[2398]	858	call abort_physic("initracer","hdo cycle issue",1)
[2312]	859	endif
	860	if (igcm_hdo_ice.eq.0) then
	861	write(,) "initracer: error !!"
	862	write(,) " cannot use hdo option without ",
	863	& "an hdo_ice tracer !"
[2398]	864	call abort_physic("initracer","hdo cycle issue",1)
[2312]	865	endif
	866	endif
	867
	868
[1617]	869	if (co2clouds) then
	870	!verify that we have co2_ice and co2 tracers
	871	if (igcm_co2 .eq. 0) then
	872	write(,) "initracer: error !!"
	873	write(,) " cannot use co2 clouds option without ",
	874	& "a co2 tracer !"
[2302]	875	call abort_physic("initracer","co2 clouds issue",1)
[2589]	876	end if
[1617]	877	if (igcm_co2_ice .eq. 0) then
	878	write(,) "initracer: error !!"
	879	write(,) " cannot use co2 clouds option without ",
	880	& "a co2_ice tracer !"
[2302]	881	call abort_physic("initracer","co2 clouds issue",1)
[2589]	882	end if
[2659]	883	if (igcm_ccnco2_number .eq. 0) then
	884	write(,) "initracer: error !!"
	885	write(,) " cannot use co2 clouds option without ",
	886	& "a ccnco2_number tracer !"
	887	call abort_physic("initracer","co2 clouds issue",1)
	888	end if
	889	if (igcm_ccnco2_mass .eq. 0) then
	890	write(,) "initracer: error !!"
	891	write(,) " cannot use co2 clouds option without ",
	892	& "a ccnco2_mass tracer !"
	893	call abort_physic("initracer","co2 clouds issue",1)
	894	end if
	895	if (co2useh2o) then
	896	if (igcm_ccnco2_h2o_number .eq. 0) then
	897	write(,) "initracer: error !!"
	898	write(,) " cannot use co2 clouds option without ",
	899	& "a ccnco2_h2o_number tracer !"
	900	call abort_physic("initracer","co2 clouds issue",1)
	901	end if
	902	if (igcm_ccnco2_h2o_mass_ice .eq. 0) then
	903	write(,) "initracer: error !!"
	904	write(,) " cannot use co2 clouds option without ",
	905	& "a ccnco2_h2o_mass_ice tracer !"
	906	call abort_physic("initracer","co2 clouds issue",1)
	907	end if
	908	if (igcm_ccnco2_h2o_mass_ccn .eq. 0) then
	909	write(,) "initracer: error !!"
	910	write(,) " cannot use co2 clouds option without ",
	911	& "a ccnco2_h2o_mass_ccn tracer !"
	912	call abort_physic("initracer","co2 clouds issue",1)
	913	end if
	914	end if
	915	if (meteo_flux) then
	916	if (igcm_ccnco2_meteor_number .eq. 0) then
	917	write(,) "initracer: error !!"
	918	write(,) " cannot use co2 clouds option without ",
	919	& "a ccnco2_meteor_number tracer !"
	920	call abort_physic("initracer","co2 clouds issue",1)
	921	end if
	922	if (igcm_ccnco2_meteor_mass .eq. 0) then
	923	write(,) "initracer: error !!"
	924	write(,) " cannot use co2 clouds option without ",
	925	& "a ccnco2_h2o_mass_ice tracer !"
	926	call abort_physic("initracer","co2 clouds issue",1)
	927	end if
	928	if (igcm_ccnco2_meteor_mass .eq. 0) then
	929	write(,) "initracer: error !!"
	930	write(,) " cannot use co2 clouds option without ",
	931	& "a ccnco2_meteor_mass tracer !"
	932	call abort_physic("initracer","co2 clouds issue",1)
	933	end if
	934	end if
[1617]	935	endif
[1974]	936
	937	if (rdstorm) then
	938	! verify that we indeed have stormdust_mass and stormdust_number tracers
	939	if (igcm_stormdust_mass.eq.0) then
	940	write(,) "initracer: error !!"
	941	write(,) " cannot use rdstorm option without ",
	942	& "a stormdust_mass tracer !"
[2302]	943	call abort_physic("initracer","rdstorm issue",1)
[1974]	944	endif
	945	if (igcm_stormdust_number.eq.0) then
	946	write(,) "initracer: error !!"
	947	write(,) " cannot use rdstorm option without ",
	948	& "a stormdust_number tracer !"
[2302]	949	call abort_physic("initracer","rdstorm issue",1)
[1974]	950	endif
	951	endif
[2199]	952
[2628]	953	if (topflows) then
[2199]	954	! verify that we indeed have topdust_mass and topdust_number tracers
	955	if (igcm_topdust_mass.eq.0) then
	956	write(,) "initracer: error !!"
[2628]	957	write(,) " cannot use topflows option without ",
[2199]	958	& "a topdust_mass tracer !"
[2628]	959	call abort_physic("initracer","topflows issue",1)
[2199]	960	endif
	961	if (igcm_topdust_number.eq.0) then
	962	write(,) "initracer: error !!"
[2628]	963	write(,) " cannot use topflows option without ",
[2199]	964	& "a topdust_number tracer !"
[2628]	965	call abort_physic("initracer","topflows issue",1)
[2199]	966	endif
	967	endif
[1974]	968
[1380]	969	if (callnlte) then ! NLTE requirements
	970	if (nltemodel.ge.1) then
	971	! check that co2, co, o and n2 tracers are available
	972	if (igcm_co2.eq.0) then
	973	write(,) "initracer: error !!"
	974	write(,) " with nltemodel>0, we need the co2 tracer!"
[2302]	975	call abort_physic("initracer","missing co2 tracer",1)
[1380]	976	endif
	977	if (igcm_co.eq.0) then
	978	write(,) "initracer: error !!"
	979	write(,) " with nltemodel>0, we need the co tracer!"
[2302]	980	call abort_physic("initracer","missing co tracer",1)
[1380]	981	endif
	982	if (igcm_o.eq.0) then
	983	write(,) "initracer: error !!"
	984	write(,) " with nltemodel>0, we need the o tracer!"
[2302]	985	call abort_physic("initracer","missing o tracer",1)
[1380]	986	endif
	987	if (igcm_n2.eq.0) then
	988	write(,) "initracer: error !!"
	989	write(,) " with nltemodel>0, we need the n2 tracer!"
[2302]	990	call abort_physic("initracer","missing n2 tracer",1)
[1380]	991	endif
	992	endif
	993	endif
	994
[358]	995	if (scavenging) then
	996	! verify that we indeed have ccn_mass and ccn_number tracers
[1617]	997	if (igcm_ccn_mass.eq.0 .and. igcm_ccnco2_mass.eq.0) then
[358]	998	write(,) "initracer: error !!"
	999	write(,) " cannot use scavenging option without ",
[1617]	1000	& "a ccn_mass or ccnco2_mass tracer !"
[2302]	1001	call abort_physic("initracer","scavenging issue",1)
[358]	1002	endif
[1617]	1003	if (igcm_ccn_number.eq.0 .and. igcm_ccnco2_number.eq.0 ) then
[358]	1004	write(,) "initracer: error !!"
	1005	write(,) " cannot use scavenging option without ",
[1617]	1006	& "a ccn_number or ccnco2_number tracer !"
[2302]	1007	call abort_physic("initracer","scavenging issue",1)
[358]	1008	endif
	1009	endif ! of if (scavenging)
	1010
[38]	1011	if (photochem .or. callthermos) then
	1012	! verify that we indeed have the chemistry tracers
	1013	if (igcm_co2.eq.0) then
	1014	write(,) "initracer: error !!"
	1015	write(,) " cannot use chemistry option without ",
	1016	& "a co2 tracer !"
[2302]	1017	call abort_physic("initracer","missing co2 tracer",1)
[38]	1018	endif
	1019	if (igcm_co.eq.0) then
	1020	write(,) "initracer: error !!"
	1021	write(,) " cannot use chemistry option without ",
	1022	& "a co tracer !"
[2302]	1023	call abort_physic("initracer","missing co tracer",1)
[38]	1024	endif
	1025	if (igcm_o.eq.0) then
	1026	write(,) "initracer: error !!"
	1027	write(,) " cannot use chemistry option without ",
	1028	& "a o tracer !"
[2302]	1029	call abort_physic("initracer","missing o tracer",1)
[38]	1030	endif
	1031	if (igcm_o1d.eq.0) then
	1032	write(,) "initracer: error !!"
	1033	write(,) " cannot use chemistry option without ",
	1034	& "a o1d tracer !"
[2302]	1035	call abort_physic("initracer","missing o1d tracer",1)
[38]	1036	endif
	1037	if (igcm_o2.eq.0) then
	1038	write(,) "initracer: error !!"
	1039	write(,) " cannot use chemistry option without ",
	1040	& "an o2 tracer !"
[2302]	1041	call abort_physic("initracer","missing o2 tracer",1)
[38]	1042	endif
	1043	if (igcm_o3.eq.0) then
	1044	write(,) "initracer: error !!"
	1045	write(,) " cannot use chemistry option without ",
	1046	& "an o3 tracer !"
[2302]	1047	call abort_physic("initracer","missing o3 tracer",1)
[38]	1048	endif
	1049	if (igcm_h.eq.0) then
	1050	write(,) "initracer: error !!"
	1051	write(,) " cannot use chemistry option without ",
	1052	& "a h tracer !"
[2302]	1053	call abort_physic("initracer","missing h tracer",1)
[38]	1054	endif
	1055	if (igcm_h2.eq.0) then
	1056	write(,) "initracer: error !!"
	1057	write(,) " cannot use chemistry option without ",
	1058	& "a h2 tracer !"
[2302]	1059	call abort_physic("initracer","missing h2 tracer",1)
[38]	1060	endif
	1061	if (igcm_oh.eq.0) then
	1062	write(,) "initracer: error !!"
	1063	write(,) " cannot use chemistry option without ",
	1064	& "an oh tracer !"
[2302]	1065	call abort_physic("initracer","missing oh tracer",1)
[38]	1066	endif
	1067	if (igcm_ho2.eq.0) then
	1068	write(,) "initracer: error !!"
	1069	write(,) " cannot use chemistry option without ",
	1070	& "a ho2 tracer !"
[2302]	1071	call abort_physic("initracer","missing ho2 tracer",1)
[1720]	1072	endif
[38]	1073	if (igcm_h2o2.eq.0) then
	1074	write(,) "initracer: error !!"
	1075	write(,) " cannot use chemistry option without ",
	1076	& "a h2o2 tracer !"
[2302]	1077	call abort_physic("initracer","missing h2o2 tracer",1)
[38]	1078	endif
	1079	if (igcm_n2.eq.0) then
	1080	write(,) "initracer: error !!"
	1081	write(,) " cannot use chemistry option without ",
	1082	& "a n2 tracer !"
[2302]	1083	call abort_physic("initracer","missing n2 tracer",1)
[38]	1084	endif
	1085	if (igcm_ar.eq.0) then
	1086	write(,) "initracer: error !!"
	1087	write(,) " cannot use chemistry option without ",
	1088	& "an ar tracer !"
[2302]	1089	call abort_physic("initracer","missing ar tracer",1)
[38]	1090	endif
	1091	endif ! of if (photochem .or. callthermos)
	1092
[2659]	1093	! Initialisation for CO2 clouds
	1094	if (co2clouds) then
	1095	radius(igcm_ccnco2_mass) = radius(igcm_dust_mass)
	1096	alpha_lift(igcm_ccnco2_mass) = 1e-30
	1097	alpha_devil(igcm_ccnco2_mass) = 1e-30
	1098	rho_q(igcm_ccnco2_mass) = rho_dust
	1099	radius(igcm_ccnco2_number) = radius(igcm_ccnco2_mass)
	1100	alpha_lift(igcm_ccnco2_number) = alpha_lift(igcm_ccnco2_mass)
	1101	alpha_devil(igcm_ccnco2_number) = alpha_devil(igcm_ccnco2_mass)
	1102	rho_q(igcm_ccnco2_number) = rho_q(igcm_ccnco2_mass)
	1103
	1104	radius(igcm_co2)=0.
	1105	alpha_lift(igcm_co2) =0.
	1106	alpha_devil(igcm_co2)=0.
	1107	radius(igcm_co2_ice)=1.e-8
	1108	rho_q(igcm_co2_ice)=rho_ice_co2
	1109	alpha_lift(igcm_co2_ice) =0.
	1110	alpha_devil(igcm_co2_ice)=0.
	1111
	1112	endif
	1113
[38]	1114	end

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

Download in other formats:

Original Format