program d_dy_read
c
c This program reads TAO/TRITON and PIRATA ascii-format density
c   files, for example d15n38w_10m.ascii. It creates an array called 
c   d, which is evenly spaced in time, and an array called iqual
c   which contains the data quality for each depth.
c
c You can easily adapt this program to your needs.
c
c Programmed by Dai McClurg, NOAA/PMEL/OCRD, August 1999
c
      implicit none
c
      integer nz, nt
      parameter(nz = 30, nt = 10000)
c
      integer k, n, m
c
      integer nblock, nk, ndep, nn, ntime, n1, n2
c
      integer kdep(nz), idep(nz)
      integer iqual(nz,nt), isrc(nz,nt), idate(nt), ihms(nt)
c
      real flag, depth(nz), d(nz,nt)
c
      character infile*80, header*132, aqual*50, frmt*160
c
c .......................................................................
c
      write(*,*) ' Enter the input density file name'
      read(*,'(a)') infile
c
      open(1,file=infile,status='old',form='formatted')
c 
c Read total number of days, depths and blocks of data.
c
      read(1,10) ntime, ndep, nblock
   10 format(49x,i5,7x,i3,8x,i3)
c
      write(*,*) ntime, ndep, nblock
c
c Read the missing data flag
c
      read(1,20) flag
   20 format(56x,f7.3)
c
      write(*,*) flag
c
c  Initialize t array to flag and iqual array to 5.
c
      do k = 1, nz
        do n = 1, nt
          d(k,n) = flag
          iqual(k,n) = 5
           isrc(k,n) = 0
        enddo
      enddo
c
c Read the data
c
      do m = 1, nblock
        read(1,30) n1, n2, nn, nk
        call blank(frmt)
        write(frmt,140) nk
        read(1,frmt) (kdep(k),k=1,nk)
        call blank(frmt)
        write(frmt,150) nk
        read(1,frmt) (idep(kdep(k)),k=1,nk)
        do k = 1, nk
          depth(kdep(k)) = real(idep(kdep(k)))
        enddo
        read(1,'(a)') header
        call blank(aqual)
        call blank(frmt)
        write(frmt,160) nk, nk
        do n = n1, n2
          read(1,frmt) idate(n), ihms(n), 
     .       (d(kdep(k),n),k=1,nk), aqual(1:nk),
     .       (isrc(kdep(k),n),k=1,nk) 
          do k = 1, nk
            if(aqual(k:k) .eq. 'C') then
              iqual(kdep(k),n) = -9 ! See AAA below
            else
              read(aqual(k:k),'(i1)') iqual(kdep(k),n)
            endif
          enddo
        enddo
      enddo
c
   30 format(50x,i6,3x,i6,1x,i6,7x,i3)
  140 format('(15x,',i3,'i7)')
  150 format('(15x,',i3,'i7)')
  160 format('(1x,i8,1x,i4,1x,',i3,'f7.3,1x,a',i3,',1x,',i3,'i1)')
c
      close(1)
c
c Write out the depth, data, and quality arrays to the 
c   standard output. 
c
      write(*,*) 'depth = ', (depth(k),k=1,ndep)
c
c For some files this statement may be too long for your max output
c   record length on your terminal. If so, comment out these lines.
c
      nk = ndep
c
      call blank(frmt)
      write(frmt,70) ndep
c
      do n = 1, ntime
        write(*,frmt) idate(n), ihms(n),(d(k,n),k=1,ndep),
     .                 (iqual(k,n),k=1,ndep), 
     .                  (isrc(k,n),k=1,ndep), n
      enddo
c
   70 format('(1x,i8,1x,i6,1x,',i3,'f7.3,1x,',i3,'i2,1x,',i3,'i2,i7)')
c
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine blank(string)
c
c blank out the string from 1 to its declared length
c
      character*(*) string
c
      integer i
c
      do i = 1, len(string)
        string(i:i) = ' '
      enddo
c
      return
      end
c
c AAA
c
c A quality of -9 means to check with freitag@pmel.noaa.gov
c   for information. There may be a technical report which
c   gives further information.