Index: /trunk/LMDZ.GENERIC/changelog.txt
===================================================================
--- /trunk/LMDZ.GENERIC/changelog.txt	(revision 3239)
+++ /trunk/LMDZ.GENERIC/changelog.txt	(revision 3240)
@@ -1873,2 +1873,3 @@
 Change mineral aerosol names in suaer_corrk.F90 to make them consistent with our standard format
 Removed the outdated 'How_to_add_an_aerosol_to_the_GCM.pdf' file, the information is now in the wiki page : https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Advanced_Use_of_the_GCM#How_to_Change_the_Aerosols_Optical_Properties
+Add 'optpropgen' routines in the SVN (code written by JB Madeleine and Franck Montmessin, to generate aerosol optical properties)
Index: /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/bhmie.f90
===================================================================
--- /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/bhmie.f90	(revision 3240)
+++ /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/bhmie.f90	(revision 3240)
@@ -0,0 +1,214 @@
+module bhmie
+    use math, only: pi
+
+    implicit none
+
+    integer, parameter :: &
+        mxnang = 1000, &  ! Note: important that MXNANG be consistent with dimension of S1 and S2 in calling routine!
+        s_size = 2 * MXNANG - 1
+
+    save
+
+    contains
+        SUBROUTINE BH_MIE(X, REFREL, NANG, S1, S2, QEXT, QSCA, QBACK, GSCA)
+            implicit none
+            ! Declare parameters:
+
+            INTEGER NMXX
+            !      PARAMETER(MXNANG=1000,NMXX=15000)
+            PARAMETER(NMXX = 150000)
+            ! Arguments:
+            INTEGER NANG
+            doubleprecision GSCA, QBACK, QEXT, QSCA, X
+            COMPLEX(kind=8) REFREL
+            COMPLEX(kind=8) S1(s_size), S2(s_size)
+            ! Local variables:
+            INTEGER J, JJ, N, NSTOP, NMX, NN
+            DOUBLE PRECISION CHI, CHI0, CHI1, DANG, DX, EN, FN, P, PSI, PSI0, PSI1, &
+                    THETA, XSTOP, YMOD
+            DOUBLE PRECISION AMU(MXNANG), pi_(MXNANG), PI0(MXNANG), PI1(MXNANG), &
+                    TAU(MXNANG)
+            COMPLEX(kind=8) AN, AN1, BN, BN1, DREFRL, XI, XI1, Y
+            COMPLEX(kind=8) D(NMXX)
+            !***********************************************************************
+            ! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine
+            !    to calculate scattering and absorption by a homogenous isotropic
+            !    sphere.
+            ! Given:
+            !    X = 2*pi*a/lambda
+            !    REFREL = (complex refr. index of sphere)/(real index of medium)
+            !    NANG = number of angles between 0 and 90 degrees
+            !           (will calculate 2*NANG-1 directions from 0 to 180 deg.)
+            !           if called with NANG<2, will set NANG=2 and will compute
+            !           scattering for theta=0,90,180.
+            ! Returns:
+            !    S1(1 - 2*NANG-1) = -i*f_22 (incid. E perp. to scatt. plane,
+            !                                scatt. E perp. to scatt. plane)
+            !    S2(1 - 2*NANG-1) = -i*f_11 (incid. E parr. to scatt. plane,
+            !                                scatt. E parr. to scatt. plane)
+            !    QEXT = C_ext/pi*a**2 = efficiency factor for extinction
+            !    QSCA = C_sca/pi*a**2 = efficiency factor for scattering
+            !    QBACK = (dC_sca/domega)/pi*a**2
+            !          = backscattering efficiency [NB: this is (1/4*pi) smaller
+            !            than the "radar backscattering efficiency"; see Bohren &
+            !            Huffman 1983 pp. 120-123]
+            !    GSCA = <cos(theta)> for scattering
+            !
+            ! Original program taken from Bohren and Huffman (1983), Appendix A
+            ! Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26
+            ! in order to compute <cos(theta)>
+            ! 91/05/07 (BTD): Modified to allow NANG=1
+            ! 91/08/15 (BTD): Corrected error (failure to initialize P)
+            ! 91/08/15 (BTD): Modified to enhance vectorizability.
+            ! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1
+            ! 91/08/15 (BTD): Changed definition of QBACK.
+            ! 92/01/08 (BTD): Converted to full double precision and double complex
+            !                 eliminated 2 unneed lines of code
+            !                 eliminated redundant variables (e.g. APSI,APSI0)
+            !                 renamed RN -> EN = double precision N
+            !                 Note that DOUBLE COMPLEX and DCMPLX are not part
+            !                 of f77 standard, so this version may not be fully
+            !                 portable.  In event that portable version is
+            !                 needed, use src/bhmie_f77.f
+            ! 93/06/01 (BTD): Changed AMAX1 to generic function MAX
+            !***********************************************************************
+            !*** Safety checks
+            IF(NANG>MXNANG)STOP '***Error: NANG > MXNANG in bhmie'
+            IF(NANG<2)NANG = 2
+
+            DX = X
+            DREFRL = REFREL
+            Y = X * DREFRL
+            YMOD = ABS(Y)
+            !
+            !*** Series expansion terminated after NSTOP terms
+            !    Logarithmic derivatives calculated from NMX on down
+            XSTOP = X + 4d0 * X ** (1d0 / 3d0) + 2d0
+            NMX = int(MAX(XSTOP, YMOD)) + 15
+            ! BTD experiment 91/1/15: add one more term to series and compare results
+            !      NMX=AMAX1(XSTOP,YMOD)+16
+            ! test: compute 7001 wavelengths between .0001 and 1000 micron
+            ! for a=1.0micron SiC grain.  When NMX increased by 1, only a single
+            ! computed number changed (out of 4*7001) and it only changed by 1/8387
+            ! conclusion: we are indeed retaining enough terms in series!
+            NSTOP = int(XSTOP)
+            !
+            IF(NMX>NMXX)THEN
+                WRITE(0, *)'Error: NMX > NMXX=', NMXX, ' for |m|x=', YMOD
+                STOP
+            ENDIF
+            !*** Require NANG.GE.1 in order to calculate scattering intensities
+            DANG = 0d0
+            IF(NANG>1)DANG = 0.5d0 * pi / DBLE(NANG - 1)
+            DO J = 1, NANG
+                THETA = DBLE(J - 1) * DANG
+                AMU(J) = COS(THETA)
+            end do
+            DO J = 1, NANG
+                PI0(J) = 0d0
+                PI1(J) = 1d0
+            end do
+            NN = 2 * NANG - 1
+            DO J = 1, NN
+                S1(J) = (0d0, 0d0)
+                S2(J) = (0d0, 0d0)
+            end do
+            !
+            !*** Logarithmic derivative D(J) calculated by downward recurrence
+            !    beginning with initial value (0.,0.) at J=NMX
+            !
+            D(NMX) = (0d0, 0d0)
+            NN = NMX - 1
+            DO N = 1, NN
+                EN = NMX - N + 1
+                D(NMX - N) = (EN / Y) - (1. / (D(NMX - N + 1) + EN / Y))
+            end do
+            !
+            !*** Riccati-Bessel functions with real argument X
+            !    calculated by upward recurrence
+            !
+            PSI0 = COS(DX)
+            PSI1 = SIN(DX)
+            CHI0 = -SIN(DX)
+            CHI1 = COS(DX)
+            XI1 = DCMPLX(PSI1, -CHI1)
+            QSCA = 0d0
+            GSCA = 0d0
+            P = -1d0
+            DO N = 1, NSTOP
+                EN = N
+                FN = (2d0 * EN + 1d0) / (EN * (EN + 1d0))
+                ! for given N, PSI  = psi_n        CHI  = chi_n
+                !              PSI1 = psi_{n-1}    CHI1 = chi_{n-1}
+                !              PSI0 = psi_{n-2}    CHI0 = chi_{n-2}
+                ! Calculate psi_n and chi_n
+                PSI = (2d0 * EN - 1d0) * PSI1 / DX - PSI0
+                CHI = (2d0 * EN - 1d0) * CHI1 / DX - CHI0
+                XI = DCMPLX(PSI, -CHI)
+                !
+                !*** Store previous values of AN and BN for use
+                !    in computation of g=<cos(theta)>
+                IF(N>1)THEN
+                    AN1 = AN
+                    BN1 = BN
+                ENDIF
+                !
+                !*** Compute AN and BN:
+                AN = (D(N) / DREFRL + EN / DX) * PSI - PSI1
+                AN = AN / ((D(N) / DREFRL + EN / DX) * XI - XI1)
+                BN = (DREFRL * D(N) + EN / DX) * PSI - PSI1
+                BN = BN / ((DREFRL * D(N) + EN / DX) * XI - XI1)
+                !
+                !*** Augment sums for Qsca and g=<cos(theta)>
+                QSCA = QSCA + (2d0 * EN + 1d0) * (ABS(AN)**2 + ABS(BN)**2)
+                GSCA = GSCA + ((2d0 * EN + 1d0) / (EN * (EN + 1d0))) * &
+                        (dble(AN) * dble(BN) + AIMAG(AN) * AIMAG(BN))
+                IF(N>1)THEN
+                    GSCA = GSCA + ((EN - 1d0) * (EN + 1d0) / EN) * &
+                            (dble(AN1) * dble(AN) + AIMAG(AN1) * AIMAG(AN) + &
+                                    dble(BN1) * dble(BN) + AIMAG(BN1) * AIMAG(BN))
+                ENDIF
+                !
+                !*** Now calculate scattering intensity pattern
+                !    First do angles from 0 to 90
+                DO J = 1, NANG
+                    JJ = 2 * NANG - J
+                    pi_(J) = PI1(J)
+                    TAU(J) = EN * AMU(J) * pi_(J) - (EN + 1d0) * PI0(J)
+                    S1(J) = S1(J) + FN * (AN * pi_(J) + BN * TAU(J))
+                    S2(J) = S2(J) + FN * (AN * TAU(J) + BN * pi_(J))
+                end do
+                !
+                !*** Now do angles greater than 90 using pi_ and TAU from
+                !    angles less than 90.
+                !    P=1 for N=1,3,...; P=-1 for N=2,4,...
+                P = -P
+                DO J = 1, NANG - 1
+                    JJ = 2 * NANG - J
+                    S1(JJ) = S1(JJ) + FN * P * (AN * pi_(J) - BN * TAU(J))
+                    S2(JJ) = S2(JJ) + FN * P * (BN * pi_(J) - AN * TAU(J))
+                end do
+                PSI0 = PSI1
+                PSI1 = PSI
+                CHI0 = CHI1
+                CHI1 = CHI
+                XI1 = DCMPLX(PSI1, -CHI1)
+                !
+                !*** Compute pi_n for next value of n
+                !    For each angle J, compute pi_n+1
+                !    from pi_ = pi_n , PI0 = pi_n-1
+                DO J = 1, NANG
+                    PI1(J) = ((2. * EN + 1.) * AMU(J) * pi_(J) - (EN + 1.) * PI0(J)) / EN
+                    PI0(J) = pi_(J)
+                end do
+            end do
+            !
+            !*** Have summed sufficient terms.
+            !    Now compute QSCA,QEXT,QBACK,and GSCA
+            GSCA = 2d0 * GSCA / QSCA
+            QSCA = (2d0 / (DX * DX)) * QSCA
+            QEXT = (4d0 / (DX * DX)) * dble(S1(1))
+            QBACK = (ABS(S1(2 * NANG - 1)) / DX)**2 / pi
+            RETURN
+        END
+end module bhmie
Index: /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/math.f90
===================================================================
--- /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/math.f90	(revision 3240)
+++ /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/math.f90	(revision 3240)
@@ -0,0 +1,1756 @@
+module math
+    ! """
+    ! Contains useful math functions.
+    ! """
+    implicit none
+
+    integer, parameter :: &
+        dp = selected_real_kind(15,300)
+
+    real(kind=dp), parameter :: &
+        pi = 3.141592653589793238460_dp
+
+    doubleprecision, parameter :: &
+        sqrtpi = 0.5641895835477563d0  ! = 1 / sqrt(pi)
+
+    doubleprecision, parameter :: &
+        prec_high = 10d0 ** (-precision(0d0)), prec_low = 10d0 ** (-precision(0.))
+
+    save
+
+    contains
+        doubleprecision function chi2(observed_data, calculated_data)
+            ! """
+            ! Evaluate the goodeness of fit of calculated data in respect to observations.
+            !
+            ! input:
+            !   observed_data: array of observed data
+            !   calculated_data: array of calculated data
+            ! """
+            implicit None
+
+            doubleprecision, dimension(:), intent(in) :: observed_data, calculated_data
+
+            chi2 = sum((observed_data(:) - calculated_data(:))**2 / calculated_data(:))
+        end function chi2
+
+
+        doubleprecision function chi2_reduced(observed_data, calculated_data, input_deviation)
+            ! """
+            ! Evaluate the goodeness of fit of calculated data in respect to observations and observational errors.
+            !
+            ! input:
+            !   observed_data: array of observed data
+            !   calculated_data: array of calculated data
+            !   input_deviation: array of observational error
+            ! """
+            implicit None
+
+            doubleprecision, dimension(:), intent(in) :: observed_data, calculated_data, input_deviation
+
+            chi2_reduced = sum((observed_data(:) - calculated_data(:))**2 / input_deviation(:)**2)
+        end function chi2_reduced
+
+
+        doubleprecision function deg2rad(angle)
+            ! """
+            ! Return an angle in radians from an angle in degrees.
+            !
+            ! input:
+            !   angle: (degree) an angle
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: angle
+
+            deg2rad = angle * (pi / 180.0_dp)
+        end function deg2rad
+
+
+        doubleprecision function ellipse_polar_form(semi_major_axis, semi_minor_axis, angle)
+            ! """
+            ! Return the ellipse polar form of an ellipse relative to its center.
+            !
+            ! inputs:
+            !   angle: (deg) angle between the semi-major axis, the center, and the point of the ellipse
+            !   semi_major_axis: semi-major axis of the ellipse
+            !   semi_minor_axis: semi-minor axis of the ellipse
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: angle, semi_major_axis, semi_minor_axis
+
+            doubleprecision :: &
+                theta
+
+            theta = deg2rad(angle)
+
+            ellipse_polar_form = semi_major_axis * semi_minor_axis / &
+                sqrt((semi_minor_axis * cos(theta))**2 + (semi_major_axis * sin(theta))**2)
+
+            return
+
+        end function ellipse_polar_form
+
+
+        doubleprecision function gaussian(x, fwhm)
+            ! """
+            ! Return the gaussian of a value, for a given Full Width Half Maximum.
+            !
+            ! inputs:
+            !   x: value for which to calculate the gaussian
+            !   fwhm: Full Width Half Maximum of the gaussian function
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: x, fwhm
+
+            doubleprecision :: &
+                sigma
+
+            sigma = fwhm / (2D0 * sqrt(2D0 * log(2D0)))
+            gaussian = 1D0 / (sigma * sqrt(2D0 * pi)) * exp(-1D0/2D0 * (x / sigma)**2D0)
+        end function gaussian
+
+
+        doubleprecision function gaussian_noise() result(n)
+            ! """
+            ! Return a random value following a standard normal distriubtion PDF.
+            !
+            ! output:
+            !   n: random gaussian noise
+            ! """
+            implicit none
+
+            doubleprecision ::&
+                r  ! random number between 0 and 1
+
+            call init_random_seed
+            call random_number(r)
+
+            n = sqrt(2D0) * erfinv(2 * r - 1)
+
+            return
+
+        end function gaussian_noise
+
+
+        doubleprecision function sec(angle)
+            ! """
+            ! Return the secant of an angle in degrees.
+            !
+            ! input:
+            !   angle: (degree) an angle
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: angle
+
+            sec = 1D0 / cos(deg2rad(angle))
+        end function sec
+
+
+        doubleprecision function sgn(value)
+            ! """
+            ! Return the sign of a value.
+            ! Not to be confused with Fortran built-in function SIGN(A, B).
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: value
+
+            if (value >= 0D0) then
+                sgn = 1D0
+            else
+                sgn = -1D0
+            end if
+
+            return
+        end function sgn
+
+
+        doubleprecision function sinc_fwhm(x, fwhm)
+            ! """
+            ! Return the sine cardinal of a value, for a given Full Width Half Maximum.
+            !
+            ! inputs:
+            !   x: value for which to calculate sinc
+            !   fwhm: Full Width Half Maximum of the sinc function
+            ! """
+            implicit none
+
+            doubleprecision,intent(in) :: x, fwhm
+
+            doubleprecision, parameter :: &
+                k = 1D0 / (2D0 * 1.89549D0)  ! sinc fwhm constant
+
+            if (x > 0D0 - tiny(0.) .and. x < 0D0 - tiny(0.)) then
+                sinc_fwhm = 1D0  ! avoid NaN at x = 0
+            else
+                sinc_fwhm = sin(x / (k * fwhm)) / (x / (k * fwhm))
+            end if
+        end function sinc_fwhm
+
+
+        function arange(start, stop, step) result(array)
+            ! """
+            ! Return evenly spaced values within a given interval.
+            !
+            ! inputs:
+            !   start: start of interval; the array starts with this value
+            !   stop: end of interval; the array does not inlude this value
+            !   step: spacing between values
+            !
+            ! outputs:
+            !   array: array of evenly spaced values
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: start, stop, step
+            doubleprecision, dimension(:), allocatable :: array
+
+            integer :: &
+                i, &  ! index
+                n     ! number of elements in array
+
+            n = ceiling((stop - start) / step)
+
+            allocate(array(n))
+
+            do i = 1, n
+                array(i) = start + (i - 1) * step
+            end do
+
+            return
+        end function arange
+
+
+        function arange_include(start, stop, step) result(array)
+            ! """
+            ! Return evenly spaced values including a given interval.
+            !
+            ! inputs:
+            !   start: start of interval; the array starts with this value
+            !   stop: end of interval; the last element of array is always greater than this value by a max of 1 step
+            !   step: spacing between values
+            !
+            ! outputs:
+            !   array: array of evenly spaced values
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: start, stop, step
+            doubleprecision, dimension(:), allocatable :: array
+
+            integer :: &
+                i, &  ! index
+                n     ! number of elements in array
+
+            n = ceiling((stop - start) / step) + 1  ! operation +1 ensures that stop is included in array
+
+            allocate(array(n))
+
+            do i = 1, n
+                array(i) = start + (i - 1) * step
+            end do
+
+            return
+        end function arange_include
+
+
+        function convolve(signal, filter) result(convolved_signal)
+            ! """
+            ! Convolve the signal by the filter using classical convolution.
+            !
+            ! inputs:
+            !   signal: the signal array
+            !   filter: the filter array
+            !
+            ! output:
+            !   convolved_signal: signal convolved by filter
+            !
+            ! source: https://fortrandev.wordpress.com/2013/04/01/fortran-convolution-algorithm/
+            ! """
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: signal, filter
+            doubleprecision, dimension(:), allocatable :: convolved_signal, convolution
+
+            integer :: &
+                i, &
+                j, &
+                k, &
+                size_filter, &
+                size_signal, &
+                start_index
+
+            size_signal = size(signal)
+            size_filter = size(filter)
+
+            start_index = floor(size_filter / 2d0)
+
+            allocate(convolved_signal(size_signal), convolution(size_signal + size_filter))
+
+            ! Last part
+            do i = size_signal, size_signal + size_filter
+                convolution(i) = 0D0
+                j = size_signal
+
+                do k = 1, size_filter
+                    convolution(i) = convolution(i) + signal(j) * filter(k)
+                    j = j - 1
+                end do
+            end do
+
+            ! Middle part
+            do i = size_filter, size_signal
+                convolution(i) = 0D0
+                j = i
+
+                do k = 1, size_filter
+                    convolution(i) = convolution(i) + signal(j) * filter(k)
+                    j = j - 1
+                end do
+            end do
+
+            ! First part
+            do i = 1, size_filter
+                convolution(i) = 0D0
+                j = i
+                k = 1
+
+                do while (j > 0)
+                    convolution(i) = convolution(i) + signal(j) * filter(k)
+                    j = j - 1
+                    k = k + 1
+                end do
+            end do
+
+            convolved_signal(:) = convolution(start_index:start_index + size_signal - 1)
+        end function convolve
+
+
+        function slide_convolve(signal, filter) result(convolved_signal)
+            ! """
+            ! Slide convolve the signal by a sliding filter using classical convolution.
+            !
+            ! inputs:
+            !   signal: the signal array
+            !   filter: the filter 2D-array
+            !
+            ! output:
+            !   convolved_signal: signal convolved by filter
+            ! """
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: signal
+            doubleprecision, dimension(:, :), intent(in) :: filter
+            doubleprecision, dimension(:), allocatable :: convolved_signal, convolution
+
+            integer :: &
+                i, &  ! index
+                j, &  ! index
+                k, &  ! index
+                size_filter, &  ! size of the filter
+                size_signal, &  ! size of the signal
+                start_index     ! start index of the result
+
+            size_signal = size(signal)
+            size_filter = size(filter, 1)
+
+            ! Check sizes
+            ! TODO [low] remove the size constraint
+            if (size_signal < size_filter) then
+                print '("ERROR: slide_convolve: filter size must be lower than signal size, &
+                        &but sizes are ", I10, "and ", I10)', size_filter, size_signal
+                print '("Catched error: ", ES15.8)', signal(size_filter)
+                stop  ! be sure the program stops here
+            end if
+
+            start_index = floor(size_filter / 2D0) + 1
+
+            allocate(convolved_signal(size_signal), convolution(size_signal + size_filter))
+
+            convolution(:) = 0D0
+
+            ! Last part
+            do i = size_signal + 1, start_index + size_signal - 1
+                j = size_signal
+
+                do k = 1, size_filter
+                    convolution(i) = convolution(i) + signal(j) * filter(k, j)
+                    j = j - 1
+                end do
+            end do
+
+            ! Middle part
+            do i = size_filter, size_signal
+                j = i
+
+                do k = 1, size_filter
+                    convolution(i) = convolution(i) + signal(j) * filter(k, j)
+                    j = j - 1
+                end do
+            end do
+
+            ! First part
+            do i = start_index + 1, size_filter - 1
+                j = i
+
+                do k = 1, i
+                    convolution(i) = convolution(i) + signal(j) * filter(k, j)
+                    j = j - 1
+                end do
+            end do
+
+            ! First index
+            i = start_index
+            convolution(i) = 0D0
+            j = i
+
+            do k = 1, i
+                convolution(i) = convolution(i) + signal(j) * filter(k, j)
+                j = j - 1
+            end do
+
+            k = i + 1
+            convolution(i) = convolution(i) + signal(1) * filter(k, 1)
+
+            convolved_signal(:) = convolution(start_index:start_index + size_signal - 1)
+        end function slide_convolve
+
+
+        function search_sorted(array, value) result(index)
+            ! """
+            ! Find the index into a sorted array such that the corresponding value is the closest to 'value'.
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: value
+            doubleprecision, dimension(:), intent(in) :: array
+
+            integer :: index
+
+            integer :: &
+                i_low, &
+                i_high, &
+                i_mid
+
+            i_low = 1
+            i_mid = 0
+            i_high = size(array)
+
+            if(value < array(1)) then
+                index = 1
+                return
+            elseif(value > array(i_high)) then
+                index = i_high
+                return
+            end if
+
+            do while(i_low <= i_high)
+                i_mid = i_low + (i_high - i_low) / 2
+
+                if(value < array(i_mid)) then
+                    i_high = i_mid - 1
+                elseif(value > array(i_mid)) then
+                    i_low = i_mid + 1
+                else
+                    index = i_mid
+                    return
+                end if
+            end do
+
+            if(array(i_low) - value < value - array(i_high)) then
+                index = i_low
+            else
+                index = i_high
+            end if
+
+            return
+        end function search_sorted
+
+
+        function erfinv(x) result(r)
+            ! """
+            ! Calculate the inverse of the erf function.
+            !
+            ! input:
+            !   x: a number between -1 and 1
+            !
+            ! output:
+            !   r: x = erf(r)
+            ! """
+            implicit none
+
+            doubleprecision, parameter :: &  ! erfinv-approximation factors
+                a0 = 0.886226899, &
+                a1 = -1.645349621, &
+                a2 = 0.914624893, &
+                a3 = -0.140543331, &
+                b0 = 1, &
+                b1 = -2.118377725, &
+                b2 = 1.442710462, &
+                b3 = -0.329097515, &
+                b4 = 0.012229801, &
+                c0 = -1.970840454, &
+                c1 = -1.62490649, &
+                c2 = 3.429567803, &
+                c3 = 1.641345311, &
+                d0 = 1, &
+                d1 = 3.543889200, &
+                d2 = 1.637067800
+
+            doubleprecision :: x, r
+
+            integer ::&
+                sign_x  ! sign of elements of array x
+
+            doubleprecision :: &
+                x2, &  ! square of x
+                y  ! intermediate value
+
+            if(x < -1 .or. x > 1) then
+                print '("ERROR: erfinv(x): x must be in [-1;1]")'
+                stop
+            end if
+
+            if (x > -tiny(0.) .and. x < tiny(0.)) then
+                r = 0
+                return
+            end if
+
+            if (x > 0) then
+                sign_x = 1
+            else
+                sign_x = -1
+                x = -x
+            end if
+
+            if (x <= 0.7) then
+                x2 = x * x
+                r = x * (((a3 * x2 + a2) * x2 + a1) * x2 + a0)
+                r = r / (((b4 * x2 + b3) * x2 + b2) * x2 + b1) * x2 + b0
+            else
+                y = sqrt(-log((1 - x) / 2))
+                r = (((c3 * y + c2) * y + c1) * y + c0)
+                r = r / ((d2 * y + d1) * y + d0)
+            end if
+
+            r = r * sign_x
+            x = x * sign_x
+
+            r = r - (erf(r) - x) / (2 / sqrt(PI) * exp (-r * r))
+            r = r - (erf(r) - x) / (2 / sqrt(PI) * exp (-r * r))
+
+            return
+
+        end function erfinv
+
+
+        function interp(x_new, x, y) result(y_new)
+            ! """
+            ! Interpolate array y of abscisse x to array y_new of abscisse x_new.
+            !
+            ! inputs:
+            !   x_new: abscisses on which y will be interpolated
+            !   x: abscisses of y
+            !   y: array to interpole
+            !
+            ! output:
+            !   y_new: interpolation of y(x) on abscisses x_new
+            ! """
+
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: x, y, x_new
+            doubleprecision, dimension(size(x_new)) :: y_new
+
+            integer ::&
+                i, &  ! index
+                j     ! index
+
+            i = 1
+            j = 1
+            y_new(:) = 0D0
+
+            if(size(x) /= size(y)) then
+                print '("ERROR: interp: x and y must have the same size")'
+                stop
+            end if
+
+            if(x(1) < x(size(x))) then  ! ascending numerical order x array
+                if(x_new(1) < x(1) .or. x_new(size(x_new)) > x(size(x))) then
+                    print '("ERROR: interp: x_new is outside x boundaries")'
+                    print *, x(1), ' < ' , x_new(1),  '--', x_new(size(x_new)), ' < ', x(size(x))
+                    stop
+                end if
+
+                do while (i <= size(x_new))
+                    if(x_new(i) <= x(j + 1)) then
+                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)
+
+                        i = i + 1
+                    else
+                        j = j + 1
+                    end if
+                end do
+            else  ! descending numerical order x array
+                if(x_new(1) > x(1) .or. x_new(size(x_new)) < x(size(x))) then
+                    print '("ERROR: interp: x_new is outside x boundaries")'
+                    print *,x(1), ' > ', x_new(1), '--', x_new(size(x_new)), ' > ', x(size(x))
+                    stop
+                end if
+
+                do while (i <= size(x_new))
+                    if(x_new(i) >= x(j + 1)) then
+                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)
+
+                        i = i + 1
+                    else
+                        j = j + 1
+                    end if
+                end do
+            end if
+
+            return
+        end function interp
+
+
+        function interp_fast(x_new, x, y) result(y_new)
+            ! """
+            ! Interpolate array y of abscisse x to array y_new of abscisse x_new.
+            ! Assumptions (no check performed):
+            !   - x_new is within x boundaries
+            !   - x and y have the same size
+            !
+            ! inputs:
+            !   x_new: abscisses on which y will be interpolated
+            !   x: abscisses of y
+            !   y: array to interpole
+            !
+            ! output:
+            !   y_new: interpolation of y(x) on abscisses x_new
+            ! """
+
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: x, y, x_new
+            doubleprecision, dimension(size(x_new)) :: y_new
+
+            integer ::&
+                i, &  ! index
+                j     ! index
+
+            i = 1
+            j = 1
+            y_new(:) = 0D0
+
+            if(x(1) < x(size(x))) then  ! ascending numerical order x array
+                do while(i <= size(x_new))
+                    if(x_new(i) <= x(j + 1)) then
+                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)
+
+                        i = i + 1
+                    else
+                        j = j + 1
+                    end if
+                end do
+            else  ! descending numerical order x array
+                do while(i <= size(x_new))
+                    if(x_new(i) >= x(j + 1)) then
+                        y_new(i) = (x_new(i) - x(j)) * (y(j + 1) - y(j)) / (x(j + 1) - x(j)) + y(j)
+
+                        i = i + 1
+                    else
+                        j = j + 1
+                    end if
+                end do
+            end if
+
+            return
+        end function interp_fast
+
+
+        function interp_ex(x_new, x, y) result(y_new)
+            ! """
+            ! Interpolate array y of abscisse x to array y_new of abscisse x_new.
+            ! The points outside x range are linearly extrapolated.
+            ! Data must be in strict increasing order.
+            !
+            ! inputs:
+            !   x_new: abscisses on which y will be interpolated
+            !   x: abscisses of y
+            !   y: array to interpole
+            !
+            ! output:
+            !   y_new: interpolation of y(x) on abscisses x_new
+            ! """
+
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: x, y, x_new
+            doubleprecision, dimension(size(x_new)) :: y_new
+
+            integer ::&
+                i, &  ! index
+                i_max, &
+                i_min, &
+                n     ! index
+
+            n = size(x)
+            i_min = 0
+            i_max = 0
+
+            ! Check order
+            if(x(1) > x(n)) then
+                if(x_new(1) < x_new(size(x_new))) then
+                    write(*, '("Error: interp_ex: x_new numerical order must be the same than x (descending)")')
+
+                    stop
+                end if
+
+                ! Find where x_new is in data range
+                do i = 1, size(x_new)
+                    if(x_new(i) <= x(1) .and. i_min == 0) then
+                        i_min = i
+                    end if
+
+                    if(x_new(i) >= x(n)) then
+                        i_max = i
+                    end if
+
+                    if(x_new(i) < x(n)) then
+                        exit
+                    end if
+                end do
+            else
+                if(x_new(1) > x_new(size(x_new))) then
+                    write(*, '("Error: interp_ex: x_new numerical order must be the same than x (ascending)")')
+
+                    stop
+                end if
+
+                ! Find where x_new is in data range
+                do i = 1, size(x_new)
+                    if(x_new(i) >= x(1) .and. i_min == 0) then
+                        i_min = i
+                    end if
+
+                    if(x_new(i) <= x(n)) then
+                        i_max = i
+                    end if
+
+                    if(x_new(i) > x(n)) then
+                        exit
+                    end if
+                end do
+            end if
+
+            ! Interpolation within data range
+            if(i_min > 0 .and. i_max > 0) then
+                y_new(i_min:i_max) = interp_fast(x_new(i_min:i_max), x, y)
+            end if
+
+            if(i_min == 0 .or. i_max == 0) then
+                if(i_min > 0) then
+                    ! Higher bound extrapolation
+                    do i = 1, size(x_new)
+                        y_new(i) = (x_new(i) - x(n - 1)) * (y(n) - y(n - 1)) / (x(n) - x(n - 1)) + y(n - 1)
+                    end do
+                else if (i_max > 0) then
+                    ! Lower bound extrapolation
+                    do i = 1, size(x_new)
+                        y_new(i) = (x_new(i) - x(1)) * (y(2) - y(1)) / (x(2) - x(1)) + y(1)
+                    end do
+                else
+                    write(*, '("Error: interp_ex: unable to determine boundary")')
+
+                    stop
+                end if
+            else
+                ! Lower bound extrapolation
+                if(i_min > 1) then
+                    do i = 1, i_min - 1
+                        y_new(i) = (x_new(i) - x(1)) * (y(2) - y(1)) / (x(2) - x(1)) + y(1)
+                    end do
+                end if
+
+                ! Higher bound extrapolation
+                if(i_max < size(x_new)) then
+                    do i = i_max + 1, size(x_new)
+                        y_new(i) = (x_new(i) - x(n - 1)) * (y(n) - y(n - 1)) / (x(n) - x(n - 1)) + y(n - 1)
+                    end do
+                end if
+            end if
+        end function interp_ex
+
+
+        function interp_ex_0d(x_new, x, y) result(y_new)
+            ! """
+            ! Interpolate array y of abscisse x to value y_new at x_new.
+            ! This function only takes a double as x_new, not an array. Return a double, not an array.
+            ! The points outside x range are linearly extrapolated.
+            ! Data must be in strict increasing order.
+            !
+            ! inputs:
+            !   x_new: value at which y will be interpolated
+            !   x: abscisses of y
+            !   y: array to interpole
+            !
+            ! output:
+            !   y_new: interpolation of y(x) at x_new
+            ! """
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: x, y
+            doubleprecision :: x_new, y_new
+
+            doubleprecision, dimension(1) :: arr_tmp
+
+            arr_tmp = interp_ex([x_new], x, y)
+            y_new = arr_tmp(1)
+        end function interp_ex_0d
+
+
+        function mean_restep(x_new, x, y) result(y_new)
+            ! """
+            ! Change the values of y of abscisse x so that the new values y_new of abcisse x_new are the mean of y
+            ! between a step of x_new centered on x_new.
+            ! Both x and x_new must be regularly spaced.
+            ! Example:
+            ! >>> x_new = [0, 3, 6]
+            ! >>> x = [0, 1, 2, 3, 4, 5, 6]
+            ! >>> y = [1, 2, 4, 3, 6, 2, 0]
+            ! >>> mean_restep(x_new, x, y)
+            ! >>> [1.50, 4.33, 1.00]
+            !
+            ! inputs:
+            !   x_new: abscisses on which y will be interpolated
+            !   x: abscisses of y
+            !   y: array to interpole
+            !
+            ! output:
+            !   y_new: interpolation of y(x) on abscisses x_new
+            ! """
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: x, y, x_new
+            doubleprecision, dimension(size(x_new)) :: y_new
+
+            integer :: &
+                i, &  ! index
+                j, &  ! index
+                n     ! number of elements
+
+            doubleprecision :: &
+                sum_x, &
+                min_x, &
+                max_x
+
+            call check_inputs()
+
+            j = 1
+
+            ! First value
+            i = 1
+            n = 0
+            sum_x = 0D0
+            max_x = x_new(i) + (x_new(i + 1) - x_new(i)) / 2D0
+            min_x = x_new(i)
+
+            do while (x(j) < max_x .and. j < size(x))
+                if (x(j) >= min_x) then
+                    sum_x = sum_x + y(j)
+                    n = n + 1
+                end if
+
+                j = j + 1
+            end do
+
+            y_new(i) = sum_x / n
+
+            ! Intermediate values
+            do i = 2, size(x_new) - 1
+                n = 0
+                sum_x = 0D0
+                max_x = x_new(i) + (x_new(i + 1) - x_new(i)) / 2D0
+                min_x = x_new(i) - (x_new(i) - x_new(i - 1)) / 2D0
+
+                do while (x(j) < max_x .and. j < size(x))
+                    if (x(j) >= min_x) then
+                        sum_x = sum_x + y(j)
+                        n = n + 1
+                    end if
+
+                    j = j + 1
+                end do
+
+                y_new(i) = sum_x / n
+            end do
+
+            ! Last value
+            i = size(x_new)
+            n = 0
+            sum_x = 0D0
+            max_x = x_new(i)
+            min_x = x_new(i) - (x_new(i) - x_new(i - 1)) / 2D0
+
+            do while (x(j) < max_x .and. j < size(x))
+                if (x(j) >= min_x) then
+                    sum_x = sum_x + y(j)
+                    n = n + 1
+                end if
+
+                j = j + 1
+            end do
+
+            if(n == 0) then
+                y_new(i) = y(j)
+            else
+                y_new(i) = sum_x / n
+            end if
+
+            return
+
+            contains
+                subroutine check_inputs()
+                    implicit none
+
+                    if (size(x) /= size(y)) then
+                        print '("ERROR: mean_restep: x and y must have the same size")'
+                        stop
+                    end if
+
+                    if (x(1) > x(size(x))) then
+                        print '("ERROR: mean_restep: x must be in increasing order")'
+                        stop
+                    end if
+
+                    if(x_new(1) < x(1) .or. x_new(size(x_new)) > x(size(x))) then
+                        print '("ERROR: interp: x_new is outside x boundaries")'
+                        print *,x_new(1),x(1),x_new(size(x_new)),x(size(x))
+                        stop
+                    end if
+                end subroutine check_inputs
+        end function mean_restep
+
+
+        function reverse_array(array) result(reversed_array)
+            ! """
+            ! Reverse a 1-D double precision array.
+            ! """
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: array
+            integer :: head, tail
+
+            doubleprecision, dimension(size(array)) :: reversed_array
+
+            head = 1
+            tail = size(array)
+            reversed_array(:) = 0d0
+
+            do while(head < tail)
+                reversed_array(tail) = array(head)
+                reversed_array(head) = array(tail)
+                head = head + 1
+                tail = tail - 1
+            end do
+
+            return
+        end function reverse_array
+
+
+        function voigt(x, y)
+            ! """
+            ! Calculate the voigt function using an algorithm written by Humlicek JQSRT, 27, 437 (1982).
+            ! Calculate the complex probability function W(z) = exp(-z^2) * erfc(-z^2) in the superior complex plan
+            ! (i.e. for y >= 0). The real part of this function is the Voigt function.
+            !
+            ! The article shows that the rational function W can be written as (Eq. 11):
+            !   W_n(z) = (-iz) * sum_{k=1}^{n/2}(c_k((-iz)^2)^{k-1}) / ((-iz)^n + sum_{k=1}^{n/2}(d_k((-iz)^2)^{k-1}))
+            ! With z = x + iy, (-iz = t). The coefficent c and d are given in Table 3 of the article.
+            ! The maximal relative error on the imaginary and real parts is < 1e-4.
+            !
+            ! inputs:
+            !   x: real part coordinate
+            !   y: imaginary part coordinate
+            !
+            ! output:
+            !   voigt: value of the voigt function at z = x + iy
+            ! """
+            implicit none
+
+            doubleprecision, dimension(:), intent(in) :: x
+            doubleprecision, intent(in) :: y
+            doubleprecision, dimension(size(x)) :: voigt
+
+            integer :: &
+                i  ! index
+
+            double precision, dimension(size(x)) :: &
+                s  ! intermediate value
+
+            complex (kind=8), dimension(size(x)) :: &
+                t, &  ! intermediate value
+                u, &  ! intermediate value
+                w     ! intermediate value
+
+            do i=1,size(x)
+                t(i) = cmplx(y, -x(i), kind=8)
+            end do
+
+            s(:) = abs(x(:)) + y
+
+            where (s >= 15D0)  ! n = 2 (region i)
+                w = t * sqrtpi / (0.5D0 + t**2)
+                voigt = dble(w * sqrtpi)
+            elsewhere (s >= 5.5D0)  ! n = 4 (region II)
+                u = t**2
+                w = t * (1.410474D0 + u * sqrtpi) / (0.75D0 + u * (3D0 + u))
+                voigt = dble(w * sqrtpi)
+            elsewhere (0.195D0 * abs(x) - 0.176D0 <= y)  ! n = 6 (region III)
+                w = (16.4955D0 + t * (20.20933D0 + t * (11.96482D0 + t * (3.778987D0 + t * 0.5642236D0)))) / &
+                         (16.4955D0 + t * (38.82363D0 + t * (39.27121D0 + t * (21.69274D0 + t * (6.699398D0 + t)))))
+                voigt = dble(w * sqrtpi)
+            elsewhere  ! n = 8 (region IV)
+                u = t**2
+                w = exp(u) - &
+                         t * (36183.31D0 - &
+                              u * (3321.9905D0 - &
+                                   u * (1540.787D0 - &
+                                        u * (219.0313D0 - u *(35.76683D0 - u * (1.320522D0 - u * sqrtpi)))))) / &
+                         (32066.6D0 - &
+                          u * (24322.84D0 - &
+                               u * (9022.228D0 - &
+                                    u * (2186.181D0 - u * (364.2191D0 - u * (61.57037D0 - u * (1.841439D0 - u)))))))
+                voigt = dble(w * sqrtpi)
+            end where
+
+            return
+        end function voigt
+
+
+        doubleprecision function voigt_from_data(x, y, v) result(voigt)
+            ! """
+            ! """
+            implicit none
+
+            doubleprecision, intent(in) :: &
+                    x, &
+                    y, &
+                    v(400, 100)
+
+            doubleprecision, parameter :: &
+                    dp = 0.0025d0, &
+                    ds = 0.01d0, &
+                    yl = 1d0 / 99d0, &
+                    xl = 399d0
+
+            integer :: &
+                    ip, &
+                    is
+
+            doubleprecision :: &
+                    a, &
+                    b, &
+                    pressure_space, &
+                    x2, &
+                    x2y2, &
+                    s, &
+                    y2, &
+                    f
+
+            if(y > 5.4d0) then
+                x2 = x * x
+                y2 = y * y
+                x2y2 = x2 + y2
+                voigt = sqrtpi * y * (1d0 + (3d0 * x2 - y2) / (2d0 * x2y2 * x2y2) + &
+                        0.75d0 * (5d0 * x2 * x2 + y2 * y2 - 10d0 * x2 * y2) / (x2y2 * x2y2 * x2y2 * x2y2)) / x2y2
+            else
+                if(y <= yl) then
+                    b = exp(-x * x)
+                    pressure_space = 1d0 / (1d0 + x)
+                    ip = int(pressure_space / dp)
+                    f = pressure_space / dp - ip
+                    if(ip >= 1) then
+                        if(ip < 400) then
+                            a = (1d0 - f) * v(ip, 1) + f * v(ip + 1, 1)
+                        else
+                            a = v(ip, 1)
+                        end if
+                    else
+                        a = f * v(ip + 1, 1)
+                    end if
+                    a = a * a
+                    f = y / yl
+                    voigt = (1d0 - f) * b + f * a
+                else
+                    if(x > xl) then
+                        s = y / (1d0 + y)
+                        is = int(s / ds)
+                        b = v(1, is) * xl / x
+                        b = b * b
+                        a = v(1, is + 1) * xl / x
+                        a = a * a
+                        f = s / ds - is
+                        voigt = (1d0 - f) * b + f * a
+                    else
+                        s = y / (1d0 + y)
+                        pressure_space = 1d0 / (1d0 + x)
+                        is = int(s / ds)
+                        ip = int(pressure_space / dp)
+                        f = pressure_space / dp - ip
+                        if(ip < 400) then
+                            b = (1d0 - f) * v(ip, is) + f * v(ip + 1, is)
+                            a = (1d0 - f) * v(ip, is + 1) + f * v(ip + 1, is + 1)
+                        else
+                            b = v(ip, is)
+                            a = v(ip, is + 1)
+                        end if
+                        b = b * b
+                        a = a * a
+                        f = s / ds - is
+                        voigt = (1d0 - f) * b + f * a
+                    end if
+                end if
+            end if
+
+            return
+        end function voigt_from_data
+
+
+        recursive subroutine fft(x)
+            ! """
+            ! Calculate the Cooley-Tukey FFT functions.
+            !
+            ! input:
+            !   x: double complex vector of size 2^n
+            !
+            ! notes:
+            !   Source: https://rosettacode.org/wiki/Fast_Fourier_transform#Fortran
+            ! """
+            implicit none
+
+            complex(kind=dp), dimension(:), intent(inout) :: x
+            complex(kind=dp) ::&
+                t
+
+            integer ::&
+                i, &  ! index
+                n     ! size of array x
+
+            complex(kind=dp), dimension(:), allocatable ::&
+                even, &  ! even-number indexed values of x
+                odd      ! odd-number indexed values of x
+
+            n = size(x)
+
+            if(n <= 1) return
+
+            allocate(odd((n+1)/2))
+            allocate(even(n/2))
+
+            ! divide
+            odd(:) = x(1:n:2)
+            even(:) = x(2:n:2)
+
+            ! conquer
+            call fft(odd)
+            call fft(even)
+
+            ! combine
+            do i = 1, n/2
+               t = exp(cmplx(0.0_dp, -2.0_dp * pi * real(i-1, dp) / real(n, dp), kind=dp)) * even(i)
+               x(i) = odd(i) + t
+               x(i+n/2) = odd(i) - t
+            end do
+
+            deallocate(odd)
+            deallocate(even)
+
+        end subroutine fft
+
+
+        recursive subroutine quicksort(array)
+            ! """
+            ! Sort double precision numbers into ascending numerical order.
+            !
+            ! input:
+            !   array: double precision array to sort
+            !
+            ! output:
+            !   array: sorted double precision array
+            !
+            ! notes:
+            !   Author: t-nissie, some tweaks by 1AdAstra1
+            !   Source: https://gist.github.com/t-nissie/479f0f16966925fa29ea
+            ! """
+            implicit none
+
+            double precision, intent(inout), dimension(:) :: array
+
+            double precision :: &
+                x, &  ! pivot point
+                tmp   ! temporary value
+
+            integer :: &
+                first = 1, & ! index of the beginning of the array
+                last, &      ! index of the end of the array
+                insertion_size_threshold = 32
+
+            integer :: &
+                i, &
+                j
+
+            last = size(array, 1)
+
+            if(last < insertion_size_threshold) then  ! use insertion sort on small arrays
+                do i = 2, last
+                    tmp = array(i)
+
+                    do j = i - 1, 1, -1
+                        if(array(j) < tmp) then
+                            exit
+                        else
+                            array(j + 1) = array(j)
+                        end if
+                    end do
+
+                    array(j + 1) = tmp
+                end do
+            else
+                x = array((first+last)/2)
+                i = first
+                j = last
+
+                do
+                    do while (array(i) < x)
+                        i = i + 1
+                    end do
+
+                    do while (x < array(j))
+                        j = j - 1
+                    end do
+
+                    if (i >= j) then
+                        exit
+                    end if
+
+                    tmp = array(i)
+                    array(i) = array(j)
+                    array(j) = tmp
+
+                    i = i + 1
+                    j = j - 1
+                end do
+
+                if (first < i - 1) call quicksort(array(first:i-1))
+                if (j + 1 < last)  call quicksort(array(j+1:last))
+            end if
+        end subroutine quicksort
+
+
+        recursive subroutine quicksort_index(array, sorted_index)
+            ! """
+            ! Sort double precision numbers into ascending numerical order.
+            !
+            ! input:
+            !   array: double precision array to sort
+            !
+            ! output:
+            !   sorted_index: sorted index of the array
+            ! """
+            implicit none
+
+            double precision, intent(inout), dimension(:) :: array
+            integer, intent(out), dimension(:) :: sorted_index
+
+            integer :: &
+                i
+
+            do i = 1, size(array)
+                sorted_index(i) = i
+            end do
+
+            call partition(array, sorted_index)
+
+            contains
+                recursive subroutine partition(array, sorted_index)
+                    implicit none
+
+                    integer, intent(inout), dimension(:) :: sorted_index
+                    double precision, intent(inout), dimension(:) :: array
+
+                    integer :: &
+                        tmp_i
+
+                    double precision :: &
+                        x, &  ! pivot point
+                        tmp   ! temporary value
+
+                    integer :: &
+                        first = 1, & ! index of the beginning of the array
+                        last, &      ! index of the end of the array
+                        insertion_size_threshold = 128
+
+                    integer :: &
+                        i, &
+                        j
+
+                    last = size(array, 1)
+
+                    if(last < insertion_size_threshold) then  ! use insertion sort on small arrays
+                        do i = 2, last
+                            tmp = array(i)
+                            tmp_i = sorted_index(i)
+
+                            do j = i - 1, 1, -1
+                                if(array(j) < tmp) then
+                                    exit
+                                else
+                                    array(j + 1) = array(j)
+                                    sorted_index(j + 1) = sorted_index(j)
+                                end if
+                            end do
+
+                                array(j + 1) = tmp
+                                sorted_index(j + 1) = tmp_i
+                        end do
+                    else
+                        x = array((first + last) / 2)
+                        i = first
+                        j = last
+
+                        do
+                            do while (array(i) < x)
+                                i = i + 1
+                            end do
+
+                            do while (x < array(j))
+                                j = j - 1
+                            end do
+
+                            if (i >= j) then
+                                exit
+                            end if
+
+                            tmp = array(i)
+                            array(i) = array(j)
+                            array(j) = tmp
+
+                            tmp_i = sorted_index(i)
+                            sorted_index(i) = sorted_index(j)
+                            sorted_index(j) = tmp_i
+
+                            i = i + 1
+                            j = j - 1
+                        end do
+
+                        if (first < i - 1) call partition(array(first:i-1), sorted_index(first:i-1))
+                        if (j + 1 < last)  call partition(array(j+1:last), sorted_index(j+1:last))
+                    end if
+                end subroutine partition
+        end subroutine quicksort_index
+
+
+        subroutine reallocate_1Ddouble(double_array, new_size)
+            ! """
+            ! Change the size of a 1D double precision array.
+            ! :param double_array: array to change the size
+            ! :param new_size: new size of the array
+            ! :return double_array: resized array
+            ! """
+            implicit none
+
+            integer, intent(in) :: new_size
+            doubleprecision, intent(inout), dimension(:), allocatable :: double_array
+
+            integer :: &
+                size_array  ! size of the array
+            doubleprecision, dimension(:), allocatable :: &
+                tmp         ! temporary array
+
+            size_array = size(double_array)
+
+            call move_alloc(double_array, tmp)
+            allocate(double_array(new_size))
+
+            double_array(:) = 0d0
+
+            size_array = min(new_size, size_array)
+
+            double_array(:size_array) = tmp(:size_array)
+        end subroutine reallocate_1Ddouble
+
+
+        subroutine reallocate_1Dinteger(integer_array, new_size)
+            ! """
+            ! Change the size of a 1D integer array.
+            ! :param integer_array: array to change the size
+            ! :param new_size: new size of the array
+            ! :return integer_array: resized array
+            ! """
+            implicit none
+
+            integer, intent(in) :: new_size
+            integer, intent(inout), dimension(:), allocatable :: integer_array
+
+            integer :: &
+                size_array  ! size of the array
+            integer, dimension(:), allocatable :: &
+                tmp         ! temporary array
+
+            size_array = size(integer_array)
+
+            call move_alloc(integer_array, tmp)
+            allocate(integer_array(new_size))
+
+            integer_array(:) = 0
+
+            size_array = min(new_size, size_array)
+
+            integer_array(:size_array) = tmp(:size_array)
+        end subroutine reallocate_1Dinteger
+
+
+        subroutine reallocate_2Ddouble(double_array, new_size_1, new_size_2)
+            ! """
+            ! Change the size of a 2D double precision array.
+            ! :param double_array: array to change the size
+            ! :param new_size_1: new size of dimension 1 of the array
+            ! :param new_size_2: new size of dimension 2 the array
+            ! :return double_array: resized array
+            ! """
+            implicit none
+
+            integer, intent(in) :: new_size_1, new_size_2
+            doubleprecision, intent(inout), dimension(:, :), allocatable :: double_array
+
+            integer :: &
+                shape_array(2)  ! shape of the array
+
+            doubleprecision, dimension(:, :), allocatable :: &
+                tmp             ! temporary array
+
+            shape_array = shape(double_array)
+
+            call move_alloc(double_array, tmp)
+            allocate(double_array(new_size_1, new_size_2))
+
+            double_array(:, :) = 0d0
+
+            shape_array(1) = min(new_size_1, shape_array(1))
+            shape_array(2) = min(new_size_2, shape_array(2))
+
+            double_array(:shape_array(1), :shape_array(2)) = tmp(:shape_array(1), :shape_array(2))
+        end subroutine reallocate_2Ddouble
+
+
+        subroutine reallocate_2Dinteger(integer_array, new_size_1, new_size_2)
+            ! """
+            ! Change the size of a 2D integer array.
+            ! :param integer_array: array to change the size
+            ! :param new_size_1: new size of dimension 1 of the array
+            ! :param new_size_2: new size of dimension 2 the array
+            ! :return integer_array: resized array
+            ! """
+            implicit none
+
+            integer, intent(in) :: new_size_1, new_size_2
+            integer, intent(inout), dimension(:, :), allocatable :: integer_array
+
+            integer :: &
+                shape_array(2)  ! shape of the array
+
+            integer, dimension(:, :), allocatable :: &
+                tmp             ! temporary array
+
+            shape_array = shape(integer_array)
+
+            call move_alloc(integer_array, tmp)
+            allocate(integer_array(new_size_1, new_size_2))
+
+            integer_array(:, :) = 0
+
+            shape_array(1) = min(new_size_1, shape_array(1))
+            shape_array(2) = min(new_size_2, shape_array(2))
+
+            integer_array(:shape_array(1), :shape_array(2)) = tmp(:shape_array(1), :shape_array(2))
+        end subroutine reallocate_2Dinteger
+
+
+        subroutine reallocate_3Ddouble(double_array, new_size_1, new_size_2, new_size_3)
+            ! """
+            ! Change the size of a 2D double precision array.
+            ! :param double_array: array to change the size
+            ! :param new_size_1: new size of dimension 1 of the array
+            ! :param new_size_2: new size of dimension 2 of the array
+            ! :param new_size_3: new size of dimension 3 of the array
+            ! :return double_array: resized array
+            ! """
+            implicit none
+
+            integer, intent(in) :: new_size_1, new_size_2, new_size_3
+            doubleprecision, intent(inout), dimension(:, :, :), allocatable :: double_array
+
+            integer :: &
+                shape_array(3)  ! shape of the array
+
+            doubleprecision, dimension(:, :, :), allocatable :: &
+                tmp             ! temporary array
+
+            shape_array = shape(double_array)
+
+            call move_alloc(double_array, tmp)
+            allocate(double_array(new_size_1, new_size_2, new_size_3))
+
+            double_array(:, :, :) = 0d0
+
+            shape_array(1) = min(new_size_1, shape_array(1))
+            shape_array(2) = min(new_size_2, shape_array(2))
+            shape_array(3) = min(new_size_3, shape_array(3))
+
+            double_array(:shape_array(1), :shape_array(2), :shape_array(3)) = &
+                tmp(:shape_array(1), :shape_array(2), :shape_array(3))
+        end subroutine reallocate_3Ddouble
+
+
+        subroutine ifft(x)
+            ! """
+            ! Calculate the Cooley-Tukey IFFT functions.
+            !
+            ! input:
+            !   x: double complex vector of size 2^n
+            ! """
+            implicit none
+
+            complex(kind=dp), dimension(:), intent(inout)  :: x
+
+            call ifft_r(x)
+
+            x = x / real(size(x), dp)
+
+            contains
+                recursive subroutine ifft_r(x)
+                    implicit none
+
+                    complex(kind=dp), dimension(:), intent(inout) :: x
+
+                    integer ::&
+                        i, &  ! index
+                        n     ! size of array x
+
+                    complex(kind=dp) ::&
+                        t  ! intermediate value
+
+                    complex(kind=dp), dimension(:), allocatable ::&
+                        even, &  ! even-number indexed values of x
+                        odd      ! odd-number indexed values of x
+
+                    n=size(x)
+
+                    if(n <= 1) return
+
+                    allocate(odd((n + 1) / 2))
+                    allocate(even(n / 2))
+
+                    ! divide
+                    odd(:) = x(1:n:2)
+                    even(:) = x(2:n:2)
+
+                    ! conquer
+                    call ifft_r(odd)
+                    call ifft_r(even)
+
+                    ! combine
+                    do i = 1, n/2
+                       t = exp(cmplx(0.0_dp, 2.0_dp * pi * real(i - 1, dp) / real(n, dp), kind=dp)) * even(i)
+                       x(i) = odd(i) + t
+                       x(i+n/2) = odd(i) - t
+                    end do
+
+                    deallocate(odd)
+                    deallocate(even)
+
+                end subroutine ifft_r
+
+        end subroutine ifft
+
+
+        subroutine init_random_seed()
+            ! """
+            ! Initialize the random seed with a varying seed in order to ensure a different random number sequence for
+            ! each invocation of the program.
+            !
+            ! notes:
+            !   Using the random number generator file takes far too long to be useful.
+            !   Source: https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
+            ! """
+            use iso_fortran_env, only: int64
+
+            implicit none
+
+            integer :: &
+                i, &      ! index
+                n, &      ! size of random seed
+                dt(8), &  ! date and time
+                pid       ! pid
+
+            integer, allocatable :: &
+                seed(:)  ! random generato seed
+
+            integer(int64) :: &
+                t  ! (s) time
+
+            call random_seed(size = n)
+
+            allocate(seed(n))
+
+            ! XOR:ing the current time and pid. The PID is useful in case one launches multiple instances of the same
+            ! program in parallel.
+            call system_clock(t)
+
+            if (t == 0) then
+                call date_and_time(values=dt)
+
+                t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
+                   + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
+                   + dt(3) * 24_int64 * 60 * 60 * 1000 &
+                   + dt(5) * 60 * 60 * 1000 &
+                   + dt(6) * 60 * 1000 + dt(7) * 1000 &
+                   + dt(8)
+            end if
+
+            pid = getpid()
+            t = ieor(t, int(pid, kind(t)))
+
+            do i = 1, n
+                seed(i) = lcg(t)
+            end do
+
+            call random_seed(put=seed)
+
+            contains
+                function lcg(s)
+                    integer :: lcg
+                    integer(int64) :: s
+
+                    if (s == 0) then
+                        s = 104729
+                    else
+                        s = mod(s, 4294967296_int64)
+                    end if
+                        s = mod(s * 279470273_int64, 4294967291_int64)
+                        lcg = int(mod(s, int(huge(0), int64)), kind(0))
+                end function lcg
+        end subroutine init_random_seed
+
+
+        subroutine matinv(a, n)
+            ! """
+            ! Inverse a 2D matrix of dimension (n,n)
+            ! The original matrix a(n,n) will be destroyed during the calculation
+            ! Method: Based on Doolittle LU factorization for Ax=b
+            ! Source: http://ww2.odu.edu/~agodunov/computing/programs/book2/Ch06/Inverse.f90
+            ! by Alex G. December 2009
+            !
+            ! inputs:
+            !   a(n,n): array of coefficients for matrix A
+            !   n: dimension
+            !
+            ! output:
+            !   c(n,n): inverse matrix of A
+            ! """
+            implicit none
+
+            integer, intent(in) :: n
+            double precision, intent(inout) :: a(n,n)
+
+            double precision :: &
+                l(n,n), &  ! lower triangular matrix
+                U(n,n), &  ! upper triangular matrix
+                b(n), &    ! intermediate value
+                c(n,n), &  ! temporary matrix
+                d(n), &    ! intermediate value
+                x(n)       ! intermediate value
+
+            double precision :: &
+                coeff  ! intermediate value
+
+            integer :: &
+                i, &  ! index
+                j, &  ! index
+                k     ! index
+
+            ! Step 0: initialization for matrices l and U and b
+            l(:, :) = 0d0
+            U(:, :) = 0d0
+            b = 0d0
+            c(:, :) = 0d0
+
+            ! Step 1: forward elimination
+            do k = 1, n - 1
+                do i = k + 1, n
+                    coeff = a(i, k) / a(k, k)
+
+                    l(i, k) = coeff
+
+                    do j = k + 1, n
+                        a(i, j) = a(i, j) - coeff * a(k, j)
+                    end do
+                end do
+            end do
+
+            ! Step 2: prepare l and U matrices
+            ! l matrix is a matrix of the elimination coefficient + the diagonal elements are 1.0
+            do i = 1, n
+                l(i, i) = 1.0
+            end do
+
+            ! U matrix is the upper triangular part of A
+            do j = 1, n
+                do i = 1, j
+                    U(i, j) = a(i, j)
+                end do
+            end do
+
+            ! Step 3: compute columns of the inverse matrix C
+            do k = 1, n
+                b(k) = 1.0
+                d(1) = b(1)
+
+                ! Step 3a: Solve Ld=b using the forward substitution
+                do i = 2, n
+                    d(i) = b(i)
+
+                    do j=1, i - 1
+                        d(i) = d(i) - l(i, j) * d(j)
+                    end do
+                end do
+
+                ! Step 3b: Solve Ux=d using the back substitution
+                x(n) = d(n) / U(n,n)
+
+                do i = n - 1, 1, -1
+                    x(i) = d(i)
+
+                    do j = n, i + 1, -1
+                        x(i) = x(i) - U(i, j) * x(j)
+                    end do
+
+                    x(i) = x(i) / u(i, i)
+                end do
+
+                ! Step 3c: fill the solutions x(n) into column k of C
+                do i = 1, n
+                    c(i, k) = x(i)
+                end do
+
+                b(k) = 0d0
+            end do
+
+            a(:, :) = c(:, :)
+        end subroutine matinv
+end module math
Index: /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/optpropgen.f90
===================================================================
--- /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/optpropgen.f90	(revision 3240)
+++ /trunk/LMDZ.GENERIC/utilities/aerosols_utilities/optpropgen.f90	(revision 3240)
@@ -0,0 +1,367 @@
+!     -------------------------------------------------------------
+!     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('')//"&
+                &<https://gitlab.obspm.fr/dblain/exorem>"
+
+            write(*, '(A)') help_str
+        end subroutine print_help
+END
