source: TOOLS/CMIP6_FORCING/AER_OPTICS/refr_ind_sulfate_MACC_v2.f @ 4226

      PROGRAM ammonium_sulfate
!--based on Tang and Munkelwitz JGR 1994 and Tang JGR 1997 
      IMPLICIT NONE 
      REAL xms  !--solute mass fraction
      REAL aw !-water activity (RH here) 
      REAL molality !-molality
      REAL rhop !-dry droplet density
      REAL nr !--real refractive index
      REAL RR, Rwater, Ramm !--molar refraction 
      PARAMETER (Rwater=3.717, Ramm=23.50) 
      REAL VV !--molal volume
      REAL Mwater, Mamm 
      PARAMETER (Mwater=18.02, Mamm=132.14)
      REAL y !--mole fraction
      PARAMETER (rhop=1.769) !--g cm-3
      REAL B0, B1, B2, B3, B4 
      PARAMETER (B0=110.65495, B1=-367.59197,B2=504.62934) 
      PARAMETER (B3=-315.43839,B4=67.70824)
      REAL C1, C2, C3, C4 
      PARAMETER (C1=-2.715E-3,C2=3.113E-5,C3=-2.336E-6,C4=1.412E-8)
      REAL A0, A1, A2, A3
      PARAMETER (A0=0.9971, A1=5.92E-3, A2=-5.036E-6, A3=1.024E-8)
c-------------------------
      INTEGER nbre_rh,IRH,wve
      INTEGER nbre_dry,ii
      PARAMETER(nbre_rh=12,nbre_dry=4)
      REAL DELTA
      REAL RH_tab(nbre_rh)
      DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
      INTEGER Nwvmax
      PARAMETER(Nwvmax=100)
c
c--wavelength in m
      REAL lambda_min, lambda_max, lambda(Nwvmax)
      PARAMETER (lambda_min=0.2E-6, lambda_max=5.0E-6)

      REAL n_rr(nbre_rh,Nwvmax)
      REAL n_ii(nbre_rh,Nwvmax)
      REAL ratio(nbre_rh)      !--diameter growth factor
      REAL rhod(nbre_rh)       !--density
      REAL n_r_exact(nbre_rh)  !--refractive index
c
      DO wve=1,Nwvmax
        lambda(wve)=
     .  lambda_min+FLOAT(wve-1)*(lambda_max-lambda_min)/FLOAT(Nwvmax-1)
      ENDDO 
c
c--xms is the varying variable as per Tang and Munkelwitz 1994
      OPEN (unit=4,file='ri_sul_v2')
      OPEN (unit=10,file='growth_sul_v2')
c
      DO IRH=1,nbre_rh
c
      IF (IRH.LE.nbre_dry) THEN 
c
      ratio(IRH)=1.0
      rhod(IRH)=rhop
      n_r_exact(IRH)=1.521
c
      ELSE 
c
      aw=rh_tab(IRH)/100.   !--water activity
      molality=B0+aw*(B1+aw*(B2+aw*(B3+aw*B4)))   !--molality
      xms=(molality*Mamm/1000.)/(molality*Mamm/1000.+1.)*100. !--solute mass fraction
c
      rhod(IRH)=A0+xms*(A1+xms*(A2+xms*A3))       !--density
      ratio(IRH)=(rhop/rhod(IRH))**(1./3.) * (xms/100.)**(-1./3.)  !--growth factor
      y=Mwater*xms/100./(Mamm*(1.-xms/100.)+xms/100.*Mwater)
      VV=1./rhod(IRH)*(Mwater+(Mamm-Mwater)*y)
      RR=Rwater+(Ramm-Rwater)*y
      nr=(2.*RR/VV+1.)/(1.-RR/VV)
      n_r_exact(IRH)=sqrt(nr)
c
      ENDIF
c
      ENDDO
c
      write(10,20)'RH', (rh_tab(IRH),IRH=1,nbre_rh)
      write(10,10)"growth factor",(ratio(IRH),IRH=1,nbre_rh)
      write(10,10)"density factor",(rhod(IRH)/rhod(1),IRH=1,nbre_rh)
      write(10,10)'rho @ RH', (rhod(IRH),IRH=1,nbre_rh)
      write(10,10)'nr @ RH', (n_r_exact(IRH),IRH=1,nbre_rh)
10    FORMAT(A14, F6.3,11(',',F6.3))
20    FORMAT(A14, F6.1,11(',',F6.1))

c
c-----Computation of n_r at different wavelengths using
c     n_r(lambda)=n_r(0.589)-0.03(lambda-0.589) from 
c     Palmer and Williams (1975) and  Kent et al. (1983).
c     Above computed n_r_exact is assumed to be at 0.589 nm
c
      DO IRH=1,nbre_rh
      DO wve=1,Nwvmax
c
      n_rr(IRH,wve)=n_r_exact(irh)-0.03*(lambda(wve)*1.E6-0.589)
      n_ii(IRH,wve)=0.0
c
      WRITE(4,*) RH_tab(IRH), lambda(wve), n_rr(IRH,wve), n_ii(IRH,wve)
c
      ENDDO
      ENDDO
c
63    FORMAT(F10.6,E13.6)
c
      END