subroutine get_seaice(xlats_ij,xlons_ij,ny,nx,sice)
!
! Abstract: get sea ice analysis at atime and target resolution (nx,ny)
!
!
! input & output
!
 real,    dimension(nx*ny), intent(in)  :: xlats_ij   ! latitudes of target grids (nx*ny)
 real,    dimension(nx*ny), intent(in)  :: xlons_ij   ! latitudes of target grids (nx*ny)
 real,    dimension(nx,ny), intent(out) :: sice       ! sea ice atime (nx,ny)
 integer, intent(in) :: nx,ny
! local declare
 integer, dimension(nx*ny)              :: mask_ij    ! mask at target grids (0 = water, 1 = land, 2 = sea ice) (nx*ny)
 real,    dimension(nx*ny)              :: sice_ij    ! sea ice  analysis at target grids (nx*ny)

 real,    allocatable, dimension(:,:)   :: sice0      ! sea ice analysis
 integer, allocatable, dimension(:,:)   :: smask      ! surface mask of sea ice
 real,    allocatable, dimension(:)     :: sxlats     ! latitudes of sea ice analysis
 real,    allocatable, dimension(:)     :: sxlons     ! longitudes of sea ice analysis

 integer :: nxs,nys,mon1,mon2,sfcflag,i,j
 character (len=6), parameter :: fin_iceanl='iceanl'  ! sea ice analysis file name
!
! get the dimensions of the sea ice analysis & allocate the related arrays
!
 call get_dim_grb(fin_iceanl,nys,nxs)
 allocate( sice0(nxs,nys),smask(nxs,nys),sxlats(nys),sxlons(nxs) )
!
! read sea ice analysis
!
  call read_seaice_grb(fin_iceanl,nys,nxs,sxlats,sxlons,sice0) 
!
! get sice (nx by ny lat/lon)
!
 if ( nx == nxs .and. ny == nys ) then
    sice(:,:) = sice0(:,:)
    write(*,'(a,4I8,2F9.3)') 'same dimensions,nx,ny,nxs,nys : ',nx,ny,nxs,nys,minval(sice),maxval(sice)
 else
    sfcflag=0
    smask=0
    mask_ij=0
    write(*,'(a,4I8)') 'different dimensions,nx,ny,nxs,nys : ',nx,ny,nxs,nys
    call lalo_to_tile(sice0,   smask,   sxlats,  sxlons,  nys, nxs, &
                      sice_ij, mask_ij, xlats_ij,xlons_ij,ny,  nx, &
                      sfcflag,0.0,0,0.0)
    write(*,'(a,2F9.3)') 'done with lalo_to_tile for sice',minval(sice_ij),maxval(sice_ij)
    sice(:,:) = reshape (sice_ij, (/nx,ny/) )
 endif
end subroutine get_seaice

subroutine get_dim_grb(file_ice,mlat,mlon)
!                .      .    .                                       .
! abstract:   get_grb_dim :  get dimensions
!   prgmmr: xu li            org: np23                date: 2019-03-13
!
!
! program history log:
!
!   input argument list:
!     file - file name of GRIB SST file
! output
!     mlat,mlon
!     call subs: getgbh
!
!
!$$$
  implicit none

! Declare passed variables and arrays
  character (len=*),  intent(in)  :: file_ice
  integer,            intent(out) :: mlat,mlon
! Declare local parameters
  integer,parameter:: lu = 22   ! FORTRAN unit number of GRIB SST file

  integer :: iret,i,j,kf,kg,k
  integer, dimension(22):: jgds,kgds
  integer, dimension(25):: jpds,kpds

!************+******************************************************************************
!
! Open analysis file (GRIB)
  write(*,*) 'get_dim_grb, file_ice : ',file_ice
  call baopenr(lu,trim(file_ice),iret)
  if (iret /= 0 ) then
     write(6,*)'get_dim_grb:  ***ERROR*** opening sea ice file',file_ice
     stop
  endif

! Define variables for read
  j=-1
  jpds=-1
  jgds=-1
  jpds(5)=91          ! sea ice variable
  jpds(6)=102         ! surface
  call getgbh(lu,0,j,jpds,jgds,kg,kf,k,kpds,kgds,iret)

  mlat = kgds(3)   ! number points on longitude circle (360)
  mlon = kgds(2)   ! number points on latitude circle (720)

  write(*,*) 'sea ice nys, nxs : ',mlat, mlon

  call baclose(lu,iret)
  if (iret /= 0 ) then
     write(6,*)'get_dim_grb:  ***ERROR*** close sea ice file'
     stop
  endif

end subroutine get_dim_grb

subroutine read_seaice_grb(file_sice,mlat,mlon,xlats,xlons,sice)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rdgrbsice                   read grib1 SST analysis
!   prgmmr: xu li            org: np23                date: 2003-04-04
!
! abstract: read sea ice analysis (GRIB format) and save it as expanded and transposed array
!
!     Subroutine rdgrbsice must be compiled with the NCEP W3 library
!     and the BACIO library.
!
!
! program history log:
!   2003-04-04  xu li, bert katz
!   2005-04-18  treadon - fill southern and northern rows of sice grid
!                         with mean of adjacent row.  This treatment is
!                         consistent with rdgesfc.f90 and rdgesig.f90
!
!   input argument list:
!     file_sice - file name of GRIB sea ice file
!     mlat      - dimension in y-direction
!     mlon      - dimension in x-direction
!
!     call subs: getgbh, getgb
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  implicit none

! Declare passed variables and arrays
  character(6),               intent(in)  :: file_sice
  integer,                    intent(in)  :: mlat,mlon
  real, dimension(mlat),      intent(out) :: xlats
  real, dimension(mlon),      intent(out) :: xlons
  real, dimension(mlon,mlat), intent(out) :: sice

! Declare local parameters
  integer,parameter:: lu_sice = 22   ! FORTRAN unit number of GRIB SST file
  real, parameter :: deg2rad = 3.141593/180.0

! Declare local variables and arrays
  logical(1), allocatable, dimension(:) ::  lb

  integer :: nlat_sice,nlon_sice,ilon,jlat
  integer :: iret,ni,nj
  integer :: mscan,kb1
  integer :: jincdir,i,iincdir,kb2,kb3,kf,kg,k,j,jf
  integer, dimension(22):: jgds,kgds
  integer, dimension(25):: jpds,kpds

  real :: x0,y0,dres
  real, allocatable, dimension(:) :: f
  
!************+******************************************************************************
!
! Open SST analysis file (GRIB)
  call baopenr(lu_sice,trim(file_sice),iret)
  if (iret /= 0 ) then
     write(6,*)'RDGRBICE:  ***ERROR*** opening Sea Ice file'
     stop
  endif
! Define SST variables for read
  j=-1
  jpds=-1
  jgds=-1
  jpds(5)=91          ! Sea Ice fraction (1=ice; 0=no ice)
  jpds(6)=102         ! surface
  call getgbh(lu_sice,0,j,jpds,jgds,kg,kf,k,kpds,kgds,iret)

  nlat_sice = kgds(3)   ! number points on longitude circle (360)
  nlon_sice = kgds(2)   ! number points on latitude circle (720)

  if ( mlat /= nlat_sice .or. mlon /= nlon_sice ) then
     write(*,*) 'STOP inconsistent dimensions : ',nlat_sice, nlon_sice,mlat,mlon
     stop
  endif
!
! get xlats and xlons
!
  dres = 180.0/real(mlat)
  y0 = 0.5*dres-90.0
  x0 = 0.5*dres

! Get lat_sst & lon_sst
  do jlat = 1, mlat
     xlats(jlat) = y0 + real(jlat-1)*dres
  enddo

  do ilon = 1, mlon
     xlons(ilon) = (x0 + real(ilon-1)*dres)
  enddo
! Allocate arrays
  allocate(lb(nlat_sice*nlon_sice))
  allocate(f(nlat_sice*nlon_sice))
  jf=nlat_sice*nlon_sice

! Read in the analysis
  call getgb(lu_sice,0,jf,j,jpds,jgds,kf,k,kpds,kgds,lb,f,iret)
  if (iret /= 0) then
     write(6,*)'RDGRBSICE:  ***ERROR*** reading sice analysis data record'
     deallocate(lb,f)
     stop
  endif

! Load dimensions and grid specs.  Check for unusual values
  ni=kgds(2)                ! 720 for 0.5 x 0.5
  nj=kgds(3)                ! 360 for 0.5 x 0.5 resolution

  mscan=kgds(11)
  kb1=ibits(mscan,7,1)   ! i scan direction
  kb2=ibits(mscan,6,1)   ! j scan direction
  kb3=ibits(mscan,5,1)   ! (i,j) or (j,i)

! Get i and j scanning directions from kb1 and kb2.
! 0 yields +1, 1 yields -1. +1 is west to east, -1 is east to west.
  iincdir = 1-kb1*2

! 0 yields -1, 1 yields +1. +1 is south to north, -1 is north to south.
  jincdir = kb2*2 - 1
  do k=1,kf

!    kb3 from scan mode indicates if i points are consecutive
!    or if j points are consecutive
     if(kb3==0)then     !  (i,j)
        i=(ni+1)*kb1+(mod(k-1,ni)+1)*iincdir
        j=(nj+1)*(1-kb2)+jincdir*((k-1)/ni+1)
     else                !  (j,i)
        j=(nj+1)*(1-kb2)+(mod(k-1,nj)+1)*jincdir
        i=(ni+1)*kb1+iincdir*((k-1)/nj+1)
     endif
     sice(i,j)=f(k)
  end do
     
  deallocate(lb,f)
  
end subroutine read_seaice_grb