program read_ascii_grd   

c   sample program to read in gridded data in ascii format 
c   on the SCSMEX/GAME grid 

c   g77 read_ascii_grd.f -o read_ascii_grd.x

      character date*8

c   xflag sets values below ground to missing value flag 
c   if equal to .true.

      logical    xflag
      parameter (xflag = .true.)

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

c   open dates file for times to read in
c   "date" is an 8 character variable with the format (yymmddhh)
c   where: yy - 98
c          mm - 05 or 06
c          dd - day of month
c          hh - hour in UTC (00, 06, 12 or 18) 

      open(1, file='dates', status='old', form='formatted')

c   call init_grid to define grid coordinates

      call init_grid

      nt = 0
    5 read(1,'(a8)',end=99) date
      print *, 'reading in data for ', date
      nt = nt + 1

c   read in ascii gridded data

      call read_ascii(date, xflag)

      go to 5

   99 print *, 'finished reading in  data for ', nt,' times'

      end

      subroutine read_ascii(date, xflag)
c
c   routine to read the gridded ascii fields
c
c   grid (glon) goes from 80E to 130E, 1 degree resolution
c   grid (glat) goes from 10S to 40N, 1 degree resolution
c   grid (gpr)  goes from 1025 to 25mb, 25mb resolution
c   note: the bottom pressure level (gpr = 1025 ) is actually 
c         at the surface        
c
c   p - surface pressure (mb)
c   z - height (m)
c   t - temperature (C)
c   q - specific humidity (g/kg)
c   u - zonal wind component (m/s)
c   v - meridional wind component (m/s)
c

      integer   ig,    jg,    kg
      parameter(ig=51, jg=51, kg=41)
      real xmis
      parameter (xmis = -999.0)

c   analyses coordinates

      real    glon(ig), glat(jg), gpr(kg)
      common /grid_coor/ glon, glat, gpr

c   gridded analyses fields

      integer          z
      real t, q, u, v, psfc

      common /grdflds/ z(ig,jg,kg), t(ig,jg,kg), q(ig,jg,kg),
     2                 u(ig,jg,kg), v(ig,jg,kg), psfc(ig,jg) 

      character  date*8, data_label*20
      logical    xflag

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

      open(30, file='p.'//date, status='unknown', form='formatted')
      open(31, file='u.'//date, status='unknown', form='formatted')
      open(32, file='v.'//date, status='unknown', form='formatted')
      open(33, file='q.'//date, status='unknown', form='formatted')
      open(34, file='t.'//date, status='unknown', form='formatted')
      open(35, file='z.'//date, status='unknown', form='formatted')

c   ****  read in surface pressure (mb) 

      read(30,'(a20)') data_label 
      do i=1,ig
         read(30,55) (psfc(i,j),j=1,jg)
      enddo
c     print *, data_label, psfc(1,1)

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

      do k=1,kg
         read(31,'(a20)') data_label 
         do i=1,ig
            read(31,50) (u(i,j,k),j=1,jg)
         enddo
c        print *, data_label, u(1,1,k)
      enddo

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

      do k=1,kg
         read(32,'(a20)') data_label 
         do i=1,ig
             read(32,50) (v(i,j,k),j=1,jg)
         enddo
c        print *, data_label, v(1,1,k)
      enddo

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

      do k=1,kg
         read(33,'(a20)') data_label 
         do i=1,ig
             read(33,60)  (q(i,j,k),j=1,jg)
         enddo
c        print *, data_label, q(1,1,k)
      enddo

c   ****  read in temperature (C) 

      do k=1,kg
         read(34,'(a20)') data_label 
         do i=1,ig
             read(34,50) (t(i,j,k),j=1,jg)
         enddo
c        print *, data_label, t(1,1,k)
      enddo

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

      do k=1,kg
         read(35,'(a20)') data_label 
         do i=1,ig
             read(35,40) (z(i,j,k), j=1,jg)
         enddo
         print *, data_label, z(1,1,k)
      enddo

      close(30)
      close(31)
      close(32)
      close(33)
      close(34)
      close(35)

c  optional test: set values below surface to missing value flag
c                 if xflag = .true.

      if(xflag) then
         do i=1,ig
         do j=1,jg
            do k=2,kg
               if(gpr(k) .gt. psfc(i,j)) then
                  z(i,j,k) = ifix(xmis)
                  t(i,j,k) = xmis
                  q(i,j,k) = xmis
                  u(i,j,k) = xmis
                  v(i,j,k) = xmis
               endif
            enddo
         enddo
         enddo
      endif

c     format statements

   40 format(2(17i6,/),17i6)
   50 format(2(17f6.1,/),17f6.1)
   55 format(2(17f7.1,/),17f7.1)
   60 format(2(17f6.2,/),17f6.2)

      return
      end

      subroutine init_grid

      integer   ig,    jg,    kg
      parameter(ig=51, jg=51, kg=41)

c   SCSMEX/GAME gridded analyses coordinates

      real    glon(ig), glat(jg), gpr(kg)
      common /grid_coor/ glon, glat, gpr

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

      gridsp = 1.
      slat   = -10.
      slon   =  80.

      do i=1,ig
         glon(i) = slon + float(i-1)*gridsp
      enddo

      do j=1,jg
         glat(j) = slat + float(j-1)*gridsp
      enddo

      dp = 25.
      do k=1,kg
         gpr(k) = 1025. - float(k-1)*dp
      enddo

      return
      end