program readera40daily

!-- read some sample 4*daily ECMWF ERA-40 analysis data
!--   SET nrec to match the number of records in the input files
!-- the data exists on 23 pressure levels (mb):
!--   1,2,3,5,7,10,20,30,50,70,100,150,200,250,300,400,500,600,700,775,850,925,1000

!-- variables and units:
!--  z (geopotential):         m^2/s^2
!--  t (temperature):          K
!--  q (specific humidity):    kg/kg
!--  w (vertical velocity):    Pa/s
!--  u (zonal velocity):       m/s
!--  v (meridional velocity)   m/s
!--
!--  psfc (surface pressure)   Pa

   implicit none
   integer, parameter :: nlon=144, nlat=73, nlev=23, nrec=124
   integer :: ilev, irec
   real, dimension(nlon,nlat,nrec,nlev) :: z, t, q, u, v, w
   real, dimension(nlon,nlat,nrec) :: psfc
   real, dimension(nlon,nlat) :: sfcgeo
   character (len=80) :: infile
   character (len=10) :: file_prefix = 'daily_4/'
   character (len=30) :: zfile='z200001.bin', tfile='t200001.bin', &
     qfile='q200001.bin', ufile='u200001.bin', vfile='v200001.bin', &
     wfile='w200001.bin', psfcfile='psfc200001.bin'
   integer :: status
   logical :: exans
      
!-- read surface geopotential (orography, which is time invariant)
   infile = trim(file_prefix) // 'sfcgeo.bin'
   write(*,*) 'infile = ',trim(infile)
   inquire(file=infile,exist=exans)
   if (.not. exans) then
     write(*,*) 'File not found: ',infile
     stop
   else
     open(unit=10,file=infile,form='unformatted',action='read')
     read(10,iostat=status) sfcgeo
     if (status < 0) then
       write(*,*) "read error on",zfile
       stop
     endif
     write(*,*) 'min/max(sfcgeo) = ',minval(sfcgeo),maxval(sfcgeo)
     close(10)
   endif

!-- read geopotential file
   infile = trim(file_prefix) // trim(zfile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_upperair(infile, z)
   write(*,*) 'min/max(z) = ',minval(z),maxval(z)

!-- read temperature file
   infile = trim(file_prefix) // trim(tfile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_upperair(infile, t)
   write(*,*) 'min/max(t) = ',minval(t),maxval(t)

!-- read specific humidity file
   infile = trim(file_prefix) // trim(qfile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_upperair(infile, q)
   write(*,*) 'min/max(q) = ',minval(q),maxval(q)

!-- read zonal wind file
   infile = trim(file_prefix) // trim(ufile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_upperair(infile, u)
   write(*,*) 'min/max(u) = ',minval(u),maxval(u)

!-- read meridional wind file
   infile = trim(file_prefix) // trim(vfile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_upperair(infile, v)
   write(*,*) 'min/max(v) = ',minval(v),maxval(v)

!-- read vertical velocity file
   infile = trim(file_prefix) // trim(wfile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_upperair(infile, w)
   write(*,*) 'min/max(w) = ',minval(w),maxval(w)

!-- read log(sfc pressure) file
   infile = trim(file_prefix) // trim(psfcfile)
   write(*,*) 'infile = ',trim(infile)
   call read_era40_sfc(infile, psfc)
   psfc = exp(psfc)
   write(*,*) 'min/max(psfc) = ',minval(psfc),maxval(psfc)

   contains
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine read_era40_upperair (infile, thedata)
  
   character (len=*) :: infile
   real, dimension(:,:,:,:) :: thedata
   real, dimension(nlon,nlat) :: data
   integer :: 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='unformatted',action='read')
     do irec = 1,nrec
       do ilev = 1,nlev
         read(10,iostat=status) data
         if (status < 0) then
           write(*,*) "read error on",zfile
           stop
         else
           thedata(:,:,irec,ilev) = data(:,:)
         endif
       enddo
     enddo 
     close(10)
   endif

   end subroutine read_era40_upperair
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine read_era40_sfc (infile, thedata)
  
   character (len=*) :: infile
   real, dimension(:,:,:) :: thedata
   real, dimension(nlon,nlat) :: data
   
   inquire(file=infile,exist=exans)
   if (.not. exans) then
     write(*,*) 'File not found: ',infile
     stop
   else
     open(unit=10,file=infile,form='unformatted',action='read')
     do irec = 1,nrec
       read(10,iostat=status) data
       if (status < 0) then
         write(*,*) "read error on",infile
         stop
       else
         thedata(:,:,irec) = data(:,:)
       endif
     enddo 
     close(10)
   endif

   end subroutine read_era40_sfc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
end program readera40daily