Context Navigation

initracer.F90 @ 2709

Last change on this file since 2709 was 2706, checked in by aslmd, 4 years ago
is_condensable has replaced is_generic, some changes in getting back tracers names
File size: 21.1 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_condensable)) allocate(is_condensable(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_condensable(:)= 0 !LT
136	rho_q(:) = 0. !need to be init here if we want to read it from modern traceur with get_tracdat
137
138
139	! Initialization: copy tracer names from dynamics
140	do iq=1,nq
141	noms(iq)=nametrac(iq)
142	write(,)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
143	enddo
144
145
146	! Identify tracers by their names: (and set corresponding values of mmol)
147	! 0. initialize tracer indexes to zero:
148	! NB: igcm_* indexes are commons in 'tracer.h'
149	do iq=1,nq
150	igcm_dustbin(iq)=0
151	enddo
152	igcm_dust_mass=0
153	igcm_dust_number=0
154	igcm_h2o_vap=0
155	igcm_h2o_ice=0
156	igcm_co2=0
157	igcm_co=0
158	igcm_o=0
159	igcm_o1d=0
160	igcm_o2=0
161	igcm_o3=0
162	igcm_h=0
163	igcm_h2=0
164	igcm_oh=0
165	igcm_ho2=0
166	igcm_h2o2=0
167	igcm_n2=0
168	igcm_n=0
169	igcm_n2d=0
170	igcm_no=0
171	igcm_no2=0
172	igcm_ar=0
173	igcm_ar_n2=0
174	igcm_co2_ice=0
175
176	igcm_ch4=0
177	igcm_ch3=0
178	igcm_ch=0
179	igcm_3ch2=0
180	igcm_1ch2=0
181	igcm_cho=0
182	igcm_ch2o=0
183	igcm_ch3o=0
184	igcm_c=0
185	igcm_c2=0
186	igcm_c2h=0
187	igcm_c2h2=0
188	igcm_c2h3=0
189	igcm_c2h4=0
190	igcm_c2h6=0
191	igcm_ch2co=0
192	igcm_ch3co=0
193	igcm_hcaer=0
194	! ! LT generic tracers
195	! igcm_generic_Fe_vap = 0
196	! igcm_generic_Fe_ice = 0
197	! igcm_generic_Cr_vap = 0
198	! igcm_generic_Cr_ice = 0
199
200
201	! 1. find dust tracers
202	count=0
203
204	! 2. find chemistry and water tracers
205	do iq=1,nq
206	if (noms(iq).eq."co2") then
207	igcm_co2=iq
208	mmol(igcm_co2)=44.
209	count=count+1
210	! write(,) 'co2: count=',count
211	endif
212	if (noms(iq).eq."co2_ice") then
213	igcm_co2_ice=iq
214	mmol(igcm_co2_ice)=44.
215	count=count+1
216	! write(,) 'co2_ice: count=',count
217	endif
218	if (noms(iq).eq."h2o_vap") then
219	igcm_h2o_vap=iq
220	mmol(igcm_h2o_vap)=18.
221	count=count+1
222	! write(,) 'h2o_vap: count=',count
223	endif
224	if (noms(iq).eq."h2o_ice") then
225	igcm_h2o_ice=iq
226	mmol(igcm_h2o_ice)=18.
227	count=count+1
228	! write(,) 'h2o_ice: count=',count
229	endif
230	if (noms(iq).eq."co") then
231	igcm_co=iq
232	mmol(igcm_co)=28.
233	count=count+1
234	endif
235	if (noms(iq).eq."o") then
236	igcm_o=iq
237	mmol(igcm_o)=16.
238	count=count+1
239	endif
240	if (noms(iq).eq."o1d") then
241	igcm_o1d=iq
242	mmol(igcm_o1d)=16.
243	count=count+1
244	endif
245	if (noms(iq).eq."o2") then
246	igcm_o2=iq
247	mmol(igcm_o2)=32.
248	count=count+1
249	endif
250	if (noms(iq).eq."o3") then
251	igcm_o3=iq
252	mmol(igcm_o3)=48.
253	count=count+1
254	endif
255	if (noms(iq).eq."h") then
256	igcm_h=iq
257	mmol(igcm_h)=1.
258	count=count+1
259	endif
260	if (noms(iq).eq."h2") then
261	igcm_h2=iq
262	mmol(igcm_h2)=2.
263	count=count+1
264	endif
265	if (noms(iq).eq."oh") then
266	igcm_oh=iq
267	mmol(igcm_oh)=17.
268	count=count+1
269	endif
270	if (noms(iq).eq."ho2") then
271	igcm_ho2=iq
272	mmol(igcm_ho2)=33.
273	count=count+1
274	endif
275	if (noms(iq).eq."h2o2") then
276	igcm_h2o2=iq
277	mmol(igcm_h2o2)=34.
278	count=count+1
279	endif
280	if (noms(iq).eq."n2") then
281	igcm_n2=iq
282	mmol(igcm_n2)=28.
283	count=count+1
284	endif
285	if (noms(iq).eq."ch4") then
286	igcm_ch4=iq
287	mmol(igcm_ch4)=16.
288	count=count+1
289	endif
290	if (noms(iq).eq."ar") then
291	igcm_ar=iq
292	mmol(igcm_ar)=40.
293	count=count+1
294	endif
295	if (noms(iq).eq."n") then
296	igcm_n=iq
297	mmol(igcm_n)=14.
298	count=count+1
299	endif
300	if (noms(iq).eq."no") then
301	igcm_no=iq
302	mmol(igcm_no)=30.
303	count=count+1
304	endif
305	if (noms(iq).eq."no2") then
306	igcm_no2=iq
307	mmol(igcm_no2)=46.
308	count=count+1
309	endif
310	if (noms(iq).eq."n2d") then
311	igcm_n2d=iq
312	mmol(igcm_n2d)=28.
313	count=count+1
314	endif
315	if (noms(iq).eq."ch3") then
316	igcm_ch3=iq
317	mmol(igcm_ch3)=15.
318	count=count+1
319	endif
320	if (noms(iq).eq."ch") then
321	igcm_ch=iq
322	mmol(igcm_ch)=13.
323	count=count+1
324	endif
325	if (noms(iq).eq."3ch2") then
326	igcm_3ch2=iq
327	mmol(igcm_3ch2)=14.
328	count=count+1
329	endif
330	if (noms(iq).eq."1ch2") then
331	igcm_1ch2=iq
332	mmol(igcm_1ch2)=14.
333	count=count+1
334	endif
335	if (noms(iq).eq."cho") then
336	igcm_cho=iq
337	mmol(igcm_cho)=29.
338	count=count+1
339	endif
340	if (noms(iq).eq."ch2o") then
341	igcm_ch2o=iq
342	mmol(igcm_ch2o)=30.
343	count=count+1
344	endif
345	if (noms(iq).eq."ch3o") then
346	igcm_ch3o=iq
347	mmol(igcm_ch3o)=31.
348	count=count+1
349	endif
350	if (noms(iq).eq."c") then
351	igcm_c=iq
352	mmol(igcm_c)=12.
353	count=count+1
354	endif
355	if (noms(iq).eq."c2") then
356	igcm_c2=iq
357	mmol(igcm_c2)=24.
358	count=count+1
359	endif
360	if (noms(iq).eq."c2h") then
361	igcm_c2h=iq
362	mmol(igcm_c2h)=25.
363	count=count+1
364	endif
365	if (noms(iq).eq."c2h2") then
366	igcm_c2h2=iq
367	mmol(igcm_c2h2)=26.
368	count=count+1
369	endif
370	if (noms(iq).eq."c2h3") then
371	igcm_c2h3=iq
372	mmol(igcm_c2h3)=27.
373	count=count+1
374	endif
375	if (noms(iq).eq."c2h4") then
376	igcm_c2h4=iq
377	mmol(igcm_c2h4)=28.
378	count=count+1
379	endif
380	if (noms(iq).eq."c2h6") then
381	igcm_c2h6=iq
382	mmol(igcm_c2h6)=30.
383	count=count+1
384	endif
385	if (noms(iq).eq."ch2co") then
386	igcm_ch2co=iq
387	mmol(igcm_ch2co)=42.
388	count=count+1
389	endif
390	if (noms(iq).eq."ch3co") then
391	igcm_ch3co=iq
392	mmol(igcm_ch3co)=43.
393	count=count+1
394	endif
395	if (noms(iq).eq."hcaer") then
396	igcm_hcaer=iq
397	mmol(igcm_hcaer)=50.
398	count=count+1
399	endif
400
401	enddo ! of do iq=1,nq
402
403	! 3. find condensable traceurs different from h2o and co2
404	do iq=1,nq
405	if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then
406	count=count+1
407	endif
408	if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then
409	count=count+1
410	endif
411
412	enddo ! of do iq=1,nq
413
414	! check that we identified all tracers:
415	if (count.ne.nq) then
416	write(,) "initracer: found only ",count," tracers"
417	write(,) " expected ",nq
418	do iq=1,count
419	write(,)' ',iq,' ',trim(noms(iq))
420	enddo
421	! stop
422	else
423	write(,) "initracer: found all expected tracers, namely:"
424	do iq=1,nq
425	write(,)' ',iq,' ',trim(noms(iq))
426	enddo
427	endif
428
429	! Get data of tracers
430	do iq=1,nqtot
431	if (is_master) read(407,'(A)') tracline
432	! call bcast(tracline)
433	call get_tracdat(iq, tracline)
434	enddo
435
436	if (is_master) close(407)
437
438	! Calculate number of species in the chemistry
439	nesp = sum(is_chim)
440	write(,) 'Number of species in the chemistry nesp = ',nesp
441
442	! Calculate number of generic tracers
443	ngt = sum(is_condensable)
444	write(,) 'Number of generic tracer is ngt = ',ngt
445	! Processing modern traceur options
446	if(moderntracdef) then
447	call ini_recombin
448	endif
449
450	!------------------------------------------------------------
451	! Initialisation tracers ....
452	!------------------------------------------------------------
453	! rho_q(1:nq)=0
454
455	rho_dust=2500. ! Mars dust density (kg.m-3)
456	rho_ice=920. ! Water ice density (kg.m-3)
457	rho_co2=1620. ! CO2 ice density (kg.m-3)
458
459	! Initialization for water vapor
460	! ------------------------------
461	if(water) then
462	radius(igcm_h2o_vap)=0.
463	Qext(igcm_h2o_vap)=0.
464	alpha_lift(igcm_h2o_vap) =0.
465	alpha_devil(igcm_h2o_vap)=0.
466	qextrhor(igcm_h2o_vap)= 0.
467
468	!! this is defined in surfdat_h.F90
469	IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
470	IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
471
472	do ig=1,ngrid
473	if (ngrid.ne.1) watercaptag(ig)=.false.
474	dryness(ig) = 1.
475	enddo
476
477
478	radius(igcm_h2o_ice)=3.e-6
479	rho_q(igcm_h2o_ice)=rho_ice
480	Qext(igcm_h2o_ice)=0.
481	! alpha_lift(igcm_h2o_ice) =0.
482	! alpha_devil(igcm_h2o_ice)=0.
483	qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice) &
484	/ (rho_ice*radius(igcm_h2o_ice))
485
486
487	end if ! (water)
488
489
490	!
491	! some extra (possibly redundant) sanity checks for tracers:
492	! ---------------------------------------------------------
493	if (water) then
494	! verify that we indeed have h2o_vap and h2o_ice tracers
495	if (igcm_h2o_vap.eq.0) then
496	write(,) "initracer: error !!"
497	write(,) " cannot use water option without ",&
498	"an h2o_vap tracer !"
499	stop
500	endif
501	if (igcm_h2o_ice.eq.0) then
502	write(,) "initracer: error !!"
503	write(,) " cannot use water option without ",&
504	"an h2o_ice tracer !"
505	stop
506	endif
507	endif
508
509
510	! Output for records:
511	! ~~~~~~~~~~~~~~~~~~
512	write(,)
513	Write(,) '******** initracer : dust transport parameters :'
514	write(,) 'alpha_lift = ', alpha_lift
515	write(,) 'alpha_devil = ', alpha_devil
516	write(,) 'radius = ', radius
517	write(,) 'Qext = ', qext
518	write(,)
519
520	contains
521
522	subroutine get_tracdat(iq,tracline)
523	!-------------------ADDING NEW OPTIONS-------------------
524	! Duplicate if sentence for adding new options
525	! Do not forget to add the new variables in tracer_h.F90
526	! Do not forget to allocate and initialize the new variables above
527	! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
528	!-------------------------------------------------------
529	use tracer_h
530	implicit none
531	integer, intent(in) :: iq ! tracer index
532	character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
533
534	read(tracline,*) noms(iq)
535	! JVO 20 : We should add a sanity check aborting when duplicates in names !
536	write(,)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
537	! option mmol
538	if (index(tracline,'mmol=' ) /= 0) then
539	read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
540	mmol(iq)
541	write(,) ' Parameter value (traceur.def) : mmol=', &
542	mmol(iq)
543	else
544	write(,) ' Parameter value (default) : mmol=', &
545	mmol(iq)
546	end if
547	! option aki
548	if (index(tracline,'aki=' ) /= 0) then
549	read(tracline(index(tracline,'aki=')+len('aki='):),*) &
550	aki(iq)
551	write(,) ' Parameter value (traceur.def) : aki=', &
552	aki(iq)
553	else
554	write(,) ' Parameter value (default) : aki=', &
555	aki(iq)
556	end if
557	! option cpi
558	if (index(tracline,'cpi=' ) /= 0) then
559	read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
560	cpi(iq)
561	write(,) ' Parameter value (traceur.def) : cpi=', &
562	cpi(iq)
563	else
564	write(,) ' Parameter value (default) : cpi=', &
565	cpi(iq)
566	end if
567	! option is_chim
568	if (index(tracline,'is_chim=' ) /= 0) then
569	read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) &
570	is_chim(iq)
571	write(,) ' Parameter value (traceur.def) : is_chim=', &
572	is_chim(iq)
573	else
574	write(,) ' Parameter value (default) : is_chim=', &
575	is_chim(iq)
576	end if
577	! option is_rad
578	if (index(tracline,'is_rad=') /= 0) then
579	read(tracline(index(tracline,'is_rad=') &
580	+len('is_rad='):),*) is_rad(iq)
581	write(,) ' Parameter value (traceur.def) : is_rad=', &
582	is_rad(iq)
583	else
584	write(,) ' Parameter value (default) : is_rad=', &
585	is_rad(iq)
586	end if
587	! option is_recomb
588	if (index(tracline,'is_recomb=') /= 0) then
589	read(tracline(index(tracline,'is_recomb=') &
590	+len('is_recomb='):),*) is_recomb(iq)
591	write(,) ' Parameter value (traceur.def) : is_recomb=', &
592	is_recomb(iq)
593	else
594	write(,) ' Parameter value (default) : is_recomb=', &
595	is_recomb(iq)
596	end if
597	! option is_recomb_qset
598	if (index(tracline,'is_recomb_qset=') /= 0) then
599	read(tracline(index(tracline,'is_recomb_qset=') &
600	+len('is_recomb_qset='):),*) is_recomb_qset(iq)
601	write(,) ' Parameter value (traceur.def) :'// &
602	' is_recomb_qset=', &
603	is_recomb_qset(iq)
604	else
605	write(,) ' Parameter value (default) :'// &
606	' is_recomb_qset=', &
607	is_recomb_qset(iq)
608	endif
609	! option is_recomb_qotf
610	if (index(tracline,'is_recomb_qotf=') /= 0) then
611	read(tracline(index(tracline,'is_recomb_qotf=') &
612	+len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
613	write(,) ' Parameter value (traceur.def) :'// &
614	' is_recomb_qotf=', &
615	is_recomb_qotf(iq)
616	else
617	write(,) ' Parameter value (default) :'// &
618	' is_recomb_qotf=', &
619	is_recomb_qotf(iq)
620	end if
621	!option is_condensable (LT)
622	if (index(tracline,'is_condensable=') /=0) then
623	read(tracline(index(tracline,'is_condensable=') &
624	+len('is_condensable='):),*) is_condensable(iq)
625	write(,) ' Parameter value (traceur.def) :'// &
626	' is_condensable=', is_condensable(iq)
627	else
628	write(,) ' Parameter value (default) :'// &
629	' is_condensable=', is_condensable(iq)
630	endif
631	!option radius
632	if (index(tracline,'radius=') .ne. 0) then
633	read(tracline(index(tracline,'radius=') &
634	+len('radius='):),*) radius(iq)
635	write(,)'Parameter value (traceur.def) :'// &
636	"radius=",radius(iq), "m"
637	else
638	write(,) ' Parameter value (default) :'// &
639	' radius=', radius(iq)," m"
640	endif
641	!option rho
642	if (index(tracline,'rho=') .ne. 0) then
643	read(tracline(index(tracline,'rho=') &
644	+len('rho='):),*) rho_q(iq)
645	write(,)'Parameter value (traceur.def) :'// &
646	"rho=",rho_q(iq)
647	else
648	write(,) ' Parameter value (default) :'// &
649	' rho=', rho_q(iq)
650	endif
651	end subroutine get_tracdat
652
653	end
654

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

Download in other formats:

Original Format