Context Navigation

initracer.F @ 3571

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

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

Download in other formats:

Original Format