c to evaluate the rate coefficient given the output of the fits
c pgmartin, august 1995
c
      real*8 t,tt
      real*8 g,gcheck
      common /gmat/ g(349,349)
      real*8 tnorm
      common /ftype/ tnorm,mm,moff
      common /error/ ierr
      real*8 tpre,xpre,encon,encont
      common /precomp/ tpre,xpre,encont
c
      call ack
c input from original fit
      tnorm=4500.
      mm=13
      moff=4
c
      call readfmat
      ierr=0
      call checkempty
      call checkdiag
      call checkinvalid
      if(ierr.eq.0) then
         write(6,*) ' elementary tests passed '
      else
         write(6,*) ' elementary tests failed '
         stop
      endif
c
      call readenergies
c
 1    continue
      write(6,*) ' '
      write(6,*) ' ===== start interactive input ====='
      write(6,*) ' temperature (-ve to stop)'
      read(5,*) t
      if(t.lt.0.) stop
      call tcheck(t)
c precompute common factors
      tpre=t/tnorm
      xpre=dlog10(tpre)
      encon=1.43876881d0
      encont=encon/t
c
      write(6,*) ' making the full gamma matrix '
      call gmake
c
      write(6,*) ' 0: write and read matrix g.dat; 1: skip'
      read(5,*) jfast   
      if(jfast.ne.1) then
         call gwrite(t)
         call gread(tt)
      endif
c     
c now test it against h2hone.f, original data, etc.
c     
      write(6,*) ' '
      write(6,*) ' v,j for initial state and final state'
      write(6,*) ' v,j for the continuum (cid) is v>=99, j=anything '
      write(6,*) ' for new start -1 0 0 0 '
 2    continue
      write(6,*) ' new pair of v,j states '
      read(5,*) ivi,iji,ivo,ijo
      if(ivi.lt.0) go to 1
      call vjcol(ivi,iji,icoli)
      call vjcol(ivo,ijo,icolo)
      if(icoli.eq.-1 .or. icolo.eq.-1) go to 2
      gcheck=g(icolo,icoli)
      if(icoli.eq.icolo) then
         write(6,326) gcheck
      else
         write(6,324) gcheck
      endif
      write(6,*) ' '
      go to 2
c
 324  format(1x,' final rate coefficient in cm3 s-1: ',1pd11.3)
 326  format(1x,' final exit rate coefficient in cm3 s-1: ',1pd11.3)
      end
c
      subroutine colvj(icol,n,j)
c  this subroutine returns quantum numbers given the column index
c  used in master equation programs
c
c  input
c  icol:  column index  1, 349
c     para:             1, 177
c     ortho:            178, 348 (171)
c     dissociation:     349
c  output
c  n      state vibrational quantum number 0, 14
c  j      state rotational quantum number 0, 38 or 1, 37
c  n=j=99, ind=1 for dissociation
c
      dimension iv(39)
c  this array is the maximum of v+1 for j=0,38
c      data iv/15,15,15,15,15,14,14,14,14,13,13,13,13,12,12,11,11,10,10,
      data iv/15,15,15,15,15,14,14,14,14,13,13,13,12,12,12,11,11,10,10,
     19,9,9,8,8,7,7,6,6,5,5,4,4,3,3,2,2,1,1,1/
      idiss=349
      ioff=177
c  compute v and j
      if(icol.ne.idiss) then
         if(icol.gt.ioff) then
c fortran indices
c ortho
            istart=2
            icheck=icol-ioff
         else
            icheck=icol
            istart=1
         endif
         do 100 i=istart,39,2
            icheck=icheck-iv(i)
            if(icheck.le.0)go to 101
 100     continue
         write(6,9999)icol
 9999    format(1x,'icol=',i3,' how did it get here?')
         stop
 101     continue
         j=i-1
         n=icheck+iv(i)-1
      else
c dissociation
         n=99
         j=99
      end if
      return
      end
c
      subroutine vjcol(n,j,icol)
c  routine using general sindex to find the column index for populations
c
c  input
c  n      state vibrational quantum number 0, 14
c  j      state rotational quantum number 0, 38 or 1, 37
c  n=j=99, ind=1 for dissociation (used in other programs)
c  but n>=99, j=anything works
c
c  output
c  icol:  column index  1, 349
c     para:             1, 177
c     ortho:            178, 348 (171)
c     dissociation:     349
c
      ioff=177
      if(n.lt.99) then
         call sindex(n,j,icol)
         if(icol.lt.0) go to 1
         if(mod(j,2).eq.1) icol=icol+ioff
      else
         icol=349
      end if
 1    continue
      return
      end
c
      subroutine sindex(n,j,indvj)
c  this subroutine returns indvj [1-178] given 
c   n     vibrational quantum number
c   j     rotational quantum number
      dimension iv(39)
c      data iv/15,15,15,15,15,14,14,14,14,13,13,13,13,12,12,11,11,10,10,
      data iv/15,15,15,15,15,14,14,14,14,13,13,13,12,12,12,11,11,10,10,
     19,9,9,8,8,7,7,6,6,5,5,4,4,3,3,2,2,1,1,1/
c specifically use >=99 * as continuum -- this combination does not get here
      if(n.gt.14 .and. n.lt.99) go to 200
      if(j.gt.38 .or. j.lt.0) go to 200
      if( (n+1).gt.iv(j+1) .or. n.lt.0) go to 200
      istart=1
      if(mod(j,2).eq.1)istart=2
      icount=1
      istop=j
c avoid the uncertain convention of a single pass or not on negative first step
      if(istart.gt.istop)go to 20
      do 10 i=istart,istop,2
      icount=icount+iv(i)
   10 continue
   20 continue
      icount=icount+n
      indvj=icount
      return
 200  continue
      write(6,9998) n,j
 9998 format(1x,'n, j = ',2i5,' illegal v,j combination')
      indvj=-1
      return
      end
c
      real*8 function fun2(nn,xc)
c evaluate analytic function of temperature
      implicit real*8(a-h,o-z)
      dimension xc(nn)
      common /precomp/ tpre,xpre,encont
c xin is log10tnorm
c      t=10.d0**x
c modified arrhenius
c a t**b e-c/t
c a +bx  +cx*x -c/t
c      fun2=xc(1) +xc(2)*x  + xc(3)*x*x -xc(4)/t
      fun2=xc(1) + (xc(2)  + xc(3)*xpre)*xpre -xc(4)/tpre
      return
      end
c
      subroutine readfmat
c read the matrices of fitting coefficients a, b, c, d
      parameter(m349=349,m4=4)
      real*4 nr,ex
      common /nrex/ nr(m4,m349,m349),ex(m4,m349,m349)
      write(6,*) ' reading fitting parameters '
      write(6,*) ' please stand by; do not adjust your set'
      open(unit=1,file='nr.dat',status='old')
      open(unit=2,file='ex.dat',status='old')
      read(1,3200) nr
      read(2,3200) ex
      close(unit=1)
      close(unit=2)
c make this more compact for distribution; requires formatted read
c the last field for d adds more accuracy for the reconstruction
 3200 format(4(3f8.3,f9.4))
c 320  format(12f9.3)
      return
      end
c
      subroutine checkempty
c
c check that empty elements have all been replaced during the fits
c
      common /error/ ierr
      parameter(m349=349,m4=4)
      real*4 nr,ex
      common /nrex/ nr(m4,m349,m349),ex(m4,m349,m349)
c
      iee=0
      do 251 j=1,m349
         do 250 i=1,m349
            call empty(nr,m4,m349,i,j,ie)
            if(ie.ne.0) iee=1
            call empty(ex,m4,m349,i,j,ie)
            if(ie.ne.0) iee=1
c skip out at first i,j occurence
            if(iee.ne.0) then
               ierr=1
               go to 253
            endif
 250     continue
 251  continue
c
c      write(6,*) ' ok: empty elements all replaced '
c
 253  continue
      return
      end
c
      subroutine empty(x,m4,m349,i,j,ie)
      dimension x(m4,m349,m349)
c there should be none of these left at the end
      parameter(zdiag=-800.,zinvalid=-900.,zempty=-600.)
c
      call nzj(zempty,x,m4,m349,i,j,ie)
      if(ie.eq.1) write(6,*) ' empty elements above not all replaced '
      return
      end
c
      subroutine checkdiag
c
c check that diagonal elements have all been inserted with zdiag
c
      common /error/ ierr
      parameter(m349=349,m4=4)
      real*4 nr,ex
      common /nrex/ nr(m4,m349,m349),ex(m4,m349,m349)
c
      iee=0
      do 251 j=1,m349
         call diag(nr,m4,m349,j,ie)
         if(ie.ne.0) iee=1
         call diag(ex,m4,m349,j,ie)
         if(ie.ne.0) iee=1
c skip out at first i,j occurence
         if(iee.ne.0) then
            ierr=1
            go to 253
         endif
 251  continue
c
c      write(6,*) ' ok: diagonal elements all replaced '
c
 253  continue
      return
      end
c
      subroutine diag(x,m4,m349,j,ie)
      dimension x(m4,m349,m349)
      parameter(zdiag=-800.,zinvalid=-900.,zempty=-600.)
c
      call yzj(zdiag,x,m4,m349,j,j,ie)
      if(ie.eq.1) write(6,*) ' diag elements above not all replaced '
      return
      end
c
      subroutine checkinvalid
c
c check that invalid transition elements have all been replaced
c
      common /error/ ierr
      parameter(m349=349,m4=4)
      real*4 nr,ex
      common /nrex/ nr(m4,m349,m349),ex(m4,m349,m349)
      parameter(zdiag=-800.,zinvalid=-900.,zempty=-600.)
c
      iee=0
      do 251 j=1,m349-1
         do 250 i=1,m349-1
            call nzj(zinvalid,ex,m4,m349,i,j,iee)
            if(iee.ne.0) then
               write(6,*) ' invalid elements in ex '
               ierr=1
               go to 253
            endif
 250     continue
 251  continue
c
c      write(6,*) ' ok: no invalid elements in ex '
c
 253  continue
c
      iee=0
      do 1251 icoli=1,m349-1
         do 1250 icolo=1,m349-1
            call colvj(icoli,ivi,iji)
            call colvj(icolo,ivo,ijo)
            iopi=mod(iji,2)+1
            iopo=mod(ijo,2)+1
c whether invalid or not depends on parity when non-reactive
            if(iopi.eq.iopo) then
               call nzj(zinvalid,nr,m4,m349,icolo,icoli,iee)
c flag a 1 when numbers agree and they should not
               if(iee.ne.0) then
                  write(6,*) ' incorrect invalid elements in nr '
                  go to 1253
               endif
            else
               call yzj(zinvalid,nr,m4,m349,icolo,icoli,iee)
c flag a 1 when numbers disagree and they should
               if(iee.ne.0) then
                  write(6,*) ' missing invalid elements in above nr '
                  ierr=1
                  go to 1253
               endif
            endif
 1250    continue
 1251 continue
c
c      write(6,*) ' ok: no incorrect or missing invalid elements in nr '
c
 1253 continue
c
c now do the continuum elements
      iee=0
      do 2250 icolo=1,m349-1
         call yzj(zinvalid,nr,m4,m349,icolo,m349,ie)
         if(ie.ne.0) iee=1
         call yzj(zinvalid,ex,m4,m349,icolo,m349,ie)
         if(ie.ne.0) iee=1
         if(iee.ne.0) then
            write(6,*) ' missing invalid elements in column 349 above'
            ierr=1
            go to 2253
         endif
 2250 continue
c
c      write(6,*) ' ok: no missing invalid elements in column 349'
c
 2253 continue
c     
      iee=0
      do 2252 icoli=1,m349-1
         call nzj(zinvalid,nr,m4,m349,m349,icoli,iee)
         if(iee.ne.0) then
            write(6,*) ' incorrect invalid elements in nr row 349'
            ierr=1
            go to 2255
         endif
 2252 continue
c
c      write(6,*) ' ok: no incorrect invalid elements in nr row 349'
c
 2255 continue
c
      iee=0
      do 2251 icoli=1,m349-1
         call yzj(zinvalid,ex,m4,m349,m349,icoli,iee)
         if(iee.ne.0) then
            write(6,*) ' missing invalid elements in ex row 349 above'
            ierr=1
            go to 2254
         endif
 2251 continue
c
c      write(6,*) ' ok: no missing invalid elements in ex row 349'
c
 2254 continue
c
      return
      end
c
      subroutine izj(ze,x,m4,m349,i,j,ie)
      dimension x(m4,m349,m349)
c flag a 1 when numbers agree
      ie=0
      eps=1.e-5
      zmax=0.
      do 1 k=1,m4
         z=abs(x(k,i,j)/ze-1.)
         if(z.gt.zmax) zmax=z
 1    continue
      if(zmax.le.eps) ie=1
      return
      end
c
      subroutine yzj(ze,x,m4,m349,i,j,ie)
      dimension x(m4,m349,m349)
c flag a 1 when numbers disagree and they should
      ie=0
      eps=1.e-5
      zmax=0.
      do 1 k=1,m4
         z=abs(x(k,i,j)/ze-1.)
         if(z.gt.zmax) zmax=z
 1    continue
      if(zmax.le.eps) then
c do nothing
      else
         write(6,320) i,j,(x(k,i,j),k=1,m4),ze
         ie=1
      endif
 320  format(2i5,12f9.3)
      return
      end
c
      subroutine nzj(ze,x,m4,m349,i,j,ie)
      dimension x(m4,m349,m349)
c flag a 1 when numbers agree and they should not
      ie=0
      eps=1.e-5
      zmax=0.
      do 1 k=1,m4
         z=abs(x(k,i,j)/ze-1.)
         if(z.gt.zmax) zmax=z
 1    continue
      if(zmax.le.eps) then
c all m4 numbers close to ze for this i,j
         write(6,320) i,j,(x(k,i,j),k=1,m4),ze
         ie=1
      endif
 320  format(2i5,12f9.3)
      return
      end
c
      subroutine readenergies
c input energies used for detailed balance
      real*8 ecol,eni,eno
      common /ener/ ecol(348),eni,eno
      open(unit=3,file='ecol.dat',status='old')
      read(3,*) ecol
      close(unit=3)
      return
      end
c

      subroutine rec(ivi,iji,ivo,ijo,gammau,gammad)
c reconstruct, do detailed balance and sum exchange and non-reactive
      implicit real*8 (a-h,o-z)
      common /ener/ ecol(348),eni,eno
      common /precomp/ tpre,xpre,encont
c
      call vjcol(ivi,iji,icoli)
      call vjcol(ivo,ijo,icolo)
      iopi=mod(iji,2)+1
      iopo=mod(ijo,2)+1
c
c to continuum
      if(icolo.eq.349) then
         call reconstruct(icolo,icoli,iopo,iopi,ratenr1,rateex1)
         gammau=ratenr1
         return
      endif
c
c otherwise only enter with eni.lt.eno, via common
c up
c look at reverse transition, since want up from down 
      call reconstruct(icoli,icolo,iopi,iopo,ratenr2,rateex2)
      gammad = ratenr2 + rateex2
c
c now do the detailed balance 
c
      gii=(2.d0*dfloat(iopi-1)+1.d0)*(2.d0*dfloat(iji)+1.d0)
      goo=(2.d0*dfloat(iopo-1)+1.d0)*(2.d0*dfloat(ijo)+1.d0)
      conv=(goo/gii)*dexp(-(eno-eni)*encont)
      ratenr2=ratenr2*conv
      rateex2=rateex2*conv
      gammau = ratenr2 + rateex2
      return
      end
c
      subroutine reconstruct(icolo,icoli,iopo,iopi,fratenr,frateex)
      implicit real*8 (a-h,o-z)
c input arrays
      parameter(m349=349,m4=4)
      real*4 nr,ex
      common /nrex/ nr(m4,m349,m349),ex(m4,m349,m349)
      dimension xc(m4)
      common /ftype/ tnorm,mm,moff
c for checking
      real*4 zzero,zzeros
      parameter(zzero=-700.,zzeros=-701.)
c
c nr
c
c screen the ones that should be zero (no valid transition, including diagonal)
c no diagonal
      if(icoli.eq.icolo) go to 6
c invalid: nothing from continuum
      if(icoli.eq.m349) go to 6
c invalid: no non-reactive to different parity
      if(iopi.ne.iopo .and. icolo.ne.m349) go to 6
c check if zero
      call izj(zzero,nr,m4,m349,icolo,icoli,ie)
c flag a 1 when numbers agree 
      call izj(zzeros,nr,m4,m349,icolo,icoli,ies)
      if(ie.eq.0 .and. ies.eq.0) then
c finally a non-zero rate coefficient to reconstruct
         do 4 k=1,m4
            xc(k)=nr(k,icolo,icoli)
 4       continue
c a +bx  +cx*x -d(1/t - 1) going back to  a +bx  +cx*x -d/t in fun2
         xc(1)=xc(1)+xc(m4)
         ratenr=fun2(m4,xc)
         fratenr=10.d0**ratenr
      else
         fratenr=0.d0
      endif
      go to 7
c
c nr invalid transitions (including diagonal) arrive here
c
 6    continue
      fratenr=0.d0
c
c ex
c
 7    continue
c screen the ones that should be zero (no valid transition and diagonal)
c no diagonal
      if(icoli.eq.icolo) go to 8
c invalid: nothing from continuum
      if(icoli.eq.m349) go to 8
c invalid: nothing to continuum
      if(icolo.eq.m349) go to 8
      call izj(zzero,ex,m4,m349,icolo,icoli,ie)
      call izj(zzeros,ex,m4,m349,icolo,icoli,ies)
      if(ie.eq.0 .and. ies.eq.0) then
c not zero
         do 5 k=1,m4
            xc(k)=ex(k,icolo,icoli)
 5       continue
         xc(1)=xc(1)+xc(m4)
         rateex=fun2(m4,xc)
         frateex=10.d0**rateex
      else
         frateex=0.d0
      endif
      return
c
c ex invalid (and diagonal) transitions arrive here
c
 8    continue
      frateex=0.d0
c
      return
      end
c
      subroutine tcheck(t)
      real*8 t,temp
      real*8 tnorm
      common /ftype/ tnorm,mm,moff
      dimension temp(22)
      data temp/1.d2,1.4d2,2.d2,3.d2,4.5d2,6.d2,7.8d2,
     1          1.d3,1.4d3,2.d3,3.d3,4.5d3,6.d3,7.8d3,
     1          1.d4,1.4d4,2.d4,3.d4,4.5d4,6.d4,7.8d4,
     1          1.d5/
c
      if(t.lt.temp(1+moff)) then
         write(6,*) ' warning: temperature', t, 'too small'
      endif
      if(t.gt.temp(mm+moff)) then
         write(6,*) ' warning: temperature', t, 'too large'
      endif
      return
      end
c
      subroutine gwrite(t)
      real*8 g,t
      common /gmat/ g(349,349)
      write(6,*) ' writing g.dat '
      open(unit=4,file='g.dat',status='unknown')
      write(4,111) t
 111  format(1x,1pe7.1,' non-reactive plus reactive with diagonal')
      write(4,3000) g
 3000 format(1p8d11.3)
      close(unit=4)
      return
      end
c
      subroutine gread(t)
      implicit real*8 (a-h,o-z)
      common /gmat/ g(349,349)
      common /ener/ ecol(348),eni,eno
      common /precomp/ tpre,xpre,encont
      open(unit=14,file='g.dat',status='old')
      read(14,*) t
      write(6,*) ' reading g.dat '
c there is no space between some numbers divided by a minus sign
c unformatted read causes a problem on the alpha but not the sgi
      read(14,3000) g
 3000 format(1p8d11.3)
      close(unit=14)
c
c warning:
c if this were a read into a master equation program,
c detailed balance (up from down) would be required 
c to maintain machine precision, since numbers not stored
c to high precision [e.g., format(1p8d10.3)]
c
      encon=1.43876881d0
      encont=encon/t
      do 1000 icoli=1,348
         call colvj(icoli,ivi,iji)
         do 2000 icolo=1,348
            if(icolo.ne.icoli) then
               eni=ecol(icoli)
               eno=ecol(icolo)
               if(eni.lt.eno) then
                  call colvj(icolo,ivo,ijo)
                  iopi=mod(iji,2)+1
                  iopo=mod(ijo,2)+1
c
c now do the detailed balance 
c
                  gii=(2.d0*dfloat(iopi-1)+1.d0)*
     $                 (2.d0*dfloat(iji)+1.d0)
                  goo=(2.d0*dfloat(iopo-1)+1.d0)*
     $                 (2.d0*dfloat(ijo)+1.d0)
                  conv=(goo/gii)*dexp(-(eno-eni)*encont)
                  g(icolo,icoli)=g(icoli,icolo)*conv
               endif
            endif
 2000    continue
 1000 continue
c also reconstruct the diagonal
      do 1001 icoli=1,349
         g(icoli,icoli)=0.d0
         gamma=0.d0
         do 2001 icolo=1,349
            gamma=gamma+g(icolo,icoli)
 2001    continue
            g(icoli,icoli)=-gamma
 1001    continue
c
      return
      end
c
      subroutine gmake
      real*8 gammau,gammad
      real*8 g
      real*8 ecol,eni,eno
      common /gmat/ g(349,349)
      common /ener/ ecol(348),eni,eno
c
      do 1000 icoli=1,348
         call colvj(icoli,ivi,iji)
         call colvj(349,ivo,ijo)
         call rec(ivi,iji,ivo,ijo,gammau,gammad)
         g(349,icoli)=gammau
         g(icoli,349)=0.d0
         do 2000 icolo=1,348
            if(icolo.eq.icoli) then
               g(icolo,icoli)=0.d0
            else
               eni=ecol(icoli)
               eno=ecol(icolo)
               if(eni.lt.eno) then
                  call colvj(icolo,ivo,ijo)
                  call rec(ivi,iji,ivo,ijo,gammau,gammad)
                  g(icolo,icoli)=gammau
                  g(icoli,icolo)=gammad
               endif
            endif
 2000    continue
 1000 continue
c construct the diagonal
      do 1001 icoli=1,349
         gamma=0.d0
         do 2001 icolo=1,349
            gamma=gamma+g(icolo,icoli)
 2001    continue
            g(icoli,icoli)=-gamma
 1001    continue
      return
      end
c
      subroutine ack
      write(6,*) ' Collisional rate coefficients for H + H2'
      write(6,*) ' P.G. Martin and M.E. Mandy 1996'
      write(6,*) ' Astrophysical Journal Letters, 455, L89'
      write(6,*) ' '
      return
      end
c
