program rd_ascii_v2
c
c   f77 rd_ascii_v2.f -o rd_ascii_v2.x
c     
c   routine to read ascii Version 2 data set of gridded data fields
c   user will have to specify path (where data resides) in 
c   subroutine rd_ascii
c
c   the data is stored using the following convention
c   i - longitude, lon(i) = 140. + (i-1)*1. where i goes from 1, 41
c   j - latitude,  lat(j) = -10. + (i-1)*1. where j goes from 1, 21
c   k - pressure level, k=1 contains the surface data
c       plev(k) = 1000. - (k-2)*25. where k goes from 2, 41

      integer    i, j, k
      integer    imax,    jmax,    kmax
      parameter (imax=41, jmax=21, kmax=41)
      real wlon, slat, gridsp, plev(kmax)

c   surface pressure reduced to sea level

      real psl(imax,jmax)

c   upper air fields

      real zec(imax,jmax,kmax)
      real z(imax,jmax,kmax), q(imax,jmax,kmax)
      real t(imax,jmax,kmax), u(imax,jmax,kmax), v(imax,jmax,kmax)
      real d(imax,jmax,kmax), w(imax,jmax,kmax)

      real lat(jmax), lon(imax)
      character  date*8

c
      parameter (npds = 480 )       
c     parameter (npds = 1 )       

      integer nday(12), iyr, imo, idy, ihr 

      data nday /31,29,31,30,31,30,31,31,30,31,30,31/

      integer iflgcord

c   *************************************************************************

      date = '92110100'
      read(date,'(4i2.2)') iyr, imo, idy, ihr
      itim = 0

c   define pressure levels for data
c   k = 1 contains the surface data

      dp = 25.
      do k=2,kmax
         plev(k) = 1000. - (k-2)*dp
      enddo

c   define lat/lon coordinates of the data

      wlon = 140.
      slat = -10.
      gridsp = 1.

c   define latitude of all points
      do j=1,jmax
         lat(j) = slat + gridsp*(j-1)
      enddo

c   define longitude of all points
      do i=1,imax
         lon(i) = wlon + gridsp*(i-1)
      enddo

   5  continue
      itim = itim + 1

c   call rd_ascii to read in ascii grids at time: date

      call rd_ascii(date, psl, zec, z, t, q, u, v, w, d, 
     2              imax, jmax, kmax)
c
c   define a new date
c
      if (itim .eq. npds) go to 100
      ihr = ihr + 6
      if (ihr .eq. 24) then
        ihr = 0
        idy = idy + 1
        if (idy .gt. nday(imo)) then
          idy = 1
          imo = imo + 1
          if (imo .gt. 12) then
            imo = 1
            iyr = iyr + 1
          endif
        endif
      endif
      write(date,'(4i2.2)') iyr, imo, idy, ihr

      go to 5   

  100 continue

      end

      subroutine rd_ascii(date, psl, zec, z, t, q, u, v, w, d, 
     2               ig, jg, kg)

      real psl(ig,jg), zec(ig,jg,kg)
      real z(ig,jg,kg), t(ig,jg,kg), q(ig,jg,kg) 
      real u(ig,jg,kg), v(ig,jg,kg), w(ig,jg,kg) 
      real d(ig,jg,kg)

      character date*8, path*44, header*50

c
c   open data files
c
      path = '/home/paulc/toga/data/grid/ascii_lsa_v2/tmp/'

      open(19, file=path//'zec.'//date, status='old', form='formatted')
      open(20, file=path//'z.'//date, status='old', form='formatted')
      open(21, file=path//'t.'//date, status='old', form='formatted')
      open(22, file=path//'q.'//date, status='old', form='formatted')
      open(23, file=path//'u.'//date, status='old', form='formatted')
      open(24, file=path//'v.'//date, status='old', form='formatted')
      open(25, file=path//'d.'//date, status='old', form='formatted')
      open(26, file=path//'w.'//date, status='old', form='formatted')
      open(27, file=path//'p.'//date, status='old', form='formatted')

c   ********** read in  sea level pressure  ***************************

      npts = ig*jg

      iu = 27
      read(iu,'(a50)') header 
      do i=1,ig
         read(iu,30) (psl(i,j),j=1,jg)
      enddo
      call fndmima(npts, psl, pmin, pmax)

c   ********** read in ECMWF height data (m)  ****************************

      npts = ig*jg*kg

      iu = 19
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,35) (zec(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, zec, zecmin, zecmax)

c   ********** read in height data (m)  **********************************

      iu = 20
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,35) (z(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, z, zmin, zmax)

c   ********** read in temp. data (C)  **********************************

      iu = 21
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,50) (t(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, t, tmin, tmax)

c   ********** read in specific humidity data (g/kg)  *******************

      iu = 22
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,65) (q(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, q, qmin, qmax)

c   ********** read in zonal wind component (m/s)  *******************

      iu = 23
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,50) (u(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, u, umin, umax)

c   ********** read in meridional wind component (m/s)  *******************

      iu = 24
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,50) (v(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, v, vmin, vmax)

c   ********** read divergence (1/s)*1.e6  ************************

      iu = 25
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,50) (d(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, d, dmin, dmax)

c   ********** read in  omega (mb/hr)  ************************

      iu = 26
      do k=1,kg
         read(iu,'(a50)') header
         do i=1,ig
            read(iu,60) (w(i,j,k),j=1,jg)
         enddo
      enddo
      call fndmima(npts, w, wmin, wmax)
      write(6,70) date, pmin, pmax, zmin, zmax, tmin, tmax, qmin, qmax,
     2                  umin, umax, vmin, vmax, dmin, dmax, wmin, wmax

c   close files

      do iu=20,27
         close(iu)
      enddo

c  format statements *********************************************

   30 format(21f6.0)
   35 format(21f7.0)
   50 format(21f6.1)
   60 format(21f6.2)
   65 format(21e10.3)
   70 format(a8,1x,2f6.0,2f7.0,14f6.1)

      return
      end

      subroutine fndmima(npts, fld, xmin, xmax)
c
c   finds min and max values in fld
c
      real fld(npts), xmin, xmax
      real       big
      parameter (big = 1.e99)

c
      xmin = +big
      xmax = -big
      do i=1,npts
         xmin = min(xmin, fld(i))
         xmax = max(xmax, fld(i))
      enddo

      return
      end