program top500rd
  !
  ! Read Compustat file and write out data for top 500 R&D spenders
  ! 
  implicit none 
  
  integer                     :: i,j,k,l,m,found,gvkey,datamonth,dataday,    &
                                 datayear,fyear,sic,costat,re

  integer, dimension(64)      :: yrvec
  integer, dimension(500,64)  :: gvkeyrd,naicsrd,sicrd
  integer, dimension(64)      :: numrd
  real*8, dimension(500,64)   :: toprd,topad,topce,topat,topsl,topem
  real*8, dimension(64)       :: sumrd,sumad,sumce,sumat,sumsl,sumem
  character(len=14)           :: naics
  character(len=1)            :: digit1,digit2,digit3,digit4
  character(len=3)            :: loc,fic
  character(len=4)            :: blank4
  character(len=7)            :: blank7
  real*8                      :: at,sale,xad,xrd,capx,capxv,emp

  open(unit=5, file='company.inp')
  open(unit=7, file='top500rd.out')
  open(unit=8, file='top500rd.dat')

  naicsrd  = -99
  sicrd    = -99
  gvkeyrd  = -99
  toprd    = -99.
  topad    = -99.
  topce    = -99.
  topat    = -99.
  topsl    = -99.
  topem    = -99.
  yrvec(1) = 1950
  do i=2,64
    yrvec(i) = yrvec(i-1)+1
  enddo
  read(5,'(a14)') naics
  sumrd    = 0.
  numrd    = 0
  sumad    = 0.
  sumce    = 0.
  sumat    = 0.
  sumsl    = 0.
  sumem    = 0.
  do i=1,399864
    read(5,'(i6,i6,i3,i5,i10,a4,a3,a14,i4,i7,a7,a3,f22.4,f22.4,f22.4,f22.4,  &
      f22.4,f22.4,f22.4)') gvkey,datamonth,dataday,datayear,fyear,blank4,    &
      loc,naics,sic,costat,blank7,fic,at,sale,xad,xrd,capx,capxv,emp

    digit1     = naics(5:5)
    digit2     = naics(6:6)
    digit3     = naics(7:7)
    digit4     = naics(8:8)
    re         = 0;
    if ((digit1=='5').and.(digit2=='3')) re=1

    if ((at>0).and.(loc=='USA').and.(fic=='USA').and.(re==0)) then

      l=0
      if (digit1 == '1') l=1000
      if (digit1 == '2') l=2000
      if (digit1 == '3') l=3000
      if (digit1 == '4') l=4000
      if (digit1 == '5') l=5000
      if (digit1 == '6') l=6000
      if (digit1 == '7') l=7000
      if (digit1 == '8') l=8000
      if (digit1 == '9') l=9000
      if (digit2 == '1') l=l+100
      if (digit2 == '2') l=l+200
      if (digit2 == '3') l=l+300
      if (digit2 == '4') l=l+400
      if (digit2 == '5') l=l+500
      if (digit2 == '6') l=l+600
      if (digit2 == '7') l=l+700
      if (digit2 == '8') l=l+800
      if (digit2 == '9') l=l+900
      if (digit3 == '1') l=l+10
      if (digit3 == '2') l=l+20
      if (digit3 == '3') l=l+30
      if (digit3 == '4') l=l+40
      if (digit3 == '5') l=l+50
      if (digit3 == '6') l=l+60
      if (digit3 == '7') l=l+70
      if (digit3 == '8') l=l+80
      if (digit3 == '9') l=l+90
      if (digit4 == '1') l=l+1
      if (digit4 == '2') l=l+2
      if (digit4 == '3') l=l+3
      if (digit4 == '4') l=l+4
      if (digit4 == '5') l=l+5
      if (digit4 == '6') l=l+6
      if (digit4 == '7') l=l+7
      if (digit4 == '8') l=l+8
      if (digit4 == '9') l=l+9

      j        = fyear-1950+1

      if (xrd >-998.) then
        sumrd(j)  = sumrd(j)+xrd
        numrd(j)  = numrd(j)+1
      endif
      if (xad >-998.)  sumad(j)  = sumad(j)+xad
      if (capxv>-998.) sumce(j)  = sumce(j)+capxv
      if (at  >-998.)  sumat(j)  = sumat(j)+at
      if (sale>-998.)  sumsl(j)  = sumsl(j)+sale
      if (emp >-998.)  sumem(j)  = sumem(j)+emp

      !
      ! Research and development
      !
      if (xrd>toprd(1,j)) then
        do k=500,2,-1
          toprd(k,j)   = toprd(k-1,j)
          topad(k,j)   = topad(k-1,j)
          topce(k,j)   = topce(k-1,j)
          topat(k,j)   = topat(k-1,j)
          topsl(k,j)   = topsl(k-1,j)
          topem(k,j)   = topem(k-1,j)
          naicsrd(k,j) = naicsrd(k-1,j)
          sicrd(k,j)   = sicrd(k-1,j)
          gvkeyrd(k,j) = gvkeyrd(k-1,j)
        enddo
        toprd(1,j)     = xrd
        topad(1,j)     = xad
        topce(1,j)     = capxv
        topat(1,j)     = at
        topsl(1,j)     = sale
        topem(1,j)     = emp
        naicsrd(1,j)   = l
        sicrd(1,j)     = sic
        gvkeyrd(1,j)   = gvkey
      else
        m              = 0
        do k=2,500
          if ((xrd<=toprd(k-1,j)).and.(xrd>toprd(k,j))) m = k
        enddo

        if (m>0) then

          do k=500,m+1,-1
            toprd(k,j)    = toprd(k-1,j)
            topad(k,j)    = topad(k-1,j)
            topce(k,j)    = topce(k-1,j)
            topat(k,j)    = topat(k-1,j)
            topsl(k,j)    = topsl(k-1,j)
            topem(k,j)    = topem(k-1,j)
            naicsrd(k,j)  = naicsrd(k-1,j)
            sicrd(k,j)    = sicrd(k-1,j)
            gvkeyrd(k,j)  = gvkeyrd(k-1,j)
          enddo
          toprd(m,j)   = xrd 
          topad(m,j)   = xad 
          topce(m,j)   = capxv
          topat(m,j)   = at
          topsl(m,j)   = sale
          topem(m,j)   = emp
          naicsrd(m,j) = l
          sicrd(m,j)   = sic
          gvkeyrd(m,j) = gvkey

        endif
      endif 

    endif

  enddo


  write(7,*) 'Statistics for the Top 500 R&D Spenders'
  write(7,*) '======================================='
  write(7,*)
  write(7,*) 'GVKEY'
  write(7,*) '====='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  write(8,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,i11))') (gvkeyrd(i,j),j=1,64)
    write(8,'(64(1x,i11))') (gvkeyrd(i,j),j=1,64)
  enddo
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'NAICS'
  write(7,*) '====='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,i11))') (naicsrd(i,j),j=1,64)
    write(8,'(64(1x,i11))') (naicsrd(i,j),j=1,64)
  enddo
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'SIC'
  write(7,*) '==='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,i11))') (sicrd(i,j),j=1,64)
    write(8,'(64(1x,i11))') (sicrd(i,j),j=1,64)
  enddo
  write(7,*)
  write(7,*)
  write(8,*)


  write(7,*) 'R&D'
  write(7,*) '==='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,f11.1))') (toprd(i,j),j=1,64)
    write(8,'(64(1x,f11.1))') (toprd(i,j),j=1,64)
  enddo
  write(7,'(64(1x,f11.1))') (sumrd(j),j=1,64)
  write(8,'(64(1x,f11.1))') (sumrd(j),j=1,64)
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'Ad Expenditures'
  write(7,*) '==============='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,f11.1))') (topad(i,j),j=1,64)
    write(8,'(64(1x,f11.1))') (topad(i,j),j=1,64)
  enddo
  write(7,'(64(1x,f11.1))') (sumad(j),j=1,64)
  write(8,'(64(1x,f11.1))') (sumad(j),j=1,64)
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'Capital Expenditures'
  write(7,*) '===================='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,f11.1))') (topce(i,j),j=1,64)
    write(8,'(64(1x,f11.1))') (topce(i,j),j=1,64)
  enddo
  write(7,'(64(1x,f11.1))') (sumce(j),j=1,64)
  write(8,'(64(1x,f11.1))') (sumce(j),j=1,64)
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'Assets'
  write(7,*) '======'
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,f11.1))') (topat(i,j),j=1,64)
    write(8,'(64(1x,f11.1))') (topat(i,j),j=1,64)
  enddo
  write(7,'(64(1x,f11.1))') (sumat(j),j=1,64)
  write(8,'(64(1x,f11.1))') (sumat(j),j=1,64)
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'Sales'
  write(7,*) '====='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,f11.1))') (topsl(i,j),j=1,64)
    write(8,'(64(1x,f11.1))') (topsl(i,j),j=1,64)
  enddo
  write(7,'(64(1x,f11.1))') (sumsl(j),j=1,64)
  write(8,'(64(1x,f11.1))') (sumsl(j),j=1,64)
  write(7,*)
  write(7,*)
  write(8,*)

  write(7,*) 'Employees'
  write(7,*) '========='
  write(7,'(64(1x,i11))') (yrvec(j),j=1,64)
  do i=1,500
    write(7,'(64(1x,f11.1))') (topem(i,j),j=1,64)
    write(8,'(64(1x,f11.1))') (topem(i,j),j=1,64)
  enddo
  write(7,'(64(1x,f11.1))') (sumem(j),j=1,64)
  write(8,'(64(1x,f11.1))') (sumem(j),j=1,64)
  write(7,*)
  write(7,*)
  write(8,*)

  write(8,'(64(1x,i11))') (numrd(j),j=1,64)

end program top500rd