program readgrd
c
c   program to read gridded ascii fields at surface and upper levels
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, elon, slat, nlat, gridsp, plev(kmax)

c   surface pressure and pressure reduced to sea level

      real ps(imax,jmax), psl(imax,jmax)

c   upper air fields

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

      real lat(jmax), lon(imax)
      character header*50, tmp*8

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.
      elon = 180.
      slat = -10.
      nlat = 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
C
      print *, 'Specify date (e.g., 92120100)...'  
      read(5,'(a8)') tmp
      print *, 'opening files for time: ', tmp

      open(20,file='u.'//tmp, status='unknown', form='formatted')
      open(21,file='v.'//tmp, status='unknown', form='formatted')
      open(22,file='q.'//tmp, status='unknown', form='formatted')
      open(23,file='t.'//tmp, status='unknown', form='formatted')
      open(24,file='z.'//tmp, status='unknown', form='formatted')
      open(25,file='sfc.'//tmp, status='unknown', form='formatted')

c   ********** read in surface fields ***********************************

c   read in surface pressure

         read(25,'(a50)') header 
         do i=1,imax
             read(25,40) (ps(i,j),j=1,jmax)
         enddo

c   read in sea level pressure (mb)

         read(25,'(a50)') header
         do i=1,imax
             read(25,40) (psl(i,j),j=1,jmax)
         enddo

c   read in surface height (m)     

         read(25,'(a50)') header
         do i=1,imax
             read(25,40) (z(i,j,1),j=1,jmax)
         enddo

c   read in surface temperature (m)     

         read(25,'(a50)') header
         do i=1,imax
             read(25,50) (t(i,j,1),j=1,jmax)
         enddo

c   read in surface specific humidity (g/kg) 

         read(25,'(a50)') header
         do i=1,imax
             read(25,60) (q(i,j,1),j=1,jmax)
         enddo

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

         read(25,'(a50)') header
         do i=1,imax
             read(25,50) (u(i,j,1),j=1,jmax)
         enddo

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

         read(25,'(a50)') header
         do i=1,imax
             read(25,50) (v(i,j,1),j=1,jmax)
         enddo

c   ********** read in upper level fields  ***********************************

c   read in zonal wind component (m/s) 

      do k=2,kmax
         read(20,'(a50)') header
         do i=1,imax
             read(20,50) (u(i,j,k),j=1,jmax)
         enddo
      enddo

c   read in meridional wind component (m/s)

      do k=2,kmax
         read(21,'(a50)') header
         do i=1,imax
             read(21,50) (v(i,j,k),j=1,jmax)
         enddo
      enddo

c   read in specific humidity (gr/kg) 

      do k=2,kmax
         read(22,'(a50)') header
         do i=1,imax
             read(22,60) (q(i,j,k),j=1,jmax)
         enddo
      enddo

c   read in temperature (C) 

      do k=2,kmax
         read(23,'(a50)') header
         do i=1,imax
             read(23,50) (t(i,j,k),j=1,jmax)
         enddo
      enddo

c   read in geopotential height (m) 

      do k=2,kmax
         read(24,'(a50)') header
         do i=1,imax
             read(24,40) (z(i,j,k),j=1,jmax)
         enddo
      enddo

   40 format(21i6)
   50 format(21f6.1)
   60 format(21f6.2)

      end