program wrascii
c
c Write ascii file of iyd,f107d,f107a,kp(8) from data obtained from
c ngdc (see http://www.ngdc.noaa.gov/ngdc.html). This file will in
c turn be read by ~/tgcminp/getndcs/wrfile which will write the 
c direct access file used by tgcm17.
c
c 9/94: Wrote this code in order to write new direct access file
c       for tgcm. Current file is /FOSTER/tgcminp/79001-92031.gpi.dat.
c       New file will be /FOSTER/tgcminp/79001-93365.gpi.dat
c
      implicit none
      integer luin,luout,mxdays,iq,nbad,imo,ida
      parameter(mxdays=1500)
      character*120 flnmin,flnmout
      real f107d(mxdays),f107a(mxdays),rkp(8,mxdays)
      integer ii,i,iyd(mxdays),iy2,kp(8),ndays,lflnmin,lflnmout
      integer lenstr,idate2iyd
      data luin/20/, luout/21/
c
      call clearstr(flnmin)
      call clearstr(flnmout)
c     write(flnmin,"('kp_ap.1992.dat')") 
c     write(flnmin,"('kp_ap.1993.dat')") 
c     write(flnmin,"('kp_ap.1994.dat')") 
c     write(flnmin,"('kp_ap.1992-1994.dat')") 
c     write(flnmin,"('kp_ap.1994-1997.dat')") 
      write(flnmin,"('kp_ap.1997-1999.dat')") 
      write(flnmout,"('/d/foster/tgcminp/ngdc/ngdc.dat')")
      lflnmin = lenstr(flnmin)
      open(luin,file=flnmin(1:lflnmin),status='OLD',err=900)
      do i=1,mxdays
        read(luin,"(3i2,6x,8i2,37x,f5.1,i1)",end=100) iy2,imo,ida,kp,
     +    f107d(i),iq
c       iyd(i) = idate2iyd(imo,ida,1900+iy2)
        iyd(i) = idate2iyd(imo,ida,iy2)
        do ii=1,8
          rkp(ii,i) = float(kp(ii)/10)+float(kp(ii)-kp(ii)/10*10)*.1
        enddo
        f107a(i) = 0.
        if (mod(i,10).eq.0) write(6,"('i=',i3,' iyd=',i5,' f107d=',
     +    f5.1,' rkp=',8f5.1)") i,iyd(i),f107d(i),(rkp(ii,i),ii=1,8)
      enddo
      write(6,"('>>> need to increase mxdays (did not reach eof',
     +  ' with mxdays=',i4,') <<<')") mxdays
      stop 'mxdays'
 100  continue
      ndays = i-1
      write(6,"('ndays=',i3)") ndays
c
c f107d == 0. if no observation was made (iq=3). Fill these in
c by simple interpolation if possible:
c
      nbad = 0
      do i=2,ndays-1
        if (f107d(i).le.0.) then
          if (f107d(i-1).le.0..or.f107d(i+1).le.0.)
     +      write(6,"('>>> warning: 2 or more bad f107d at i=',
     +        i3)") i
          f107d(i) = .5*(f107d(i-1)+f107d(i+1))
          nbad = nbad+1
        endif
      enddo
      if (nbad.gt.0) write(6,"('>>> nbad f107d = ',i3)") nbad
c
c Find as many f107a (81-day ave) values as possible:
c (leave rest zero)
c
      if (ndays.ge.81) then
        do i=41,ndays
          if (i+40.gt.ndays) then
            ii = i-1
            goto 101
          endif
          f107a(i) = 0.
          do ii=i-40,i+40
            f107a(i) = f107a(i) + f107d(ii)
          enddo
          f107a(i) = f107a(i)/81.
        enddo
 101    continue
        write(6,"('Calculated 81-day ave f107a for days ',
     +    '41 to ',i3)") ii
      endif
c
c Write new file:
c
      lflnmout = lenstr(flnmout)
      open(luout,file=flnmout(1:lflnmout),status='NEW',err=901)
      write(6,"('Opened new file ',a)") flnmout(1:lflnmout)
      do i=1,ndays
        write(luout,"(i5,2f8.1,8f5.1)") iyd(i),f107d(i),f107a(i),
     +    (rkp(ii,i),ii=1,8)
      enddo
      close(luout)
      close(luin)
      write(6,"('done')")
      stop 'done'
 900  write(6,"('>>> Error opening input file ',a)") flnmin(1:lflnmin)
      stop 'flnmin'
 901  write(6,"('>>> Error opening new file ',a)") flnmout(1:lflnmout)
      stop 'flnmout'
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function lenstr(str)
      character*(*) str
c
c Return index to last non-blank char in str
c
      length = len(str)
      do i=length,1,-1
        if (str(i:i).ne.' ') then
          lenstr = i
          return
        endif
      enddo
      lenstr = 0
      return
      end
c
      subroutine clearstr(str)
c 
c Set given string to all blanks 
c 
      character*(*) str
      length = len(str)
      do i=1,length
        str(i:i) = ' '
      enddo
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function idate2iyd(imo,ida,iyr)
      dimension ndmon(12)
c                 J  F  M  A  M  J  J  A  S  O  N  D
      data ndmon/31,28,31,30,31,30,31,31,30,31,30,31/
c
c given iyr (2-digit), imo, ida, convert to iyd (yyddd):
c
      if (ichdate(imo,ida,iyr).ne.0) then
        write(6,"('>>> idate2iyd: bad idate iyr,imo,ida=',3i3)")
     +    iyr,imo,ida
        idate2iyd = 0
        return
      endif
      ndmon(2) = 28
      if (mod(iyr,4).eq.0) ndmon(2) = 29
      iday = 0
      do i=1,imo-1
        iday = iday+ndmon(i)
      enddo
      idate2iyd = iyr*1000+iday+ida
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function ichdate(im,id,iy)
      dimension ndmon(12)
c                 J  F  M  A  M  J  J  A  S  O  N  D
      data ndmon/31,28,31,30,31,30,31,31,30,31,30,31/
c
c Check validity of date im/id/iy
c
      ichdate = 0
      if (im.lt.1.or.im.gt.12) then
        write(6,"('>>> ichdate: bad month=',i3)") im
        ichdate = 1
        return
      endif
      if (iy.lt.1.or.iy.gt.1999) then
        write(6,"('>>> ichdate: bad year=',i3)") iy
        ichdate = 1
        return
      endif
      ndmon(2) = 28
      if (mod(iy,4).eq.0) ndmon(2) = 29
      if (id.lt.1.or.id.gt.ndmon(im)) then
        write(6 ,"('>>> ichdate: bad day=',i4,' (month=',i2,')')")
     +    id,im
        ichdate = 1
        return
      endif
      return
      end