
      program wrfile
c
c Read ascii data file with iyd,f107d,f107a,rkp(8) (file written
c   by ~/tgcminp/getndcs), and write unformatted direct access file
c   for later use by time-gcm
c There are 4779 records for days 79001 to 92031
c
c 9/28/94: Updating with ngdc data (see ~/tgcminp/ngdc) to 94080
c          There are 5559 records for days 79001 to 94080
c          Wrote /FOSTER/tgcminp/79001-94080.dat
c 9/28/94: Added option for calculating missing f107a values
c
      parameter(mxrec=10000)
      parameter(irec=11*8)
      dimension f107d(mxrec),f107a(mxrec),rkp(8,mxrec),iyd(mxrec)
      character*56 filein,fileout
      integer :: iy2,imo,ida,kp(8)
      data luin/10/, luout/11/
      data filein
     | /'ngdc.1979001-2001031.dat                                '/
!    + /'79001-99111.dat                                         '/
!    + /'79001-97212.dat                                         '/
c    + /'79001-94080.dat                                         '/
c    + /'79001-92031.dat                                         '/
c    + /'/d/foster/tgcminp/getndcs/79001-92031.dat               '/
      data fileout
     | /'1979001-2001031.bin                                     '/
!    + /'79001-99111.bin                                         '/
!    + /'79001-97212.bin                                         '/
c    + /'79001-94080.bin                                         '/
c    + /'79001-92031.bin                                         '/
c    + /'/d/foster/tgcminp/getndcs/wrfile/79001-92031.bin        '/
c
      lenin = lenstr(filein)
      lenout = lenstr(fileout)
      open(luin,file=filein(1:lenin),status='OLD',err=800)
      open(luout,file=fileout(1:lenout),access='DIRECT',
     +  form='UNFORMATTED',recl=irec,err=801)
c
c Read ascii file:
c
      nrec = 0
      do i=1,mxrec
!       read(luin,"(i5,2f8.1,8f5.1)",end=100) iyd(i),f107d(i),
!    +    f107a(i),(rkp(ii,i),ii=1,8)
        read(luin,"(3i2,6x,8i2,37x,f5.1)",end=100) iy2,imo,ida,kp,
     +    f107d(i)
        iyd(i) = idate2iyd(imo,ida,iy2,i)
        do ii=1,8
          rkp(ii,i) = float(kp(ii)/10)+float(kp(ii)-kp(ii)/10*10)*.1
        enddo

!       write(6,"('i=',i4,' iyd=',i8,' f107d=',f6.2)") i,iyd(i),
!    |    f107d(i)

        f107a(i) = 0.
        nrec = nrec+1
      enddo
      write(6,"('>>> need to increase mxrec (did not reach eof',
     +  ' with mxrec=',i6,') <<<')") mxrec
      stop 'mxrec'
 100  continue
      write(6,"('EOF: nrec=',i6)") nrec
c
c Calculate 81-day ave if necessary and possible:
c
      nf107afix = 0
      nf107abad = 0
      do i=1,nrec
        if (f107a(i).le.0.) then
          if (i.ge.41.and.i+40.le.nrec) then
            nf107afix = nf107afix+1
            f107a(i) = 0.
            do ii=i-40,i+40
              f107a(i) = f107a(i) + f107d(ii)
            enddo
            f107a(i) = f107a(i)/81.
          else
            nf107abad = nf107abad+1
          endif
        endif
c
c Write to new direct access file:
c
        write(luout,rec=i) iyd(i),f107d(i),f107a(i),(rkp(ii,i),ii=1,8)

!       write(6,"('write: i=',i4,' iyd=',i8,' f107d,a=',2f8.2)")
!    |    i,iyd(i),f107d(i),f107a(i)

      enddo
      if (nf107afix.gt.0) write(6,"('nf107afix = ',i6)") nf107afix
      if (nf107abad.gt.0) write(6,"('>>> nf107abad = ',i6)") nf107abad
      stop 'done'
 800  continue
      write(6,"('Error opening file ',a)") filein(1:lenin)
      stop 'filein'
 801  continue
      write(6,"('Error opening file ',a)") fileout(1:lenout)
      stop 'fileout'
      stop
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
!-----------------------------------------------------------------------
      integer function idate2iyd(imo,ida,iyr,irec)
      implicit none
!
! 2/16/00: Return 7-digit iyd, i.e., 4-digit leading year followed
!          by 3-digit day, e.g., 1999001,2000365
!
! Args:
      integer,intent(in) :: imo,ida,iyr,irec
!
! Local:
      integer :: i,ndmon(12),iday,iyear
      integer,external :: ichdate
!
!                 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/
!
! given iyr (2-digit), imo, ida, convert to iyd (yyddd):
!
      if (ichdate(imo,ida,iyr).ne.0) then
        write(6,"('>>> idate2iyd: bad idate iyr,imo,ida=',3i3,
     |    ' irec=',i6)") iyr,imo,ida,irec
        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
      if (iyr < 50) then ! assume 21st century
        iyear = 2000+iyr
      else               ! assume 20th century
        iyear = 1900+iyr
      endif
      idate2iyd = iyear*1000+iday+ida
      return
      end
!-----------------------------------------------------------------------
      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.0.or.iy.gt.2050) 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 
