!
      module geodata_module

      implicit none

      type geo_type
         integer :: ntimes                      ! number of data points
         integer(kind=8),allocatable :: nsec(:) ! seconds from the beginning of a year
         real, allocatable :: value(:)          ! geo data value
      end type geo_type
      type(geo_type) :: 
     |     geo_f107,                ! f107
     |     geo_f107a,               ! f107a
     |     geo_power,               ! hemisphere power
     |     geo_ctpoten,             ! cross tail potential
     |     geo_byimf                ! IMF By component

      contains
!-----------------------------------------------------------------------
      subroutine rd_geodata(geo_data)

! Args
      character(len=80),intent(in) :: geo_data
! locals
      integer :: i,n,istat,seq
      integer :: nx=3
      real,allocatable :: time(:,:)
!
! external
      integer(kind=8),external :: mtime_to_nsec

      open(12,file=trim(geo_data),status='old')
!
! F107
!     
      read(12,"(i7)") geo_f107%ntimes
      n =geo_f107%ntimes

      if(allocated(time)) deallocate(time)
      allocate(time(nx,n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_f107%nsec)) deallocate(geo_f107%nsec)
      allocate(geo_f107%nsec(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_f107%value)) deallocate(geo_f107%value)
      allocate(geo_f107%value(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      do i=1, n
          read(12,"(i7,4f6.1)") seq,time(1:3,i),geo_f107%value(i)
          ! f107 must be positive
          if (geo_f107%value(i) <0.) then
             call shutdown('input f107 is negative')
          endif
          geo_f107%nsec(i) = mtime_to_nsec(int(time(1:3,i))) 
      enddo 
!
! F107a
!
      read(12,"(i7)") geo_f107a%ntimes
      n =geo_f107a%ntimes

      if(allocated(time)) deallocate(time)
      allocate(time(nx,n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_f107a%nsec)) deallocate(geo_f107a%nsec)
      allocate(geo_f107a%nsec(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_f107a%value)) deallocate(geo_f107a%value)
      allocate(geo_f107a%value(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      do i=1, n
          read(12,"(i7,4f6.1)") seq,time(1:3,i),geo_f107a%value(i)
          ! f107a must be positive
          if (geo_f107a%value(i) <0.) then
             call shutdown('input f107a is negative')
          endif
          geo_f107a%nsec(i) = mtime_to_nsec(int(time(1:3,i))) 
      enddo 
!
! Hemispheric Power
!
      read(12,"(i7)") geo_power%ntimes
      n =geo_power%ntimes

      if(allocated(time)) deallocate(time)
      allocate(time(nx,n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_power%nsec)) deallocate(geo_power%nsec)
      allocate(geo_power%nsec(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_power%value)) deallocate(geo_power%value)
      allocate(geo_power%value(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      do i=1, n
          read(12,"(i7,4f6.1)") seq,time(1:3,i),geo_power%value(i)
          ! power must be positive
          if (geo_power%value(i) <0.) then
             call shutdown('input hemisphere power is negative')
          endif
          geo_power%nsec(i) = mtime_to_nsec(int(time(1:3,i))) 
      enddo 
!
! Cross Polar cap Potential
!
      read(12,"(i7)") geo_ctpoten%ntimes
      n =geo_ctpoten%ntimes

      if(allocated(time)) deallocate(time)
      allocate(time(nx,n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_ctpoten%nsec)) deallocate(geo_ctpoten%nsec)
      allocate(geo_ctpoten%nsec(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_ctpoten%value)) deallocate(geo_ctpoten%value)
      allocate(geo_ctpoten%value(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      do i=1, n
          read(12,"(i7,4f6.1)") seq,time(1:3,i),geo_ctpoten%value(i)
          ! ctpoten must be positive
          if (geo_ctpoten%value(i) <0.) then
             call shutdown('input cross tail potential is negative')
          endif
          geo_ctpoten%nsec(i) = mtime_to_nsec(int(time(1:3,i))) 
      enddo 
!
! IMF By component
!
      read(12,"(i7)") geo_byimf%ntimes
      n =geo_byimf%ntimes

      if(allocated(time)) deallocate(time)
      allocate(time(nx,n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_byimf%nsec)) deallocate(geo_byimf%nsec)
      allocate(geo_byimf%nsec(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      if(allocated(geo_byimf%value)) deallocate(geo_byimf%value)
      allocate(geo_byimf%value(n),stat=istat)
      if (istat/= 0) print*,'rd_geodata: fail to allocate'

      do i=1, n
          read(12,"(i7,4f6.1)") seq,time(1:3,i),geo_byimf%value(i)
          geo_byimf%nsec(i) = mtime_to_nsec(int(time(1:3,i))) 
      enddo 
!
! Close data file
!
      close(12)
      write(6,*)'--------------------------------------'
      write(6,*)'   --  Finished geo_cond reading  --   '
      write(6,*)'--------------------------------------'
   
      end subroutine rd_geodata
!-----------------------------------------------------------------------
      subroutine get_geodata(nsec,f107,f107a,ctpoten,power,byimf) 
!
! Args:
!     
      integer(kind=8),intent(in) :: nsec   
      real,intent(out) :: f107,f107a,ctpoten,power,byimf
! locals
!
      integer :: ind,ind1,ind2

!
! external
!     
      integer,external::long_bsearch
      real,external :: finterp_bigints

!
! f107
!
      ind=long_bsearch(geo_f107%nsec(:),1,geo_f107%ntimes,nsec)
      if (ind ==-1) then
         call shutdown('There is no f107 data for the required time')
      endif

      ind1=ind
      ind2=ind+1
      if (ind2 > geo_f107%ntimes) then
         f107 = geo_f107%value(ind1)
      else
         f107 = finterp_bigints(geo_f107%value(ind1),
     |          geo_f107%value(ind2),geo_f107%nsec(ind1),
     |          geo_f107%nsec(ind2),nsec)
      endif
!
! f107a
!
      ind=long_bsearch(geo_f107a%nsec(:),1,geo_f107a%ntimes,nsec)
      if (ind ==-1) then
         call shutdown('There is no f107a data for the required time')
      endif

      ind1=ind
      ind2=ind+1
      if (ind2 > geo_f107a%ntimes) then
         f107a = geo_f107a%value(ind1)
      else
         f107a = finterp_bigints(geo_f107a%value(ind1),
     |          geo_f107a%value(ind2),geo_f107a%nsec(ind1),
     |          geo_f107a%nsec(ind2),nsec)
      endif
!
! power
!
      ind=long_bsearch(geo_power%nsec(:),1,geo_power%ntimes,nsec)
      if (ind ==-1) then
         call shutdown('There is no power data for the required time')
      endif

      ind1=ind
      ind2=ind+1
      if (ind2 > geo_power%ntimes) then
         power = geo_power%value(ind1)
      else
         power = finterp_bigints(geo_power%value(ind1),
     |          geo_power%value(ind2),geo_power%nsec(ind1),
     |          geo_power%nsec(ind2),nsec)
      endif
!
! cross tail potential
!
      ind=long_bsearch(geo_ctpoten%nsec(:),1,geo_ctpoten%ntimes,nsec)
      if (ind ==-1) then
         call shutdown('There is no ctpoten data for the required time')
      endif

      ind1=ind
      ind2=ind+1
      if (ind2 > geo_ctpoten%ntimes) then
         ctpoten = geo_ctpoten%value(ind1)
      else
         ctpoten = finterp_bigints(geo_ctpoten%value(ind1),
     |          geo_ctpoten%value(ind2),geo_ctpoten%nsec(ind1),
     |          geo_ctpoten%nsec(ind2),nsec)
      endif
!
! byimf
!
      ind=long_bsearch(geo_byimf%nsec(:),1,geo_byimf%ntimes,nsec)
      if (ind ==-1) then
         call shutdown('There is no byimf data for the required time')
      endif

      ind1=ind
      ind2=ind+1
      if (ind2 > geo_byimf%ntimes) then
         byimf = geo_byimf%value(ind1)
      else
         byimf = finterp_bigints(geo_byimf%value(ind1),
     |          geo_byimf%value(ind2),geo_byimf%nsec(ind1),
     |          geo_byimf%nsec(ind2),nsec)
      endif
      end subroutine get_geodata
!-----------------------------------------------------------------------------
 
      end module geodata_module