Andy,
Here's some code I had around. Note that at the time I wrote this I
didn't need to read in grids that have two GEMPAK levels or times
associated with them, so the code below is perhaps too simple, but it
should give you some idea of what calls to the underlying GEMPAK library
are necessary. Some source lines have wrapped.
Cheers,
Steve
module gempak
implicit none
real, parameter :: RMISSD = -9999.
integer, parameter :: LLNNAV = 256, LLNANL = 128, LLGDHD = 128
integer, parameter :: GFNAMELEN = 256, GTLEN = 20, GPLEN = 12,
PROJLEN = 4
type GemFileT
real, dimension(LLNNAV) :: navBlock = 0
real, dimension(LLNANL) :: analBlock = 0
integer :: gemID = 10, maxGrids = 1000
logical :: writeAllowed = .false.
character(GFNAMELEN) :: filename = "gemfile.gem"
end type GemFileT
contains !
=======================================================================================
! init_gem initializes the GEMPAK interface.
subroutine init_gem()
integer :: status
call in_bdta(status)
call handle_gerr(status, "ID_BDTA")
call gd_init(status)
call handle_gerr(status, "GD_INIT")
call dg_intl(status)
call handle_gerr(status, "DG_INTL")
end subroutine init_gem
!
===============================================================================================
! open_gemfile opens a GEMPAK file for further processing. The
arguments are:
! gemFileName - The name of the GEMPAK file, including the path
! readOnly - Is the file going to be read only?
! gemFile - The GEMPAK file handle
subroutine open_gemfile(gemFileName, readOnly, gemFile, stat)
character(*), intent(in) :: gemFileName
logical, intent(in) :: readOnly
type(GemFileT), intent(out) :: gemFile
integer, optional, intent(out) :: stat
integer :: status
if (len(gemFileName) > len(gemFile%filename)) &
print *, "WARNING: Filename truncated in open_gemfile!"
gemFile%filename = gemFileName
gemFile%writeAllowed = .not. readOnly
call gd_open(gemFile%filename, gemFile%writeAllowed, LLNANL,
LLNNAV, gemFile%gemID, &
gemFile%analBlock, gemFile%navBlock, gemFile%maxGrids, status)
call handle_gerr(status, "GD_OPEN")
if (present(stat)) stat = status
end subroutine open_gemfile
!
===============================================================================================
! create_gemfile creates a new GEMPAK file for further processing.
The arguments are:
! gemFileName - The name of the GEMPAK file, including the path
! gemFileToCopy - The GEMPAK file handle of the file being copied
! maxGrids - The maximum number of grids allowed in the new file
! gemFile - The GEMPAK file handle
subroutine create_gemfile(gemFileName, gemFileToCopy, maxGrids, gemFile)
character(*), intent(in) :: gemFileName
type(GemFileT), intent(in) :: gemFileToCopy
integer, intent(in) :: maxGrids
type(GemFileT), intent(out) :: gemFile
integer :: ihdrsz = 2, status
! Copy GEMPAK file metadata
gemFile%navBlock = gemFileToCopy%navBlock
gemFile%analBlock = gemFileToCopy%analBlock
gemFile%maxGrids = maxGrids
gemFile%writeAllowed = .true.
gemFile%filename = gemFileName
! Create the file
call gd_cref(gemFile%filename, LLNNAV, gemFile%navBlock, LLNANL,
gemFile%analBlock, ihdrsz, &
gemFile%maxGrids, gemFile%gemID, status)
call handle_gerr(status, "GD_CREF")
end subroutine create_gemfile
!
===============================================================================================
! read_gemfile reads the specified grid from the GEMPAK file. The
arguments are:
! gemFile - The GEMPAK file handle
! gemTime - The GEMPAK time
! level - The GEMPAK level
! coord - The GEMPAK vertical coordinate ID
! parm - The GEMPAK grid name
! grid - The actual gridded data
! stat - The GEMPAK status indicator
subroutine read_gemfile(gemFile, gemTime, level, coord, parm, grid, stat)
type(GemFileT), intent(in) :: gemFile
character(*), intent(in) :: gemTime, parm
integer, intent(in) :: level, coord
real, dimension(:,:), intent(out) :: grid
integer, optional, intent(out) :: stat
integer, dimension(LLGDHD) :: ighdr
integer, dimension(2) :: levels
integer :: igx, igy, status
character(GTLEN), dimension(2) :: gdattm
gdattm = (/ gemTime, repeat(" ", GTLEN) /)
levels = (/ level, -1 /)
igx = size(grid,1)
igy = size(grid,2)
call gd_rdat(gemFile%gemID, gdattm, levels, coord, parm, grid, igx,
igy, ighdr, status)
call handle_gerr(status, "GD_RDAT")
if (present(stat)) stat = status
end subroutine read_gemfile
!
===============================================================================================
! write_gemfile writes the specified grid to the GEMPAK file. The
arguments are:
! gemFile - The GEMPAK file handle
! gemTime - The GEMPAK time
! level - The GEMPAK level
! coord - The GEMPAK vertical coordinate ID
! parm - The GEMPAK grid name
! grid - The actual gridded data
! stat - The GEMPAK status indicator
subroutine write_gemfile(gemFile, gemTime, level, coord, parm, grid,
overwrite, stat)
type(GemFileT), intent(in) :: gemFile
character(*), intent(in) :: gemTime, parm
integer, intent(in) :: level, coord
real, dimension(:,:), intent(in) :: grid
logical, intent(in) :: overwrite
integer, optional, intent(out) :: stat
integer, dimension(2) :: ighdr = 0, levels
integer :: igx, igy, status
character(GTLEN), dimension(2) :: gdattm
gdattm = (/ gemTime, repeat(" ", GTLEN) /)
levels = (/ level, -1 /)
igx = size(grid,1)
igy = size(grid,2)
call gd_wdat(gemFile%gemID, grid, igx, igy, ighdr, gdattm, levels,
coord, parm, overwrite, &
status)
call handle_gerr(status, "GD_RDAT")
if (present(stat)) stat = status
end subroutine write_gemfile
!
===============================================================================================
! close_gemfile closes the specified GEMPAK file. The arguments are:
! gemFile - The GEMPAK file handle
subroutine close_gemfile(gemFile)
type(GemFileT), intent(in) :: gemFile
integer :: status
call gd_clos(gemFile%gemID, status)
call handle_gerr(status, "GD_CLOS")
end subroutine close_gemfile
!
===============================================================================================
! handle_gerr examines the result code from a gemlib routine. If
necessary, it then prints an
! error message and aborts the program.
subroutine handle_gerr(status, routine)
integer, intent(in) :: status
character(len=*), intent(in) :: routine
if (status /= 0) write (*,"(a,a,i3)") routine, ": status = ", status
end subroutine handle_gerr
end module gempak
On 11/24/2014 09:42 AM, Andrew Penny - NOAA Affiliate wrote:
Hi all,
I'm new to GEMPAK and am wondering if someone has a simple fortran
routine to read in data from a GEMPAK data file?
Thanks,
Andy
_______________________________________________
gembud mailing list
gembud@xxxxxxxxxxxxxxxx
For list information or to unsubscribe, visit:
http://www.unidata.ucar.edu/mailing_lists/
--
Steve Decker, Teaching Instructor
Department of Environmental Sciences Phone: (848) 932-5750
Rutgers University Fax: (732) 932-8644
14 College Farm Rd Email: decker@xxxxxxxxxxxxxxxxxx
New Brunswick, NJ 08901-8551
