! ------------------------------------------------------------- ! MIE CODE designed for the LMD/GCM ! This code generates the look-up tables used by to GCM to ! compute the cloud scattering parameters in each grid box. ! ! Mie code: Bohren and Huffman (1983), modified by B.T.Draine ! Interface and integration: F. Montmessin ! Implementation into the LMD/GCM: J.-B. Madeleine 08W38 ! ------------------------------------------------------------- PROGRAM optpropgen use bhmie, only: bh_mie, s_size use math, only: pi IMPLICIT NONE ! ------------------------------------------------------------- ! PARAMETERS USED BY THE MIE CODE BEFORE CONVOLUTION | ! param_mie | ! ------------------------------------------------------------- character(len=*), parameter :: & optpropgen_version = '1.0.1' integer, parameter :: n_args = 2, n_opt_args = 2, file_name_size = 255 character(len=file_name_size) :: file_in, file_out character(len=file_name_size), dimension(n_args) :: args character(len=file_name_size), dimension(n_opt_args) :: opt_args doubleprecision :: & rmin, & ! Min. and max radius of the rmax ! interval in METERS INTEGER nsize, nsun ! Number of integration subintervals ! Useful when activice is true. ! Parameters for water ice: PARAMETER (nsize = 10000, nsun = 207) ! 127 VIS 81 IR ! ------------------------------------------------------------- doubleprecision QextMono(nsize, nsun), & QscatMono(nsize, nsun), gMono(nsize, nsun) doubleprecision CextMono(nsize, nsun), CscatMono(nsize, nsun) DOUBLE PRECISION radiustab(nsize) ! VARIABLES USED BY THE BHMIE CODE doubleprecision x(nsun) doubleprecision ni(nsun) doubleprecision nr(nsun) doubleprecision xf COMPLEX(kind=8) refrel INTEGER NANG COMPLEX(kind=8) S1(s_size), S2(s_size) doubleprecision QBACK ! VARIABLES USED FOR THE OUTPUT FILES CHARACTER(LEN = 5) radius_id ! ------------------------------------------------------------- ! PARAMETERS USED BY THE CONVOLUTION | ! param_conv | ! ------------------------------------------------------------- INTEGER nreff PARAMETER (nreff = 30) ! ------------------------------------------------------------- doubleprecision reff(nreff), nueff doubleprecision logvratreff DOUBLE PRECISION vratreff doubleprecision r_g, sigma_g DOUBLE PRECISION dfi DOUBLE PRECISION dfi_tmp(nsize + 1) doubleprecision radiusint(nsize + 1) doubleprecision Qext_tmp(nsun), omeg_tmp(nsun), g_tmp(nsun), pirr_tmp(nsun) doubleprecision Qscat_tmp(nsun) doubleprecision Cext_tmp(nsun), Cscat_tmp(nsun) doubleprecision :: Qextint(nsun, nreff) doubleprecision :: Qscatint(nsun, nreff) doubleprecision :: omegint(nsun, nreff) doubleprecision :: gint(nsun, nreff) ! Local variables doubleprecision Qext_mono(nsun, nsize), & Omeg_mono(nsun, nsize), & g_mono(nsun, nsize) DOUBLE PRECISION :: vrat INTEGER i, j, l call get_optpropgen_command_arguments(args, opt_args) write(*, '("Optpropgen ", A, /, "____", /)') optpropgen_version call main() contains subroutine get_optpropgen_command_arguments(args, opt_args) ! """ ! Retrieve the txt2h5 arguments that were passed on the command line when it was invoked. ! :param args: array storing the arguments ! """ implicit none character(len=file_name_size), dimension(:), intent(inout) :: args, opt_args character(len=file_name_size) :: arg integer :: i, j, n_args, n_opt_args j = 0 n_args = size(args) n_opt_args = size(opt_args) args(:) = '' opt_args(:) = '' do i = 1, command_argument_count() call get_command_argument(i, arg) if(trim(arg) == '-v' .or. trim(arg) == '--version') then write(*, '("Optropgen ", A)') optpropgen_version stop elseif(trim(arg) == '-h' .or. trim(arg) == '--help') then call print_help() stop else j = j + 1 end if if (len_trim(arg) == 0) then if (j < n_args) then write(*, '("Error: txt2h5 takes ", I0, " arguments but ", I0, " were given")') & n_args, j stop end if exit end if if(j > n_args) then opt_args(j - n_args) = trim(arg) if (j > n_args + n_opt_args) then write(*, '("Error: txt2h5 takes ", I0, " arguments but ", I0, " were given")') & n_args + n_opt_args, j stop end if else if (j > 0) then args(j) = trim(arg) end if end do end subroutine get_optpropgen_command_arguments subroutine main ! Initialization file_in = trim(args(1)) file_out = trim(args(2)) if(opt_args(1) /= '') then read(opt_args(1), *) rmin else rmin = 1d-7 end if if(opt_args(2) /= '') then read(opt_args(2), *) rmax else rmax = 1d-4 end if ! Main radiustab(1) = rmin radiustab(nsize) = rmax vrat = log(rmax / rmin) / dble(nsize - 1) * 3d0 vrat = exp(vrat) do i = 2, nsize - 1 radiustab(i) = radiustab(i - 1) * vrat**(1. / 3.) enddo write(*, '("Calculating for ", I0," particle radii between ", ES7.1, " and ", ES7.1, " m...")') & nsize, radiustab(1), radiustab(nsize) ! Number of angles between 0 and 90 degrees NANG = 10 ! ------------------------------------------------------------- ! OPTICAL INDICES: FILENAME | ! param_file | ! ------------------------------------------------------------- write(*, '("Reading relative refractive indices in file ''", A,"''")') & trim(file_in) OPEN(1, file=trim(file_in)) ! ------------------------------------------------------------- do l = 1, nsun ! do l = nsun, 1,-1 ! for Fe & Mg2SiO4 read(1, *) x(l), nr(l), ni(l) enddo close(1) write(*, '("Running BHMie...")') DO i = 1, nsize do l = 1, nsun refrel = cmplx(nr(l), ni(l), kind=8) xf = 2d0 * pi * radiustab(i) / x(l) call bh_mie(xf, refrel, NANG, S1, S2, & QextMono(i, l), QscatMono(i, l), QBACK, gMono(i, l)) CextMono(i, l) = QextMono(i, l) * pi * radiustab(i) * radiustab(i) CscatMono(i, l) = QscatMono(i, l) * pi * radiustab(i) * radiustab(i) Qext_mono(l, i) = QextMono(i, l) Omeg_mono(l, i) = QscatMono(i, l) / QextMono(i, l) g_mono(l, i) = gMono(i, l) enddo ENDDO ! ------------------------------------------------------------- ! PARAMETERS USED BY THE CONVOLUTION | ! param_conv | ! ------------------------------------------------------------- nueff = 0.3 ! Effective variance of the log-normal distr. !! reff(1) = rmin ! Minimum effective radius IF (nreff > 1) THEN reff(nreff) = rmax ! Maximum effective radius ! ------------------------------------------------------------- logvratreff = log(reff(nreff) / reff(1)) / float(nreff - 1) * 3. vratreff = exp(logvratreff) do i = 2, nreff - 1 reff(i) = reff(i - 1) * vratreff**(1. / 3.) enddo ENDIF ! Integration radius and effective variance radiusint(1) = 1.e-9 DO i = 2, nsize radiusint(i) = ((2. * vrat) / (vrat + 1.))**(1. / 3.) * & radiustab(i - 1) ENDDO radiusint(nsize + 1) = 1.e-2 ! Integration write(*, '("Integrating...")') DO j = 1, nreff sigma_g = log(1. + nueff) ! r_g and sigma_g are defined in r_g = exp(2.5 * sigma_g) ! [hansen_1974], "Light scattering in sigma_g = sqrt(sigma_g) ! planetary atmospheres", r_g = reff(j) / r_g ! Space Science Reviews 16 527-610. Qext_tmp(:) = 0. Qscat_tmp(:) = 0. omeg_tmp(:) = 0. g_tmp(:) = 0. pirr_tmp(:) = 0. Cext_tmp(:) = 0. Cscat_tmp(:) = 0. dfi_tmp(:) = log(radiusint(:) / r_g) / sqrt(2.) / sigma_g DO i = 1, nsize dfi = 0.5 * (erf(dfi_tmp(i + 1)) - erf(dfi_tmp(i))) DO l = 1, nsun Cext_tmp(l) = Cext_tmp(l) + CextMono(i, l) * dfi Cscat_tmp(l) = Cscat_tmp(l) + & CscatMono(i, l) * dfi Qext_tmp(l) = Qext_tmp(l) + & QextMono(i, l) * pi * radiustab(i) * radiustab(i) * dfi Qscat_tmp(l) = Qscat_tmp(l) + & QscatMono(i, l) * pi * radiustab(i) * radiustab(i) * dfi g_tmp(l) = g_tmp(l) + & gMono(i, l) * pi * radiustab(i) * radiustab(i) * & QscatMono(i, l) * dfi pirr_tmp(l) = pirr_tmp(l) + & pi * radiustab(i) * radiustab(i) * dfi ENDDO ENDDO DO l = 1, nsun Qextint(l, j) = Qext_tmp(l) / pirr_tmp(l) Qscatint(l, j) = Qscat_tmp(l) / pirr_tmp(l) omegint(l, j) = Qscat_tmp(l) / Qext_tmp(l) gint(l, j) = g_tmp(l) / Qscatint(l, j) / pirr_tmp(l) ENDDO ENDDO ! Writing the LMD/GCM output file ! ------------------------------- ! ------------------------------------------------------------- ! OUTPUT: FILENAME | ! param_output | ! ------------------------------------------------------------- write(*, '("Writting optical constants in file ''", A,"''")') & trim(file_out) OPEN(FILE=trim(file_out), UNIT = 60, & FORM = 'formatted', STATUS = 'replace') ! ------------------------------------------------------------- WRITE(UNIT = 60, FMT = '(a31)') '# Number of wavelengths (nwvl):' WRITE(UNIT = 60, FMT = '(i4)') nsun WRITE(UNIT = 60, FMT = '(a27)') '# Number of radius (nsize):' WRITE(UNIT = 60, FMT = '(i5)') nreff WRITE(UNIT = 60, FMT = '(a24)') '# Wavelength axis (wvl):' WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') x WRITE(UNIT = 60, FMT = '(a30)') '# Particle size axis (radius):' WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') reff WRITE(UNIT = 60, FMT = '(a29)') '# Extinction coef. Qext (ep):' DO j = 1, nreff WRITE(UNIT = radius_id, FMT = '(I5)') j WRITE(UNIT = 60, FMT = '(a21)') '# Radius number ' // radius_id WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') Qextint(:, j) ENDDO WRITE(UNIT = 60, FMT = '(a28)') '# Single Scat Albedo (omeg):' DO j = 1, nreff WRITE(UNIT = radius_id, FMT = '(I5)') j WRITE(UNIT = 60, FMT = '(a21)') '# Radius number ' // radius_id WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') omegint(:, j) ENDDO WRITE(UNIT = 60, FMT = '(a29)') '# Asymmetry Factor (gfactor):' DO j = 1, nreff WRITE(UNIT = radius_id, FMT = '(I5)') j WRITE(UNIT = 60, FMT = '(a21)') '# Radius number ' // radius_id WRITE(UNIT = 60, FMT = '(5(1x,e12.6))') gint(:, j) ENDDO CLOSE(60) ! ------------------------------------------------------------- write(*, '("Done")') end subroutine main subroutine print_help() implicit none character(len=*), parameter :: help_str = & "Usage: ./optpropgen.exe [options] path/to/refractive_index_file.dat path/to/output_file.ocst.txt & &[min_radius max_radius]"& //new_line('')//"& &Example: ./optpropgen.exe ../data/refractive_indices/optind_ice.dat & &../data/refractive_indices/optind_ice.dat 1e-7 1e-4"//new_line('')//"& &Options: "//new_line('')//"& & --help Display this information."//new_line('')//"& & --version Display software version information."//new_line('')//"& &"//new_line('')//"& &Refractive index ASCII (.dat) files must have 3 columns in the format:"//new_line('')//"& &wavelength(m) real_ref_rel imaginary_ref_rel"//new_line('')//"& &Where ref_rel is the complex refractive index of a sphere over the real index of the medium& &"//new_line('')//"& &Gitlab of the software:"//new_line('')//"& &" write(*, '(A)') help_str end subroutine print_help END