Program Gasprint

* Reads the Gastherm data file, and writes the reactions in a more readable form,
*   and writes log K's and regression coefficients for intermediate temperatures.
*   Derived from program Solprint, a similar program for the Soltherm data file.
*   Contributors to Solprint include M.Reed, N.Spycher, J.DeBraal, J.Palandri
*   August ** 2023 - J.Palandri

      implicit none

      real*8 coef(9), weight, wt(50)

      integer inp1, iout1, eofflag
      integer bstot, ispec(9), i, itot, jj, l, ns

      character*8  blank, dat(10), hyd, name, cmpname(50)
      character*80 dum

      common /stuff/ coef, wt, ispec, cmpname

      inp1    =  1
      iout1   = 10
      eofflag =  0
      ns      =  0

      blank = '        '
      hyd   = 'H2      '
      name  = '12345678'

      weight = 0.0

      open (inp1,  file = 'gastherm.dat', status = 'old')
      open (iout1, file = 'GasthermPrint.txt', status = 'unknown')

      write (iout1,'(//,a,a)') 'For each record below, lines 1 and 2',
     +      ' show the balanced reaction.'
      write (iout1,'(a,a)') 'Line 3: Log K''s at  25  100  200  300 ',
     +      '400 500 600 700 800 deg C.'
      write (iout1,'(a,a)') 'Line 4: Log K''s at 900 1000 1100 1200',
     +      ' deg C (all at one bar.)'
      write (iout1,'(a,a)') 'Line 5: Regression coefficients ',
     +      '(a,b,c,d,e) for log K''s as f(T).'
      write (iout1,'(a,/,a)') '     Where T is in Kelvins and ',
     +      '     f(T) = a + b/T + c*T + d/(T**2) + e*log(T).'

      write (*,'(//2x,a)')
     +   'Input the number of basis component gases in Gastherm >>  '
      read *, bstot
      write (iout1, '(//2x,a17,/)') 'Component Species'

*Skip all the entries before the component species; read and write them
      do while (dum(1:8) .ne. hyd)
        read (inp1, '(a80)', iostat = eofflag) dum
        if (eofflag.lt.0) then
          write (*,'(//2x,A58)')
     +      'Stop error: EOF read before the H2 basis component found'
          stop
        endif
      end do
      backspace (inp1)

      do 10 i = 1, bstot
        read  (inp1, '(a8,4x,f10.4)') cmpname(i), wt(i)
        write (iout1, '(5x, a8, 4x, f10.4)') cmpname(i), wt(i)
   10 continue

*Read the data for derived gases, and print it
      write (iout1, '(//2x,a15,/)') 'Derived Gases'
      do while (name.ne.blank)
        read (inp1,'(a8,8x,i1,11(f6.3,i2),f5.3, i2)') name,
     +       itot, (coef(i), ispec(i),i = 1, itot)
        if (name.ne.blank) then
          ns = ns + 1
          itot = itot - 1
          do 20 i = 1, itot
            jj = i+1
            coef(i) = coef(jj)
            ispec(i) = ispec(jj)

*If H2 is in the derived gas but is not in the I=1 position, flag it.
            if(i.gt.1.and.ispec(i).eq.1) then
              write(iout1,663) name
              write(*,663) name
  663              format('  WARNING: Species ',a8, ' has H2 in a'
     +             ' position beyond I = 2.')
            endif
   20     continue
          call print (name, itot, weight)
          do 30 l = 1,3
            read (inp1, '(8x, 10a8)') (dat(i), i = 1, 10)
            write (iout1, '(4x, 10a8)') (dat(i), i = 1, 10)
   30     continue
          write (iout1,*)
        endif
      end do
      write (*,500) ns
      write (iout1,500) ns
  500       format ('     There are ', I4,' derived gases.',//)
      read (inp1,*)

*Read data for gases and minerals, skipping inappropriate data, and write it.
      write (iout1, '(/2x,a,/)') 'Minerals, Solids, and Liquids'

      ns = 0

      name  = '12345678'

      do 40 i = 1, 9
        coef(i) = 0.0
        ispec(i) = 0
   40 continue
      do while (name.ne.blank)
        read (inp1, 520) name, weight, itot, (coef(i),
     +       ispec(i), i = 1, itot)
  520   format (a8, 2x, f7.2, i2, 1x, 9(f6.3, i2))
        if (name.ne.blank) then
          ns = ns + 1
          call print (name, itot, weight)
          do 100 l = 1, 3
            read (inp1, '(8x, 10a8)') (dat(i), i = 1, 10)
            write (iout1, '(4x, 10a8)') (dat(i), i = 1, 10)
  100     continue
          write (iout1, '(2x)')
        end if
      end do
      write (*,530) ns
      write (iout1,530) ns
  530       format ('     There are ', I4,' Minerals, solids, and ',
     +             'liquids.',//)
      end


      Subroutine Print (name, itot, weight)
* =====================================================================|
* Prints the reaction in readable form

      implicit none

      real*8 coef(9), prodcf(9), reacf(9), weight, wt(50), wtm

      integer i, inp1, iout1, iprod, ireac, ispec(9), itot, pp, tright

      character *  3 char1(9), char2(9)
      character *  8 name, prodnm(9), reacnm(9), cmpname(50)
      character * 80 blank, buffer

      external tright

      common /stuff/ coef, wt, ispec, cmpname

      data char1 /'   ', 8*' + '/
      data char2 /' = ', 8*' + '/
      data blank /' '/

      parameter (inp1 = 1, iout1 = 10)
      ireac = 1
      iprod = 0
      reacf(1) = 1.0
      reacnm(1) = name
      wtm = 0.0

      do 40 i = 1, itot
        wtm   = wtm + coef(i) * wt(ispec(i))
        if (coef(i)) 10, 10, 20
   10     continue
          ireac = ireac + 1
          reacf(ireac) = -1.0 * coef(i)
          reacnm(ireac) = cmpname(ispec(i))
          go to 30

   20     iprod = iprod + 1
        prodcf(iprod) = coef(i)
        prodnm(iprod) = cmpname(ispec(i))
   30   continue
   40 continue

      pp = 1
      buffer = ' '

      do 43 i = 1, ireac
        if (i .ne. 1) then
          write (buffer(pp:), '(a3)') char1(i)
          pp = pp + 3
        end if

        if (tright(reacf(i)) .eq. 1) then
          write (buffer(pp:), '(f5.1)') reacf(i)
          pp = pp + 5
        else if (tright(reacf(i)) .eq. 2) then
          write (buffer(pp:), '(f6.2)') reacf(i)
          pp = pp + 6
        else
          write (buffer(pp:), '(f7.3)') reacf(i)
          pp = pp + 7
        end if

        write (buffer(pp:), '(1x, a8)') reacnm(i)
        pp = pp + 9

        if (pp + 20 .ge. 80) then
          write (iout1, '(a80)') buffer
          pp = 2
          buffer = ' '
        end if

   43 continue

      if (buffer .ne. blank) then
        write (iout1, '(a80)') buffer
      end if

      pp = 5
      buffer = ' '

      do 45 i = 1, iprod
        write (buffer(pp:), '(a3)') char2(i)
        pp = pp + 3

        if (tright(prodcf(i)) .eq. 1) then
          write (buffer(pp:), '(f5.1)') prodcf(i)
          pp = pp + 5
        else if (tright(prodcf(i)) .eq. 2) then
          write (buffer(pp:), '(f6.2)') prodcf(i)
          pp = pp + 6
        else
          write (buffer(pp:), '(f7.3)') prodcf(i)
          pp = pp + 7
        end if

        write (buffer(pp:), '(1x, a8)') prodnm(i)
        pp = pp + 9

        if (pp + 20 .ge. 80) then
          write (iout1, '(a80)') buffer
          pp = 5
          buffer = ' '
        end if

   45 continue

      if (buffer .ne. blank) then
        write (iout1, '(a80)') buffer
      end if

      wtm = wtm - weight

   50 continue
      if(weight.eq.0.) return
      if(abs(wtm).lt.0.1) return

      write (iout1, 415) wtm
  415 format(/5x,' *** warning: mol.weight balance is: ',f8.3/)

      write (*,416) wtm, name
  416 format(/5x,'*** warning: mol.weight balance is: ',f8.3,' for ',a8)

      return
      end


      Integer Function Tright (stuff)
* =====================================================================|
* For each reaction coefficient, finds and returns the appropriate count
*   of significant figures to write after the decimal,
*   to greatly reduce the number of zeros shown.

      real * 8 cpumin, fractn, number, stuff
      parameter (cpumin = 1.0d-35)
      intrinsic dint, dabs

      number = stuff

      fractn = dint(number) - number
      if (dabs(fractn) .le. cpumin) then
        tright = 1
        return
      end if

      number = number * 10
      fractn = dint(number) - number
      if (dabs(fractn) .le. cpumin) then
        tright = 1
        return
      end if

      number = number * 10
      fractn = dint(number) - number
      if (dabs(fractn) .le. cpumin) then
        tright = 2
        return
      end if

      tright = 3
      return
      end