program FITSlist

  implicit none

  integer :: status, bitpix, naxis, naxes(2)
  ! status ... FITS status (0=no error)
  ! naxis .. number of axes in image (we require =2)
  ! naxes .. dimension of the image (2-element array) 

  integer :: i,j,blocksize,pcount,gcount,minvalue
  logical :: extend, simple,anyf
  ! required by cFITSIO

  real, dimension(:,:), allocatable :: d
  ! data matrix

  character(len=666) :: name = 'image.fits'
  ! name .. fill with name of the image to open

  status = 0
  call ftopen(25,name,0,blocksize,status)
  call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
  allocate(d(naxes(1),naxes(2)))
  call ftg2de(25,1,minvalue,naxes(1),naxes(1),naxes(2),d,anyf,status)
  call ftclos(25,status)

  ! print value of a random pixel
  write(*,*) '# d(1,1)=',d(1,1)

  ! print last tree values of first row
  write(*,*) '# d(1,-10:)=',d(1,size(d,2)-3:)
  
  ! print of a submatrix with indexes
  do i = 1,10
     do j = 50,60
        write(*,*) i,j,d(i,j)
     end do
  end do

end program FITSlist
