NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
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 = 128integer, 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 indicatorsubroutine 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
gembud archives: