Context Navigation

initracer.F90 @ 2694

Last change on this file since 2694 was 2693, checked in by aslmd, 3 years ago
Added radius and rho option in modern traceur
File size: 20.5 KB

Line
1	SUBROUTINE initracer(ngrid,nq,nametrac)
2
3	use surfdat_h, ONLY: dryness, watercaptag
4	USE tracer_h
5	USE callkeys_mod, only: water
6	USE recombin_corrk_mod, ONLY: ini_recombin
7	USE mod_phys_lmdz_para, only: is_master, bcast
8	IMPLICIT NONE
9	!=======================================================================
10	! subject:
11	! --------
12	! Initialization related to tracer
13	! (transported dust, water, chemical species, ice...)
14	!
15	! Name of the tracer
16	!
17	! Test of dimension :
18	! Initialize COMMON tracer in tracer.h, using tracer names provided
19	! by the argument nametrac
20	!
21	! author: F.Forget
22	! ------
23	! Ehouarn Millour (oct. 2008) identify tracers by their names
24	! Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
25	!=======================================================================
26
27	integer,intent(in) :: ngrid,nq
28	character(len=30),intent(in) :: nametrac(nq) ! name of the tracer from dynamics
29
30	character(len=500) :: tracline ! to read traceur.def lines
31	! character(len=30) :: txt ! to store some text
32	integer iq,ig,count,ierr
33
34	!-----------------------------------------------------------------------
35	! radius(nq) ! aerosol particle radius (m)
36	! rho_q(nq) ! tracer densities (kg.m-3)
37	! qext(nq) ! Single Scat. Extinction coeff at 0.67 um
38	! alpha_lift(nq) ! saltation vertical flux/horiz flux ratio (m-1)
39	! alpha_devil(nq) ! lifting coeeficient by dust devil
40	! rho_dust ! Mars dust density
41	! rho_ice ! Water ice density
42	! doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio
43	! varian ! Characteristic variance of log-normal distribution
44	!-----------------------------------------------------------------------
45
46	if (is_master) then ! only the master proc/thread needs do this
47
48	moderntracdef=.false. ! For modern traceur.def (default false, old type)
49
50	open(407, form = 'formatted', status = 'old', &
51	file = 'traceur.def', iostat=ierr)
52	if (ierr /=0) then
53	! call abort_physic('initracer',&
54	! 'Problem in opening traceur.def',1)
55	print*,'Problem in opening traceur.def'
56	stop
57	end if
58	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
59	READ(407,'(A)') tracline
60	IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
61	READ(tracline,*) nqtot ! Try standard traceur.def
62	ELSE
63	moderntracdef = .true.
64	DO
65	READ(407,'(A)',iostat=ierr) tracline
66	IF (ierr==0) THEN
67	IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
68	READ(tracline,*) nqtot
69	! Temporary not implemented solution
70	if (nqtot/=nq) then
71	! call abort_physic('initracer','Different number of '// &
72	! 'tracers in dynamics and physics not managed yet',1)
73	print*,'!= nbr oftracers in dynamics and physics not managed yet'
74	stop
75	endif
76	EXIT
77	ENDIF
78	ELSE ! If pb, or if reached EOF without having found nqtot
79	! call abort_physic('initracer','Unable to read numbers '// &
80	! 'of tracers in traceur.def',1)
81	print*,"unable to read numbers of tracer in tracer.def"
82	stop
83	ENDIF
84	ENDDO
85	ENDIF ! if modern or standard traceur.def
86
87	endif ! of if (is_master)
88
89	! share the information with other procs/threads (if any)
90	CALL bcast(nqtot)
91	CALL bcast(moderntracdef)
92
93	!! -----------------------------------------------------------------------
94	!! For the moment number of tracers in dynamics and physics are equal
95	nqtot=nq
96	!! we allocate once for all arrays in common in tracer_h.F90
97	!! (supposedly those are not used before call to initracer)
98	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
99	IF (.NOT.ALLOCATED(mmol)) ALLOCATE(mmol(nq))
100	IF (.NOT.ALLOCATED(aki)) ALLOCATE(aki(nqtot))
101	IF (.NOT.ALLOCATED(cpi)) ALLOCATE(cpi(nqtot))
102	IF (.NOT.ALLOCATED(radius)) ALLOCATE(radius(nq))
103	IF (.NOT.ALLOCATED(rho_q)) ALLOCATE(rho_q(nq))
104	IF (.NOT.ALLOCATED(qext)) ALLOCATE(qext(nq))
105	IF (.NOT.ALLOCATED(alpha_lift)) ALLOCATE(alpha_lift(nq))
106	IF (.NOT.ALLOCATED(alpha_devil)) ALLOCATE(alpha_devil(nq))
107	IF (.NOT.ALLOCATED(qextrhor)) ALLOCATE(qextrhor(nq))
108	IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))
109	IF (.NOT.ALLOCATED(is_chim)) ALLOCATE(is_chim(nqtot))
110	IF (.NOT.ALLOCATED(is_rad)) ALLOCATE(is_rad(nqtot))
111	IF (.NOT.ALLOCATED(is_recomb)) ALLOCATE(is_recomb(nqtot))
112	IF (.NOT.ALLOCATED(is_recomb_qset)) THEN
113	ALLOCATE(is_recomb_qset(nqtot))
114	ENDIF
115	IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN
116	ALLOCATE(is_recomb_qotf(nqtot))
117	ENDIF
118	IF (.NOT. allocated(is_generic)) allocate(is_generic(nq)) !LT
119	!! initialization
120	alpha_lift(:) = 0.
121	alpha_devil(:) = 0.
122	mmol(:) = 0.
123	aki(:) = 0.
124	cpi(:) = 0.
125	is_chim(:) = 0
126	is_rad(:) = 0
127	is_recomb(:) = 0
128	is_recomb_qset(:) = 0
129	is_recomb_qotf(:) = 0
130
131	! Added by JVO 2017 : these arrays are handled later
132	! -> initialization is the least we can do, please !!!
133	radius(:)=0.
134	qext(:)=0.
135	is_generic(:)= 0 !LT
136
137
138	! Initialization: copy tracer names from dynamics
139	do iq=1,nq
140	noms(iq)=nametrac(iq)
141	write(,)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
142	enddo
143
144
145	! Identify tracers by their names: (and set corresponding values of mmol)
146	! 0. initialize tracer indexes to zero:
147	! NB: igcm_* indexes are commons in 'tracer.h'
148	do iq=1,nq
149	igcm_dustbin(iq)=0
150	enddo
151	igcm_dust_mass=0
152	igcm_dust_number=0
153	igcm_h2o_vap=0
154	igcm_h2o_ice=0
155	igcm_co2=0
156	igcm_co=0
157	igcm_o=0
158	igcm_o1d=0
159	igcm_o2=0
160	igcm_o3=0
161	igcm_h=0
162	igcm_h2=0
163	igcm_oh=0
164	igcm_ho2=0
165	igcm_h2o2=0
166	igcm_n2=0
167	igcm_n=0
168	igcm_n2d=0
169	igcm_no=0
170	igcm_no2=0
171	igcm_ar=0
172	igcm_ar_n2=0
173	igcm_co2_ice=0
174
175	igcm_ch4=0
176	igcm_ch3=0
177	igcm_ch=0
178	igcm_3ch2=0
179	igcm_1ch2=0
180	igcm_cho=0
181	igcm_ch2o=0
182	igcm_ch3o=0
183	igcm_c=0
184	igcm_c2=0
185	igcm_c2h=0
186	igcm_c2h2=0
187	igcm_c2h3=0
188	igcm_c2h4=0
189	igcm_c2h6=0
190	igcm_ch2co=0
191	igcm_ch3co=0
192	igcm_hcaer=0
193	! ! LT generic tracers
194	! igcm_generic_Fe_vap = 0
195	! igcm_generic_Fe_ice = 0
196	! igcm_generic_Cr_vap = 0
197	! igcm_generic_Cr_ice = 0
198
199
200	! 1. find dust tracers
201	count=0
202
203	! 2. find chemistry and water tracers
204	do iq=1,nq
205	if (noms(iq).eq."co2") then
206	igcm_co2=iq
207	mmol(igcm_co2)=44.
208	count=count+1
209	! write(,) 'co2: count=',count
210	endif
211	if (noms(iq).eq."co2_ice") then
212	igcm_co2_ice=iq
213	mmol(igcm_co2_ice)=44.
214	count=count+1
215	! write(,) 'co2_ice: count=',count
216	endif
217	if (noms(iq).eq."h2o_vap") then
218	igcm_h2o_vap=iq
219	mmol(igcm_h2o_vap)=18.
220	count=count+1
221	! write(,) 'h2o_vap: count=',count
222	endif
223	if (noms(iq).eq."h2o_ice") then
224	igcm_h2o_ice=iq
225	mmol(igcm_h2o_ice)=18.
226	count=count+1
227	! write(,) 'h2o_ice: count=',count
228	endif
229	if (noms(iq).eq."co") then
230	igcm_co=iq
231	mmol(igcm_co)=28.
232	count=count+1
233	endif
234	if (noms(iq).eq."o") then
235	igcm_o=iq
236	mmol(igcm_o)=16.
237	count=count+1
238	endif
239	if (noms(iq).eq."o1d") then
240	igcm_o1d=iq
241	mmol(igcm_o1d)=16.
242	count=count+1
243	endif
244	if (noms(iq).eq."o2") then
245	igcm_o2=iq
246	mmol(igcm_o2)=32.
247	count=count+1
248	endif
249	if (noms(iq).eq."o3") then
250	igcm_o3=iq
251	mmol(igcm_o3)=48.
252	count=count+1
253	endif
254	if (noms(iq).eq."h") then
255	igcm_h=iq
256	mmol(igcm_h)=1.
257	count=count+1
258	endif
259	if (noms(iq).eq."h2") then
260	igcm_h2=iq
261	mmol(igcm_h2)=2.
262	count=count+1
263	endif
264	if (noms(iq).eq."oh") then
265	igcm_oh=iq
266	mmol(igcm_oh)=17.
267	count=count+1
268	endif
269	if (noms(iq).eq."ho2") then
270	igcm_ho2=iq
271	mmol(igcm_ho2)=33.
272	count=count+1
273	endif
274	if (noms(iq).eq."h2o2") then
275	igcm_h2o2=iq
276	mmol(igcm_h2o2)=34.
277	count=count+1
278	endif
279	if (noms(iq).eq."n2") then
280	igcm_n2=iq
281	mmol(igcm_n2)=28.
282	count=count+1
283	endif
284	if (noms(iq).eq."ch4") then
285	igcm_ch4=iq
286	mmol(igcm_ch4)=16.
287	count=count+1
288	endif
289	if (noms(iq).eq."ar") then
290	igcm_ar=iq
291	mmol(igcm_ar)=40.
292	count=count+1
293	endif
294	if (noms(iq).eq."n") then
295	igcm_n=iq
296	mmol(igcm_n)=14.
297	count=count+1
298	endif
299	if (noms(iq).eq."no") then
300	igcm_no=iq
301	mmol(igcm_no)=30.
302	count=count+1
303	endif
304	if (noms(iq).eq."no2") then
305	igcm_no2=iq
306	mmol(igcm_no2)=46.
307	count=count+1
308	endif
309	if (noms(iq).eq."n2d") then
310	igcm_n2d=iq
311	mmol(igcm_n2d)=28.
312	count=count+1
313	endif
314	if (noms(iq).eq."ch3") then
315	igcm_ch3=iq
316	mmol(igcm_ch3)=15.
317	count=count+1
318	endif
319	if (noms(iq).eq."ch") then
320	igcm_ch=iq
321	mmol(igcm_ch)=13.
322	count=count+1
323	endif
324	if (noms(iq).eq."3ch2") then
325	igcm_3ch2=iq
326	mmol(igcm_3ch2)=14.
327	count=count+1
328	endif
329	if (noms(iq).eq."1ch2") then
330	igcm_1ch2=iq
331	mmol(igcm_1ch2)=14.
332	count=count+1
333	endif
334	if (noms(iq).eq."cho") then
335	igcm_cho=iq
336	mmol(igcm_cho)=29.
337	count=count+1
338	endif
339	if (noms(iq).eq."ch2o") then
340	igcm_ch2o=iq
341	mmol(igcm_ch2o)=30.
342	count=count+1
343	endif
344	if (noms(iq).eq."ch3o") then
345	igcm_ch3o=iq
346	mmol(igcm_ch3o)=31.
347	count=count+1
348	endif
349	if (noms(iq).eq."c") then
350	igcm_c=iq
351	mmol(igcm_c)=12.
352	count=count+1
353	endif
354	if (noms(iq).eq."c2") then
355	igcm_c2=iq
356	mmol(igcm_c2)=24.
357	count=count+1
358	endif
359	if (noms(iq).eq."c2h") then
360	igcm_c2h=iq
361	mmol(igcm_c2h)=25.
362	count=count+1
363	endif
364	if (noms(iq).eq."c2h2") then
365	igcm_c2h2=iq
366	mmol(igcm_c2h2)=26.
367	count=count+1
368	endif
369	if (noms(iq).eq."c2h3") then
370	igcm_c2h3=iq
371	mmol(igcm_c2h3)=27.
372	count=count+1
373	endif
374	if (noms(iq).eq."c2h4") then
375	igcm_c2h4=iq
376	mmol(igcm_c2h4)=28.
377	count=count+1
378	endif
379	if (noms(iq).eq."c2h6") then
380	igcm_c2h6=iq
381	mmol(igcm_c2h6)=30.
382	count=count+1
383	endif
384	if (noms(iq).eq."ch2co") then
385	igcm_ch2co=iq
386	mmol(igcm_ch2co)=42.
387	count=count+1
388	endif
389	if (noms(iq).eq."ch3co") then
390	igcm_ch3co=iq
391	mmol(igcm_ch3co)=43.
392	count=count+1
393	endif
394	if (noms(iq).eq."hcaer") then
395	igcm_hcaer=iq
396	mmol(igcm_hcaer)=50.
397	count=count+1
398	endif
399
400	enddo ! of do iq=1,nq
401
402	! check that we identified all tracers:
403	if (count.ne.nq) then
404	write(,) "initracer: found only ",count," tracers"
405	write(,) " expected ",nq
406	do iq=1,count
407	write(,)' ',iq,' ',trim(noms(iq))
408	enddo
409	! stop
410	else
411	write(,) "initracer: found all expected tracers, namely:"
412	do iq=1,nq
413	write(,)' ',iq,' ',trim(noms(iq))
414	enddo
415	endif
416
417	! Get data of tracers
418	do iq=1,nqtot
419	if (is_master) read(407,'(A)') tracline
420	! call bcast(tracline)
421	call get_tracdat(iq, tracline)
422	enddo
423
424	if (is_master) close(407)
425
426	! Calculate number of species in the chemistry
427	nesp = sum(is_chim)
428	write(,) 'Number of species in the chemistry nesp = ',nesp
429
430	! Calculate number of generic tracers
431	ngt = sum(is_generic)
432	write(,) 'Number of generic tracer is ngt = ',ngt
433	! Processing modern traceur options
434	if(moderntracdef) then
435	call ini_recombin
436	endif
437
438	!------------------------------------------------------------
439	! Initialisation tracers ....
440	!------------------------------------------------------------
441	rho_q(1:nq)=0
442
443	rho_dust=2500. ! Mars dust density (kg.m-3)
444	rho_ice=920. ! Water ice density (kg.m-3)
445	rho_co2=1620. ! CO2 ice density (kg.m-3)
446
447	! Initialization for water vapor
448	! ------------------------------
449	if(water) then
450	radius(igcm_h2o_vap)=0.
451	Qext(igcm_h2o_vap)=0.
452	alpha_lift(igcm_h2o_vap) =0.
453	alpha_devil(igcm_h2o_vap)=0.
454	qextrhor(igcm_h2o_vap)= 0.
455
456	!! this is defined in surfdat_h.F90
457	IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
458	IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
459
460	do ig=1,ngrid
461	if (ngrid.ne.1) watercaptag(ig)=.false.
462	dryness(ig) = 1.
463	enddo
464
465
466	radius(igcm_h2o_ice)=3.e-6
467	rho_q(igcm_h2o_ice)=rho_ice
468	Qext(igcm_h2o_ice)=0.
469	! alpha_lift(igcm_h2o_ice) =0.
470	! alpha_devil(igcm_h2o_ice)=0.
471	qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice) &
472	/ (rho_ice*radius(igcm_h2o_ice))
473
474
475	end if ! (water)
476
477
478	!
479	! some extra (possibly redundant) sanity checks for tracers:
480	! ---------------------------------------------------------
481	if (water) then
482	! verify that we indeed have h2o_vap and h2o_ice tracers
483	if (igcm_h2o_vap.eq.0) then
484	write(,) "initracer: error !!"
485	write(,) " cannot use water option without ",&
486	"an h2o_vap tracer !"
487	stop
488	endif
489	if (igcm_h2o_ice.eq.0) then
490	write(,) "initracer: error !!"
491	write(,) " cannot use water option without ",&
492	"an h2o_ice tracer !"
493	stop
494	endif
495	endif
496
497
498	! Output for records:
499	! ~~~~~~~~~~~~~~~~~~
500	write(,)
501	Write(,) '******** initracer : dust transport parameters :'
502	write(,) 'alpha_lift = ', alpha_lift
503	write(,) 'alpha_devil = ', alpha_devil
504	write(,) 'radius = ', radius
505	write(,) 'Qext = ', qext
506	write(,)
507
508	contains
509
510	subroutine get_tracdat(iq,tracline)
511	!-------------------ADDING NEW OPTIONS-------------------
512	! Duplicate if sentence for adding new options
513	! Do not forget to add the new variables in tracer_h.F90
514	! Do not forget to allocate and initialize the new variables above
515	! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
516	!-------------------------------------------------------
517	use tracer_h
518	implicit none
519	integer, intent(in) :: iq ! tracer index
520	character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
521
522	read(tracline,*) noms(iq)
523	! JVO 20 : We should add a sanity check aborting when duplicates in names !
524	write(,)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
525	! option mmol
526	if (index(tracline,'mmol=' ) /= 0) then
527	read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
528	mmol(iq)
529	write(,) ' Parameter value (traceur.def) : mmol=', &
530	mmol(iq)
531	else
532	write(,) ' Parameter value (default) : mmol=', &
533	mmol(iq)
534	end if
535	! option aki
536	if (index(tracline,'aki=' ) /= 0) then
537	read(tracline(index(tracline,'aki=')+len('aki='):),*) &
538	aki(iq)
539	write(,) ' Parameter value (traceur.def) : aki=', &
540	aki(iq)
541	else
542	write(,) ' Parameter value (default) : aki=', &
543	aki(iq)
544	end if
545	! option cpi
546	if (index(tracline,'cpi=' ) /= 0) then
547	read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
548	cpi(iq)
549	write(,) ' Parameter value (traceur.def) : cpi=', &
550	cpi(iq)
551	else
552	write(,) ' Parameter value (default) : cpi=', &
553	cpi(iq)
554	end if
555	! option is_chim
556	if (index(tracline,'is_chim=' ) /= 0) then
557	read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) &
558	is_chim(iq)
559	write(,) ' Parameter value (traceur.def) : is_chim=', &
560	is_chim(iq)
561	else
562	write(,) ' Parameter value (default) : is_chim=', &
563	is_chim(iq)
564	end if
565	! option is_rad
566	if (index(tracline,'is_rad=') /= 0) then
567	read(tracline(index(tracline,'is_rad=') &
568	+len('is_rad='):),*) is_rad(iq)
569	write(,) ' Parameter value (traceur.def) : is_rad=', &
570	is_rad(iq)
571	else
572	write(,) ' Parameter value (default) : is_rad=', &
573	is_rad(iq)
574	end if
575	! option is_recomb
576	if (index(tracline,'is_recomb=') /= 0) then
577	read(tracline(index(tracline,'is_recomb=') &
578	+len('is_recomb='):),*) is_recomb(iq)
579	write(,) ' Parameter value (traceur.def) : is_recomb=', &
580	is_recomb(iq)
581	else
582	write(,) ' Parameter value (default) : is_recomb=', &
583	is_recomb(iq)
584	end if
585	! option is_recomb_qset
586	if (index(tracline,'is_recomb_qset=') /= 0) then
587	read(tracline(index(tracline,'is_recomb_qset=') &
588	+len('is_recomb_qset='):),*) is_recomb_qset(iq)
589	write(,) ' Parameter value (traceur.def) :'// &
590	' is_recomb_qset=', &
591	is_recomb_qset(iq)
592	else
593	write(,) ' Parameter value (default) :'// &
594	' is_recomb_qset=', &
595	is_recomb_qset(iq)
596	endif
597	! option is_recomb_qotf
598	if (index(tracline,'is_recomb_qotf=') /= 0) then
599	read(tracline(index(tracline,'is_recomb_qotf=') &
600	+len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
601	write(,) ' Parameter value (traceur.def) :'// &
602	' is_recomb_qotf=', &
603	is_recomb_qotf(iq)
604	else
605	write(,) ' Parameter value (default) :'// &
606	' is_recomb_qotf=', &
607	is_recomb_qotf(iq)
608	end if
609	!option is_generic (LT)
610	if (index(tracline,'is_generic=') /=0) then
611	read(tracline(index(tracline,'is_generic=') &
612	+len('is_generic='):),*) is_generic(iq)
613	write(,) ' Parameter value (traceur.def) :'// &
614	' is_generic=', is_generic(iq)
615	else
616	write(,) ' Parameter value (default) :'// &
617	' is_generic=', is_generic(iq)
618	endif
619	!option radius
620	if (index(tracline,'radius=') .ne. 0) then
621	read(tracline(index(tracline,'radius=') &
622	+len('radius='):),*) radius(iq)
623	write(,)'Parameter value (traceur.def) :'// &
624	"radius=",radius(iq), "m"
625	else
626	write(,) ' Parameter value (default) :'// &
627	' radius=', radius(iq)," m"
628	endif
629	!option rho
630	if (index(tracline,'rho=') .ne. 0) then
631	read(tracline(index(tracline,'rho=') &
632	+len('rho='):),*) rho_q(iq)
633	write(,)'Parameter value (traceur.def) :'// &
634	"rho=",rho_q(iq)
635	else
636	write(,) ' Parameter value (default) :'// &
637	' rho=', rho_q(iq)
638	endif
639	end subroutine get_tracdat
640
641	end
642

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

Download in other formats:

Original Format