Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

optpropgen.f90 @ 3585

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

Rev	Line
[3240]	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