Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

optpropgen.f90 @ 3537

Last change on this file since 3537 was 3240, checked in by mturbet, 10 months ago
add optpropgen aerosol optical property code in the svn
File size: 14.4 KB

Line
1	! -------------------------------------------------------------
2	! MIE CODE designed for the LMD/GCM
3	! This code generates the look-up tables used by to GCM to
4	! compute the cloud scattering parameters in each grid box.
5	!
6	! Mie code: Bohren and Huffman (1983), modified by B.T.Draine
7	! Interface and integration: F. Montmessin
8	! Implementation into the LMD/GCM: J.-B. Madeleine 08W38
9	! -------------------------------------------------------------
10
11	PROGRAM optpropgen
12	use bhmie, only: bh_mie, s_size
13	use math, only: pi
14
15	IMPLICIT NONE
16
17	! -------------------------------------------------------------
18	! PARAMETERS USED BY THE MIE CODE BEFORE CONVOLUTION \|
19	! param_mie \|
20	! -------------------------------------------------------------
21	character(len=*), parameter :: &
22	optpropgen_version = '1.0.1'
23
24	integer, parameter :: n_args = 2, n_opt_args = 2, file_name_size = 255
25	character(len=file_name_size) :: file_in, file_out
26	character(len=file_name_size), dimension(n_args) :: args
27	character(len=file_name_size), dimension(n_opt_args) :: opt_args
28	doubleprecision :: &
29	rmin, & ! Min. and max radius of the
30	rmax ! interval in METERS
31
32	INTEGER nsize, nsun ! Number of integration subintervals
33	! Useful when activice is true.
34	! Parameters for water ice:
35	PARAMETER (nsize = 10000, nsun = 207) ! 127 VIS 81 IR
36
37	! -------------------------------------------------------------
38
39	doubleprecision QextMono(nsize, nsun), &
40	QscatMono(nsize, nsun), gMono(nsize, nsun)
41	doubleprecision CextMono(nsize, nsun), CscatMono(nsize, nsun)
42	DOUBLE PRECISION radiustab(nsize)
43
44	! VARIABLES USED BY THE BHMIE CODE
45
46	doubleprecision x(nsun)
47	doubleprecision ni(nsun)
48	doubleprecision nr(nsun)
49	doubleprecision xf
50	COMPLEX(kind=8) refrel
51	INTEGER NANG
52	COMPLEX(kind=8) S1(s_size), S2(s_size)
53
54	doubleprecision QBACK
55
56	! VARIABLES USED FOR THE OUTPUT FILES
57
58	CHARACTER(LEN = 5) radius_id
59
60	! -------------------------------------------------------------
61	! PARAMETERS USED BY THE CONVOLUTION \|
62	! param_conv \|
63	! -------------------------------------------------------------
64	INTEGER nreff
65	PARAMETER (nreff = 30)
66	! -------------------------------------------------------------
67	doubleprecision reff(nreff), nueff
68	doubleprecision logvratreff
69	DOUBLE PRECISION vratreff
70	doubleprecision r_g, sigma_g
71	DOUBLE PRECISION dfi
72	DOUBLE PRECISION dfi_tmp(nsize + 1)
73	doubleprecision radiusint(nsize + 1)
74	doubleprecision Qext_tmp(nsun), omeg_tmp(nsun), g_tmp(nsun), pirr_tmp(nsun)
75	doubleprecision Qscat_tmp(nsun)
76	doubleprecision Cext_tmp(nsun), Cscat_tmp(nsun)
77
78	doubleprecision :: Qextint(nsun, nreff)
79	doubleprecision :: Qscatint(nsun, nreff)
80	doubleprecision :: omegint(nsun, nreff)
81	doubleprecision :: gint(nsun, nreff)
82
83	! Local variables
84
85	doubleprecision Qext_mono(nsun, nsize), &
86	Omeg_mono(nsun, nsize), &
87	g_mono(nsun, nsize)
88	DOUBLE PRECISION :: vrat
89
90	INTEGER i, j, l
91
92	call get_optpropgen_command_arguments(args, opt_args)
93
94	write(*, '("Optpropgen ", A, /, "____", /)') optpropgen_version
95
96	call main()
97
98	contains
99	subroutine get_optpropgen_command_arguments(args, opt_args)
100	! """
101	! Retrieve the txt2h5 arguments that were passed on the command line when it was invoked.
102	! :param args: array storing the arguments
103	! """
104	implicit none
105
106	character(len=file_name_size), dimension(:), intent(inout) :: args, opt_args
107
108	character(len=file_name_size) :: arg
109	integer :: i, j, n_args, n_opt_args
110
111	j = 0
112	n_args = size(args)
113	n_opt_args = size(opt_args)
114
115	args(:) = ''
116	opt_args(:) = ''
117
118	do i = 1, command_argument_count()
119	call get_command_argument(i, arg)
120
121	if(trim(arg) == '-v' .or. trim(arg) == '--version') then
122	write(*, '("Optropgen ", A)') optpropgen_version
123
124	stop
125	elseif(trim(arg) == '-h' .or. trim(arg) == '--help') then
126	call print_help()
127
128	stop
129	else
130	j = j + 1
131	end if
132
133	if (len_trim(arg) == 0) then
134	if (j < n_args) then
135	write(*, '("Error: txt2h5 takes ", I0, " arguments but ", I0, " were given")') &
136	n_args, j
137	stop
138	end if
139
140	exit
141	end if
142
143	if(j > n_args) then
144	opt_args(j - n_args) = trim(arg)
145
146	if (j > n_args + n_opt_args) then
147	write(*, '("Error: txt2h5 takes ", I0, " arguments but ", I0, " were given")') &
148	n_args + n_opt_args, j
149
150	stop
151	end if
152	else if (j > 0) then
153	args(j) = trim(arg)
154	end if
155	end do
156	end subroutine get_optpropgen_command_arguments
157
158	subroutine main
159	! Initialization
160	file_in = trim(args(1))
161	file_out = trim(args(2))
162
163	if(opt_args(1) /= '') then
164	read(opt_args(1), *) rmin
165	else
166	rmin = 1d-7
167	end if
168
169	if(opt_args(2) /= '') then
170	read(opt_args(2), *) rmax
171	else
172	rmax = 1d-4
173	end if
174
175	! Main
176	radiustab(1) = rmin
177	radiustab(nsize) = rmax
178
179	vrat = log(rmax / rmin) / dble(nsize - 1) * 3d0
180	vrat = exp(vrat)
181
182	do i = 2, nsize - 1
183	radiustab(i) = radiustab(i - 1) * vrat**(1. / 3.)
184	enddo
185
186	write(*, '("Calculating for ", I0," particle radii between ", ES7.1, " and ", ES7.1, " m...")') &
187	nsize, radiustab(1), radiustab(nsize)
188
189	! Number of angles between 0 and 90 degrees
190	NANG = 10
191
192	! -------------------------------------------------------------
193	! OPTICAL INDICES: FILENAME \|
194	! param_file \|
195	! -------------------------------------------------------------
196	write(*, '("Reading relative refractive indices in file ''", A,"''")') &
197	trim(file_in)
198
199	OPEN(1, file=trim(file_in))
200
201	! -------------------------------------------------------------
202	do l = 1, nsun
203	! do l = nsun, 1,-1 ! for Fe & Mg2SiO4
204	read(1, *) x(l), nr(l), ni(l)
205	enddo
206	close(1)
207
208	write(*, '("Running BHMie...")')
209	DO i = 1, nsize
210	do l = 1, nsun
211	refrel = cmplx(nr(l), ni(l), kind=8)
212	xf = 2d0 * pi * radiustab(i) / x(l)
213	call bh_mie(xf, refrel, NANG, S1, S2, &
214	QextMono(i, l), QscatMono(i, l), QBACK, gMono(i, l))
215	CextMono(i, l) = QextMono(i, l) * pi * radiustab(i) * radiustab(i)
216	CscatMono(i, l) = QscatMono(i, l) * pi * radiustab(i) * radiustab(i)
217	Qext_mono(l, i) = QextMono(i, l)
218	Omeg_mono(l, i) = QscatMono(i, l) / QextMono(i, l)
219	g_mono(l, i) = gMono(i, l)
220	enddo
221	ENDDO
222
223	! -------------------------------------------------------------
224	! PARAMETERS USED BY THE CONVOLUTION \|
225	! param_conv \|
226	! -------------------------------------------------------------
227	nueff = 0.3 ! Effective variance of the log-normal distr. !!
228	reff(1) = rmin ! Minimum effective radius
229
230	IF (nreff > 1) THEN
231	reff(nreff) = rmax ! Maximum effective radius
232	! -------------------------------------------------------------
233	logvratreff = log(reff(nreff) / reff(1)) / float(nreff - 1) * 3.
234	vratreff = exp(logvratreff)
235	do i = 2, nreff - 1
236	reff(i) = reff(i - 1) * vratreff**(1. / 3.)
237	enddo
238	ENDIF
239
240	! Integration radius and effective variance
241
242	radiusint(1) = 1.e-9
243	DO i = 2, nsize
244	radiusint(i) = ((2. * vrat) / (vrat + 1.))*(1. / 3.) &
245	radiustab(i - 1)
246	ENDDO
247	radiusint(nsize + 1) = 1.e-2
248
249	! Integration
250	write(*, '("Integrating...")')
251	DO j = 1, nreff
252	sigma_g = log(1. + nueff) ! r_g and sigma_g are defined in
253	r_g = exp(2.5 * sigma_g) ! [hansen_1974], "Light scattering in
254	sigma_g = sqrt(sigma_g) ! planetary atmospheres",
255	r_g = reff(j) / r_g ! Space Science Reviews 16 527-610.
256
257	Qext_tmp(:) = 0.
258	Qscat_tmp(:) = 0.
259	omeg_tmp(:) = 0.
260	g_tmp(:) = 0.
261	pirr_tmp(:) = 0.
262	Cext_tmp(:) = 0.
263	Cscat_tmp(:) = 0.
264
265	dfi_tmp(:) = log(radiusint(:) / r_g) / sqrt(2.) / sigma_g
266	DO i = 1, nsize
267	dfi = 0.5 * (erf(dfi_tmp(i + 1)) - erf(dfi_tmp(i)))
268
269	DO l = 1, nsun
270	Cext_tmp(l) = Cext_tmp(l) + CextMono(i, l) * dfi
271	Cscat_tmp(l) = Cscat_tmp(l) + &
272	CscatMono(i, l) * dfi
273	Qext_tmp(l) = Qext_tmp(l) + &
274	QextMono(i, l) * pi * radiustab(i) * radiustab(i) * dfi
275	Qscat_tmp(l) = Qscat_tmp(l) + &
276	QscatMono(i, l) * pi * radiustab(i) * radiustab(i) * dfi
277	g_tmp(l) = g_tmp(l) + &
278	gMono(i, l) * pi * radiustab(i) * radiustab(i) * &
279	QscatMono(i, l) * dfi
280	pirr_tmp(l) = pirr_tmp(l) + &
281	pi * radiustab(i) * radiustab(i) * dfi
282	ENDDO
283	ENDDO
284
285	DO l = 1, nsun
286	Qextint(l, j) = Qext_tmp(l) / pirr_tmp(l)
287	Qscatint(l, j) = Qscat_tmp(l) / pirr_tmp(l)
288	omegint(l, j) = Qscat_tmp(l) / Qext_tmp(l)
289	gint(l, j) = g_tmp(l) / Qscatint(l, j) / pirr_tmp(l)
290	ENDDO
291
292	ENDDO
293
294	! Writing the LMD/GCM output file
295	! -------------------------------
296
297	! -------------------------------------------------------------
298	! OUTPUT: FILENAME \|
299	! param_output \|
300	! -------------------------------------------------------------
301	write(*, '("Writting optical constants in file ''", A,"''")') &
302	trim(file_out)
303	OPEN(FILE=trim(file_out), UNIT = 60, &
304	FORM = 'formatted', STATUS = 'replace')
305	! -------------------------------------------------------------
306
307	WRITE(UNIT = 60, FMT = '(a31)') '# Number of wavelengths (nwvl):'
308	WRITE(UNIT = 60, FMT = '(i4)') nsun
309
310	WRITE(UNIT = 60, FMT = '(a27)') '# Number of radius (nsize):'
311	WRITE(UNIT = 60, FMT = '(i5)') nreff
312
313	WRITE(UNIT = 60, FMT = '(a24)') '# Wavelength axis (wvl):'
314	WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') x
315
316	WRITE(UNIT = 60, FMT = '(a30)') '# Particle size axis (radius):'
317	WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') reff
318
319	WRITE(UNIT = 60, FMT = '(a29)') '# Extinction coef. Qext (ep):'
320	DO j = 1, nreff
321	WRITE(UNIT = radius_id, FMT = '(I5)') j
322	WRITE(UNIT = 60, FMT = '(a21)') '# Radius number ' // radius_id
323	WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') Qextint(:, j)
324	ENDDO
325
326	WRITE(UNIT = 60, FMT = '(a28)') '# Single Scat Albedo (omeg):'
327	DO j = 1, nreff
328	WRITE(UNIT = radius_id, FMT = '(I5)') j
329	WRITE(UNIT = 60, FMT = '(a21)') '# Radius number ' // radius_id
330	WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') omegint(:, j)
331	ENDDO
332
333	WRITE(UNIT = 60, FMT = '(a29)') '# Asymmetry Factor (gfactor):'
334	DO j = 1, nreff
335	WRITE(UNIT = radius_id, FMT = '(I5)') j
336	WRITE(UNIT = 60, FMT = '(a21)') '# Radius number ' // radius_id
337	WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') gint(:, j)
338	ENDDO
339
340	CLOSE(60)
341	! -------------------------------------------------------------
342	write(*, '("Done")')
343	end subroutine main
344
345	subroutine print_help()
346	implicit none
347
348	character(len=*), parameter :: help_str = &
349	"Usage: ./optpropgen.exe [options] path/to/refractive_index_file.dat path/to/output_file.ocst.txt &
350	&[min_radius max_radius]"&
351	//new_line('')//"&
352	&Example: ./optpropgen.exe ../data/refractive_indices/optind_ice.dat &
353	&../data/refractive_indices/optind_ice.dat 1e-7 1e-4"//new_line('')//"&
354	&Options: "//new_line('')//"&
355	& --help Display this information."//new_line('')//"&
356	& --version Display software version information."//new_line('')//"&
357	&"//new_line('')//"&
358	&Refractive index ASCII (.dat) files must have 3 columns in the format:"//new_line('')//"&
359	&wavelength(m) real_ref_rel imaginary_ref_rel"//new_line('')//"&
360	&Where ref_rel is the complex refractive index of a sphere over the real index of the medium&
361	&"//new_line('')//"&
362	&Gitlab of the software:"//new_line('')//"&
363	&<https://gitlab.obspm.fr/dblain/exorem>"
364
365	write(*, '(A)') help_str
366	end subroutine print_help
367	END

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

Download in other formats:

Original Format