!
      program mkncgpi
!
! This program reads ascii files containing NOAA flux and kp data,
!   and writes a netcdf file for use by tgcm. Yearly ascii data files
!   are obtained from 
!
! ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/KP_AP/
!
! These yearly files are catted together to make the input file 
!   for this code.
!
      implicit none
      integer luin,mxdays,iq,nbad,imo,ida
      parameter(mxdays=10000)
      character*120 flnmin
      real f107d(mxdays),f107a(mxdays),rkp(8,mxdays),ftmp
      integer ii,i,iyd(mxdays),iy2,kp(8),ndays
      integer,external :: idate2iyd
      data luin/20/
c
      flnmin = " "
!     flnmin = "ngdc.1979-2000.dat"
!     flnmin = "ngdc.1979001-2001031.dat"
      flnmin = "ngdc.1979001-2002059.dat"
      open(luin,file=trim(flnmin),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
        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
        f107a(i) = 0.
!       if (mod(i,10).eq.0) write(6,"('i=',i5,' iyd=',i8,' 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=',i5,') <<<')") mxdays
      stop 'mxdays'
 100  continue
      ndays = i-1
      write(6,"('Finished read: ndays=',i5)") 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=',
     +        i5)") i
          f107d(i) = .5*(f107d(i-1)+f107d(i+1))
          nbad = nbad+1
        endif
      enddo
      if (nbad.gt.0) write(6,"('>>> nbad f107d = ',i5)") 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 ',i5)") ii
!
! Replace zeroes at begining and end with closest f107a:
! Front end:
        do i=1,ndays
          if (f107a(i) /= 0.) then
            ftmp = f107a(i)
            exit
          endif
        enddo
        do i=1,ndays
          if (f107a(i) /= 0.) exit
          f107a(i) = ftmp
        enddo
!
! Back end:
        do i=ndays,1,-1
          if (f107a(i) /= 0.) then
            ftmp = f107a(i)
            exit
          endif
        enddo
        do i=ndays,1,-1
          if (f107a(i) /= 0.) exit
          f107a(i) = ftmp
        enddo
!
!       write(6,"('f107a=',/(8f8.2))") (f107a(i),i=1,ndays)
      endif
!
      close(luin)
!
! Write netcdf file:
!     write(6,"('f107d=',/,(6e12.4))") f107d(1:ndays)
      call wrnc(iyd,f107d,f107a,rkp,ndays)
!
      write(6,"('done')")
      stop 'done'
 900  write(6,"('>>> Error opening input file ',a)") trim(flnmin)
      stop 'flnmin'
      end program mkncgpi
!-----------------------------------------------------------------------
      subroutine wrnc(iyd,f107d,f107a,rkp,ndays)
      implicit none
      include "netcdf.inc"
!
! Args:
      integer,intent(in) :: ndays
      real,intent(in) :: f107d(ndays),f107a(ndays),rkp(8,ndays)
      integer,intent(in) :: iyd(ndays)
!
! Local:
      integer :: istat,ncid,ids1(1),ivar1(1),ids2(2)
      integer :: id_ndays,id_kp,idv_iyd,idv_f107d,idv_f107a,idv_kp
      character(len=32) :: ncfile
      character(len=80) :: char80,msspath
      character(len=120) :: char120
!
! msspath will be a global attribute, but this code does not
! do the dispose to mss:
!     write(msspath,"('/TGCM/data/gpi_1979001-2001031.nc')")
      write(msspath,"('/TGCM/data/gpi_1979001-2002059.nc')")
!
! Open new netcdf dataset:
      ncfile = " "
      ncfile = "gpi.nc"
      istat = nf_create(ncfile,NF_CLOBBER,ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'opening new netcdf dataset')
!
! Define ndays dimension:
      istat = nf_def_dim(ncid,"ndays",ndays,id_ndays)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining ndays dimension')
!
! Define kp dimension:
      istat = nf_def_dim(ncid,"nkp",8,id_kp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining kp dimension')
!
! Define variables iyd, f107d, f107a, kp:
!
! iyd:
      ids1(1) = id_ndays
      istat = nf_def_var(ncid,"year_day",NF_INT,1,ids1,idv_iyd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining variable yearday')
      write(char80,"('4-digit year followed by 3-digit day')")
      istat = nf_put_att_text(ncid,idv_iyd,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining long_name of variable year_day')
! f107d:
      istat = nf_def_var(ncid,"f107d",NF_DOUBLE,1,ids1,idv_f107d)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining variable f107d')
      write(char80,"('daily 10.7 cm solar flux')")
      istat = nf_put_att_text(ncid,idv_f107d,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining long_name of variable f107d')
! f107a:
      istat = nf_def_var(ncid,"f107a",NF_DOUBLE,1,ids1,idv_f107a)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining variable f107a')
      write(char80,"('81-day average 10.7 cm solar flux')")
      istat = nf_put_att_text(ncid,idv_f107a,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining long_name of variable f107a')
! kp:
      ids2(1) = id_kp
      ids2(2) = id_ndays
      istat = nf_def_var(ncid,"kp",NF_DOUBLE,2,ids2,idv_kp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining variable kp')
      write(char80,"('3-hourly kp index')")
      istat = nf_put_att_text(ncid,idv_kp,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining long_name of variable kp')
!
! Global file attributes:
!
! Title:
      write(char80,"('Geophysical Indices, obtained from NGDC')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"title",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute title')
!
! Beginning yearday:
! (make it the 3rd day, so model interpolation can use 1st 2 days)
      ivar1(1) = iyd(3)
      istat = nf_put_att_int(ncid,NF_GLOBAL,"yearday_beg",NF_INT,1,
     |  ivar1)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute yearday_beg')
!
! Ending yearday:
! (make it day ngpidays-2, so model interpolation can use last 2 days)
      ivar1(1) = iyd(ndays-2)
      istat = nf_put_att_int(ncid,NF_GLOBAL,"yearday_end",NF_INT,1,
     |  ivar1)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute yearday_end')
!
! mss path (this code does not dispose to mss):
      istat = nf_put_att_text(ncid,NF_GLOBAL,"ncar_mss_path",
     |  len_trim(msspath),trim(msspath))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute ncar_mss_path')
!
! NGDC URL where ascii data was acquired:
      write(char80,"('ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA',
     |  '/INDICES/KP_AP/')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"data_source_url",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute data_source_url')
!
! Local path to program that wrote the netcdf file:
      write(char80,"('/home/foster/tgcminp/nc_gpi/mkncgpi.f')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"hao_file_write_source",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute hao_file_write_source')
!
! Description:
      write(char120,"('Yearly ascii data files obtained from ',
     |  'data_source_url; nc file written by hao_file_write_source.')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"info",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining global attribute info')
!
! Exit define mode:
      istat = nf_enddef(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error return from nf_enddef')
!
! Give values to variables:
!
! iyd:
      istat = nf_put_var_int(ncid,idv_iyd,iyd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error giving values to variable iyd')
! f107d:
      istat = nf_put_var_double(ncid,idv_f107d,f107d)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error giving values to variable f107d')
! f107a:
      istat = nf_put_var_double(ncid,idv_f107a,f107a)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error giving values to variable f107a')
! kp:
      istat = nf_put_var_double(ncid,idv_kp,rkp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error giving values to variable kp')

      istat = nf_close(ncid)
      end subroutine wrnc
!-----------------------------------------------------------------------
      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 
!-------------------------------------------------------------------
      subroutine handle_ncerr(istat,msg)
!
! Handle a netcdf lib error:
!
      integer,intent(in) :: istat
      character(len=*),intent(in) :: msg
      include "netcdf.inc"
!
      write(6,"(/72('-'))")
      write(6,"('>>> Error from netcdf library:')")
      write(6,"(a)") trim(msg)
      write(6,"('istat=',i5)") istat
      write(6,"(a)") nf_strerror(istat)
      write(6,"(72('-')/)")
      return
      end subroutine handle_ncerr