module energylevs
        ! this is the module which stores the energy level quantum numbers
        ! and energies
      implicit none
      save
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: enlevs
      integer, allocatable, dimension(:) :: ij, isym4, inn
      integer, allocatable, dimension(:,:) :: iV, KJ
      integer  nlevs

      end module energylevs

!------------------------------

      program spect4

!  this program is designed to handle the output from programs dipole
!  (of the TRIATOM suite) and dipole3 (of the DVR3D suite)
!  for several runs and turn it into data sorted on the transition
!  frequency wif. From this, and input from triatom/dvr3d and rotlev runs,
!  the partition functions and integrated absorption coefficients at
!  the required temperatures.

!  the program structure consists essentially of four main parts:

!  1) the main program, spect, determines the amount of memory needed
!  for the job in hand, and calls subroutine gtmain, a machine
!  dependent routine, which allocates the memory needed. spect also
!  calls the other three main working routines.

!  2) subroutine sortsp. this routine sorts out data from program
!  dipole on ascending transition frequency and may print this out
!  if required.

!  the main input file is "itra", and contains the output from the
!  necessary dipole runs. by default itra is attached to stream 13.
!  itra is read in the main program to determine the size of the
!  problem. the main program also creates a scratch file, item (16),
!  to cut down the amount of data sorted in sortsp.

!  subroutine sortsp creates a scratch file ispe for calculations of
!  spectral intensities, if required. if this file has been saved from
!  a previous run, sortsp may be skipped. ispe defaults to stream 15.

!  3) subroutine pfcalc. this routine calculates the partition function
!  from energy levels calculated by programs triatom and rotlevd at the
!  temperature required for each run.

!  data necessary to calculate the partition function is given on
!  file "ilev". this is defaulted to stream 14.
!  the first eigenenergy on ilev must be the ground state for j=0,
!  (q=0). this is the zero energy of the problem. ilev is read in the
!  main program to determine the size of the problem and in subroutine
!  pfcal! to access the data.

!  4) subroutine spectm. this routine calculates boltzman weighted
!  intensities at the required temperatures and in the frequency ranges
!  determined by the input parameters given by input data lines 4-n.

!  spectm reads the sorted data from sortsp off scratch file ispe.
!  if required it writes out data for a suitable graphics package to
!  data stream iplot (default 20). one option is to produce a spectrum
!  of gaussian profiles by calling subroutine profil. users may wish to
!  substitute this routine with other profile types.
!  lines can also be written to a linelist file, stream ilist (default
!  36), if zlist is .true.


!  input data is supplied to the program on stream 5, as follows:

!  1) namelist prt: four logical parameters:

!                   zout(false) this should be true if the sorted
!                               line strengths are to be written
!                               to the lineprinter. zout is set to
!                               true automatically if zspe is false.

!                   zsort(true) if false subroutine sortsp is skipped.

!                   zspe (true) if false the program stops after sortsp.
!                               units of ispe are atomic units.

!                   zpfun(true) calculates the partition function
!                               from energy levels supplied from
!                               DVR3DRJZ and ROTLEV3/3B.
!                               if zpfun false, the partition function
!                               is set to q read in below. 
!                   ZBAND(false) Set this to true if a specific band is required
!                                Must then provide upper and lower
!                                vibrational quantum numbers.
!
!                   the default values of itra(13), ilev(14), ispe(15),
!                   and item(16) may also be reset on this namelist
!                   if required.

!                   gz(0.0)     absolute energy of the ground state
!                               in cm-1 (not used if zpfun=.true.).

!   If zsort=.true. the range of data to be sorted may be set by setting
!                   any or all of wsmin, wsmax, emin, emax, jmax, smin.

!                   emin(-1.d27) the minimum value of e" for which data
!                                is to be written out.

!                   emax(+1.d27) the maximum value of e" for which data
!                                is to be written out. if emax is not
!                                specified all values of e" above emin
!                                are considered.

!                   jmax(-1)    the maximum value of j" for which
!                               transition frequencies, et! are to be
!                               written out in subroutine sortspe.
!                               if jmax is not specified, all j" values
!                               are considered.

!                   smin(0.0d0) The minimum value of the square of the 
!                               transition dipole (D**2).
!                               
!                   smax        The maximum value of the square of the
!                               transition dipole (D**2)        
!
!                   wsmax       The maximum frequency of a transition (cm-1).
!
!                   wsmin       The minimum frequency of a transition (cm-1)

!  judicious use of the parameters on this namelist can enable certain types
!  of transitions to be separated - e.g. fundamental bands, hot bands etc.

!  2) a title, format 9a8.

!  3) ge, go. format 2f10.0   nuclear degeneracy factors for homonuclear
!             diatomic molecules for the j even and j odd states.
!             defaults of 1.0d0 supplied if left totally blank.
!             if the molecule does not contain a homonuclear diatomic
!             part, ge and go are ignored.

!  4) temp, xmin, wmin, wmax, dwl, q. format 6f10.0

!      temp       temperature in k;
!      xmin       the lowest relative intensity to be printed out;
!      wmin and wmax define the transition frequency range in cm-1.
!      dwl        the width of the gaussian profile for the lines (see below).
!      q          is the partition function (only used if ZPFUN=F).
!                 use of default (q=1.0) means the absolute absortpion 
!                  coefficients are not correct. 

!  5) namelist spe. this namelist is read in subroutine spectm.
!     emin1, emax1, emin2, emax2, tinte and jmax control the print out
!     of spectral intensities.

!     tinte set the value below which the intensities are printed out.

!     zemit (.false.) gives emissivities if .true.

!     in addition, there is a parameter zplot(false)  which creates an
!     ascii file suitable for use with a graphics package. this can
!     consist of stick spectra or a gaussian broadened spectrum if
!     zprof (.false.) is .true. at npoint points, it is in:

!     1) frequency in wavenumbers against intensity if zfreq (.true.)
!        is .true.

!     2) wavelength in microns against intensity if zfreq .false.

!     3) either of the above against einstein a-coefficient times
!        spin weighting in zeinst (.false.) is .true.
!        zeinst (.true.) sets zemit to .true.
!     in addition if zene is .true the energy levels and the rotational
!     levels are printed in the output file iplot.

!  if zprof is .true., stream idat takes data from spectm to profil.

!  The program makes use of the dpsort freeware routines, found at
!  http://netlib2.cs.utk.edu/
! 
!  Updated to f90 to use dynamic memory allocation and to preselect transitions
!  by GJH & JT 2001.

      use energylevs
      implicit double precision(a-h,o-y), logical(z)
!      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: enlevs
!      integer, allocatable, dimension(:) :: ij, isym4, inn
!      integer, allocatable, dimension(:,:) :: VV, KJ
      character(len=8) title(9)
      namelist/prt/ zout, zsort, zspe, zpfun, itra, ilev, ispe, item, &
                    wsmax,wsmin, emin, emax, jmax, smin, gz, zband
      common/logic/ zout, zsort, zspe, zpfun, zembed
      common/base/ ibase1,ibase2
      common/timing/itime0
!      common/energylevs/ enlevs,ij,isym4,inn, vv,kj, nenlevs

!      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ee1, ee2, ss
      data autocm/ 2.19474624d+05/, autode/ 2.5417662d0/
!     detosc converts from s(f-i) in debye**2 to seconds**(-1)
      data detosc/ 3.136186d-07/
      nlevs=221097 ! this is the total number of energy levels

!     set time zero for timer
      call SYSTEM_CLOCK(itime0,irate2,imax2)

!     preset namelist defaults

      zout= .false.
      zsort= .true.
      zspe= .true.
      zpfun= .true.
      zband= .false.
      itra= 13
      ilev= 14
      ispe= 15
      item= 16
      wsmin= 0.0d0
      wsmax= 1.0d6
      smin=0.0d0
      jmax=500
      emin= -1.0d27
      emax= +1.0d27
      gz=0.0d0

! NOTE IDIA IS HARD CODED HERE
      idia=-2


      write(6,*) "Program SPECTRA (version of January 2004)"

!     read in logical variables and title
      open(unit=5,file='spec.job',form='formatted')
      read(5,prt)
      read(5,100) title
100   format(9a8)
      write(6,202) title
202   format(/5x,9a8/)
      if (zband) read(5,*) ivu1, ivu2, ivu3, ivl1, ivl2, ivl3

! allocate memory of the energy levels
       allocate(enlevs(nlevs),ij(nlevs),isym4(nlevs),inn(nlevs))
       allocate(kj(nlevs,3), iv(nlevs,3))

       call read_energies

!     determine energy zero from levels file if available
      if (zpfun) then
       gz=enlevs(1)
      endif
      write(6,203) gz
203   format(5x,'Ground state energy set to',f14.5,' cm-1')
!     determine number of transitions to be processed and create item

      maxr=0
      nr= 0
      nmax1= 0
      nmax2= 0
      open(unit=ispe,form='unformatted')
      if(zsort) then
         write(6,1000) wsmin,wsmax,emin,emax,smin,jmax
 1000    format(/5x,'Sorting requested, list preselected using:', &
                /5x,'Minimum frequency,   WSMIN =',d12.5,' cm-1', &
                /5x,'Maximum frequency,   WSMAX =',d12.5,' cm-1', &
                /5x,'Minimum energy,       EMIN =',d12.5,' cm-1', &
                /5x,'Maximum energy,       EMAX =',d12.5,' cm-1', &
                /5x,'Minimum linestrength, SMIN =',d12.5,' D**2', &
                /5x,'Maximum rotations     JMAX =',i5)

         wsmi=wsmin/autocm
         wsma=wsmax/autocm
         emi=(emin+GZ)/autocm
         ema=(emax+GZ)/autocm
         smi=smin/autode**2

! begining of 1st select phase
!
! This phase reads in the data line by line and calculates 
! S, frequncy and other parameters then if criterior are met
! outputs to item.
!
      open(unit=itra,form='formatted')
      open(unit=item,form='unformatted')
      rewind itra

      maxr=0
      nr=0
10    read(itra,102,end=90) index1,index2,AA
         maxr=maxr+1
         EE1=enlevs(index1)
         EE2=enlevs(index2)

         if(zband) then
          if(iv(index1,1) .eq. ivu1 .and. iv(index1,2) .eq. ivu2 .and. &
             iv(index1,3) .eq. ivu3 .and. iv(index2,1) .eq. ivl1 .and. &
             iv(index2,2) .eq. ivl2 .and. iv(index2,3) .eq. ivl3) goto 15
 
          if(iv(index1,1) .eq. ivl1 .and. iv(index1,2) .eq. ivl2 .and. &
             iv(index1,3) .eq. ivl3 .and. iv(index2,1) .eq. ivu1 .and. &
             iv(index2,2) .eq. ivu2 .and. iv(index2,3) .eq. ivu3) goto 15
          goto 10
         endif

15       if (ee1.ge.emin .and. ee1.le.emax) then
          if (ee2.ge.emin .and. ee2.le.emax) then
           w= ee2 - ee1
           if (abs(w).ge.wsmin .and. abs(w).le.wsmax) then
            call isym_reverse(isym4(index1), ij(index1), ipar1, kmin1)
            call isym_reverse(isym4(index2), ij(index2), ipar2, kmin2)

           if (w.ge.0.0d0) then 
            ss=(dble(2*ij(index2) + 1)*aa)/abs(w*w*w*detosc)
            else
             ss=-(dble(2*ij(index1) + 1)*aa)/(w*w*w*detosc)
            endif
!--
            if (ss.ge.smin) then
             ss=ss/autode**2
             ee1=ee1/autocm
             ee2=ee2/autocm
             w=w/autocm
             nr=nr+1

             if (w.ge.0.0d0) then 
              write(item) ipar1,ij(index2),kmin2,inn(index2), &
                          ij(index1),kmin1,inn(index1),index2,index1, &
                          ee2,ee1,w,ss,1,1
             else
              write(item) ipar1,ij(index1),kmin1,inn(index1),&
                          ij(index2),kmin2,inn(index2),index1,index2,  &
                          ee1,ee2,-w,ss,1,1
             endif
            endif
           endif
          endif
         endif
       goto 10

90     close(itra)

         write(6,1010) maxr,itra,nr,item
 1010 format(/,i10,' transition records read from unit ITRA =',i3, &
             /,i10,' selected and      written to unit ITEM =',I3)
         if (nr.le.0) then
            write(6,900) 
  900       format(/'No lines selected: stop')
            stop
         endif
      else
         rewind ispe
         read(ispe)
13       read(ispe,end=93)
         nr= nr+1
         goto 13
93       continue
         write(6,1020) nr,ispe
 1020    format(/5x,'List already sorted:', &
                /i10,' transition records found on unit ISPE =',i3) 
      endif

      write (6,*) "call spmain"
      call spmain(ilev,ispe,nr,nmax1,nmax2,item,idia,gz)

      deallocate(enlevs,ij,isym4,inn, iV,kj)
      stop
102   format(2i7,1p,e10.3)
      end

! ------------------------------------------------------------------
      subroutine read_energies!(E,ij,isym4,inn,nen)
      use energylevs
      implicit double precision(a-h,o-y), logical(z)
!      common/energylevs/ enlevs,ij,isym4,inn, vvkj, nen
!      dimension E(NeN), ij(nen), inn(nen), isym4(nen),vv(nen,3),kj(nen,3)
      dimension kjt(3), ivt(3)

      ien=80
      icount=0

      open(unit=ien,form='formatted')

      do 10 i=1,nlevs
       read(ien,101,end=20) index,ijt, isymt, innt, Etemp,(iVt(j),j=1,3),&
                            (KJt(j),j=1,3)
       ij(index)=ijt
       isym4(index)=isymt
       inn(index)=innt
       enlevs(index)=Etemp

       do 30 j=1,3
        kj(index,j)=kjt(j)
        iv(index,j)=ivt(j)
30     continue

!       kj(index,1)=ij(index)

       icount=icount+1
10    continue
20    close(ien)

      write(6,*) icount," Energy levels read from stream ", ien

      return
101   format(i6,i3,i2,i5,f14.6,3i4,3i3)
      end
! ------------------------------------------------------------------
      subroutine isym_reverse(isym,j,ipart,kmin)
      implicit double precision(a-h,o-y), logical(z)

      if(isym .eq. 1) then
       ipart=0
       kmin=mod(J+1,2)
       return
      endif

      if(isym .eq. 2) then
        ipart=0
        kmin=mod(J,2)
       return
      endif

      if(isym .eq. 3) then
        ipart=1
        kmin=mod(j,2)
       return
      endif

      if(isym .eq. 4) then
        ipart=1
        kmin=mod(j+1,2)
       return
      endif

      write(6,*) "isym_reverse failed"
      write(6,*) "isym not between 1 and 4"
      stop

      return
      end
!-------------------------------------------------------------------

      subroutine spmain(ilev,ispe,nr,nmax1,nmax2,item,idia,gz)

!     this is the effective main program of spectra 
      implicit double precision(a-h,o-y), logical(z)
      common/logic/ zout, zsort, zspe, zpfun, zembed
      
      data dwl/0.0d0/,x1/1.0d0/,x0/0.0d0/

!     sort out the spectrum on frequencies

      if (zsort) call sortsp(nr, item, ispe)
      if (.not.zspe) stop
      rewind ispe
      read(ispe) zembed

!     read in the nuclear spin factors


      read(5,101) ge, go
101   format(6f10.0)
      if(ge.le.x0 .and. go.le.x0) then
        ge= x1
        go= x1
      endif
      write(6,205) ge, go
205   format(/5x,'Even spin factor = ',f6.3,'  odd spin factor = ',f6.3)

!     read in the conditions for spectra required.
! 
      read(5,101, end=92) temp, xmin, wmin, wmax, dwl, q

!     print out run conditions

      write(6,206) temp
206   format(/5x,'Temperature set to: ',f8.2,' K'//)
      if (xmin.ne.x0) then
         write(6,207) xmin
      else
         write(6,208)
      endif
207   format(5x,'Minimum relative intensity required = ',f8.6)
208   format(5x,'All transitions printed out')
      if(wmax.ne.x0) then
         write(6,209) wmin, wmax, dwl
      else
         write(6,210) wmin, dwl
      endif
209   format(5x,'Frequency range from ',f10.3,' cm-1 to ',f10.3,' cm-1',/ &
             5x,'profile half width',f10.3/)
210   format(5x,'Frequency range from ',f10.3,' cm-1 to maximum',/ &
             5x,'profile half width',f10.3/)

!     calculate the partition function

      if (zpfun) then
         nlev=max(nmax1,nmax2)
!         call pfcalc(temp,emax,qerr,ge,go,nlev,q,ilev,idia,gz)
         call pfcalc(temp,ge,go,q,idia,gz)
         write(6,211) q
211      format(/5x,'Partition function Q =',d12.5)
      else
         if (q .le. x0) q=x1
         write(6,2061) q
2061    format(/5x,' Partition function set to Q =',d12.5)
      endif

!     calculate integrated absorption coefficients
      jdia= abs(idia)
      call spectm(xmin, wmin, wmax, dwl, &
                  ge, go, temp, q, nr, jdia, ispe, gz)

92    continue
       write(6,212)
212   format(/10x,'Job completed successfully')
      call timer
!      stop
      return
      end

! ------------------------------------------------------------------

      subroutine sortsp(nr, item, ispe)

!     this subroutine sorts the dipole data on ascending frequencies.
!     data printed out is in cm-1, debye**2 and sec-1.
!     data written to ispe for spectm is in atomic units.

      implicit double precision(a-h,o-y), logical(z)
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: da
      integer, allocatable, dimension(:,:) :: ia
      integer, allocatable, dimension(:) :: iperm


      common/logic/ zout, zsort, zspe, zpfun, zembed

      fmem = (64*nr)/1048576.0d0
      write(6,300) fmem
300   format(/5x,'Memory required to sort transitions',f10.3,' MB')
      allocate(ia(nr,9), da(nr,4), iperm(nr))

      if(.not.zspe) zout= .true.
      write(ispe) zembed
      ir= 0

!     read in of data from item.
!     input is in atomic units.

      rewind item
10    read(item,end=92,err=998) ipar1, j2,kmin2,ie2,j1,kmin1,ie1, &
                        index2,index1, e2, e1, w, s,ibase2,ibase1

      ir= ir + 1

      ia(ir,1)= ipar1
      ia(ir,2)= j2
      ia(ir,3)= kmin2
      ia(ir,4)= ie2 + ibase2 - 1
      ia(ir,5)= j1
      ia(ir,6)= kmin1
      ia(ir,7)= ie1 + ibase1 - 1
      ia(ir,8)=index2
      ia(ir,9)=index1

      da(ir,1)= e2
      da(ir,2)= e1
      da(ir,3)= w
      da(ir,4)= s
      goto 10

92    continue

!     now sort on frequencies

      ier = 0

      call dpsort(da(1,3), nr, iperm, 1, ier)
      if(ier .ne. 0) then
         write(6,*) " ERROR returned by DPSORT", ier
         stop
      endif

      do 698 ir=1,9
 
      call ipperm(ia(1,ir), nr, iperm, ier)
      if(ier .ne. 0) write(6,*) " ERROR returned by ipperm ", ier
      if(ir .le. 4) then
         call dpperm(da(1,ir), nr, iperm, ier)
         if(ier .ne. 0) write(6,*) " ERROR returned by dpperm ", ier
      endif   

  698 continue

!     write out transitions for subroutine spectm to ispe.
!     output to ispe is in atomic units.

!     print out results to lineprinter if requested.
!     output to lineprinter is in wavenumbers and debye.

      if (zout) write(6,202)
202   format(//'ipar   j2  p2  i2   j1  p1  i1 index2 index1       e2        e1',9x, &
         'freq       s(f-i)'/)
      i20= 0
      do 4 ir= 1,nr
      write(ispe) (ia(ir,ic), ic=1,9), (da(ir,ic), ic=1,4)
      if (zout) write(6,203) (ia(ir,ic), ic=1,9), (da(ir,ic), ic=1,4)
4     continue
203   format(i3,2x,3(1x,i3),1x,3(1x,i3),2i7,3(2x,f9.3),2(3x,e10.3))

      close(unit=item)
      deallocate(ia, da, iperm)

      write(6,204)
204   format(//10x,'Sort completed successfully')
      call timer
      return
998   write(6,898) ir+1
898   format(/,5x,'Failure reading from unit item, record',i6,/)
!      stop
      return
       end

! --------------------------------------------------------------

      subroutine pfcalc(temp,ge,go,q,idia,gz)

      use energylevs
      implicit double precision(a-h,o-y), logical(z)

!     program constants
!     autocm converts atomic units to wavenumbers
      data x0/0.0d0/,x1/1.0d0/,&
           hc/ 1.9864476d-16/, &
           bk/ 1.3806581d-16/, & 
           autocm/ 2.19474624d+05/

      tk= bk*temp
      c1= hc/tk
      q=x0

!     calculate the partition coefficient

      do 3 i=1,nlevs
      grot=dble(2*ij(i)+1)
      if(abs(idia) .eq. 2) then
         if(isym4(i).le.2) then
            gg= ge*grot
         else
            gg= go*grot
         endif
      else
         gg=grot
      endif

      q= q + gg*exp(-(enlevs(i)-gz)*c1)
3     continue


      if(q.le.x0) q= x1

      write(6,*) "Partition function calculated to be ", q

      return
      end


! ------------------------------------------------------------

      subroutine spectm(xmin, wmin, wmax, dwl, &
                        ge, go, temp, q, nr, jdia, ispe, gz)

!     this is the main working subroutine which calculates the
!     required integrated absorption coefficients.

!     this routine is designed to calculate integrated absorption
!     coefficients using the data produced by subroutine sortsp.
!     it uses the formula:

!           4.162034*10(-19)*gg*w*(exp(-hcw"/kt) - exp(-hcw'/kt))*s(f-i)
!     i(w)= ------------------------------------------------------------
!                                         q

!     where:
!           gg  is the nuclear spin weighting
!               times the degeneracy weighting
!           w   is the transition frequency
!           w"  is the energy of the lower level in wavenumbers
!           w'  is the energy of the upper level in wavenumbers
!           t   is the temperature at which the spectrum is required
!        s(f-i) is the line strength in debye**2 for all magnetic
!               components of the line.
!           q   is the partition function.

!     this gives the intensity in cm./molecule.

!     in many cases the partition function will not be supplied, and it
!     is thus set at unity. this must be taken into account when using
!     integrated absorption coefficients to calculate expected spectral
!     intensities at given column densities.

!     relative absorption coefficients are given relative to the maximum
!     set to unity.

!     if zemit is true, the emissivity, rather than the integrated
!     absoption coefficient, is calculated from:

!               (2j'+1)*gg*hcw*exp(-hcw'/kt)*aif
!         j(w)= --------------------------------
!                             4pi*q

!
! New Namelist logical switches:
!    Only one of the following switches may be set to true.
!
!    z90     Outputs to ilist in fort.90 style, ie. index1, index2, Einstein A.
!    zassign Outputs to ilist, El, frequency intensity and Einstein A along with energy level 
!            Assignments. 
!
!
      use energylevs
      implicit double precision(a-h,o-y), logical(z)
      double precision, allocatable, dimension(:,:) :: a
      integer, allocatable, dimension(:,:) :: iqnum
      namelist /spe/ emin1,emax1,jmax,zplot,zemit,iplot,zfreq,zeinst, &
                     emin2,emax2,zprof,idat,zene,tinte,zlist,ilist, &
                     zdop,z90,zassign,prthr,npoints, xmolm
      common/logic/ zout, zsort, zspe, zpfun, zembed
      common/base/ ibase1,ibase2
!     program constants
      data idat/19/,iplot/20/,ilist/36/, &
           hc/ 1.9864476d-16/, &
           pi/ 3.1415927d0/, & 
           bk/ 1.3806581d-16/, & 
!     autocm converts atomic units to wavenumbers
           autocm/ 2.19474624d+05/, &
!     autode converts atomic units to debye
           autode/ 2.5417662d0/, & 
!     detosc converts from s(f-i) in debye**2 to seconds**(-1)
           detosc/ 3.136186d-07/, & 
           conv  / 4.162034d-19/ 

!     detocg converts from debye to c.g.s. units.

      fmem =  (76*nr)/1048576.0d0
      write(6,300) fmem
300   format(/5x,'Memory required to generate spectrum',f10.3,' MB'/)

      allocate(iqnum(nr,9), a(nr,6))

200   format(///)

!     preset namelist parameters

      emin1= -1.0d27
      emax1= +1.0d27
      emin2= -1.0d27
      emax2= +1.0d27
      jmax= -1
      tinte= +1.0d-15
      zplot= .false.
      zlist= .false.
      zprof= .false.
      npoints=3000
      zemit= .false.
      zfreq= .true.
      zeinst= .false.
      zene= .false.
      zdop= .false.
      z90= .false.
      zassign= .false.
      prthr = 0.1d0
      xmolm = 18.0d0 !molecular mass default (H20).

!     read in namelist parameters

      read(5,spe)
      if(zeinst) zemit=.true.
      write(6,*) "SPECTM output requests:"
      if (zplot) then
         write(6,1010) iplot
 1010    format(5x,'Data for plotting written to stream IPLOT =',i3)
         open(unit=iplot,form='formatted')
         if (zprof) then
            open(unit=idat,form='formatted')
            write(6,1020) npoints,idat
 1020    format(5x,'Full line profile requested at',i5,' points,',&
                 ' scatch  IDAT =',i3)
         endif
      else
         write(6,1030)
 1030    format(5x,'Plot file not requested')
      endif
      if (zlist) then
         write(6,1040) ilist
 1040    format(5x,'Full line list    written to stream ILIST =',i3)
         open(unit=ilist,form='formatted')
      else
         write(6,1050)
 1050    format(5x,'Full line list not requested')
      endif
      if (zlist) write(6,1060) prthr
 1060 format(/5x,'Linelist printed using relative threshold PRTHR =', &
        d9.2)

!     determine the spectral range required

      nskip= 0
      ntrans= 0
      nread= 0
      rewind ispe
      read(ispe)
10    read(ispe,end=90,err=99) ipar, j2, k2, ie2, j1, k1, ie1, index2,index1,e2, e1, w
      nread= nread+1
      w= w*autocm
      if(w.lt.wmin) then
         nskip= nskip + 1
         goto 10
      endif
      if(wmax.gt.0.0) then
         if(w.le.wmax) then
            ntrans= ntrans + 1
         else
            goto 90
         endif
      else
         ntrans= ntrans + 1
      endif
      goto 10
90    continue

!     skip unrequired transitions

      rewind ispe
      read(ispe)
      do 1 ir=1,nskip
       read(ispe)
1     continue

!     read in required data from ispe:

      do 2 ir=1,ntrans
      read(ispe) (iqnum(ir,ic),ic=1,9), (a(ir,ic),ic=1,4)
      a(ir,1)= (a(ir,1)*autocm)-gz
      a(ir,2)= (a(ir,2)*autocm)-gz
      a(ir,3)= a(ir,3)*autocm
      a(ir,4)= a(ir,4)*autode*autode
2     continue

!     start the calculation

      c1= hc/(temp*bk)
      if (zemit) goto 40

!     calculate integrated absorption coefficients

      c2= conv/q

     xmax= 0.0d0
      xint= 0.0d0
      if(jdia.eq.2) then
         do 3 ir=1,ntrans
         if(iqnum(ir,1).eq.0.or.iqnum(ir,1).eq.2) then
            gg= ge
         else
            gg= go
         endif
         acur= c2*gg*a(ir,3)*a(ir,4)* &
                  (exp(-a(ir,2)*c1) - exp(-a(ir,1)*c1))
         a(ir,5)= acur
         xint= xint + acur
         if(acur.gt.xmax) xmax= acur
3        continue
      else
         do 4 ir=1,ntrans
         acur= c2*a(ir,3)*a(ir,4)* &
                  (exp(-a(ir,2)*c1) - exp(-a(ir,1)*c1))
         a(ir,5)= acur
         xint= xint + acur

         if(acur.gt.xmax) xmax= acur
4        continue
      endif
      goto 41

!     calculate emissivity coefficients

40    c2= hc*detosc/(4.0d0*pi*q)
      xmax= 0.0d0
      xint= 0.0d0
      if(jdia.eq.2) then
         do 43 ir=1,ntrans
         w= a(ir,3)
         if(iqnum(ir,1).eq.0.or.iqnum(ir,1).eq.2) then
            gg= ge
         else
            gg= go
         endif
         if(.not. zeinst) then
            acur= w*w*w*w*a(ir,4)*c2*gg* &
                  exp(-a(ir,1)*c1)
         else
            acur= w*w*w*gg*a(ir,4)*detosc
         endif
         a(ir,5)= acur
         xint= xint + acur
         if(acur.gt.xmax) xmax= acur
43       continue
      else
         do 44 ir=1,ntrans
         w= a(ir,3)
         if(.not. zeinst) then
            acur= w*w*w*w*a(ir,4)*c2* &
                  exp(-a(ir,1)*c1)
         else
            acur= w*w*w*a(ir,4)*detosc
         endif
         a(ir,5)= acur
         xint= xint + acur
         if(acur.gt.xmax) xmax= acur
44       continue
      endif
41    continue

!     relative absorption or emissivity coefficients

      inc= 0
      neg= 0
      xmax1= 0.0d0
      do 5 ir= 1,ntrans
      if(jmax.gt.-1.and.iqnum(ir,5).gt.jmax) goto 5
      if(a(ir,2).lt.emin1) goto 5
      if(a(ir,2).gt.emax1) goto 5
      if(a(ir,1).lt.emin2) goto 5
      if(a(ir,1).gt.emax2) goto 5
         if(a(ir,5).gt.xmax1) xmax1= a(ir,5)
5     continue

      do 6 ir= 1,ntrans
      a(ir,6)= a(ir,5)/xmax1
      if(a(ir,6).lt.xmin) goto 7
      if(jmax.gt.-1.and.iqnum(ir,5).gt.jmax) goto 7
      if(a(ir,2).lt.emin1) goto 7
      if(a(ir,2).gt.emax1) goto 7
      if(a(ir,1).lt.emin2) goto 7
      if(a(ir,1).gt.emax2) goto 7
         inc= inc + 1
         goto 6
7     continue
         neg= neg + 1
6     continue

      if(jmax.eq.-1) then
         jm=500
      else
         jm= jmax
      endif
      if(.not.zemit) then
         write(6,2021) xmax,xint,emin1,emax1,emin2,emax2,jm,ntrans, &
                     inc,neg
      else
         write(6,2022) xmax,xint,emin1,emax1,emin2,emax2,jm,ntrans, &
                      inc,neg
      endif
2021  format(/5x,'Maximum absorption coefficient = ',e15.6,/, &
             10x,'units are cm./molecule.'//, &
             5x,'Integrated band absorption =',e15.6,//, & 
             5x,'Transitions with lower state energy = ',e13.6,' cm-1', &
             ' to ',e13.6,' cm-1 considered.',/, & 
             5x,'transitions with upper state energy = ',e13.6,' cm-1', &
             ' to ',e13.6,' cm-1 considered.',/, & 
             5x,'Maximum allowed value of J" =',i4,//, &
             i11,' transitions within spectral range',/, &
             i11,' transitions included',/, &
             i11,' transitions neglected')
2022  format(/5x,'Maximum emission coefficient = ',e15.6,/, &
             10x,'units are erg/molecule/sr.'//, &
             5x,'Integrated band emission =',e15.6,//, &
             5x,'Transitions with lower state energy = ',e13.6,' cm-1', &
             ' to ',e13.6,' cm-1 considered.',/, &
             5x,'Transitions with upper state energy = ',e13.6,' cm-1', &
             ' to ',e13.6,' cm-1 considered.',/, &
             5x,'maximum allowed value of j" =',i4,//, &
             i11,' transitions within spectral range',/, &
             i11,' transitions included',/, &
             i11,' transitions neglected')

!     final output

      write(6,*)
      i20= 0
      write(6,203)
203   format(/'ipar  j2 p2  i2  j1 p1  i1  index2 index1   e2         e1   ', &
       '        freq          s(f-i)       abs i(w)    rel i(w)      a(if) '/)
204   format(i3,2x,2i3,i4,1x,2i3,i4,2i7,3f13.6,1x,1p,e16.8,3e11.3)
      do 8 ir=1,ntrans
      if(jmax.gt.-1.and.iqnum(ir,5).gt.jmax) goto 8
      if(a(ir,2).lt.emin1) goto 8
      if(a(ir,2).gt.emax1) goto 8
      if(a(ir,1).lt.emin2) goto 8
      if(a(ir,1).gt.emax2) goto 8
      if(a(ir,6).lt.xmin) goto 8
         w= a(ir,3)
         aa= w*w*w*a(ir,4)*detosc/dble(2*iqnum(ir,2) + 1)
         iqnum(ir,3) = abs(iqnum(ir,3) - 1)
         iqnum(ir,6) = abs(iqnum(ir,6) - 1)
         if (zlist) then
          if(z90) then
           write(ilist,208) iqnum(ir,8), iqnum(ir,9),aa
          else
           if(zassign) then
            ind2=iqnum(ir,8)
            ind1=iqnum(ir,9)
            write(ilist,209) isym4(ind2),iqnum(ir,4),(iv(ind2,ic),ic=1,3), &
              (kj(ind2,ic),ic=1,3),isym4(ind1),iqnum(ir,7),(iv(ind1,ic),ic=1,3),&
              (kj(ind1,ic),ic=1,3),ind2,ind1, &
              a(ir,2),a(ir,3),a(ir,5),aa
           else
            write(ilist,204) (iqnum(ir,ic),ic=1,9),(a(ir,ic),ic=1,6),aa
           endif
          endif
         endif
         if(a(ir,6).ge.prthr) then
            write(6,204) (iqnum(ir,ic),ic=1,9),(a(ir,ic),ic=1,6),aa
            i20= i20 + 1
         endif
         if(zplot.and..not.zprof) then
            if( a(ir,6) .ge.tinte) then
              if ( zfreq ) then
                  xval= a(ir,3)
               else
                  xval= 10000.0/a(ir,3)
               endif
               if ( .not.zene ) then
                  write(iplot,207) xval, a(ir,6)
               else
                  write(iplot,206) xval, a(ir,6),a(ir,1),a(ir,2), &
                  iqnum(ir,2),iqnum(ir,5)
               endif
            endif
         endif
         if(i20.eq.20) then
            write(6,200)
            i20= 0
         endif
8     continue

      if(zplot.and.zprof) then
         do 9 ir=1,ntrans
         if( zfreq ) then
            a(ir,6)= a(ir,3)
         else
            a(ir,6)= 10000.0/a(ir,3)
         endif
         write(idat,*) a(ir,6), a(ir,5)
9        continue
         fmin= a(1,6)
         fmax= a(ntrans,6)
         call profil(idat,dwl,temp,iplot,fmin,fmax,npoints,&
                         xmolm,zdop,zfreq)
!         call profil(idat,dwl,iplot,fmin,fmax,npoints)
      endif
206   format(1x,f11.5,2x,f13.10,2x,f12.4,2x,f12.4,2x,i3,2x,i3)
207   format(1x,f11.5,2x,f20.10)
      write(6,200)
      write(6,205)
205   format(10x,'Spectrum completed successfully')
      return
99    write(6,*) nread, j2,k2,ie2,j1,k1,ie1,e2,e1,w

      deallocate(iqnum, a)

!      stop
      return
208   format(2i7,1p,e10.3)
209   format(i2,i5,6i3,2x,i2,i5,6i3,2x,2i7,2x,2f11.4,1p,2e10.3)
      end

! ----------------------------------------------------------

      subroutine getrow(row,nrow,iunit)

      implicit double precision (a-h,o-y)
      dimension row(nrow)
      read(iunit,err=999) row
      return
999   write(6,899) iunit
899   format(/,5x,'fell over in getrow reading unit = ',i6,/)
      stop
      end

! ------------------------------------------------------------

      subroutine outrow(row,nrow,iunit)

      implicit double precision (a-h,o-y)
      dimension row(nrow)
      write(iunit) row
      return
      end
      subroutine timer
!     prints current cpu time usage                                 #030

      implicit double precision (a-h,o-y)
      common/timing/itime0
      write(6,10)
      call SYSTEM_CLOCK(itime2,irate2,imax2)
      itime=(itime2-itime0)/irate2
      write(6,1)itime
 1    format(/i10,' secs CPU time used'/)
      return
10    format(/)
      end


! --------------------------------------------------------

      subroutine profil(idat,dwl,temp,iplot,fmin,fmax,npoints,&
                            xmolm,zdop,zfreq)

! This routine calculates a Gaussian or thermal Doppler line profile
! and applies it to the previously calculated integrated absorption 
! coefficient or the emissivity. The frequency dependent total 
! absorption coefficient and emissivity for the molecule are 
! calculated at each frequency or wavelength point. 
!
! Output is to the unit idat.
! Units of the line profiled absorption coefficient are [cm**2 /molecule]
! Units of the line profiled emissivity are [ergs cm /molecule /str].
!   If units of [ergs /micron /molecule /str] are desired, the user
!   must multiply by (10000.0 / (lambda[microns])**2)
!
! ARGUMENTS
! idat     I/O unit, on which frequency/intensity data is stored.
! dwl      Gaussian half width as specified by the user
! temp     Temperature
! ipolt    I/O unit, on which data is to be output.
! fmin     minimum of frequency range   
! fmax     maximum of frequency range
!
! Input arguments from namelist SPE.
!
! npoints  Number of frequency points, default is 3000.
! Zdop     logical parameter, if true then a thermal Doppler
!          line half width is used. The default is false, in which
!          case the user specified value of the half width at half
!          maximum of a Gaussian is used.
! Zfreq    Default .true. units of wavenumber, and calculation 
!          proceeds as discussed below. If .false. units of
!          wavelength in microns. For consistency of intensity units
!          the line profiling is performed only in wavenumbers (cm**-1).
!          But intensity data is output as a function of wavelength.
! xmolm    Mean molecular mass of the species, default is 18.0d0 (H20).
!
! OTHER INPUT
!
!Input from unit idat
!
! w        Wavenumber (cm**-1) or wavelength (microns) of line.
! x        Integrated absorption coefficient (cm / molecule) or
!          emissivity (ergs /s /str) of line.
!
! OUTPUT to unit iplot
!
! wl(npoints)  Wavenumber (cm**-1) or wavelength (microns) of point.
! f(npoints)   Absorption cross section (cm**2 / molecule) or
!              emissivity (ergs cm /s /str) at the point.
!
! The Doppler profile is:
!
!         1        ( m*c*c   )      ( -m*c*c             )
!  DLP = --- * sqrt(---------) * exp(----------- * Dv*Dv )
!         v        ( 2*pi*k*T)      (  2*k*T*v*v         )
!
! Where: v is the central line frequency or wavenumber
!        dv is the displacement from the line centre
!        m is the mass of particle, T is temperature.
!
! The Doppler half width at half maximum is thus:
!
!               ( 2*log(2)*k*T ) 
!  HW = v * sqrt( ------------ )
!               (     m*c*c    )
!
! Note: log(2) is the natural log of 2, ie FORTRAN notation.
!
! the profile can now be written:
!
!          1        (  log(2)  )      ( -Dv*Dv*log(2) )
!  DLP = ---- * sqrt( -------- ) * exp( ------------- )
!         HW        (   pi     )      (  HW*HW        )
!
! the factor sqrt(log(2)) is often dropped from the definition of
! the Doppler half width, so simplifying the all the expressions.
!
! Here, however, we use the rigorous definition of 
! half width at half maximum of the Gaussian.

      implicit double precision (a-h,o-y), logical(z)
      double precision Temp, xmolm, xln2, rtln2, rtln2pi, xlambda 

      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: f, wl

      allocate(f(npoints), wl(npoints))

      write(6,*)
      write(6,*) "LINE PROFILING"

      if(zdop) then
        write(6,*) "Thermal Doppler line profiling requested."
        write(6,*) "Molecular mass set to", xmolm
      else
        write(6,*) "Gaussian profiling requested."
        write(6,*) "Half width at half maximum set to ", dwl, " cm**-1"
      endif
      if (zfreq) then
        write(6,*) "Wavenumber units of cm**-1 selected"
      else
        write(6,*) "Wavelength units of microns selected"
      endif

      df= (fmax - fmin)/dble(npoints-1)
! some constants.
      pi= acos(-1.0d0)
      xln2=log(2.0d0)
      rtln2=sqrt(xln2)
      rtln2pi=rtln2/sqrt(pi)

!     zero the intensity array

      do 1 i= 1,npoints
      f(i)= 0.0d0
      wl(i)= fmin + dble(i-1)*df
1     continue
!     read in data from idat

      rewind idat
10    read(idat,*,end=11) w, x

      if(zfreq .eqv. .false.) then
        xlambda=w
        w=10000.0d0/w
      endif

      ! set the half width at half maximum [cm**-1]
      if(zdop) then
       hw=dop_half(temp,w,xmolm)
      else
       hw=dwl
      endif

      fac= rtln2pi/hw  ! Normalisation factor [cm]

      do 12 i=1,npoints
         ! set displacement from line centre in [cm**-1]
      if(zfreq) then
       delw= wl(i) - w
      else
       delw= (10000.0d0/wl(i))-w
      endif

      if(abs(delw).gt.10.0d0*hw) goto 12
      xarg= delw/hw
      xx= x*exp(-xarg*xarg*xln2)  ! Guassian line profile, ....
      f(i)= f(i) + (xx*fac)       ! normalise and bin. 
12    continue
      goto 10

!     output the spectrum

11    continue
      do 20 i=1,npoints

      write(iplot,*) wl(i), f(i)
20    continue

      deallocate(f, wl)

      return
      end

! ----------------------------------------------------------

      function dop_half(T,nu_cm, molm)

! This function calculates the Gaussian thermal Doppler half 
! width at half maximum.

! INPUT:
!        T:        temperature
!        nu_cm:    Wavenumber (cm-1).
!        molm:     Molecular mass
! OUTPUT:
!        dop_half: doppler half width (cm-1).

      implicit none

      double precision T, Nu_cm,molm, dop_half,  C

      data C /3.5682149d-7/

      dop_half=C*sqrt(T/Molm)*nu_cm

      return
      end

! -------------------------------------------------------------

      function qint_H2O(T)
! This function calculates the ro-vibrational H2O partition funtion from the fit of
! Vidler and Tennyson.
      implicit none
      double precision qint_H2O, T, AZ_mol(7)

      double precision logZ, logtheta
      integer i

      data (AZ_mol(i),i=1,7) /-14.087469157d0,37.924324853d0, &
                                   -42.681797873d0,25.330244851d0, &
                                    -8.108512629d0,1.331068717d0, &
                                    -0.087298105d0/
      logtheta=log10(T)
      logZ=0.0d0

      do 10 i=1, 7
       logZ=logZ+(az_mol(i)*(logtheta**(i-1)))
 10   continue

      qint_H2O=10.0d0**(logZ)

      return

      end

! ------------------------------------------------------