Context Navigation

initracer.F @ 1009

Last change on this file since 1009 was 1005, checked in by emillour, 11 years ago

Mars GCM:

Bug fix: when running with photochemistry, ccns did not sediment! Fixed initracer.F. Also added that callsedim/newsedim use updated temperatures.
Converted run0 and run_mcd scripts to bash.

File size: 21.0 KB

Line
1	SUBROUTINE initracer(qsurf,co2ice)
2
3	IMPLICIT NONE
4	c=======================================================================
5	c subject:
6	c --------
7	c Initialization related to tracer
8	c (transported dust, water, chemical species, ice...)
9	c
10	c Name of the tracer
11	c
12	c Test of dimension :
13	c Initialize COMMON tracer in tracer.h, using tracer names provided
14	c by the dynamics in "advtrac.h"
15	c
16	c
17	c author: F.Forget
18	c ------
19	c Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003)
20	c Ehouarn Millour (oct. 2008) identify tracers by their names
21	c=======================================================================
22
23
24	#include "dimensions.h"
25	#include "dimphys.h"
26	#include "comcstfi.h"
27	#include "callkeys.h"
28	#include "tracer.h"
29	#include "advtrac.h"
30	#include "comgeomfi.h"
31
32	#include "surfdat.h"
33
34	real qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
35	real co2ice(ngridmx) ! co2 ice mass on surface (e.g. kg.m-2)
36	integer iq,ig,count
37	real r0_lift , reff_lift, nueff_lift
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.
42	character(len=20) :: txt ! to store some text
43
44	c-----------------------------------------------------------------------
45	c radius(nqmx) ! aerosol particle radius (m)
46	c rho_q(nqmx) ! tracer densities (kg.m-3)
47	c alpha_lift(nqmx) ! saltation vertical flux/horiz flux ratio (m-1)
48	c alpha_devil(nqmx) ! lifting coeeficient by dust devil
49	c rho_dust ! Mars dust density
50	c rho_ice ! Water ice density
51	c nuice_ref ! Effective variance nueff of the
52	c ! water-ice size distributions
53	c doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio
54	c varian ! Characteristic variance of log-normal distribution
55	c-----------------------------------------------------------------------
56
57	! Initialization: get tracer names from the dynamics and check if we are
58	! using 'old' tracer convention ('q01',q02',...)
59	! or new convention (full tracer names)
60	! check if tracers have 'old' names
61	count=0
62	do iq=1,nqmx
63	txt=" "
64	write(txt,'(a1,i2.2)') 'q',iq
65	if (txt.eq.tnom(iq)) then
66	count=count+1
67	endif
68	enddo ! of do iq=1,nqmx
69
70	if (count.eq.nqmx) then
71	write(,) "initracer: tracers seem to follow old naming ",
72	& "convention (q01,q02,...)"
73	write(,) "you should run newstart to rename them"
74	stop
75	endif
76
77	! copy tracer names from dynamics
78	do iq=1,nqmx
79	noms(iq)=tnom(iq)
80	enddo
81
82	c------------------------------------------------------------
83	c NAME and molar mass of the tracer
84	c------------------------------------------------------------
85
86	! Identify tracers by their names: (and set corresponding values of mmol)
87	! 0. initialize tracer indexes to zero:
88	do iq=1,nqmx
89	igcm_dustbin(iq)=0
90	enddo
91	igcm_dust_mass=0
92	igcm_dust_number=0
93	igcm_ccn_mass=0
94	igcm_ccn_number=0
95	igcm_dust_submicron=0
96	igcm_h2o_vap=0
97	igcm_h2o_ice=0
98	igcm_co2=0
99	igcm_co=0
100	igcm_o=0
101	igcm_o1d=0
102	igcm_o2=0
103	igcm_o3=0
104	igcm_h=0
105	igcm_h2=0
106	igcm_oh=0
107	igcm_ho2=0
108	igcm_h2o2=0
109	igcm_ch4=0
110	igcm_n2=0
111	igcm_ar=0
112	igcm_ar_n2=0
113	igcm_n=0
114	igcm_no=0
115	igcm_no2=0
116	igcm_n2d=0
117	igcm_co2plus=0
118	igcm_oplus=0
119	igcm_o2plus=0
120	igcm_coplus=0
121	igcm_cplus=0
122	igcm_nplus=0
123	igcm_noplus=0
124	igcm_n2plus=0
125	igcm_hplus=0
126	igcm_hco2plus=0
127	igcm_elec=0
128
129	! 1. find dust tracers
130	count=0
131	if (dustbin.gt.0) then
132	do iq=1,nqmx
133	txt=" "
134	write(txt,'(a4,i2.2)')'dust',count+1
135	if (noms(iq).eq.txt) then
136	count=count+1
137	igcm_dustbin(count)=iq
138	mmol(iq)=100.
139	endif
140	enddo !do iq=1,nqmx
141	endif ! of if (dustbin.gt.0)
142	if (doubleq) then
143	do iq=1,nqmx
144	if (noms(iq).eq."dust_mass") then
145	igcm_dust_mass=iq
146	count=count+1
147	endif
148	if (noms(iq).eq."dust_number") then
149	igcm_dust_number=iq
150	count=count+1
151	endif
152	enddo
153	endif ! of if (doubleq)
154	if (microphys) then
155	do iq=1,nqmx
156	if (noms(iq).eq."ccn_mass") then
157	igcm_ccn_mass=iq
158	count=count+1
159	endif
160	if (noms(iq).eq."ccn_number") then
161	igcm_ccn_number=iq
162	count=count+1
163	endif
164	enddo
165	endif ! of if (microphys)
166	if (submicron) then
167	do iq=1,nqmx
168	if (noms(iq).eq."dust_submicron") then
169	igcm_dust_submicron=iq
170	mmol(iq)=100.
171	count=count+1
172	endif
173	enddo
174	endif ! of if (submicron)
175	! 2. find chemistry and water tracers
176	do iq=1,nqmx
177	if (noms(iq).eq."co2") then
178	igcm_co2=iq
179	mmol(igcm_co2)=44.
180	count=count+1
181	endif
182	if (noms(iq).eq."co") then
183	igcm_co=iq
184	mmol(igcm_co)=28.
185	count=count+1
186	endif
187	if (noms(iq).eq."o") then
188	igcm_o=iq
189	mmol(igcm_o)=16.
190	count=count+1
191	endif
192	if (noms(iq).eq."o1d") then
193	igcm_o1d=iq
194	mmol(igcm_o1d)=16.
195	count=count+1
196	endif
197	if (noms(iq).eq."o2") then
198	igcm_o2=iq
199	mmol(igcm_o2)=32.
200	count=count+1
201	endif
202	if (noms(iq).eq."o3") then
203	igcm_o3=iq
204	mmol(igcm_o3)=48.
205	count=count+1
206	endif
207	if (noms(iq).eq."h") then
208	igcm_h=iq
209	mmol(igcm_h)=1.
210	count=count+1
211	endif
212	if (noms(iq).eq."h2") then
213	igcm_h2=iq
214	mmol(igcm_h2)=2.
215	count=count+1
216	endif
217	if (noms(iq).eq."oh") then
218	igcm_oh=iq
219	mmol(igcm_oh)=17.
220	count=count+1
221	endif
222	if (noms(iq).eq."ho2") then
223	igcm_ho2=iq
224	mmol(igcm_ho2)=33.
225	count=count+1
226	endif
227	if (noms(iq).eq."h2o2") then
228	igcm_h2o2=iq
229	mmol(igcm_h2o2)=34.
230	count=count+1
231	endif
232	if (noms(iq).eq."n2") then
233	igcm_n2=iq
234	mmol(igcm_n2)=28.
235	count=count+1
236	endif
237	if (noms(iq).eq."ch4") then
238	igcm_ch4=iq
239	mmol(igcm_ch4)=16.
240	count=count+1
241	endif
242	if (noms(iq).eq."ar") then
243	igcm_ar=iq
244	mmol(igcm_ar)=40.
245	count=count+1
246	endif
247	if (noms(iq).eq."n") then
248	igcm_n=iq
249	mmol(igcm_n)=14.
250	count=count+1
251	endif
252	if (noms(iq).eq."no") then
253	igcm_no=iq
254	mmol(igcm_no)=30.
255	count=count+1
256	endif
257	if (noms(iq).eq."no2") then
258	igcm_no2=iq
259	mmol(igcm_no2)=46.
260	count=count+1
261	endif
262	if (noms(iq).eq."n2d") then
263	igcm_n2d=iq
264	mmol(igcm_n2d)=28.
265	count=count+1
266	endif
267	if (noms(iq).eq."co2plus") then
268	igcm_co2plus=iq
269	mmol(igcm_co2plus)=44.
270	count=count+1
271	endif
272	if (noms(iq).eq."oplus") then
273	igcm_oplus=iq
274	mmol(igcm_oplus)=16.
275	count=count+1
276	endif
277	if (noms(iq).eq."o2plus") then
278	igcm_o2plus=iq
279	mmol(igcm_o2plus)=32.
280	count=count+1
281	endif
282	if (noms(iq).eq."coplus") then
283	igcm_coplus=iq
284	mmol(igcm_coplus)=28.
285	count=count+1
286	endif
287	if (noms(iq).eq."cplus") then
288	igcm_cplus=iq
289	mmol(igcm_cplus)=12.
290	count=count+1
291	endif
292	if (noms(iq).eq."nplus") then
293	igcm_nplus=iq
294	mmol(igcm_nplus)=14.
295	count=count+1
296	endif
297	if (noms(iq).eq."noplus") then
298	igcm_noplus=iq
299	mmol(igcm_noplus)=30.
300	count=count+1
301	endif
302	if (noms(iq).eq."n2plus") then
303	igcm_n2plus=iq
304	mmol(igcm_n2plus)=28.
305	count=count+1
306	endif
307	if (noms(iq).eq."hplus") then
308	igcm_hplus=iq
309	mmol(igcm_hplus)=1.
310	count=count+1
311	endif
312	if (noms(iq).eq."hco2plus") then
313	igcm_hco2plus=iq
314	mmol(igcm_hco2plus)=45.
315	count=count+1
316	endif
317	if (noms(iq).eq."elec") then
318	igcm_elec=iq
319	mmol(igcm_elec)=1./1822.89
320	count=count+1
321	endif
322	if (noms(iq).eq."h2o_vap") then
323	igcm_h2o_vap=iq
324	mmol(igcm_h2o_vap)=18.
325	count=count+1
326	endif
327	if (noms(iq).eq."h2o_ice") then
328	igcm_h2o_ice=iq
329	mmol(igcm_h2o_ice)=18.
330	count=count+1
331	endif
332	! Other stuff: e.g. for simulations using co2 + neutral gaz
333	if (noms(iq).eq."Ar_N2") then
334	igcm_ar_n2=iq
335	mmol(igcm_ar_n2)=30.
336	count=count+1
337	endif
338
339	enddo ! of do iq=1,nqmx
340
341	! check that we identified all tracers:
342	if (count.ne.nqmx) then
343	write(,) "initracer: found only ",count," tracers"
344	write(,) " expected ",nqmx
345	do iq=1,count
346	write(,)' ',iq,' ',trim(noms(iq))
347	enddo
348	stop
349	else
350	write(,) "initracer: found all expected tracers, namely:"
351	do iq=1,nqmx
352	write(,)' ',iq,' ',trim(noms(iq))
353	enddo
354	endif
355
356	! if water cycle but iceparty=.false., there will nevertheless be
357	! water ice at the surface (iceparty is not used anymore, but this
358	! part is still relevant, as we want to stay compatible with the
359	! older versions).
360	if (water.and.(igcm_h2o_ice.eq.0)) then
361	igcm_h2o_ice=igcm_h2o_vap ! so that qsurf(i_h2o_ice) is identified
362	! even though there is no q(i_h2o_ice)
363	else
364	! surface ice qsurf(i_h2o_ice) was loaded twice by phyetat0,
365	! as qsurf(i_h2o_vap) & as qsurf(i_h2o_ice), so to be clean:
366	if (igcm_h2o_vap.ne.0) then
367	qsurf(1:ngridmx,igcm_h2o_vap)=0
368	endif
369	endif
370
371	c------------------------------------------------------------
372	c Initialize tracers .... (in tracer.h)
373	c------------------------------------------------------------
374	! start by setting everything to (default) zero
375	rho_q(1:nqmx)=0 ! tracer density (kg.m-3)
376	radius(1:nqmx)=0. ! tracer particle radius (m)
377	alpha_lift(1:nqmx) =0. ! tracer saltation vertical flux/horiz flux ratio (m-1)
378	alpha_devil(1:nqmx)=0. ! tracer lifting coefficient by dust devils
379
380
381	! some reference values
382	rho_dust=2500. ! Mars dust density (kg.m-3)
383	rho_ice=920. ! Water ice density (kg.m-3)
384	nuice_ref=0.1 ! Effective variance nueff of the
385	! water-ice size distribution
386	!!!nuice_sed=0.45 ! Sedimentation effective variance
387	! of the water-ice size distribution
388
389	if (doubleq) then
390	c "doubleq" technique
391	c -------------------
392	c (transport of mass and number mixing ratio)
393	c iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
394
395	if( (nqmx.lt.2).or.(water.and.(nqmx.lt.4)) ) then
396	write(,)'initracer: nqmx is too low : nqmx=', nqmx
397	write(,)'water= ',water,' doubleq= ',doubleq
398	end if
399
400	nueff_lift = 0.5
401	varian=sqrt(log(1.+nueff_lift))
402
403	rho_q(igcm_dust_mass)=rho_dust
404	rho_q(igcm_dust_number)=rho_dust
405
406	c Intermediate calcul for computing geometric mean radius r0
407	c as a function of mass and number mixing ratio Q and N
408	c (r0 = (r3n_q * Q/ N)^(1/3))
409	r3n_q = exp(-4.5varian2)(3./4.)/(pi*rho_dust)
410
411	c Intermediate calcul for computing effective radius reff
412	c from geometric mean radius r0
413	c (reff = ref_r0 * r0)
414	ref_r0 = exp(2.5varian*2)
415
416	c lifted dust :
417	c '''''''''''
418	reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m)
419	alpha_devil(igcm_dust_mass)=9.e-9 ! dust devil lift mass coeff
420	c alpha_lift(igcm_dust_mass)=3.0e-15 ! Lifted mass coeff
421	alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
422
423	r0_lift = reff_lift/ref_r0
424	alpha_devil(igcm_dust_number)=r3n_q*
425	& alpha_devil(igcm_dust_mass)/r0_lift**3
426	alpha_lift(igcm_dust_number)=r3n_q*
427	& alpha_lift(igcm_dust_mass)/r0_lift**3
428
429	radius(igcm_dust_mass) = reff_lift
430	radius(igcm_dust_number) = reff_lift
431
432	write(,) "initracer: doubleq_param reff_lift:", reff_lift
433	write(,) "initracer: doubleq_param nueff_lift:", nueff_lift
434	write(,) "initracer: doubleq_param alpha_lift:",
435	& alpha_lift(igcm_dust_mass)
436	else
437
438	! initialize varian, which may be used (e.g. by surfacearea)
439	! even with conrath dust
440	nueff_lift = 0.5
441	varian=sqrt(log(1.+nueff_lift))
442
443	if (dustbin.gt.1) then
444	print*,'initracer: STOP!',
445	$ ' properties of dust need to be set in initracer !!!'
446	stop
447
448	else if (dustbin.eq.1) then
449
450	c This will be used for 1 dust particle size:
451	c ------------------------------------------
452	radius(igcm_dustbin(1))=3.e-6
453	alpha_lift(igcm_dustbin(1))=0.0e-6
454	alpha_devil(igcm_dustbin(1))=7.65e-9
455	rho_q(igcm_dustbin(1))=rho_dust
456
457	endif
458	end if ! (doubleq)
459
460
461	c Scavenging of dust particles by H2O clouds:
462	c ------------------------------------------
463	c Initialize the two tracers used for the CCNs
464	if (water.AND.doubleq.AND.scavenging) then
465	radius(igcm_ccn_mass) = radius(igcm_dust_mass)
466	alpha_lift(igcm_ccn_mass) = 1e-30
467	alpha_devil(igcm_ccn_mass) = 1e-30
468	rho_q(igcm_ccn_mass) = rho_dust
469
470	radius(igcm_ccn_number) = radius(igcm_ccn_mass)
471	alpha_lift(igcm_ccn_number) = alpha_lift(igcm_ccn_mass)
472	alpha_devil(igcm_ccn_number) = alpha_devil(igcm_ccn_mass)
473	rho_q(igcm_ccn_number) = rho_q(igcm_ccn_mass)
474	endif ! of if (water.AND.doubleq.AND.scavenging)
475
476	c Submicron dust mode:
477	c --------------------
478
479	if (submicron) then
480	radius(igcm_dust_submicron)=0.1e-6
481	rho_q(igcm_dust_submicron)=rho_dust
482	if (doubleq) then
483	c If doubleq is also active, we use the population ratio:
484	alpha_lift(igcm_dust_submicron) =
485	& alpha_lift(igcm_dust_number)popratio
486	& rho_q(igcm_dust_submicron)4./3.pi*
487	& radius(igcm_dust_submicron)**3.
488	alpha_devil(igcm_dust_submicron)=1.e-30
489	else
490	alpha_lift(igcm_dust_submicron)=1e-6
491	alpha_devil(igcm_dust_submicron)=1.e-30
492	endif ! (doubleq)
493	end if ! (submicron)
494
495	c Initialization for water vapor
496	c ------------------------------
497	if(water) then
498	radius(igcm_h2o_vap)=0.
499	alpha_lift(igcm_h2o_vap) =0.
500	alpha_devil(igcm_h2o_vap)=0.
501	if(water.and.(nqmx.ge.2)) then
502	radius(igcm_h2o_ice)=3.e-6
503	rho_q(igcm_h2o_ice)=rho_ice
504	alpha_lift(igcm_h2o_ice) =0.
505	alpha_devil(igcm_h2o_ice)=0.
506	elseif(water.and.(nqmx.lt.2)) then
507	write(,) 'nqmx is too low : nqmx=', nqmx
508	write(,) 'water= ',water
509	endif
510
511	end if ! (water)
512
513	c Output for records:
514	c ~~~~~~~~~~~~~~~~~~
515	write(,)
516	Write(,) '******** initracer : dust transport parameters :'
517	write(,) 'alpha_lift = ', alpha_lift
518	write(,) 'alpha_devil = ', alpha_devil
519	write(,) 'radius = ', radius
520	if(doubleq) then
521	write(,) 'reff_lift (um) = ', reff_lift
522	write(,) 'size distribution variance = ', varian
523	write(,) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
524	end if
525
526	!
527	! some extra (possibly redundant) sanity checks for tracers:
528	! ---------------------------------------------------------
529
530	if (doubleq) then
531	! verify that we indeed have dust_mass and dust_number tracers
532	if (igcm_dust_mass.eq.0) then
533	write(,) "initracer: error !!"
534	write(,) " cannot use doubleq option without ",
535	& "a dust_mass tracer !"
536	stop
537	endif
538	if (igcm_dust_number.eq.0) then
539	write(,) "initracer: error !!"
540	write(,) " cannot use doubleq option without ",
541	& "a dust_number tracer !"
542	stop
543	endif
544	endif
545
546	if ((.not.doubleq).and.(dustbin.gt.0)) then
547	! verify that we indeed have 'dustbin' dust tracers
548	count=0
549	do iq=1,dustbin
550	if (igcm_dustbin(iq).ne.0) then
551	count=count+1
552	endif
553	enddo
554	if (count.ne.dustbin) then
555	write(,) "initracer: error !!"
556	write(,) " dusbin is set to ",dustbin,
557	& " but we only have the following dust tracers:"
558	do iq=1,count
559	write(,)" ",trim(noms(igcm_dustbin(iq)))
560	enddo
561	stop
562	endif
563	endif
564
565	if (water) then
566	! verify that we indeed have h2o_vap and h2o_ice tracers
567	if (igcm_h2o_vap.eq.0) then
568	write(,) "initracer: error !!"
569	write(,) " cannot use water option without ",
570	& "an h2o_vap tracer !"
571	stop
572	endif
573	if (igcm_h2o_ice.eq.0) then
574	write(,) "initracer: error !!"
575	write(,) " cannot use water option without ",
576	& "an h2o_ice tracer !"
577	stop
578	endif
579	endif
580
581	if (scavenging) then
582	! verify that we indeed have ccn_mass and ccn_number tracers
583	if (igcm_ccn_mass.eq.0) then
584	write(,) "initracer: error !!"
585	write(,) " cannot use scavenging option without ",
586	& "a ccn_mass tracer !"
587	stop
588	endif
589	if (igcm_ccn_number.eq.0) then
590	write(,) "initracer: error !!"
591	write(,) " cannot use scavenging option without ",
592	& "a ccn_number tracer !"
593	stop
594	endif
595	endif ! of if (scavenging)
596
597	if (photochem .or. callthermos) then
598	! verify that we indeed have the chemistry tracers
599	if (igcm_co2.eq.0) then
600	write(,) "initracer: error !!"
601	write(,) " cannot use chemistry option without ",
602	& "a co2 tracer !"
603	stop
604	endif
605	if (igcm_co.eq.0) then
606	write(,) "initracer: error !!"
607	write(,) " cannot use chemistry option without ",
608	& "a co tracer !"
609	stop
610	endif
611	if (igcm_o.eq.0) then
612	write(,) "initracer: error !!"
613	write(,) " cannot use chemistry option without ",
614	& "a o tracer !"
615	stop
616	endif
617	if (igcm_o1d.eq.0) then
618	write(,) "initracer: error !!"
619	write(,) " cannot use chemistry option without ",
620	& "a o1d tracer !"
621	stop
622	endif
623	if (igcm_o2.eq.0) then
624	write(,) "initracer: error !!"
625	write(,) " cannot use chemistry option without ",
626	& "an o2 tracer !"
627	stop
628	endif
629	if (igcm_o3.eq.0) then
630	write(,) "initracer: error !!"
631	write(,) " cannot use chemistry option without ",
632	& "an o3 tracer !"
633	stop
634	endif
635	if (igcm_h.eq.0) then
636	write(,) "initracer: error !!"
637	write(,) " cannot use chemistry option without ",
638	& "a h tracer !"
639	stop
640	endif
641	if (igcm_h2.eq.0) then
642	write(,) "initracer: error !!"
643	write(,) " cannot use chemistry option without ",
644	& "a h2 tracer !"
645	stop
646	endif
647	if (igcm_oh.eq.0) then
648	write(,) "initracer: error !!"
649	write(,) " cannot use chemistry option without ",
650	& "an oh tracer !"
651	stop
652	endif
653	if (igcm_ho2.eq.0) then
654	write(,) "initracer: error !!"
655	write(,) " cannot use chemistry option without ",
656	& "a ho2 tracer !"
657	stop
658	endif
659	if (igcm_h2o2.eq.0) then
660	write(,) "initracer: error !!"
661	write(,) " cannot use chemistry option without ",
662	& "a h2o2 tracer !"
663	stop
664	endif
665	if (igcm_n2.eq.0) then
666	write(,) "initracer: error !!"
667	write(,) " cannot use chemistry option without ",
668	& "a n2 tracer !"
669	stop
670	endif
671	if (igcm_ar.eq.0) then
672	write(,) "initracer: error !!"
673	write(,) " cannot use chemistry option without ",
674	& "an ar tracer !"
675	stop
676	endif
677	endif ! of if (photochem .or. callthermos)
678
679	end

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

Download in other formats:

Original Format