model_grid.F90

Go to the documentation of this file.

 module model_grid

 use program_setup, only         : model_lsmask_file, &
                                   model_lon_file, &
                                   model_lat_file, &
                                   gfs_lpl_file

 integer                        :: grid_id_mdl
 integer                        :: imdl
 integer                        :: jmdl
 integer                        :: ijmdl
 integer, allocatable           :: ipts_mdl(:)
 integer, allocatable           :: jpts_mdl(:)

 integer                        :: kgds_mdl(200)
 integer, allocatable           :: lonsperlat_mdl (:)
 
 logical                        :: thinned

 real, allocatable              :: lats_mdl(:)
 real                           :: lat11
 real                           :: latlast
 real                           :: lon11
 real                           :: lonlast
 real, allocatable              :: lons_mdl(:)
 real, allocatable              :: lsmask_mdl(:,:)
 real, allocatable              :: lsmask_mdl_sav(:,:)
 real                           :: resol_mdl

 contains
 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


 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