Context Navigation

initracer.F90 @ 3959

Last change on this file since 3959 was 3951, checked in by debatzbr, 7 weeks ago
Pluto PCM: Add variables, indices, and flags related to microphysical clouds BBT
File size: 20.7 KB

Line
1	SUBROUTINE initracer(ngrid,nq)
2
3	use surfdat_h, ONLY: dryness
4	USE tracer_h
5	USE callkeys_mod, only: optichaze,nb_monomer,haze,fractal,fasthaze,rad_haze,&
6	callmufi,callmuclouds
7	USE recombin_corrk_mod, ONLY: ini_recombin
8	USE mod_phys_lmdz_para, only: is_master, bcast
9	use aerosol_mod, only: iaero_haze,i_haze
10	IMPLICIT NONE
11	!=======================================================================
12	! subject:
13	! --------
14	! Initialization related to tracer
15	! (transported dust, water, chemical species, ice...)
16	!
17	! Name of the tracer
18	!
19	! Test of dimension :
20	! Initialize COMMON tracer in tracer.h, using tracer names provided
21	! by the argument nametrac
22	!
23	! author: F.Forget
24	! ------
25	! Ehouarn Millour (oct. 2008): identify tracers by their names
26	! Y. Jaziri & J. Vatant d'Ollone (2020) : modern traceur.def
27	! B. de Batz de Trenquelléon (2024): specific microphysical tracers
28	!=======================================================================
29
30	integer,intent(in) :: ngrid,nq
31
32	character(len=500) :: tracline ! to read traceur.def lines
33	integer :: blank !to store the index of 1st blank when reading tracers names
34	integer iq,ig,count,ierr
35	real r0_lift , reff_lift, rho_haze
36	integer nqhaze(nq) ! to store haze tracers
37	integer i, ia, block, j
38	character(len=20) :: txt ! to store some text
39	character(LEN=20) :: tracername ! to temporarily store text
40	character(LEN=20) :: str
41
42	!-----------------------------------------------------------------------
43	! radius(nq) ! aerosol particle radius (m)
44	! rho_q(nq) ! tracer densities (kg.m-3)
45	! qext(nq) ! Single Scat. Extinction coeff at 0.67 um
46	! rho_dust ! Mars dust density
47	! rho_ice ! Water ice density
48	! doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio
49	! varian ! Characteristic variance of log-normal distribution
50	!-----------------------------------------------------------------------
51
52	if (is_master) then ! only the master proc/thread needs do this
53
54	moderntracdef=.false. ! For modern traceur.def (default false, old type)
55
56	open(407, form = 'formatted', status = 'old', &
57	file = 'traceur.def', iostat=ierr)
58	if (ierr /=0) then
59	! call abort_physic('initracer',&
60	! 'Problem in opening traceur.def',1)
61	print*,'Problem in opening traceur.def'
62	stop
63	end if
64	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
65	READ(407,'(A)') tracline
66	IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
67	READ(tracline,*) nqtot ! Try standard traceur.def
68	ELSE
69	moderntracdef = .true.
70	DO
71	READ(407,'(A)',iostat=ierr) tracline
72	IF (ierr==0) THEN
73	IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
74	READ(tracline,*) nqtot
75	! Temporary not implemented solution
76	if (nqtot/=nq) then
77	! call abort_physic('initracer','Different number of '// &
78	! 'tracers in dynamics and physics not managed yet',1)
79	print*,'!= nbr oftracers in dynamics and physics not managed yet'
80	stop
81	endif
82	EXIT
83	ENDIF
84	ELSE ! If pb, or if reached EOF without having found nqtot
85	! call abort_physic('initracer','Unable to read numbers '// &
86	! 'of tracers in traceur.def',1)
87	print*,"unable to read numbers of tracer in tracer.def"
88	stop
89	ENDIF
90	ENDDO
91	ENDIF ! if modern or standard traceur.def
92
93	endif ! of if (is_master)
94
95	! share the information with other procs/threads (if any)
96	CALL bcast(nqtot)
97	CALL bcast(moderntracdef)
98
99	!! -----------------------------------------------------------------------
100	!! For the moment number of tracers in dynamics and physics are equal
101	nqtot=nq
102	!! we allocate once for all arrays in common in tracer_h.F90
103	!! (supposedly those are not used before call to initracer)
104	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
105	IF (.NOT.ALLOCATED(mmol)) ALLOCATE(mmol(nq))
106	IF (.NOT.ALLOCATED(aki)) ALLOCATE(aki(nqtot))
107	IF (.NOT.ALLOCATED(cpi)) ALLOCATE(cpi(nqtot))
108	IF (.NOT.ALLOCATED(radius)) ALLOCATE(radius(nq))
109	IF (.NOT.ALLOCATED(rho_q)) ALLOCATE(rho_q(nq))
110	IF (.NOT.ALLOCATED(qext)) ALLOCATE(qext(nq))
111	IF (.NOT.ALLOCATED(qextrhor)) ALLOCATE(qextrhor(nq))
112	! IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))
113	IF (.NOT.ALLOCATED(is_rad)) ALLOCATE(is_rad(nqtot))
114	IF (.NOT.ALLOCATED(is_recomb)) ALLOCATE(is_recomb(nqtot))
115	IF (.NOT.ALLOCATED(is_recomb_qset)) THEN
116	ALLOCATE(is_recomb_qset(nqtot))
117	ENDIF
118	IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN
119	ALLOCATE(is_recomb_qotf(nqtot))
120	ENDIF
121
122	!! initialization
123	mmol(:) = 0.
124	aki(:) = 0.
125	cpi(:) = 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
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: Read tracers names from traceur.def
140	do iq=1,nq
141	if (is_master) read(407,'(A)') tracline
142	call bcast(tracline)
143	blank = index(tracline,' ') ! Find position of 1st blank in tracline
144	noms(iq) = tracline(1:blank) !ensure that in Modern-trac case, noms is actually the name and not all properties
145	enddo
146
147	! Identify tracers by their names: (and set corresponding values of mmol)
148	! 0. initialize tracer indexes to zero:
149	! NB: igcm_* indexes are commons in 'tracer.h'
150	igcm_n2=0
151	igcm_ch4_gas=0
152	igcm_ch4_ice=0
153	igcm_co_gas=0
154	igcm_co_ice=0
155	igcm_C2H2_mugas=0
156	igcm_C2H6_mugas=0
157	igcm_C4H2_mugas=0
158	igcm_C6H6_mugas=0
159	igcm_HCN_mugas=0
160	igcm_prec_haze=0
161
162	nqhaze(:)=0
163	i=0
164	DO iq=1,nq
165	txt=noms(iq)
166	IF (txt(1:4).eq."haze") THEN
167	i=i+1
168	nqhaze(i)=iq
169	ENDIF
170	ENDDO
171	if ((haze.or.fasthaze).and.i==0) then
172	print*, 'Haze active but no haze tracer in traceur.def'
173	stop
174	endif
175
176	igcm_haze=0
177	igcm_haze10=0
178	igcm_haze30=0
179	igcm_haze50=0
180	igcm_haze100=0
181
182	! Eddy diffusion tracers
183	igcm_eddy1e6=0
184	igcm_eddy1e7=0
185	igcm_eddy5e7=0
186	igcm_eddy1e8=0
187	igcm_eddy5e8=0
188	write(,) 'initracer: noms() ', noms
189	rho_n2=1000 ! n2 ice
190	rho_ch4_ice=520. ! rho ch4 ice
191	rho_co_ice=520. ! rho ch4 ice
192	rho_haze=800. ! haze
193
194	! 1. find dust tracers
195	count=0
196
197	! 2. find chemistry and water tracers
198	do iq=1,nq
199	if (noms(iq).eq."n2") then
200	igcm_n2=iq
201	mmol(igcm_n2)=28.
202	count=count+1
203	write(,) 'Tracer ',count,' = n2'
204	endif
205	if (noms(iq).eq."ch4_gas") then
206	igcm_ch4_gas=iq
207	mmol(igcm_ch4_gas)=16.
208	count=count+1
209	write(,) 'Tracer ',count,' = ch4 gas'
210	endif
211	if (noms(iq).eq."ch4_ice") then
212	igcm_ch4_ice=iq
213	mmol(igcm_ch4_ice)=16.
214	radius(iq)=3.e-6
215	rho_q(iq)=rho_ch4_ice
216	count=count+1
217	write(,) 'Tracer ',count,' = ch4 ice'
218	endif
219	if (noms(iq).eq."co_gas") then
220	igcm_co_gas=iq
221	mmol(igcm_co_gas)=28.
222	count=count+1
223	write(,) 'Tracer ',count,' = co gas'
224	endif
225	if (noms(iq).eq."co_ice") then
226	igcm_co_ice=iq
227	mmol(igcm_co_ice)=28.
228	radius(iq)=3.e-6
229	rho_q(iq)=rho_co_ice
230	count=count+1
231	write(,) 'Tracer ',count,' = co ice'
232	endif
233	! Microphysics related tracers
234	if (noms(iq).eq."C2H2_mugas") then
235	igcm_C2H2_mugas=iq
236	mmol(igcm_C2H2_mugas)=26.04
237	count=count+1
238	write(,) 'Tracer ',count,' = C2H2 mugas'
239	endif
240	if (noms(iq).eq."C2H6_mugas") then
241	igcm_C2H6_mugas=iq
242	mmol(igcm_C2H6_mugas)=30.07
243	count=count+1
244	write(,) 'Tracer ',count,' = C2H6 mugas'
245	endif
246	if (noms(iq).eq."C4H2_mugas") then
247	igcm_C4H2_mugas=iq
248	mmol(igcm_C4H2_mugas)=50.05
249	count=count+1
250	write(,) 'Tracer ',count,' = C4H2 mugas'
251	endif
252	if (noms(iq).eq."C6H6_mugas") then
253	igcm_C6H6_mugas=iq
254	mmol(igcm_C6H6_mugas)=78.11
255	count=count+1
256	write(,) 'Tracer ',count,' = C6H6 mugas'
257	endif
258	if (noms(iq).eq."HCN_mugas") then
259	igcm_HCN_mugas=iq
260	mmol(igcm_HCN_mugas)=27.03
261	count=count+1
262	write(,) 'Tracer ',count,' = HCN mugas'
263	endif
264	! Haze tracers
265	if (noms(iq).eq."prec_haze") then
266	igcm_prec_haze=iq
267	count=count+1
268	write(,) 'Tracer ',count,' = prec haze'
269	endif
270	if (noms(iq).eq."haze") then
271	igcm_haze=iq
272	count=count+1
273	radius(iq)=rad_haze
274	rho_q(iq)=rho_haze
275	write(,) 'Tracer ',count,' = haze'
276	endif
277	if (noms(iq).eq."haze10") then
278	igcm_haze10=iq
279	count=count+1
280	radius(iq)=10.e-9
281	rho_q(iq)=rho_haze
282	write(,) 'Tracer ',count,' = haze10'
283	endif
284	if (noms(iq).eq."haze30") then
285	igcm_haze30=iq
286	count=count+1
287	radius(iq)=30.e-9
288	rho_q(iq)=rho_haze
289	write(,) 'Tracer ',count,' = haze30'
290	endif
291	if (noms(iq).eq."haze50") then
292	igcm_haze50=iq
293	count=count+1
294	radius(iq)=50.e-9
295	rho_q(iq)=rho_haze
296	write(,) 'Tracer ',count,' = haze50'
297	endif
298	if (noms(iq).eq."haze100") then
299	igcm_haze100=iq
300	count=count+1
301	radius(iq)=100.e-9
302	rho_q(iq)=rho_haze
303	write(,) 'Tracer ',count,' = haze100'
304	endif
305	! Eddy diffusion tracers
306	if (noms(iq).eq."eddy1e6") then
307	igcm_eddy1e6=iq
308	count=count+1
309	write(,) 'Tracer ',count,' = eddy1e6'
310	endif
311	if (noms(iq).eq."eddy1e7") then
312	igcm_eddy1e7=iq
313	count=count+1
314	write(,) 'Tracer ',count,' = eddy1e7'
315	endif
316	if (noms(iq).eq."eddy5e7") then
317	igcm_eddy5e7=iq
318	count=count+1
319	write(,) 'Tracer ',count,' = eddy5e7'
320	endif
321	if (noms(iq).eq."eddy1e8") then
322	igcm_eddy1e8=iq
323	count=count+1
324	write(,) 'Tracer ',count,' = eddy1e8'
325	endif
326	if (noms(iq).eq."eddy5e8") then
327	igcm_eddy5e8=iq
328	count=count+1
329	write(,) 'Tracer ',count,' = eddy5e8'
330	endif
331	enddo ! of do iq=1,nq
332
333	! 3. Find microphysics tracers:
334	! By convention they all have the prefix "mu_" (case sensitive !)
335	nmicro = 0
336	nmicro_ices = 0
337	IF (callmufi) THEN
338	DO iq=1,nq
339	str = noms(iq)
340	IF (str(1:3) == "mu_") THEN
341	nmicro = nmicro+1
342	count = count+1
343	ENDIF
344	ENDDO
345
346	! Checking the expected number of microphysical tracers:
347	! No cloud: nmicro = 4 aer.
348	! Clouds: nmicro = 4 aer. + 2 ccn + 1(+) ices
349	IF (callmuclouds) THEN
350	IF (nmicro < 7) THEN
351	WRITE(,) "initracer:error:"," Inconsistent number of microphysical tracers"
352	WRITE(,) "expected at least 7 tracers (clouds: on),", nmicro, " given"
353	CALL abort
354	ELSE
355	nmicro_ices = nmicro - 6
356	ENDIF
357
358	ELSE
359	IF (nmicro < 4) THEN
360	WRITE(,) "initracer:error:"," Inconsistent number of microphysical tracers"
361	WRITE(,) "expected at least 4 tracers (clouds: off),", nmicro, " given"
362	CALL abort
363	ELSE IF (nmicro > 4) THEN
364	WRITE(,) "!!! WARNING !!! initracer: I was expecting only four tracers, you gave me more."
365	CALL abort
366	ENDIF
367	ENDIF ! end of callmuclouds
368
369	! Microphysics indexes share the same values than original tracname.
370	IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
371	j = 1
372	DO i=1,nq
373	str = noms(i)
374	IF (str(1:3) == "mu_") THEN
375	micro_indx(j) = i
376	j=j+1
377	ENDIF
378	ENDDO
379
380	! Cloud related indexes initialize in inimufi subroutine.
381	IF (.NOT.ALLOCATED(micro_ice_indx)) ALLOCATE(micro_ice_indx(nmicro_ices))
382	IF (.NOT.ALLOCATED(micro_gas_indx)) ALLOCATE(micro_gas_indx(nmicro_ices))
383
384	ELSE
385	IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
386	IF (.NOT.ALLOCATED(micro_ice_indx)) ALLOCATE(micro_ice_indx(nmicro_ices))
387	IF (.NOT.ALLOCATED(micro_gas_indx)) ALLOCATE(micro_gas_indx(nmicro_ices))
388
389	ENDIF ! end of callmufi
390
391	! Check that we identified all tracers:
392	if (count.ne.nq) then
393	write(,) "initracer: found only ",count," tracers"
394	write(,) " expected ",nq
395	do iq=1,count
396	write(,)' ',iq,' ',trim(noms(iq))
397	enddo
398	else
399	write(,) "initracer: found all expected tracers, namely:"
400	do iq=1,nq
401	write(,)' ',iq,' ',trim(noms(iq))
402	enddo
403	endif
404
405	! Get data of tracers. Need to rewind traceur.def first
406	if (is_master) then
407	rewind(407)
408	do
409	read(407,'(A)') tracline
410	if (index(tracline,"#") .ne. 1) then
411	exit
412	endif
413	enddo
414	endif
415	do iq=1,nqtot
416	if (is_master) read(407,'(A)') tracline
417	call bcast(tracline)
418	call get_tracdat(iq, tracline)
419	enddo
420
421	if (is_master) close(407)
422
423	! Processing modern traceur options
424	if(moderntracdef) then
425	call ini_recombin
426	endif
427
428	!------------------------------------------------------------
429	! Initialisation tracers ....
430	!------------------------------------------------------------
431	! rho_q(1:nq)=0
432
433	rho_n2=1000. ! N2 ice density (kg.m-3)
434
435	lw_co=274000.
436	lw_ch4=586700.
437	lw_n2=2.5e5
438	write(,) "lw_n2 = ", lw_n2
439
440	if (callmufi) then
441	if (optichaze) then
442	if (callmuclouds) then
443	iaero_haze = 3
444	else
445	iaero_haze = 2
446	endif ! end of callmuclouds
447	endif ! end optichaze
448
449	elseif (haze) then
450	! the sedimentation radius remains radius(igcm_haze)
451	if (fractal) then
452	nmono=nb_monomer
453	else
454	nmono=1
455	endif ! end fractal
456
457	ia = 0
458	if (optichaze) then
459	ia = ia + 1
460	iaero_haze = ia
461	write(,) '--- number of haze aerosol = ', iaero_haze
462	block = 0 ! Only one type of haze is active : the first one set in traceur.def
463	do iq=1,nq
464	tracername=noms(iq)
465	write(,) "--> tracername ",iq,'/',nq,' = ',tracername
466	if (tracername(1:4).eq."haze".and.block.eq.0) then
467	i_haze=iq
468	block=1
469	write(,) "i_haze=",i_haze
470	write(,) "Careful: if you set many haze traceurs in &
471	traceur.def,only ",tracername," will be radiatively active &
472	(first one in traceur.def)"
473	endif
474	enddo
475	endif ! end optichaze
476	endif ! end callmufi or haze
477
478	! Output for records:
479	! ~~~~~~~~~~~~~~~~~~
480	write(,)
481	Write(,) '******** initracer : dust transport parameters :'
482	write(,) 'radius = ', radius
483	write(,) 'Qext = ', qext
484	write(,)
485
486	contains
487
488	subroutine get_tracdat(iq,tracline)
489	!-------------------ADDING NEW OPTIONS-------------------
490	! Duplicate if sentence for adding new options
491	! Do not forget to add the new variables in tracer_h.F90
492	! Do not forget to allocate and initialize the new variables above
493	! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
494	!-------------------------------------------------------
495	use tracer_h
496	implicit none
497	integer, intent(in) :: iq ! tracer index
498	character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
499
500	read(tracline,*) noms(iq)
501	! JVO 20 : We should add a sanity check aborting when duplicates in names !
502	write(,)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
503	! option mmol
504	if (index(tracline,'mmol=' ) /= 0) then
505	read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
506	mmol(iq)
507	write(,) ' Parameter value (traceur.def) : mmol=', &
508	mmol(iq)
509	else
510	write(,) ' Parameter value (default) : mmol=', &
511	mmol(iq)
512	end if
513	! option aki
514	if (index(tracline,'aki=' ) /= 0) then
515	read(tracline(index(tracline,'aki=')+len('aki='):),*) &
516	aki(iq)
517	write(,) ' Parameter value (traceur.def) : aki=', &
518	aki(iq)
519	else
520	write(,) ' Parameter value (default) : aki=', &
521	aki(iq)
522	end if
523	! option cpi
524	if (index(tracline,'cpi=' ) /= 0) then
525	read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
526	cpi(iq)
527	write(,) ' Parameter value (traceur.def) : cpi=', &
528	cpi(iq)
529	else
530	write(,) ' Parameter value (default) : cpi=', &
531	cpi(iq)
532	end if
533	! option is_rad
534	if (index(tracline,'is_rad=') /= 0) then
535	read(tracline(index(tracline,'is_rad=') &
536	+len('is_rad='):),*) is_rad(iq)
537	write(,) ' Parameter value (traceur.def) : is_rad=', &
538	is_rad(iq)
539	else
540	write(,) ' Parameter value (default) : is_rad=', &
541	is_rad(iq)
542	end if
543	! option is_recomb
544	if (index(tracline,'is_recomb=') /= 0) then
545	read(tracline(index(tracline,'is_recomb=') &
546	+len('is_recomb='):),*) is_recomb(iq)
547	write(,) ' Parameter value (traceur.def) : is_recomb=', &
548	is_recomb(iq)
549	else
550	write(,) ' Parameter value (default) : is_recomb=', &
551	is_recomb(iq)
552	end if
553	! option is_recomb_qset
554	if (index(tracline,'is_recomb_qset=') /= 0) then
555	read(tracline(index(tracline,'is_recomb_qset=') &
556	+len('is_recomb_qset='):),*) is_recomb_qset(iq)
557	write(,) ' Parameter value (traceur.def) :'// &
558	' is_recomb_qset=', &
559	is_recomb_qset(iq)
560	else
561	write(,) ' Parameter value (default) :'// &
562	' is_recomb_qset=', &
563	is_recomb_qset(iq)
564	endif
565	! option is_recomb_qotf
566	if (index(tracline,'is_recomb_qotf=') /= 0) then
567	read(tracline(index(tracline,'is_recomb_qotf=') &
568	+len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
569	write(,) ' Parameter value (traceur.def) :'// &
570	' is_recomb_qotf=', &
571	is_recomb_qotf(iq)
572	else
573	write(,) ' Parameter value (default) :'// &
574	' is_recomb_qotf=', &
575	is_recomb_qotf(iq)
576	end if
577	!option radius
578	if (index(tracline,'radius=') .ne. 0) then
579	read(tracline(index(tracline,'radius=') &
580	+len('radius='):),*) radius(iq)
581	write(,)'Parameter value (traceur.def) :'// &
582	"radius=",radius(iq), "m"
583	else
584	write(,) ' Parameter value (default) :'// &
585	' radius=', radius(iq)," m"
586	endif
587	!option rho
588	if (index(tracline,'rho=') .ne. 0) then
589	read(tracline(index(tracline,'rho=') &
590	+len('rho='):),*) rho_q(iq)
591	write(,)'Parameter value (traceur.def) :'// &
592	"rho=",rho_q(iq)
593	else
594	write(,) ' Parameter value (default) :'// &
595	' rho=', rho_q(iq)
596	endif
597	end subroutine get_tracdat
598
599	end
600

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

Download in other formats:

Original Format