        PROGRAM mie
        IMPLICIT NONE
C
C-------Mie computations for a size distribution 
C       of homogeneous spheres.
c
C==========================================================
C--Ref : Toon and Ackerman, Applied Optics, 1981
C        Stephens, CSIRO, 1979
C to be compiled with double precision option
C AUTHOR: Olivier Boucher
C
C-------SIZE distribution properties----------------
C--sigma_g : geometric standard deviation 
C--r_0     : geometric number mean radius (um)/modal radius
C--Ntot    : total concentration in m-3 
c
        REAL rmin_dry, rmax_dry,r_0_dry   !----integral bounds in  m
        REAL rmin, rmax    !----integral bounds in  m
        REAL sigma_g, r_0, Ntot
        PARAMETER (rmin_dry=0.002E-6, rmax_dry=20.0E-6)
        PARAMETER (r_0_dry=0.0355E-6,sigma_g=2.0)
        PARAMETER (Ntot=1.0)
        INTEGER Nwv, Nwvmax
        REAL min_wv, increment_wv
        PARAMETER (Nwvmax=38, min_wv=0.3E-6, increment_wv=0.1E-6)
c
        REAL volume,surface,rgv,rgn,rho
        PARAMETER (rho=1.769E3)  !--dry density kg/m3
c
c-------------------------------------
c
        COMPLEX m          !----refractive index m=n_r-i*n_i
        INTEGER Nmax,Nstart !--number of iterations
        COMPLEX k2, k3, z1, z2
        COMPLEX u1,u5,u6,u8
        COMPLEX a(1:21000), b(1:21000)
        COMPLEX I 
        INTEGER n  !--loop index
        REAL pi, nnn
        REAL nombre,rayon,dry_mass
        COMPLEX nn
        REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
        REAL x  !--size parameter
        REAL r  !--radius
        REAL sigma_sca, sigma_ext, sigma_abs
        REAL omegatot,  gtot !--averaged parameters
        COMPLEX ksiz2(-1:21000), psiz2(1:21000)
        COMPLEX nu1z1(1:21010), nu1z2(1:21010)
        COMPLEX nu3z2(0:21000)
        REAL number, deltar
        INTEGER bin, Nbin
        PARAMETER (Nbin=2000)
c
c---------- RH growth parameters----------------
c
        INTEGER rh_int, irh
        PARAMETER(rh_int=12)
        REAL RH_tab(rh_int)
        DATA RH_tab/0.,10.,20,30.,40.,50.,60.,70.,80.,85.,90.,95./
        REAL rh_growth(rh_int)
c
        DATA rh_growth/1.000,1.00,1.00,1.0,1.169,
     .     1.220,1.282,1.363,1.485,1.581,1.732,2.085/
c
c----------refractive index are RH dependent and wavelength independent
c
        REAL ref_ind(rh_int)
        DATA ref_ind/1.53,1.53,1.53,1.53,1.44,1.43,1.42,
     .               1.41,1.40,1.39,1.38,1.36/
c
c------optical references
c
        REAL lambda_int(1:Nwvmax)
        REAL n_r(1:Nwvmax)
        REAL n_i(1:Nwvmax)
c
        REAL final_k(1:Nwvmax)
        REAL final_g(1:Nwvmax)
        REAL final_w(1:Nwvmax)

c
c-------RH dependent RI from file
c
        DO Nwv=1, Nwvmax
        lambda_int(Nwv)=min_wv+FLOAT(Nwv-1)*increment_wv
        ENDDO
c
        DO irh=1,rh_int        !---------LOOP OVER RH
c
        r_0=r_0_dry*rh_growth(irh)
        rmin=rmin_dry*rh_growth(irh)
        rmax=rmax_dry*rh_growth(irh)
c
c----------Into the MIE code-----------
c
        DO Nwv=1, Nwvmax

        n_r(Nwv)=ref_ind(irh)
        n_i(Nwv)=0.00
c
        m=CMPLX(n_r(Nwv),-n_i(Nwv))
c
        pi=4.*ATAN(1.)
c
        I=CMPLX(0.,1.)
c
          sigma_sca=0.0
          sigma_ext=0.0
          sigma_abs=0.0
          gtot=0.0
          omegatot=0.0
          dry_mass = 0.0
          volume=0.0
          surface=0.0
          nombre= 0.0
          rayon=0.0
          rgv=0.0
          rgn=0.0
c
        DO bin=0, Nbin !---loop on size bins
 
        r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
        x=2.*pi*r/lambda_int(Nwv)
        deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin))

        number=0
        number=Ntot/SQRT(2.*pi)/log(sigma_g)*
     .              exp(-0.5*(log(r/r_0)/log(sigma_g))**2)
        nombre=nombre+  number*deltar     !--m-3
        rayon =rayon +r*number*deltar
        dry_mass=dry_mass+4./3.*pi*((r/rh_growth(irh))**3)*number*
     .                     deltar*rho*1.E3          !--g/m3
        volume=volume+4./3.*pi*(r**3)*number*deltar
        surface=surface+4.*pi*r**2*number*deltar
        rgv=rgv+log(r)*r**3.*number*deltar 
        rgn=rgn+log(r)*number*deltar

        k2=m
        k3=CMPLX(1.0,0.0)

        z2=CMPLX(x,0.0)
        z1=m*z2
 
        IF (0.0.LE.x.AND.x.LE.8.) THEN
        Nmax=INT(x+4*x**(1./3.)+1.)+2
        ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
        ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
        Nmax=INT(x+4*x**(1./3.)+2.)+1
        ELSE
        WRITE(10,*) 'x out of bound, x=', x
        STOP
        ENDIF

        Nstart=Nmax+10

C-----------loop for nu1z1, nu1z2

       nu1z1(Nstart)=CMPLX(0.0,0.0)
       nu1z2(Nstart)=CMPLX(0.0,0.0)
       DO n=Nstart-1, 1 , -1
       nn=CMPLX(FLOAT(n),0.0)
       nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) 
       nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
       ENDDO

C------------loop for nu3z2

       nu3z2(0)=-I
       DO n=1, Nmax
       nn=CMPLX(FLOAT(n),0.0)
       nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) 
       ENDDO

C-----------loop for psiz2 and ksiz2 (z2)
       ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
       ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
       DO n=1,Nmax
       nn=CMPLX(FLOAT(n),0.0)
       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
       ENDDO

C-----------loop for a(n) and b(n)
    
       DO n=1, Nmax
        u1=k3*nu1z1(n) - k2*nu1z2(n)
        u5=k3*nu1z1(n) - k2*nu3z2(n)
        u6=k2*nu1z1(n) - k3*nu1z2(n)
        u8=k2*nu1z1(n) - k3*nu3z2(n)
        a(n)=psiz2(n)/ksiz2(n) * u1/u5
        b(n)=psiz2(n)/ksiz2(n) * u6/u8
       ENDDO

C-----------------final loop--------------
        Q_ext=0.0
        Q_sca=0.0
        g=0.0
        DO n=Nmax-1,1,-1
          nnn=FLOAT(n)
          Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) 
          Q_sca=Q_sca+ (2.*nnn+1.) * 
     .               REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
          g=g + nnn*(nnn+2.)/(nnn+1.) * 
     .       REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  + 
     .          (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
        ENDDO
        Q_ext=2./x**2 * Q_ext
        Q_sca=2./x**2 * Q_sca
        Q_abs=Q_ext-Q_sca
        IF (AIMAG(m).EQ.0.0) Q_abs=0.0
        omega=Q_sca/Q_ext
        g=g*4./x**2/Q_sca
c
        sigma_sca=sigma_sca+r**2*Q_sca*number*deltar
        sigma_abs=sigma_abs+r**2*Q_abs*number*deltar
        sigma_ext=sigma_ext+r**2*Q_ext*number*deltar
        omegatot=omegatot+r**2*Q_ext*omega*number*deltar
        gtot    =gtot+r**2*Q_sca*g*number*deltar 
c
        ENDDO   !---bin
C------------------------------------------------------------------

        rgv=exp(rgv/(volume*3./4./pi))
        rgn=exp(rgn/nombre)
        sigma_sca=pi*sigma_sca
        sigma_abs=pi*sigma_abs
        sigma_ext=pi*sigma_ext
        gtot=pi*gtot/sigma_sca
        omegatot=pi*omegatot/sigma_ext
        rayon=rayon/nombre
c
        final_g(Nwv)=gtot
        final_w(Nwv)=omegatot
        final_k(Nwv)=sigma_ext/dry_mass*(132./96.)
c
c--------------writing to output--------------
c
        WRITE(*,*) 'wavelength=',lambda_int(Nwv),' rh=',RH_tab(irh) 
        WRITE(*,*) 'k w g=',final_k(Nwv), final_w(Nwv), final_g(Nwv)
c
        ENDDO  !--loop on wavelength
c
        ENDDO          !-------Loop on rh_int
c        
        STOP
        END
