CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine refr_ind(temperature,frequency,refractive_index)

c	A subroutine lifted from STB's monodispersion.F code for calculating the
c	refractive index of ice from the tabulated data of Warren (1984).

	implicit none
	
	double precision temperature, frequency
      complex refractive_index
      complex warren_ri(11,3)
      real ri_real,ri_imag,ri_lo_temp_real,ri_lo_temp_imag,
     &  ri_hi_temp_real,ri_hi_temp_imag
	integer frequency_band
Cf2py intent(in) temperature,frequency
Cf2py intent(out) refractive_index
c     Set up the fixed grid of refractive index data for ice from
c     Warren's data. The first column holds the fixed frequencies, and
c     the first row holds the fixed temperatures.

      data warren_ri/(0.0,0.0),(7.531e11,0.0),(6.711e11,0.0),
     &  (5.981e11,0.0),(5.332e11,0.0),(4.751e11,0.0),(3.774e11,0.0),
     &  (2.998e11,0.0),(2.381e11,0.0),(1.199e11,0.0),(5.996e10,0.0),
     &  (-20.0,0.0),(1.7860,-1.25e-2),(1.7843,-1.15e-2),
     &  (1.7832,-1.06e-2),(1.7825,-9.77e-3),(1.7820,-9.01e-3),
     &  (1.7816,-7.66e-3),(1.7814,-6.52e-3),(1.7816,-5.54e-3),
     &  (1.7822,-3.42e-3),(1.7831,-2.1e-3),(-60.0,0.0),
     &  (1.7860,-7.34e-3),(1.7843,-6.4e-3),(1.7832,-5.6e-3),
     &  (1.7825,-5.0e-3),(1.7820,-4.52e-3),(1.7815,-3.68e-3),
     &  (1.7807,-2.99e-3),(1.7801,-2.49e-3),(1.7789,-1.55e-3),
     &  (1.7779,-9.61e-4)/

      
c       Calculate the refractive index by interpolation from a fixed
c       grid. 

        if((temperature.gt.-20.0d0).or.(temperature.lt.-100.0d0))then
          print *,'Temperature (',temperature,') out of range.'
          goto 55500
        endif
        if((frequency.gt.7.5d11).or.(frequency.lt.6.0d10))
     &    then
          print *,'Frequency (',frequency,') out of range.'
          goto 55500
        endif

        frequency_band=2
100     continue
          if(frequency.gt.dble(warren_ri(frequency_band+1,1)))
     &      goto 200
          frequency_band=frequency_band+1
          goto 100
200     continue

c       Interpolate logarithmically. The imaginary part is separated out
c       for this, because I am not sure whether logs and exps of complex
c       numbers will work OK. The real part can be linearly
c       interpolated. Interpolate firstly with respect to frequency.

        ri_hi_temp_real=real(warren_ri(frequency_band,2))
     &    +(frequency-real(warren_ri(frequency_band,1)))
     &    *(real(warren_ri(frequency_band+1,2))
     &      -real(warren_ri(frequency_band,2)))
     &    /(real(warren_ri(frequency_band+1,1))
     &      -real(warren_ri(frequency_band,1)))

        ri_hi_temp_imag=aimag(warren_ri(frequency_band,2))
     &    *exp( log(frequency/real(warren_ri(frequency_band,1)))
     &    *log(aimag(warren_ri(frequency_band+1,2))
     &      /aimag(warren_ri(frequency_band,2)))
     &    /log(real(warren_ri(frequency_band+1,1))
     &      /real(warren_ri(frequency_band,1)))
     &    )

        ri_lo_temp_real=real(warren_ri(frequency_band,3))
     &    +(frequency-real(warren_ri(frequency_band,1)))   
     &    *(real(warren_ri(frequency_band+1,3))
     &      -real(warren_ri(frequency_band,3)))
     &    /(real(warren_ri(frequency_band+1,1))
     &      -real(warren_ri(frequency_band,1)))
        
        ri_lo_temp_imag=aimag(warren_ri(frequency_band,3))
     &    *exp( log(frequency/real(warren_ri(frequency_band,1)))
     &    *log(aimag(warren_ri(frequency_band+1,3))
     &      /aimag(warren_ri(frequency_band,3)))
     &    /log(real(warren_ri(frequency_band+1,1))
     &      /real(warren_ri(frequency_band,1)))
     &    )

c       Now the interpolation w.r.t. temperature. If the temperature is
c       below -60C this will be an extrapolation, but the code is the
c       same.

        ri_real=ri_hi_temp_real
     &    +(temperature-real(warren_ri(1,2)))
     &    *(ri_lo_temp_real-ri_hi_temp_real)
     &    /(real(warren_ri(1,3))-real(warren_ri(1,2)))

        ri_imag=ri_hi_temp_imag
     &    *exp( log(temperature/real(warren_ri(1,2)))
     &    *log(ri_lo_temp_imag/ri_hi_temp_imag)
     &    /log(real(warren_ri(1,3))/real(warren_ri(1,2)))
     &    )

        refractive_index=cmplx(ri_real,ri_imag)

!        print *,'Refractive Index (Frequency =',frequency,') = '
!     & ,refractive_index
55500	continue	  
	return

	end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!