program readncep

!-- Read some sample 4*daily NCEP isentropic level reanalysis data from an ASCII file.
!-- There are five fields contained in the file:
!--   1) The pressure at 11 isentropic levels [Pa]
!--   2) The Montgomery streamfunction at 11 isentropic levels [m^2/s^2]
!--   3) The meridional wind at 11 isentropic levels [m/s]
!--   4) The surface pressure [Pa] 
!--   5) The surface geopotential (i.e., topography) [m]
!-- The upper-air data exists on 11 constant-theta levels (K):
!--   270, 280, 290, 300, 315, 330, 350, 400, 450, 550, 650
!-- Missing data (where the theta surfaces run underground) is denoted 
!--   by the value 1.0e+36.

   implicit none
   integer, parameter :: nlon=144, nlat=73, nlev=11
   real, dimension(nlon,nlat,nlev) :: pres, mntsf, vgrd
   real, dimension(nlon,nlat) :: data, psfc, zsfc
   character (len=80) :: infile='ncep.data', hdr
   integer :: ilev, status
   logical :: exans
   
   inquire(file=infile,exist=exans)
   if (.not. exans) then
     write(*,*) 'File not found: ',infile
     stop
   else
     open(unit=10,file=infile,form='formatted',action='read',status='old')
!-- read the pressure data
     do ilev = 1,nlev
       read(10,'(a80)',iostat=status) hdr 
       read(10,'(5e15.8)',iostat=status) data
       if (status < 0) then
         write(*,*) "read error on",infile
         stop
       else
         pres(:,:,ilev) = data(:,:)
       endif
       write(*,*) ilev,minval(pres(:,:,ilev)),maxval(pres(:,:,ilev))
     enddo 
!-- read the streamfunction data
     do ilev = 1,nlev
       read(10,'(a80)',iostat=status) hdr 
       read(10,'(5e15.8)',iostat=status) data
       if (status < 0) then
         write(*,*) "read error on",infile
         stop
       else
         mntsf(:,:,ilev) = data(:,:)
       endif
       write(*,*) ilev,minval(mntsf(:,:,ilev)),maxval(mntsf(:,:,ilev))
     enddo 
!-- read the v-wind data
     do ilev = 1,nlev
       read(10,'(a80)',iostat=status) hdr 
       read(10,'(5e15.8)',iostat=status) data
       if (status < 0) then
         write(*,*) "read error on",infile
         stop
       else
         vgrd(:,:,ilev) = data(:,:)
       endif
       write(*,*) ilev,minval(vgrd(:,:,ilev)),maxval(vgrd(:,:,ilev))
     enddo 
     read(10,'(a80)',iostat=status) hdr 
     read(10,'(5e15.8)',iostat=status) psfc
     write(*,*) minval(psfc),maxval(psfc)
     read(10,'(a80)',iostat=status) hdr 
     read(10,'(5e15.8)',iostat=status) zsfc
     write(*,*) minval(zsfc),maxval(zsfc)
     close(10)
   endif

end program readncep