program fire
      implicit none
      real*8 wmin,wmax,c,l,v0,lexc,lambda,epsc,depsc,epscp,depscp,
     1T,hw,dd,v,ee,iroa,ft,iroaicp,iroadcpi,iroadcpii
      integer*4 np,ne,nd,nt,i,nl,ntpm,ntpp,npc,ntc
      real*8,allocatable::eps(:),deps(:),xe(:),xd(:),w(:),i0(:),doc(:),
     1ram(:),roa(:),x(:),dcpi(:),ramdcpi(:),roaicp(:),roadcpi(:),
     1ramdcpII(:),roadcpii(:),dcpII(:)
      logical LG,lwr
      write(6,600)
600   format(
     1' "Ring of fire" solvent ROA',/,/,
     1'     Input: EPS.PRN  ECD    dye, l mol-1 cm-1',/,
     1'            DEPS.PRN UV-ABS dye, l mol-1 cm-1',/,
     1'            ROA.TAB         solvent',/,
     1'            IRR-IRL.PRN IRR+IRL.PRN alternative to ROA.TAB',/,
     1'            FIRE.OPT        options',/,/,
     1'    Output: ram,prn         Raman solvent',/,
     1'            ramdcpi,prn    ',/,
     1'            ramdcpii,prn    ',/,
     1'            roa,prn         ROA solvent',/,
     1'            roaicp,prn   ',/,
     1'            roadcpi,prn   ',/,
     1'            roadcpii,prn   ',/)
      call readopt(wmin,wmax,c,l,np,v0,lexc,T,hw,LG,lwr)
      ne=nl('EPS.PRN')
      nd=nl('DEPS.PRN')
      nt=nl('ROA.TAB')-4
      ntpm=nl('IRR-IRL.PRN')
      ntpp=nl('IRR+IRL.PRN')
      allocate(eps(ne),xe(ne),deps(nd),xd(nd))
      call rs(ne, eps,xe, 'EPS.PRN')
      call rs(nd,deps,xd,'DEPS.PRN')

c     excitation:
c     closest point and value in el spectra:
      call rc( epsc,lexc,ne, eps,xe)
      call rc(depsc,lexc,nd,deps,xd)

      if(nt.gt.0)then
       write(6,6002)nt
6002   format(i5,' parameters in ROA.TAB')
       ntc=nt
       npc=np
       allocate(w(nt),i0(nt),doc(nt),dcpi(nt),dcpII(nt))
       call rt(nt,w,i0,dcpi,dcpII,doc,'ROA.TAB')
      else
       if(ntpm.gt.0)then
        write(6,6001)ntpm,ntpp
6001    format(2i5,' points in IRR-IRL, IRR+IRL')
        if(ntpm.ne.ntpp)
     1  call report('IRR-IRL,IRR+IRL different number of points')
        npc=ntpm
        ntc=npc
        allocate(w(ntc),i0(ntc),doc(ntc),dcpi(ntc),dcpII(ntc))
        dcpii=0.0d0
        dcpi=0.0d0
        call rtprn(ntc,w,i0,doc)
       else
        call report('no DOC data')
       endif
      endif
      allocate(x(npc),ram(npc),roa(npc),
     1ramdcpi(npc),roadcpi(npc),roaicp(npc),ramdcpII(npc),roadcpii(npc))
      call inisp(wmax,wmin,npc,x,roa,roaicp,roadcpi,roadcpii,ramdcpi,
     1ramdcpii,ram,LG,dd,hw)

      if(lwr)write(6,610)c,L
 610  format(' c = ',e10.3,' M, L -',e10.3,' cm')

      do 2 i=1,ntc
c     absolute wavenumber:

      v=v0-w(i)
c     wavelength
      lambda=1.0d7/v
c     closest point and value in el spectra:
      call rc( epscp,lambda,ne, eps,xe)
      call rc(depscp,lambda,nd,deps,xd)
      ee=(epsc+epscp)*c*l
      if(ee.gt.30.0d0)then
       ee=0.0d0
      else
       ee=exp(-ee)
      endif
      ft=(1.0d0-ee*(c*l*(epscp+epsc)+1.0d0))  /
     1 (2.0d0*(epscp+epsc)*(1.0d0-ee))
c     SCP
c           (de' +DOC de) (1 - exp(-(e+e')cL) (cL(e'+e) +1)
c           ----------------------------------------------- I0
c                     2(e+e')(1- exp(-(e+e')cL))
      iroa=(depscp+doc(i)*depsc)*ft*i0(i)
      if(lwr)write(6,609)i,w(i),lambda,epsc,epscp,depsc,depscp,doc(i),
     1ft,i0(i)
609   format(i6,f8.2,' cm-1',f7.2,' nm, e:',e11.3,' ep:',
     1e11.3,' de:',e11.3,' dep:',e11.3,' DOC:',e11.3,' ft:',
     1e11.3,' I:',e11.3)
c     ICP
c           (de +DOC de') (1 - exp(-(e+e')cL) (cL(e'+e) +1)
c           ----------------------------------------------- I0
c                     2(e+e')(1- exp(-(e+e')cL))
      iroaicp=(depsc+doc(i)*depscp)*ft*i0(i)
c     DCPI
c           (de' +    de) (1 - exp(-(e+e')cL) (cL(e'+e) +1)
c           ----------------------------------------------- I0
c                     2(e+e')(1- exp(-(e+e')cL)
      iroadcpi=(depscp+depsc      )*ft*dcpi(i)
c     DCPI
c           (de' -    de) (1 - exp(-(e+e')cL) (cL(e'+e) +1)
c           ----------------------------------------------- I0
c                     2(e+e')(1- exp(-(e+e')cL)
      iroadcpii=(depscp-depsc     )*ft*dcpii(i)
      if(nt.gt.0)then
c      spectral convolution:
       call spr(iroa     ,w(i),dd,x,np,roa,     T,LG)
       call spr(iroaicp  ,w(i),dd,x,np,roaicp  ,T,LG)
       call spr(iroadcpi ,w(i),dd,x,np,roadcpi ,T,LG)
       call spr(iroadcpii,w(i),dd,x,np,roadcpii,T,LG)
       call spr(dcpi(i)  ,w(i),dd,x,np,ramdcpi ,T,LG)
       call spr(dcpii(i) ,w(i),dd,x,np,ramdcpii,T,LG)
       call spr(i0(i)    ,w(i),dd,x,np,ram     ,T,LG)
      else
c      directly this point:
       roa(i)=iroa
       x(i)=w(i)
      endif
2     continue
      
      call ws('roa.prn',npc,x,roa)
      if(nt.gt.0)then
       call ws('ram.prn',npc,x,ram)
       call ws('ramdcpi.prn',npc,x,ramdcpi)
       call ws('ramdcpii.prn',npc,x,ramdcpii)
       call ws('roaicp.prn',npc,x,roaicp)
       call ws('roadcpi.prn',npc,x,roadcpi)
       call ws('roadcpii.prn',npc,x,roadcpii)
       endif

      end

      subroutine spr(i0,wo,dd,x,np,s,temp,LG)
      implicit none
      integer*4 np,i
      real*8 i0,wo,x(*),s(*),temp,dd,const,sfc,sf
      logical LG
      do 4 i=1,np
      sfc=sf(dd,LG,x(i),wo)
      const=i0/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
4     s(i)=s(i)+sfc*const 
      return
      end

c     ==================================================================
      function sf(d,g,x,x0)
      implicit none
      logical g
      real*8 pi,spi,sf,d,x,x0,dd
      pi=4.0d0*atan(1.0d0)
      spi=dsqrt(pi)
      dd=((x-x0)/d)**2
      if(g)then
       if(dd.lt.20.0d0)then
        sf=exp(-dd)/d/spi
       else
        sf=0.0d0
       endif
      else
       if(dd.lt.1.0d10)then
        sf=1.0d0/d/(dd+1.0d0)/pi
       else
        sf=0.0d0
       endif
      endif
      return
      end

      subroutine  rc(ec,l,ne,e,xe)
      implicit none
      integer*4 ne,i
      real*8 ec,l,e(*),xe(*)
      do 1 i=1,ne-1
      if((xe(i).le.l.and.xe(i+1).ge.l).or.
     1   (xe(i).ge.l.and.xe(i+1).le.l))then
       ec=e(i)
       return
      endif
 1    continue
      write(6,*)l
      call report('no point found')
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine ws(f,n,x,s)
      implicit none
      integer*4 n,i
      real*8 x(*),s(*)
      character*(*) f
      open(9,file=f)
      do 1 i=1,n
 1    write(9,900)x(i),s(i)
 900  format(f12.2,e14.4)
      close(9)
      write(6,*)f//' written'
      return
      end


      subroutine rt(n,w,r,p,II,d,f)
c     read TAB file
      implicit none
      integer*4 n,i,j
      real*8 w(*),r(*),d(*),p(*),II(*)
      character*(*) f
      open(9,file=f)
      do 3 i=1,3
  3   read(9,*)
      do 1 i=1,n
  1   read(9,*)w(i),w(i),(r(i),j=1,3),(d(i),j=1,6),p(i),II(i),II(i)
      close(9)
      return
      end 

      subroutine rtprn(n,w,r,d)
c     read PRN files
      implicit none
      integer*4 n,i
      real*8 w(*),r(*),d(*)
      open(91,file='IRR-IRL.PRN',status='old')
      open(92,file='IRR+IRL.PRN',status='old')
      open(93,file='DOC.PRN')
      do 1 i=1,n
      read(91,*)w(i),d(i)
      read(92,*)w(i),r(i)
      if(dabs(r(i)).gt.dabs(d(i)*1.0d-9))d(i)=d(i)/r(i)
1     write(93,930)w(i),d(i)
930   format(2e13.4)
      close(91)
      close(92)
      close(93)
      return
      end 

      subroutine rs(n,s,x,f)
      implicit none
      real*8 s(*),x(*)
      character*(*) f
      integer*4 n,i
      open(9,file=f)
      do 1 i=1,n
  1   read(9,*)x(i),s(i)
      close(9)
      return
      end 

      subroutine readopt(wmin,wmax,c,l,np,v0,lexc,T,hw,LG,lwr)
      implicit none
      real*8 wmin,wmax,c,l,v0,lexc,hw,T
      integer*4 np
      logical LG,lwr
      character*4 key
      np=3901
      wmin=100.0d0
      wmax=4000.0d0
      l=0.10d0
      lexc=532.0d0
      T=293.0d0
      hw=10.0d0
      LG=.false.
      lwr=.false.
      c=1.0d-4
c     np  number of points
c     wmin, wmax ... limit for the spectra
c     lexc  excitation wave length / nm
c     v0 wavenumber in cm-1
c     T /K temperature
c     LG Gaussian peaks
c     lwr extensive output
c     hw half width/cm-1
c     concentration mol/l
      open(9,file='FIRE.OPT')
1     read(9,90,end=99,err=99)key
  90  format(a4)
      if(key(1:2).eq.'NP')read(9,*)np
      if(key.eq.'WMIN')read(9,*)wmin
      if(key.eq.'WMAX')read(9,*)wmax
      if(key.eq.'PATH')read(9,*)l
      if(key.eq.'LEXC')read(9,*)lexc
      if(key.eq.'TEMP')read(9,*)T
      if(key.eq.'FWHH')read(9,*)hw
      if(key.eq.'GAUS')read(9,*)LG
      if(key.eq.'WRIT')read(9,*)lwr
      if(key.eq.'CONC')read(9,*)c
      goto 1
  99  close(9)
      v0=1.0d7/lexc
      return
      end

      function nl(f)
c     number of lines in a file
      implicit none
      character*(*) f
      integer*4 nl,n
      open(9,file=f)
      n=0
1     read(9,*,end=99,err=99)
      n=n+1
      goto 1
99    close(9)
      nl=n
      return
      end


      subroutine inisp(wmax,wmin,np,x,roa,roaicp,roadcpi,roadcpii,
     1ramdcpi,ramdcpii,ram,LG,dd,hw)
      implicit none
      integer*4 np,i
      real*8 wmax,wmin,dw,wx,x(*),roa(*),roaicp(*),roadcpi(*),
     1roadcpii(*),ramdcpi(*),ramdcpii(*),ram(*),dd,hw
      logical LG
      dw=(wmax-wmin)/dble(np-1)
      wx=wmin-dw
      do 1 i=1,np
c     wavenumber:
      wx=wx+dw
      x(i)=wx
      roa(i)=0.0d0
      roaicp(i)=0.0d0
      roadcpi(i)=0.0d0
      roadcpii(i)=0.0d0
      ramdcpi(i)=0.0d0
      ramdcpII(i)=0.0d0
1     ram(i)=0.0d0
      if(LG)then
       dd=hw/2.0d0/dsqrt(log(2.0d0))
      else
       dd=hw/2.0d0
      endif
      return
      end