Context Navigation

initracer.F @ 3118

Last change on this file since 3118 was 3027, checked in by jbclement, 16 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

Line
1	SUBROUTINE initracer(ngrid,nq,qsurf)
2
3	use tracer_mod
4	use comcstfi_h, only: pi
5	use dust_param_mod, only: doubleq, submicron, dustbin
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 :
16	c Initialize tracer related data in tracer_mod, using tracer names provided
17	c by the dynamics in "infotrac"
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
27	include "callkeys.h"
28
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
33	integer iq,ig,count
34	real r0_lift , reff_lift, nueff_lift
35	real r0_storm,reff_storm
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.
40	character(len=30) :: txt ! to store some text
41	character(len=30), dimension(nq) :: idtracers
42
43	c-----------------------------------------------------------------------
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
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
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
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:
65	igcm_dustbin(1:nq)=0
66	igcm_co2_ice=0
67	igcm_ccnco2_mass=0
68	igcm_ccnco2_number=0
69	igcm_ccnco2_meteor_mass=0
70	igcm_ccnco2_meteor_number=0
71	igcm_ccnco2_h2o_mass_ice=0
72	igcm_ccnco2_h2o_mass_ccn=0
73	igcm_ccnco2_h2o_number=0
74	igcm_dust_mass=0
75	igcm_dust_number=0
76	igcm_ccn_mass=0
77	igcm_ccn_number=0
78	igcm_dust_submicron=0
79	igcm_h2o_vap=0
80	igcm_h2o_ice=0
81	igcm_hdo_vap=0
82	igcm_hdo_ice=0
83	igcm_stormdust_mass=0
84	igcm_stormdust_number=0
85	igcm_topdust_mass=0
86	igcm_topdust_number=0
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
94	igcm_d=0
95	igcm_hd=0
96	igcm_h2=0
97	igcm_od=0
98	igcm_do2=0
99	igcm_hdo2=0
100	igcm_oh=0
101	igcm_ho2=0
102	igcm_h2o2=0
103	igcm_ch4=0
104	igcm_n2=0
105	igcm_ar=0
106	igcm_ar_n2=0
107	igcm_n=0
108	igcm_no=0
109	igcm_no2=0
110	igcm_n2d=0
111	igcm_he=0
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
121	igcm_hco2plus=0
122	igcm_hcoplus=0
123	igcm_h2oplus=0
124	igcm_h3oplus=0
125	igcm_ohplus=0
126	igcm_elec=0
127
128	! 1. find dust tracers
129	count=0
130	if (dustbin.gt.0) then
131	do iq=1,nq
132	txt=" "
133	write(txt,'(a4,i2.2)')'dust',count+1
134	if (noms(iq).eq.txt) then
135	count=count+1
136	idtracers(count) = txt
137	igcm_dustbin(count)=iq
138	mmol(iq)=100.
139	endif
140	enddo !do iq=1,nq
141	endif ! of if (dustbin.gt.0)
142	if (doubleq) then
143	do iq=1,nq
144	if (noms(iq).eq."dust_mass") then
145	igcm_dust_mass=iq
146	count=count+1
147	idtracers(count) = "dust_mass"
148	endif
149	if (noms(iq).eq."dust_number") then
150	igcm_dust_number=iq
151	count=count+1
152	idtracers(count) = "dust_number"
153	endif
154	enddo
155	endif ! of if (doubleq)
156	if (microphys) then
157	do iq=1,nq
158	if (noms(iq).eq."ccn_mass") then
159	igcm_ccn_mass=iq
160	count=count+1
161	idtracers(count) = "ccn_mass"
162	endif
163	if (noms(iq).eq."ccn_number") then
164	igcm_ccn_number=iq
165	count=count+1
166	idtracers(count) = "ccn_number"
167	endif
168	enddo
169	endif ! of if (microphys)
170	if (submicron) then
171	do iq=1,nq
172	if (noms(iq).eq."dust_submicron") then
173	igcm_dust_submicron=iq
174	mmol(iq)=100.
175	count=count+1
176	idtracers(count) = "dust_submicron"
177	endif
178	enddo
179	endif ! of if (submicron)
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
185	idtracers(count) = "stormdust_mass"
186	endif
187	if (noms(iq).eq."stormdust_number") then
188	igcm_stormdust_number=iq
189	count=count+1
190	idtracers(count) = "stormdust_number"
191	endif
192	enddo
193	endif ! of if (rdstorm)
194	if (topflows) then
195	do iq=1,nq
196	if (noms(iq).eq."topdust_mass") then
197	igcm_topdust_mass=iq
198	count=count+1
199	idtracers(count) = "topdust_mass"
200	endif
201	if (noms(iq).eq."topdust_number") then
202	igcm_topdust_number=iq
203	count=count+1
204	idtracers(count) = "topdust_number"
205	endif
206	enddo
207	endif ! of if (topflows)
208	! 2. find chemistry and water tracers
209	do iq=1,nq
210	if (noms(iq).eq."co2") then
211	igcm_co2=iq
212	mmol(igcm_co2)=44.
213	count=count+1
214	idtracers(count) = "co2"
215	endif
216	if (noms(iq).eq."co") then
217	igcm_co=iq
218	mmol(igcm_co)=28.
219	count=count+1
220	idtracers(count) = "co"
221	endif
222	if (noms(iq).eq."o") then
223	igcm_o=iq
224	mmol(igcm_o)=16.
225	count=count+1
226	idtracers(count) = "o"
227	endif
228	if (noms(iq).eq."o1d") then
229	igcm_o1d=iq
230	mmol(igcm_o1d)=16.
231	count=count+1
232	idtracers(count) = "o1d"
233	endif
234	if (noms(iq).eq."o2") then
235	igcm_o2=iq
236	mmol(igcm_o2)=32.
237	count=count+1
238	idtracers(count) = "o2"
239	endif
240	if (noms(iq).eq."o3") then
241	igcm_o3=iq
242	mmol(igcm_o3)=48.
243	count=count+1
244	idtracers(count) = "o3"
245	endif
246	if (noms(iq).eq."h") then
247	igcm_h=iq
248	mmol(igcm_h)=1.
249	count=count+1
250	idtracers(count) = "h"
251	endif
252	if (noms(iq).eq."h2") then
253	igcm_h2=iq
254	mmol(igcm_h2)=2.
255	count=count+1
256	idtracers(count) = "h2"
257	endif
258	if (noms(iq).eq."oh") then
259	igcm_oh=iq
260	mmol(igcm_oh)=17.
261	count=count+1
262	idtracers(count) = "oh"
263	endif
264	if (noms(iq).eq."ho2") then
265	igcm_ho2=iq
266	mmol(igcm_ho2)=33.
267	count=count+1
268	idtracers(count) = "ho2"
269	endif
270	if (noms(iq).eq."h2o2") then
271	igcm_h2o2=iq
272	mmol(igcm_h2o2)=34.
273	count=count+1
274	idtracers(count) = "h2o2"
275	endif
276	if (noms(iq).eq."n2") then
277	igcm_n2=iq
278	mmol(igcm_n2)=28.
279	count=count+1
280	idtracers(count) = "n2"
281	endif
282	if (noms(iq).eq."ch4") then
283	igcm_ch4=iq
284	mmol(igcm_ch4)=16.
285	count=count+1
286	idtracers(count) = "ch4"
287	endif
288	if (noms(iq).eq."ar") then
289	igcm_ar=iq
290	mmol(igcm_ar)=40.
291	count=count+1
292	idtracers(count) = "ar"
293	endif
294	if (noms(iq).eq."n") then
295	igcm_n=iq
296	mmol(igcm_n)=14.
297	count=count+1
298	idtracers(count) = "n"
299	endif
300	if (noms(iq).eq."no") then
301	igcm_no=iq
302	mmol(igcm_no)=30.
303	count=count+1
304	idtracers(count) = "no"
305	endif
306	if (noms(iq).eq."no2") then
307	igcm_no2=iq
308	mmol(igcm_no2)=46.
309	count=count+1
310	idtracers(count) = "no2"
311	endif
312	if (noms(iq).eq."n2d") then
313	igcm_n2d=iq
314	mmol(igcm_n2d)=28.
315	count=count+1
316	idtracers(count) = "n2d"
317	endif
318	if (noms(iq).eq."he") then
319	igcm_he=iq
320	mmol(igcm_he)=4.
321	count=count+1
322	idtracers(count) = "he"
323	endif
324	if (noms(iq).eq."co2plus") then
325	igcm_co2plus=iq
326	mmol(igcm_co2plus)=44.
327	count=count+1
328	idtracers(count) = "co2plus"
329	endif
330	if (noms(iq).eq."oplus") then
331	igcm_oplus=iq
332	mmol(igcm_oplus)=16.
333	count=count+1
334	idtracers(count) = "oplus"
335	endif
336	if (noms(iq).eq."o2plus") then
337	igcm_o2plus=iq
338	mmol(igcm_o2plus)=32.
339	count=count+1
340	idtracers(count) = "o2plus"
341	endif
342	if (noms(iq).eq."coplus") then
343	igcm_coplus=iq
344	mmol(igcm_coplus)=28.
345	count=count+1
346	idtracers(count) = "coplus"
347	endif
348	if (noms(iq).eq."cplus") then
349	igcm_cplus=iq
350	mmol(igcm_cplus)=12.
351	count=count+1
352	idtracers(count) = "cplus"
353	endif
354	if (noms(iq).eq."nplus") then
355	igcm_nplus=iq
356	mmol(igcm_nplus)=14.
357	count=count+1
358	idtracers(count) = "nplus"
359	endif
360	if (noms(iq).eq."noplus") then
361	igcm_noplus=iq
362	mmol(igcm_noplus)=30.
363	count=count+1
364	idtracers(count) = "noplus"
365	endif
366	if (noms(iq).eq."n2plus") then
367	igcm_n2plus=iq
368	mmol(igcm_n2plus)=28.
369	count=count+1
370	idtracers(count) = "n2plus"
371	endif
372	if (noms(iq).eq."hplus") then
373	igcm_hplus=iq
374	mmol(igcm_hplus)=1.
375	count=count+1
376	idtracers(count) = "hplus"
377	endif
378	if (noms(iq).eq."hco2plus") then
379	igcm_hco2plus=iq
380	mmol(igcm_hco2plus)=45.
381	count=count+1
382	idtracers(count) = "hco2plus"
383	endif
384	if (noms(iq).eq."hcoplus") then
385	igcm_hcoplus=iq
386	mmol(igcm_hcoplus)=29.
387	count=count+1
388	idtracers(count) = "hcoplus"
389	endif
390	if (noms(iq).eq."h2oplus") then
391	igcm_h2oplus=iq
392	mmol(igcm_h2oplus)=18.
393	count=count+1
394	idtracers(count) = "h2oplus"
395	endif
396	if (noms(iq).eq."h3oplus") then
397	igcm_h3oplus=iq
398	mmol(igcm_h3oplus)=19.
399	count=count+1
400	idtracers(count) = "h3oplus"
401	endif
402	if (noms(iq).eq."ohplus") then
403	igcm_ohplus=iq
404	mmol(igcm_ohplus)=17.
405	count=count+1
406	idtracers(count) = "ohplus"
407	endif
408	if (noms(iq).eq."elec") then
409	igcm_elec=iq
410	mmol(igcm_elec)=1./1822.89
411	count=count+1
412	idtracers(count) = "elec"
413	endif
414	if (noms(iq).eq."h2o_vap") then
415	igcm_h2o_vap=iq
416	mmol(igcm_h2o_vap)=18.
417	count=count+1
418	idtracers(count) = "h2o_vap"
419	endif
420	if (noms(iq).eq."hdo_vap") then
421	igcm_hdo_vap=iq
422	mmol(igcm_hdo_vap)=19.
423	count=count+1
424	idtracers(count) = "hdo_vap"
425	endif
426	if (noms(iq).eq."od") then
427	igcm_od=iq
428	mmol(igcm_od)=18.
429	count=count+1
430	idtracers(count) = "od"
431	endif
432	if (noms(iq).eq."d") then
433	igcm_d=iq
434	mmol(igcm_d)=2.
435	count=count+1
436	idtracers(count) = "d"
437	endif
438	if (noms(iq).eq."hd") then
439	igcm_hd=iq
440	mmol(igcm_hd)=3.
441	count=count+1
442	idtracers(count) = "hd"
443	endif
444	if (noms(iq).eq."do2") then
445	igcm_do2=iq
446	mmol(igcm_do2)=34.
447	count=count+1
448	idtracers(count) = "do2"
449	endif
450	if (noms(iq).eq."hdo2") then
451	igcm_hdo2=iq
452	mmol(igcm_hdo2)=35.
453	count=count+1
454	idtracers(count) = "hdo2"
455	endif
456	if (noms(iq).eq."co2_ice") then
457	igcm_co2_ice=iq
458	mmol(igcm_co2_ice)=44.
459	count=count+1
460	idtracers(count) = "co2_ice"
461	endif
462	if (noms(iq).eq."h2o_ice") then
463	igcm_h2o_ice=iq
464	mmol(igcm_h2o_ice)=18.
465	count=count+1
466	idtracers(count) = "h2o_ice"
467	endif
468	if (noms(iq).eq."hdo_ice") then
469	igcm_hdo_ice=iq
470	mmol(igcm_hdo_ice)=19.
471	count=count+1
472	idtracers(count) = "hdo_ice"
473	endif
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
479	idtracers(count) = "Ar_N2"
480	endif
481	if (co2clouds) then
482	if (noms(iq).eq."ccnco2_mass") then
483	igcm_ccnco2_mass=iq
484	count=count+1
485	idtracers(count) = "ccnco2_mass"
486	endif
487	if (noms(iq).eq."ccnco2_number") then
488	igcm_ccnco2_number=iq
489	count=count+1
490	idtracers(count) = "ccnco2_number"
491	endif
492	if (meteo_flux) then
493	if (noms(iq).eq."ccnco2_meteor_mass") then
494	igcm_ccnco2_meteor_mass=iq
495	count=count+1
496	idtracers(count) = "ccnco2_meteor_mass"
497	endif
498	if (noms(iq).eq."ccnco2_meteor_number") then
499	igcm_ccnco2_meteor_number=iq
500	count=count+1
501	idtracers(count) = "ccnco2_meteor_number"
502	endif
503	end if
504	if (co2useh2o) then
505	if (noms(iq).eq."ccnco2_h2o_number") then
506	igcm_ccnco2_h2o_number=iq
507	count=count+1
508	idtracers(count) = "ccnco2_h2o_number"
509	endif
510	if (noms(iq).eq."ccnco2_h2o_mass_ice") then
511	igcm_ccnco2_h2o_mass_ice=iq
512	count=count+1
513	idtracers(count) = "ccnco2_h2o_mass_ice"
514	endif
515	if (noms(iq).eq."ccnco2_h2o_mass_ccn") then
516	igcm_ccnco2_h2o_mass_ccn=iq
517	count=count+1
518	idtracers(count) = "ccnco2_h2o_mass_ccn"
519	endif
520	end if
521	endif
522	enddo ! of do iq=1,nq
523
524	! check that we identified all tracers:
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
532	call abort_physic("initracer","tracer mismatch",1)
533	else
534	write(,) "initracer: found all expected tracers, namely:"
535	do iq = 1,nq
536	write(,) ' ',iq,' ',trim(noms(iq))
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
551	qsurf(1:ngrid,igcm_h2o_vap)=0
552	endif
553	endif
554
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
559	call abort_physic("initracer",
560	& "Tracer H2O must be initialized before HDO",1)
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"
565	call abort_physic("initracer","adapt your tracer.def",1)
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"
570	call abort_physic("initracer","adapt your tracer.def",1)
571	endif
572	endif
573
574	c------------------------------------------------------------
575	c Initialize tracers .... (in tracer_mod)
576	c------------------------------------------------------------
577	! start by setting everything to (default) zero
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
582
583
584	! some reference values
585	rho_dust=2500. ! Mars dust density (kg.m-3)
586	rho_ice=920. ! Water ice density (kg.m-3)
587	rho_ice_co2=1650.
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
595	if( (nq.lt.2).or.(water.and.(nq.lt.4))
596	* .or.(hdo.and.(nq.lt.6) )) then
597	write(,)'initracer: nq is too low : nq=', nq
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
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
629	alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
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
635	#endif
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)
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 ----------------------------------------------------------------------
674	!c mountain top dust flows scheme
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
677	if (topflows) then
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
682	end if !(topflows)
683	!c ----------------------------------------------------------------------
684
685	else
686
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
692	if (dustbin.gt.1) then
693	print*,'initracer: STOP!',
694	$ ' properties of dust need to be set in initracer !!!'
695	call abort_physic("initracer","dustbin properties issue",1)
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
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
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.
750	if(water.and.(nq.ge.2)) then
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.
755	elseif(water.and.(nq.lt.2)) then
756	write(,) 'nq is too low : nq=', nq
757	write(,) 'water= ',water
758	endif
759
760	end if ! (water)
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
780
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 !"
805	call abort_physic("initracer","doubleq issue",1)
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 !"
811	call abort_physic("initracer","doubleq issue",1)
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 !!"
825	write(,) " dustbin is set to ",dustbin,
826	& " but we only have the following dust tracers:"
827	do iq=1,count
828	write(,)" ",trim(noms(igcm_dustbin(iq)))
829	enddo
830	call abort_physic("initracer","dustbin issue",1)
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 !"
840	call abort_physic("initracer","water cycle issue",1)
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 !"
846	call abort_physic("initracer","water cycle issue",1)
847	endif
848	endif
849
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 !"
856	call abort_physic("initracer","hdo cycle issue",1)
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 !"
862	call abort_physic("initracer","hdo cycle issue",1)
863	endif
864	endif
865
866
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 !"
873	call abort_physic("initracer","co2 clouds issue",1)
874	end if
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 !"
879	call abort_physic("initracer","co2 clouds issue",1)
880	end if
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
933	endif
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 !"
941	call abort_physic("initracer","rdstorm issue",1)
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 !"
947	call abort_physic("initracer","rdstorm issue",1)
948	endif
949	endif
950
951	if (topflows) then
952	! verify that we indeed have topdust_mass and topdust_number tracers
953	if (igcm_topdust_mass.eq.0) then
954	write(,) "initracer: error !!"
955	write(,) " cannot use topflows option without ",
956	& "a topdust_mass tracer !"
957	call abort_physic("initracer","topflows issue",1)
958	endif
959	if (igcm_topdust_number.eq.0) then
960	write(,) "initracer: error !!"
961	write(,) " cannot use topflows option without ",
962	& "a topdust_number tracer !"
963	call abort_physic("initracer","topflows issue",1)
964	endif
965	endif
966
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!"
973	call abort_physic("initracer","missing co2 tracer",1)
974	endif
975	if (igcm_co.eq.0) then
976	write(,) "initracer: error !!"
977	write(,) " with nltemodel>0, we need the co tracer!"
978	call abort_physic("initracer","missing co tracer",1)
979	endif
980	if (igcm_o.eq.0) then
981	write(,) "initracer: error !!"
982	write(,) " with nltemodel>0, we need the o tracer!"
983	call abort_physic("initracer","missing o tracer",1)
984	endif
985	if (igcm_n2.eq.0) then
986	write(,) "initracer: error !!"
987	write(,) " with nltemodel>0, we need the n2 tracer!"
988	call abort_physic("initracer","missing n2 tracer",1)
989	endif
990	endif
991	endif
992
993	if (scavenging) then
994	! verify that we indeed have ccn_mass and ccn_number tracers
995	if (igcm_ccn_mass.eq.0 .and. igcm_ccnco2_mass.eq.0) then
996	write(,) "initracer: error !!"
997	write(,) " cannot use scavenging option without ",
998	& "a ccn_mass or ccnco2_mass tracer !"
999	call abort_physic("initracer","scavenging issue",1)
1000	endif
1001	if (igcm_ccn_number.eq.0 .and. igcm_ccnco2_number.eq.0 ) then
1002	write(,) "initracer: error !!"
1003	write(,) " cannot use scavenging option without ",
1004	& "a ccn_number or ccnco2_number tracer !"
1005	call abort_physic("initracer","scavenging issue",1)
1006	endif
1007	endif ! of if (scavenging)
1008
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 !"
1015	call abort_physic("initracer","missing co2 tracer",1)
1016	endif
1017	if (igcm_co.eq.0) then
1018	write(,) "initracer: error !!"
1019	write(,) " cannot use chemistry option without ",
1020	& "a co tracer !"
1021	call abort_physic("initracer","missing co tracer",1)
1022	endif
1023	if (igcm_o.eq.0) then
1024	write(,) "initracer: error !!"
1025	write(,) " cannot use chemistry option without ",
1026	& "a o tracer !"
1027	call abort_physic("initracer","missing o tracer",1)
1028	endif
1029	if (igcm_o1d.eq.0) then
1030	write(,) "initracer: error !!"
1031	write(,) " cannot use chemistry option without ",
1032	& "a o1d tracer !"
1033	call abort_physic("initracer","missing o1d tracer",1)
1034	endif
1035	if (igcm_o2.eq.0) then
1036	write(,) "initracer: error !!"
1037	write(,) " cannot use chemistry option without ",
1038	& "an o2 tracer !"
1039	call abort_physic("initracer","missing o2 tracer",1)
1040	endif
1041	if (igcm_o3.eq.0) then
1042	write(,) "initracer: error !!"
1043	write(,) " cannot use chemistry option without ",
1044	& "an o3 tracer !"
1045	call abort_physic("initracer","missing o3 tracer",1)
1046	endif
1047	if (igcm_h.eq.0) then
1048	write(,) "initracer: error !!"
1049	write(,) " cannot use chemistry option without ",
1050	& "a h tracer !"
1051	call abort_physic("initracer","missing h tracer",1)
1052	endif
1053	if (igcm_h2.eq.0) then
1054	write(,) "initracer: error !!"
1055	write(,) " cannot use chemistry option without ",
1056	& "a h2 tracer !"
1057	call abort_physic("initracer","missing h2 tracer",1)
1058	endif
1059	if (igcm_oh.eq.0) then
1060	write(,) "initracer: error !!"
1061	write(,) " cannot use chemistry option without ",
1062	& "an oh tracer !"
1063	call abort_physic("initracer","missing oh tracer",1)
1064	endif
1065	if (igcm_ho2.eq.0) then
1066	write(,) "initracer: error !!"
1067	write(,) " cannot use chemistry option without ",
1068	& "a ho2 tracer !"
1069	call abort_physic("initracer","missing ho2 tracer",1)
1070	endif
1071	if (igcm_h2o2.eq.0) then
1072	write(,) "initracer: error !!"
1073	write(,) " cannot use chemistry option without ",
1074	& "a h2o2 tracer !"
1075	call abort_physic("initracer","missing h2o2 tracer",1)
1076	endif
1077	if (igcm_n2.eq.0) then
1078	write(,) "initracer: error !!"
1079	write(,) " cannot use chemistry option without ",
1080	& "a n2 tracer !"
1081	call abort_physic("initracer","missing n2 tracer",1)
1082	endif
1083	if (igcm_ar.eq.0) then
1084	write(,) "initracer: error !!"
1085	write(,) " cannot use chemistry option without ",
1086	& "an ar tracer !"
1087	call abort_physic("initracer","missing ar tracer",1)
1088	endif
1089	endif ! of if (photochem .or. callthermos)
1090
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
1112	end

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

Download in other formats:

Original Format