Re: [gembud] read GEMPAK data with fortran

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