! --read T1534 NEMS GFS forecast files
! --to compile: ifort -mcmodel medium -shared-intel -convert big_endian -FR read_atmfh.f
! --fanglin.yang@noaa.gov   NWS/NOAA/NCEP/EMC

integer, parameter :: nlat=1536, nlon=3072, nlev=64
integer, parameter :: nvar3d=6,nvar2d=2,nrec=nvar3d*nlev+nvar2d

character(255) filename

!--forecast files, north-to-south, surface-to-top
real clwmr(nlon,nlat,nlev)            
real hgtsfc(nlon,nlat)         
real o3mr(nlon,nlat,nlev)            
real pressfc(nlon,nlat)         
real spfh(nlon,nlat,nlev)            
real tmp(nlon,nlat,nlev)            
real ugrd(nlon,nlat,nlev)       
real vgrd(nlon,nlat,nlev)            

call getarg(1,filename)
open (1, file=trim(filename),form="unformatted",access="sequential",status="unknown")

!--determine total lines then the lines to skip
nn=1
100 continue
read(1,end=200)
nn=nn+1
goto 100
200 print*, 'total lines =',nn
nskip=nn-nrec-1

!---rewind then skip header lines
rewind 1
do k=1,nskip
 read(1)
enddo

 do k=1,nlev
  read(1) ((clwmr(i,j,k),i=1,nlon),j=1,nlat)           
  print *,'clwmr=',clwmr(1,1,k),maxval(clwmr(:,:,k)),minval(clwmr(:,:,k))
 enddo

 read(1) ((hgtsfc(i,j),i=1,nlon),j=1,nlat)         
 print *,'hgt=',hgtsfc(1,1),maxval(hgtsfc(:,:)),minval(hgtsfc(:,:))

 do k=1,nlev
  read(1) ((o3mr(i,j,k),i=1,nlon),j=1,nlat)           
  print *,'o3mr=',o3mr(1,1,k),maxval(o3mr(:,:,k)),minval(o3mr(:,:,k))
 enddo

 read(1) ((pressfc(i,j),i=1,nlon),j=1,nlat)        
 print *,'psfc=',pressfc(1,1),maxval(pressfc(:,:)),minval(pressfc(:,:))

 do k=1,nlev
  read(1) ((spfh(i,j,k),i=1,nlon),j=1,nlat)           
  print *,'spfh=',spfh(1,1,k),maxval(spfh(:,:,k)),minval(spfh(:,:,k))
 enddo

 do k=1,nlev
  read(1) ((tmp(i,j,k),i=1,nlon),j=1,nlat)           
  print *,'tmp=',tmp(1,1,k),maxval(tmp(:,:,k)),minval(tmp(:,:,k))
 enddo

 do k=1,nlev
  read(1) ((ugrd(i,j,k),i=1,nlon),j=1,nlat)      
  print *,'ugrd=',ugrd(1,1,k),maxval(ugrd(:,:,k)),minval(ugrd(:,:,k))
 enddo

 do k=1,nlev
  read(1) ((vgrd(i,j,k),i=1,nlon),j=1,nlat)           
  print *,'vgrd=',ugrd(1,1,k),maxval(vgrd(:,:,k)),minval(vgrd(:,:,k))
 enddo



end