!> @file
!! @brief Read in data defining the model grid.
!! @author George Gayno org: w/np2 @date 2005-Dec-16
!> Read in data defining the model grid.
!!
!! program history log:
!! - 2005-dec-16 gayno - initial version
!! - 2007-nov-30 gayno - improved method for thinning gfs grids.
!! added nam b-grids.
!! - 2014-sep-29 gayno - add option to read model lat, lon and
!! landmask data in grib2.
!!
!! variable definitions:
!! - lonsperlat - for global grids, the number of i points
!! in each row (decrease toward pole)
!!
!! @author George Gayno org: w/np2 @date 2005-Dec-16
module model_grid
use program_setup, only : model_lsmask_file, &
model_lon_file, &
model_lat_file, &
gfs_lpl_file
integer :: grid_id_mdl !< grib id of model grid, 4-gaussian, 203-egrid
integer :: imdl !< i-dimension of model grid
integer :: jmdl !< j-dimension of model grid
integer :: ijmdl !< total number of model land points
integer, allocatable :: ipts_mdl(:) !< i index of point on full grid
integer, allocatable :: jpts_mdl(:) !< j index of point on full grid
integer :: kgds_mdl(200) !< holds grib gds info of model grid
integer, allocatable :: lonsperlat_mdl (:) !< Number of longitudes (i-points)
!! for each latitude (row). Used
!! for global thinned (reduced) grids.
logical :: thinned !< When true, global grids will run thinned
!! (number of i points decrease toward pole)
real, allocatable :: lats_mdl(:) !< Latitudes of model grid points
real :: lat11 !< Corner point latitude (1,1) of model grid.
real :: latlast !< Corner point latitude (imdl,jmdl) of model grid.
real :: lon11 !< Corner point longitude (1,1) of model grid.
real :: lonlast !< Corner point longitude (imdl,jmdl) of model grid.
real, allocatable :: lons_mdl(:) !< longitudes of model grid points
real, allocatable :: lsmask_mdl(:,:) !< land mask of model grid (0 - non land, 1-land) for global
!! grids run thinned, will contain a modified version of the original
!! mask that has land at all points encompassed by a thinned point
real, allocatable :: lsmask_mdl_sav(:,:) !< saved copy of land mask of model grid (0 - non land, 1-land)
!! only used for global thinned grids.
real :: resol_mdl !< approximate model resolution in km.
contains
!> Read mdl grid.
!!
!! program history log:
!! 2005-dec-16 gayno - initial version
!! 2007-nov-30 gayno - Improved method for thinning gfs grids.
!! Added nam b-grids.
!! 2014-sep-29 gayno - Add option to read lat,lon and mask
!! data in grib2.
!! files:
!! inputs:
!! - model latitudes (grib 1 or grib 2)
!! - model longitudes (grib 1 or grib 2)
!! - model landmask (grib 1 or grib 2)
!! - number pts per row, gfs grid (the "lonsperlat" file, ascii)
!! condition codes: all fatal
!! 76 - bad open/read gfs "lonsperlat" file
!! 79 - unrecognized model grid
!! 80 - bad open model latitude file
!! 81 - bad read of model latitude grib 1 header
!! 82 - bad read of model latitude data
!! 83 - bad open model longitude file
!! 82 - bad read of model longitude data
!! 85 - bad open model landmask file
!! 86 - bad read of model landmask data
!! 90 - model latitude file not grib 1 or grib 2
!! 91 - model longitude file not grib 1 or grib 2
!! 92 - model landmask file not grib 1 or grib 2
!!
!! @author George Gayno org: w/np2 @date 2005-dec-16
subroutine read_mdl_grid_info
use grib_mod ! grib 2 library
implicit none
character*256 :: fngrib
integer :: i, j, ij, jj
integer :: ii, iii, istart, iend, imid
integer :: iret
integer :: isgrib
integer, parameter :: iunit = 14 ! unit of input grib file
integer, parameter :: iunit2 = 34 ! unit of input lonsperlat file
integer :: jgds(200)
integer :: jpds(200)
integer :: jdisc, jgdtn, jpdtn, k
integer :: jids(200), jgdt(200), jpdt(200)
integer :: lskip
integer, parameter :: lugi = 0 ! unit of grib index file - not used
integer :: kgds(200)
integer :: kpds(200)
integer :: message_num
integer :: numbytes
integer :: numpts
logical*1, allocatable :: lbms(:)
logical :: unpack
real :: gridis, gridie, fraction, x1, r
real, allocatable :: lats_mdl_temp (:,:)
real, allocatable :: lons_mdl_temp (:,:)
type(gribfield) :: gfld
print*,"- READ MODEL GRID INFORMATION"
!-----------------------------------------------------------------------
! read latitudes on the model grid. first check if file is grib1 or 2.
!-----------------------------------------------------------------------
fngrib = model_lat_file
call grib_check(fngrib, isgrib)
if (isgrib==0) then
print*,'- FATAL ERROR: MODEL LAT FILE MUST BE GRIB1 OR GRIB2 FORMAT'
call w3tage('SNOW2MDL')
call errexit(90)
end if
print*,"- OPEN MODEL LAT FILE ", trim(fngrib)
call baopenr (iunit, fngrib, iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD OPEN, IRET IS ', iret
call w3tage('SNOW2MDL')
call errexit(80)
end if
if (isgrib==1) then ! grib 1 file
!-----------------------------------------------------------------------
! tell degribber to search for latitudes
!-----------------------------------------------------------------------
lskip = -1 ! read beginning of file
jgds = -1
jpds = -1
jpds(5) = 176 ! latitude
kgds = -1
kpds = -1
print*,"- GET GRIB HEADER"
call getgbh(iunit, lugi, lskip, jpds, jgds, numbytes, &
numpts, message_num, kpds, kgds, iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD READ OF GRIB HEADER. IRET IS ', iret
call w3tage('SNOW2MDL')
call errexit(81)
end if
!-----------------------------------------------------------------------
! save gds for gribbing the interpolated data later. also required
! by ncep ipolates library.
!-----------------------------------------------------------------------
kgds_mdl = kgds
!-----------------------------------------------------------------------
! get model grid specs from header. model resolution (km) is used
! to determine the interpolation method.
!-----------------------------------------------------------------------
grid_id_mdl = kpds(3) ! grib 1 grid id number. sect 1, oct 7
if (kgds(1) == 4) then ! gaussian grid
imdl = kgds(2) ! i-dimension of model grid
jmdl = kgds(3) ! j-dimension of model grid
resol_mdl = float(kgds(9)) / 1000.0 * 111.0
else if (kgds(1) == 203) then ! e-grid
imdl = kgds(2) ! i-dimension of model grid
jmdl = kgds(3) ! j-dimension of model grid
resol_mdl = sqrt( (float(kgds(9)) / 1000.0)**2 + &
(float(kgds(10)) / 1000.0)**2 )
resol_mdl = resol_mdl * 111.0
else if (kgds(1) == 205) then ! b-grid
imdl = kgds(2) ! i-dimension of model grid
jmdl = kgds(3) ! j-dimension of model grid
resol_mdl = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
* 0.5 * 111.0
else
print*,'- FATAL ERROR: UNRECOGNIZED MODEL GRID.'
call w3tage('SNOW2MDL')
call errexit(79)
end if
allocate(lats_mdl_temp(imdl,jmdl))
allocate(lbms(imdl*jmdl))
print*,"- DEGRIB DATA"
call getgb(iunit, lugi, (imdl*jmdl), lskip, jpds, jgds, &
numpts, message_num, kpds, kgds, lbms, lats_mdl_temp, iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD DEGRIB OF FILE. IRET IS ',iret
call w3tage('SNOW2MDL')
call errexit(82)
end if
deallocate(lbms)
lat11 = lats_mdl_temp(1,1)
latlast = lats_mdl_temp(imdl,jmdl)
elseif (isgrib==2) then ! grib 2 file
j = 0 ! search at beginning of file
jdisc = 0 ! search for discipline; 0 - meteorological products
jpdtn = -1 ! search for any product definition template number
jgdtn = -1 ! search for any grid definition template number
jids = -9999 ! array of values in identification section, set to wildcard
jgdt = -9999 ! array of values in grid definition template 3.m
jpdt = -9999 ! array of values in product definition template 4.n
jpdt(1) = 191 ! search for parameter category - misc
jpdt(2) = 1 ! search for parameter number - latitude
unpack = .true. ! unpack data
call grib2_null(gfld)
print*,"- DEGRIB DATA"
call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
unpack, k, gfld, iret)
if (iret /=0) then
print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret
call w3tage('SNOW2MDL')
call errexit(82)
endif
!-----------------------------------------------------------------------
! create the grib 1 gds array from the g2 gdt array. the grib 1
! gds info is used by ipolates and for gribbing the final snow analysis.
!-----------------------------------------------------------------------
call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds_mdl, &
imdl, jmdl, resol_mdl)
grid_id_mdl = 255 ! grib1 grid id number. n/a for grib2.
! set to 'missing'.
allocate(lats_mdl_temp(imdl,jmdl))
lats_mdl_temp = reshape (gfld%fld , (/imdl,jmdl/) )
lat11 = lats_mdl_temp(1,1)
latlast = lats_mdl_temp(imdl,jmdl)
call grib2_free(gfld)
endif ! grib1 or grib2?
call baclose(iunit,iret)
!-----------------------------------------------------------------------
! read longitudes on the model grid.
!-----------------------------------------------------------------------
fngrib = model_lon_file
call grib_check(fngrib, isgrib)
if (isgrib==0) then
print*,'- FATAL ERROR: MODEL LON FILE MUST BE GRIB1 OR GRIB2 FORMAT'
call w3tage('SNOW2MDL')
call errexit(91)
end if
print*,"- OPEN MODEL LON FILE ", trim(fngrib)
call baopenr (iunit, fngrib, iret)
if (iret /= 0) then
print*,"- FATAL ERROR: BAD OPEN. IRET IS ", iret
call w3tage('SNOW2MDL')
call errexit(83)
end if
if (isgrib==1) then ! grib 1 file
lskip = -1
kgds = -1
kpds = -1
jgds = -1
jpds = -1
jpds(5) = 177 ! longitude
allocate(lons_mdl_temp(imdl,jmdl))
allocate(lbms(imdl*jmdl))
print*,"- DEGRIB DATA"
call getgb(iunit, lugi, (imdl*jmdl), lskip, jpds, jgds, &
numpts, message_num, kpds, kgds, lbms, lons_mdl_temp, iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD DEGRIB OF DATA. IRET IS ',iret
call w3tage('SNOW2MDL')
call errexit(84)
end if
deallocate(lbms)
lon11 = lons_mdl_temp(1,1)
lonlast = lons_mdl_temp(imdl,jmdl)
elseif (isgrib==2) then ! grib2
j = 0 ! search at beginning of file
jdisc = 0 ! search for discipline; 0 - meteorological products
jpdtn = -1 ! search for any product definition template number
jgdtn = -1 ! search for any grid definition template number
jids = -9999 ! array of values in identification section, set to wildcard
jgdt = -9999 ! array of values in grid definition template 3.m
jpdt = -9999 ! array of values in product definition template 4.n
jpdt(1) = 191 ! search for parameter category - misc
jpdt(2) = 2 ! search for parameter number - longitude
unpack = .true. ! unpack data
call grib2_null(gfld)
print*,"- DEGRIB DATA"
call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
unpack, k, gfld, iret)
if (iret /=0) then
print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret
call w3tage('SNOW2MDL')
call errexit(84)
endif
allocate(lons_mdl_temp(imdl,jmdl))
lons_mdl_temp = reshape (gfld%fld , (/imdl,jmdl/) )
lon11 = lons_mdl_temp(1,1)
lonlast = lons_mdl_temp(imdl,jmdl)
call grib2_free(gfld)
endif ! grib1 or grib?
call baclose(iunit, iret)
!-----------------------------------------------------------------------
! read model land/sea mask.
!-----------------------------------------------------------------------
fngrib = model_lsmask_file
call grib_check(fngrib, isgrib)
if (isgrib==0) then
print*,'- FATAL ERROR: MODEL LANDMASK FILE MUST BE GRIB1 OR GRIB2 FORMAT'
call w3tage('SNOW2MDL')
call errexit(92)
end if
print*,"- OPEN MODEL LANDMASK FILE ", trim(fngrib)
call baopenr (iunit, fngrib, iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD OPEN OF FILE. IRET IS ', iret
call w3tage('SNOW2MDL')
call errexit(85)
end if
if (isgrib==1) then ! grib 1 file
lskip = -1
kgds = -1
kpds = -1
jpds = -1
jgds = -1
jpds(5) = 81 ! land-sea mask
allocate(lsmask_mdl(imdl,jmdl))
allocate(lbms(imdl*jmdl))
print*,"- DEGRIB DATA"
call getgb(iunit, lugi, (imdl*jmdl), lskip, jpds, jgds, &
numpts, message_num, kpds, kgds, lbms, lsmask_mdl, iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD DEGRIB OF DATA. IRET IS ',iret
call w3tage('SNOW2MDL')
call errexit(86)
end if
deallocate (lbms)
elseif (isgrib==2) then ! grib2
j = 0 ! search at beginning of file
jdisc = 2 ! search for discipline; 2 - land-sfc products
jpdtn = -1 ! search for any product definition template number
jgdtn = -1 ! search for any grid definition template number
jids = -9999 ! array of values in identification section, set to wildcard
jgdt = -9999 ! array of values in grid definition template 3.m
jpdt = -9999 ! array of values in product definition template 4.n
jpdt(1) = 0 ! search for parameter category - veg_biomass
jpdt(2) = 0 ! search for parameter number - landcover
unpack = .true. ! unpack data
call grib2_null(gfld)
print*,"- DEGRIB DATA"
call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
unpack, k, gfld, iret)
if (iret /=0) then
print*,'- FATAL ERROR: BAD DEGRIB OF FILE, IRET IS ', iret
call w3tage('SNOW2MDL')
call errexit(86)
endif
allocate(lsmask_mdl(imdl,jmdl))
lsmask_mdl = reshape (gfld%fld , (/imdl,jmdl/) )
call grib2_free(gfld)
endif
call baclose(iunit,iret)
!-----------------------------------------------------------------------
! global model runs on a thinned grid (# grid points decreases
! towards the poles). if thinned logical is set, and this is a
! gaussian grid, modify the land/sea mask to account for the
! fact that delta x increases toward the poles.
!-----------------------------------------------------------------------
thinned=.false.
if (kgds(1) == 4 .and. (len_trim(gfs_lpl_file) > 0)) then
thinned=.true.
print*,"- RUNNING A THINNED GRID"
allocate (lonsperlat_mdl(jmdl/2))
print*,"- OPEN/READ GFS LONSPERLAT FILE: ",trim(gfs_lpl_file)
open (iunit2, file=trim(gfs_lpl_file), iostat=iret)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD OPEN OF LONSPERLAT FILE. ABORT. IRET: ', iret
call w3tage('SNOW2MDL')
call errexit(76)
endif
read (iunit2,*,iostat=iret) numpts, lonsperlat_mdl
close(iunit2)
if (iret /= 0) then
print*,'- FATAL ERROR: BAD READ OF LONSPERLAT FILE. ABORT. IRET: ', iret
call w3tage('SNOW2MDL')
call errexit(76)
endif
if (numpts /= (jmdl/2)) then
print*,'- FATAL ERROR: WRONG DIMENSIION IN LONSPERLAT FILE. ABORT.'
call w3tage('SNOW2MDL')
call errexit(76)
endif
allocate (lsmask_mdl_sav(imdl,jmdl))
lsmask_mdl_sav = lsmask_mdl
lsmask_mdl = 0.0 ! this will identify land points to be processed by
! the ipolates routines.
!-----------------------------------------------------------------------
! loop over every point on the thinned grid. calculate the start/end
! bounds with respect to the full grid in the 'i' direction. if
! the thinned point contains land, set all full grid points within
! the bounds to be land. this modified mask will identify the
! points to be processed by ipolates. after the call to ipolates,
! the thinned points will be set to a linear weighting of the full points
! located within the thinned point.
!-----------------------------------------------------------------------
do j = 1, jmdl
jj = j
if (j > jmdl/2) jj = jmdl - j + 1
r = float(imdl)/ float(lonsperlat_mdl(jj))
do i = 1, lonsperlat_mdl(jj)
x1=float(i-1)*r
imid = nint(x1+1.0) ! for this thinned grid point, this is
! the nearest 'i' index on the full grid.
if (lsmask_mdl_sav(imid,j) > 0.0) then
gridis = x1+1.0-r/2.
istart = nint(gridis)
gridie = x1+1.0+r/2.
iend = nint(gridie)
do ii = istart, iend
if (ii == istart) then
fraction = 0.5 - (gridis - float(istart))
if (fraction < 0.0001) cycle
endif
if (ii == iend) then
fraction = 0.5 + (gridie - float(iend))
if (fraction < 0.0001) cycle
endif
iii = ii
if (iii < 1) iii = imdl + iii
lsmask_mdl(iii,j) = lsmask_mdl_sav(imid,j)
enddo
endif
enddo
enddo
end if
!-----------------------------------------------------------------------
! program only worries about land points. save i/j coordinate
! with respect to 2-d grid.
!-----------------------------------------------------------------------
ij = 0
do j = 1, jmdl
do i = 1, imdl
if (lsmask_mdl(i,j) > 0.0) then
ij = ij+1
end if
enddo
enddo
ijmdl = ij
if (ijmdl == 0) then ! grid has only water points, dont run
print*,' '
print*,'- MODEL GRID ONLY HAS WATER POINTS, DONT CREATE SNOW FILE.'
print*,'- NORMAL TERMINATION.'
call w3tage('SNOW2MDL')
call errexit(0)
endif
allocate (lats_mdl(ijmdl))
allocate (lons_mdl(ijmdl))
allocate (ipts_mdl(ijmdl))
allocate (jpts_mdl(ijmdl))
ij = 0
do j = 1, jmdl
do i = 1, imdl
if (lsmask_mdl(i,j) > 0.0) then
ij = ij+1
lats_mdl(ij) = lats_mdl_temp(i,j)
lons_mdl(ij) = lons_mdl_temp(i,j)
ipts_mdl(ij) = i
jpts_mdl(ij) = j
end if
enddo
enddo
deallocate (lats_mdl_temp, lons_mdl_temp)
return
end subroutine read_mdl_grid_info
!> Clean up allocatable arrays.
!!
!! This deallocate this module's allocatable array.
!!
!! program history log:
!! 2005-dec-16 gayno - initial version
!!
!! @author George Gayno org: w/np2 @date Dec 16, 2005
subroutine model_grid_cleanup
implicit none
if (allocated(lsmask_mdl)) deallocate(lsmask_mdl)
if (allocated(lats_mdl)) deallocate(lats_mdl)
if (allocated(lons_mdl)) deallocate(lons_mdl)
if (allocated(lonsperlat_mdl)) deallocate(lonsperlat_mdl)
if (allocated(ipts_mdl)) deallocate(ipts_mdl)
if (allocated(jpts_mdl)) deallocate(jpts_mdl)
return
end subroutine model_grid_cleanup
end module model_grid