Context Navigation

initracer.F90 @ 4113

Last change on this file since 4113 was 4040, checked in by tbertrand, 2 months ago
PLUTO PCM: Implementing paleo mode (to be continued) TB
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
190	! Densities :
191	rho_n2=1000. ! n2 ice
192	rho_ch4_ice=520. ! rho ch4 ice
193	rho_co_ice=520. ! rho CO ice
194	rho_haze=800. ! haze
195
196	! 1. find dust tracers
197	count=0
198
199	! 2. find chemistry and water tracers
200	do iq=1,nq
201	if (noms(iq).eq."n2") then
202	igcm_n2=iq
203	mmol(igcm_n2)=28.
204	rho_q(iq)=rho_n2
205	count=count+1
206	write(,) 'Tracer ',count,' = n2'
207	endif
208	if (noms(iq).eq."ch4_gas") then
209	igcm_ch4_gas=iq
210	mmol(igcm_ch4_gas)=16.
211	count=count+1
212	write(,) 'Tracer ',count,' = ch4 gas'
213	endif
214	if (noms(iq).eq."ch4_ice") then
215	igcm_ch4_ice=iq
216	mmol(igcm_ch4_ice)=16.
217	radius(iq)=3.e-6
218	rho_q(iq)=rho_ch4_ice
219	count=count+1
220	write(,) 'Tracer ',count,' = ch4 ice'
221	endif
222	if (noms(iq).eq."co_gas") then
223	igcm_co_gas=iq
224	mmol(igcm_co_gas)=28.
225	count=count+1
226	write(,) 'Tracer ',count,' = co gas'
227	endif
228	if (noms(iq).eq."co_ice") then
229	igcm_co_ice=iq
230	mmol(igcm_co_ice)=28.
231	radius(iq)=3.e-6
232	rho_q(iq)=rho_co_ice
233	count=count+1
234	write(,) 'Tracer ',count,' = co ice'
235	endif
236	! Microphysics related tracers
237	if (noms(iq).eq."C2H2_mugas") then
238	igcm_C2H2_mugas=iq
239	mmol(igcm_C2H2_mugas)=26.04
240	count=count+1
241	write(,) 'Tracer ',count,' = C2H2 mugas'
242	endif
243	if (noms(iq).eq."C2H6_mugas") then
244	igcm_C2H6_mugas=iq
245	mmol(igcm_C2H6_mugas)=30.07
246	count=count+1
247	write(,) 'Tracer ',count,' = C2H6 mugas'
248	endif
249	if (noms(iq).eq."C4H2_mugas") then
250	igcm_C4H2_mugas=iq
251	mmol(igcm_C4H2_mugas)=50.05
252	count=count+1
253	write(,) 'Tracer ',count,' = C4H2 mugas'
254	endif
255	if (noms(iq).eq."C6H6_mugas") then
256	igcm_C6H6_mugas=iq
257	mmol(igcm_C6H6_mugas)=78.11
258	count=count+1
259	write(,) 'Tracer ',count,' = C6H6 mugas'
260	endif
261	if (noms(iq).eq."HCN_mugas") then
262	igcm_HCN_mugas=iq
263	mmol(igcm_HCN_mugas)=27.03
264	count=count+1
265	write(,) 'Tracer ',count,' = HCN mugas'
266	endif
267	! Haze tracers
268	if (noms(iq).eq."prec_haze") then
269	igcm_prec_haze=iq
270	count=count+1
271	write(,) 'Tracer ',count,' = prec haze'
272	endif
273	if (noms(iq).eq."haze") then
274	igcm_haze=iq
275	count=count+1
276	radius(iq)=rad_haze
277	rho_q(iq)=rho_haze
278	write(,) 'Tracer ',count,' = haze'
279	endif
280	if (noms(iq).eq."haze10") then
281	igcm_haze10=iq
282	count=count+1
283	radius(iq)=10.e-9
284	rho_q(iq)=rho_haze
285	write(,) 'Tracer ',count,' = haze10'
286	endif
287	if (noms(iq).eq."haze30") then
288	igcm_haze30=iq
289	count=count+1
290	radius(iq)=30.e-9
291	rho_q(iq)=rho_haze
292	write(,) 'Tracer ',count,' = haze30'
293	endif
294	if (noms(iq).eq."haze50") then
295	igcm_haze50=iq
296	count=count+1
297	radius(iq)=50.e-9
298	rho_q(iq)=rho_haze
299	write(,) 'Tracer ',count,' = haze50'
300	endif
301	if (noms(iq).eq."haze100") then
302	igcm_haze100=iq
303	count=count+1
304	radius(iq)=100.e-9
305	rho_q(iq)=rho_haze
306	write(,) 'Tracer ',count,' = haze100'
307	endif
308	! Eddy diffusion tracers
309	if (noms(iq).eq."eddy1e6") then
310	igcm_eddy1e6=iq
311	count=count+1
312	write(,) 'Tracer ',count,' = eddy1e6'
313	endif
314	if (noms(iq).eq."eddy1e7") then
315	igcm_eddy1e7=iq
316	count=count+1
317	write(,) 'Tracer ',count,' = eddy1e7'
318	endif
319	if (noms(iq).eq."eddy5e7") then
320	igcm_eddy5e7=iq
321	count=count+1
322	write(,) 'Tracer ',count,' = eddy5e7'
323	endif
324	if (noms(iq).eq."eddy1e8") then
325	igcm_eddy1e8=iq
326	count=count+1
327	write(,) 'Tracer ',count,' = eddy1e8'
328	endif
329	if (noms(iq).eq."eddy5e8") then
330	igcm_eddy5e8=iq
331	count=count+1
332	write(,) 'Tracer ',count,' = eddy5e8'
333	endif
334	enddo ! of do iq=1,nq
335
336	! 3. Find microphysics tracers:
337	! By convention they all have the prefix "mu_" (case sensitive !)
338	nmicro = 0
339	nmicro_ices = 0
340	IF (callmufi) THEN
341	DO iq=1,nq
342	str = noms(iq)
343	IF (str(1:3) == "mu_") THEN
344	nmicro = nmicro+1
345	count = count+1
346	ENDIF
347	ENDDO
348
349	! Checking the expected number of microphysical tracers:
350	! No cloud: nmicro = 4 aer.
351	! Clouds: nmicro = 4 aer. + 2 ccn + 1(+) ices
352	IF (callmuclouds) THEN
353	IF (nmicro < 7) THEN
354	WRITE(,) "initracer:error:"," Inconsistent number of microphysical tracers"
355	WRITE(,) "expected at least 7 tracers (clouds: on),", nmicro, " given"
356	CALL abort
357	ELSE
358	nmicro_ices = nmicro - 6
359	ENDIF
360
361	ELSE
362	IF (nmicro < 4) THEN
363	WRITE(,) "initracer:error:"," Inconsistent number of microphysical tracers"
364	WRITE(,) "expected at least 4 tracers (clouds: off),", nmicro, " given"
365	CALL abort
366	ELSE IF (nmicro > 4) THEN
367	WRITE(,) "!!! WARNING !!! initracer: I was expecting only four tracers, you gave me more."
368	CALL abort
369	ENDIF
370	ENDIF ! end of callmuclouds
371
372	! Microphysics indexes share the same values than original tracname.
373	IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
374	j = 1
375	DO i=1,nq
376	str = noms(i)
377	IF (str(1:3) == "mu_") THEN
378	micro_indx(j) = i
379	j=j+1
380	ENDIF
381	ENDDO
382
383	! Cloud related indexes initialize in inimufi subroutine.
384	IF (.NOT.ALLOCATED(micro_ice_indx)) ALLOCATE(micro_ice_indx(nmicro_ices))
385	IF (.NOT.ALLOCATED(micro_gas_indx)) ALLOCATE(micro_gas_indx(nmicro_ices))
386
387	ELSE
388	IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro))
389	IF (.NOT.ALLOCATED(micro_ice_indx)) ALLOCATE(micro_ice_indx(nmicro_ices))
390	IF (.NOT.ALLOCATED(micro_gas_indx)) ALLOCATE(micro_gas_indx(nmicro_ices))
391
392	ENDIF ! end of callmufi
393
394	! Check that we identified all tracers:
395	if (count.ne.nq) then
396	write(,) "initracer: found only ",count," tracers"
397	write(,) " expected ",nq
398	do iq=1,count
399	write(,)' ',iq,' ',trim(noms(iq))
400	enddo
401	else
402	write(,) "initracer: found all expected tracers, namely:"
403	do iq=1,nq
404	write(,)' ',iq,' ',trim(noms(iq))
405	enddo
406	endif
407
408	! Get data of tracers. Need to rewind traceur.def first
409	if (is_master) then
410	rewind(407)
411	do
412	read(407,'(A)') tracline
413	if (index(tracline,"#") .ne. 1) then
414	exit
415	endif
416	enddo
417	endif
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	! Processing modern traceur options
427	if(moderntracdef) then
428	call ini_recombin
429	endif
430
431	!------------------------------------------------------------
432	! Initialisation tracers ....
433	!------------------------------------------------------------
434	! rho_q(1:nq)=0
435
436	rho_n2=1000. ! N2 ice density (kg.m-3)
437
438	lw_co=274000.
439	lw_ch4=586700.
440	lw_n2=2.5e5
441	write(,) "lw_n2 = ", lw_n2
442
443	if (callmufi) then
444	if (optichaze) then
445	if (callmuclouds) then
446	iaero_haze = 3
447	else
448	iaero_haze = 2
449	endif ! end of callmuclouds
450	endif ! end optichaze
451
452	elseif (haze) then
453	! the sedimentation radius remains radius(igcm_haze)
454	if (fractal) then
455	nmono=nb_monomer
456	else
457	nmono=1
458	endif ! end fractal
459
460	ia = 0
461	if (optichaze) then
462	ia = ia + 1
463	iaero_haze = ia
464	write(,) '--- number of haze aerosol = ', iaero_haze
465	block = 0 ! Only one type of haze is active : the first one set in traceur.def
466	do iq=1,nq
467	tracername=noms(iq)
468	write(,) "--> tracername ",iq,'/',nq,' = ',tracername
469	if (tracername(1:4).eq."haze".and.block.eq.0) then
470	i_haze=iq
471	block=1
472	write(,) "i_haze=",i_haze
473	write(,) "Careful: if you set many haze traceurs in &
474	traceur.def,only ",tracername," will be radiatively active &
475	(first one in traceur.def)"
476	endif
477	enddo
478	endif ! end optichaze
479	endif ! end callmufi or haze
480
481	! Output for records:
482	! ~~~~~~~~~~~~~~~~~~
483	write(,)
484	Write(,) '******** initracer : dust transport parameters :'
485	write(,) 'radius = ', radius
486	write(,) 'Qext = ', qext
487	write(,)
488
489	contains
490
491	subroutine get_tracdat(iq,tracline)
492	!-------------------ADDING NEW OPTIONS-------------------
493	! Duplicate if sentence for adding new options
494	! Do not forget to add the new variables in tracer_h.F90
495	! Do not forget to allocate and initialize the new variables above
496	! Please update list of options in "LMDZ.GENERIC/deftank/traceur.def.modern"
497	!-------------------------------------------------------
498	use tracer_h
499	implicit none
500	integer, intent(in) :: iq ! tracer index
501	character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
502
503	read(tracline,*) noms(iq)
504	! JVO 20 : We should add a sanity check aborting when duplicates in names !
505	write(,)"initracer: iq=",iq,"noms(iq)=",trim(noms(iq))
506	! option mmol
507	if (index(tracline,'mmol=' ) /= 0) then
508	read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
509	mmol(iq)
510	write(,) ' Parameter value (traceur.def) : mmol=', &
511	mmol(iq)
512	else
513	write(,) ' Parameter value (default) : mmol=', &
514	mmol(iq)
515	end if
516	! option aki
517	if (index(tracline,'aki=' ) /= 0) then
518	read(tracline(index(tracline,'aki=')+len('aki='):),*) &
519	aki(iq)
520	write(,) ' Parameter value (traceur.def) : aki=', &
521	aki(iq)
522	else
523	write(,) ' Parameter value (default) : aki=', &
524	aki(iq)
525	end if
526	! option cpi
527	if (index(tracline,'cpi=' ) /= 0) then
528	read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
529	cpi(iq)
530	write(,) ' Parameter value (traceur.def) : cpi=', &
531	cpi(iq)
532	else
533	write(,) ' Parameter value (default) : cpi=', &
534	cpi(iq)
535	end if
536	! option is_rad
537	if (index(tracline,'is_rad=') /= 0) then
538	read(tracline(index(tracline,'is_rad=') &
539	+len('is_rad='):),*) is_rad(iq)
540	write(,) ' Parameter value (traceur.def) : is_rad=', &
541	is_rad(iq)
542	else
543	write(,) ' Parameter value (default) : is_rad=', &
544	is_rad(iq)
545	end if
546	! option is_recomb
547	if (index(tracline,'is_recomb=') /= 0) then
548	read(tracline(index(tracline,'is_recomb=') &
549	+len('is_recomb='):),*) is_recomb(iq)
550	write(,) ' Parameter value (traceur.def) : is_recomb=', &
551	is_recomb(iq)
552	else
553	write(,) ' Parameter value (default) : is_recomb=', &
554	is_recomb(iq)
555	end if
556	! option is_recomb_qset
557	if (index(tracline,'is_recomb_qset=') /= 0) then
558	read(tracline(index(tracline,'is_recomb_qset=') &
559	+len('is_recomb_qset='):),*) is_recomb_qset(iq)
560	write(,) ' Parameter value (traceur.def) :'// &
561	' is_recomb_qset=', &
562	is_recomb_qset(iq)
563	else
564	write(,) ' Parameter value (default) :'// &
565	' is_recomb_qset=', &
566	is_recomb_qset(iq)
567	endif
568	! option is_recomb_qotf
569	if (index(tracline,'is_recomb_qotf=') /= 0) then
570	read(tracline(index(tracline,'is_recomb_qotf=') &
571	+len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
572	write(,) ' Parameter value (traceur.def) :'// &
573	' is_recomb_qotf=', &
574	is_recomb_qotf(iq)
575	else
576	write(,) ' Parameter value (default) :'// &
577	' is_recomb_qotf=', &
578	is_recomb_qotf(iq)
579	end if
580	!option radius
581	if (index(tracline,'radius=') .ne. 0) then
582	read(tracline(index(tracline,'radius=') &
583	+len('radius='):),*) radius(iq)
584	write(,)'Parameter value (traceur.def) :'// &
585	"radius=",radius(iq), "m"
586	else
587	write(,) ' Parameter value (default) :'// &
588	' radius=', radius(iq)," m"
589	endif
590	!option rho
591	if (index(tracline,'rho=') .ne. 0) then
592	read(tracline(index(tracline,'rho=') &
593	+len('rho='):),*) rho_q(iq)
594	write(,)'Parameter value (traceur.def) :'// &
595	"rho=",rho_q(iq)
596	else
597	write(,) ' Parameter value (default) :'// &
598	' rho=', rho_q(iq)
599	endif
600	end subroutine get_tracdat
601
602	end
603

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

Download in other formats:

Original Format