subroutine da_detsurtyp (snow,ice,landsea, vegtyp, soiltyp, is, ie, js, je, &
   i, j, dx, dy, dxm, dym, isflg,ob_vegtyp,ob_soiltyp, seap, icep, lndp, snop)

   !---------------------------------------------------------------------------
   ! Purpose: determine surface type at obs locations.
   !
   ! METHOD: using the background information 
   !      1. get the background landmask/snow/sea-ice
   !      2. calculate percentage of sea/ice/land/snow
   !      3. determine surface type at obs location
   !      4. find nearest grid vegtype and soil type 
   !
   !  HISTORY: 11/25/2008 - Use input surface type percentage      Tom Auligne
   !
   !---------------------------------------------------------------------------

   implicit none

   integer, intent(in)  :: is, ie, js, je
   integer, intent(in)  :: i, j
   real   , intent(in)  :: dx, dxm, dy, dym
   real   , intent(in)  :: snow(is:ie,js:je)    ! Input variable
   real   , intent(in)  :: ice(is:ie,js:je)     ! Input variable
   real   , intent(in)  :: landsea(is:ie,js:je) ! Input variable
   integer, intent(in)  :: vegtyp(is:ie,js:je)  ! Input variable
   integer, intent(in)  :: soiltyp(is:ie,js:je) ! Input variable
   integer, intent(out) :: isflg                ! Output variable 
   real,    intent(out) :: ob_vegtyp            ! Output variable
   real,    intent(out) :: ob_soiltyp            ! Output variable
   real*8,  intent(out) :: seap, icep, lndp, snop ! percentage of surface type
   
   !     isflg    - surface flag
   !                0 sea
   !                1 sea ice
   !                2 land
   !                3 snow
   !                4 mixed predominately sea
   !                5 mixed predominately sea ice
   !                6 mixed predominately land
   !                7 mixed predominately snow

   !  local variables
!   integer   ::  n, xbflag(4)   ! surface type at xb location
!                               ! 0:sea  1:sea-ice 2:land  3:snow
   real      ::  w(4),minw      ! weight at 4 xb locations

   if (trace_use) call da_trace_entry("da_detsurtyp")    

   ! 1.0 determine surface type of xb at 4 location around obs
   !-------------------------------------------------------
!   if ( nint(landsea(i  ,j  )) == 0 ) xbflag(1) = 0 ! sea
!   if ( nint(landsea(i+1,j  )) == 0 ) xbflag(2) = 0 ! sea
!   if ( nint(landsea(i  ,j+1)) == 0 ) xbflag(3) = 0 ! sea
!   if ( nint(landsea(i+1,j+1)) == 0 ) xbflag(4) = 0 ! sea
!
!   if ( nint(landsea(i  ,j  )) == 1 ) xbflag(1) = 2 ! land/snow/sea-ice
!   if ( nint(landsea(i+1,j  )) == 1 ) xbflag(2) = 2 ! land/snow/sea-ice
!   if ( nint(landsea(i  ,j+1)) == 1 ) xbflag(3) = 2 ! land/snow/sea-ice
!   if ( nint(landsea(i+1,j+1)) == 1 ) xbflag(4) = 2 ! land/snow/sea-ice
!
!   if ( nint(snow(i  ,j  )) == 1 ) xbflag(1) = 3 ! snow
!   if ( nint(snow(i+1,j  )) == 1 ) xbflag(2) = 3 ! snow
!   if ( nint(snow(i  ,j+1)) == 1 ) xbflag(3) = 3 ! snow
!   if ( nint(snow(i+1,j+1)) == 1 ) xbflag(4) = 3 ! snow
!
!   if ( nint(ice(i  ,j  )) == 1 ) xbflag(1) = 1 ! sea-ice
!   if ( nint(ice(i+1,j  )) == 1 ) xbflag(2) = 1 ! sea-ice
!   if ( nint(ice(i  ,j+1)) == 1 ) xbflag(3) = 1 ! sea-ice
!   if ( nint(ice(i+1,j+1)) == 1 ) xbflag(4) = 1 ! sea-ice 

   ! 2.0 determine surface type percentage at obs location
   !------------------------------------------------------
   !  (i,j+1) 
   !    -----------------(i+1,j+1)
   !    |   w2          |
   !    |dym        w1  |
   !    |      obs      |
   !    |------- *      |
   !    |dy  w4  |  w3  |
   !    |   dx   | dxm  |
   !    |----------------
   !   (i,j)            (i+1,j)
   !
   !--------------------------------------------------------

!   seap = 0.0
!   icep = 0.0
!   snop = 0.0
!   lndp = 0.0
   w(1) = dym*dxm   ! weight for point (i,j)
   w(2) = dym*dx    ! weight for point (i+1,j)
   w(3) = dy *dxm   ! weight for point (i,j+1)
   w(4) = dy *dx    ! weight for point (i+1,j+1)

!   do n = 1, 4
!      if (xbflag(n) == 0) seap = seap+w(n)
!      if (xbflag(n) == 1) icep = icep+w(n)
!      if (xbflag(n) == 2) lndp = lndp+w(n)
!      if (xbflag(n) == 3) snop = snop+w(n)
!   end do
   
   lndp = 1.0/SUM(w) * ( w(1)*landsea(i,j)         + w(2)*landsea(i+1,j)          + &
                         w(3)*landsea(i,j+1)       + w(4)*landsea(i+1,j+1) )
   
   icep = 1.0/SUM(w) * ( w(1)*ice(i,j)              + w(2)*ice(i+1,j)              + &
                         w(3)*ice(i,j+1)            + w(4)*ice(i+1,j+1) )
   
   snop = 1.0/SUM(w) * ( w(1)*MIN(snow(i,j),   1.0) + w(2)*MIN(snow(i+1,j),   1.0) + &
                         w(3)*MIN(snow(i,j+1), 1.0) + w(4)*MIN(snow(i+1,j+1), 1.0) )
   
   if (icep > 0.0) then
      snop = MIN(snop, 1.0 - icep)
      seap = 1.0 - icep - snop
   else
      seap = 1.0 - lndp
   end if   
   lndp = MAX(1.0_8-seap-icep-snop, 0.0_8)
   
   ! fo2d   = dym*(dxm*fm2d(i,j  ) + dx*fm2d(i+1,j  )) &
   !         + dy *(dxm*fm2d(i,j+1) + dx*fm2d(i+1,j+1))

   ! 3.0 determine final surface flag at obs location
   !-----------------------------------------
   if (seap >= 0.99) isflg = 0
   if (icep >= 0.99) isflg = 1
   if (lndp >= 0.99) isflg = 2
   if (snop >= 0.99) isflg = 3
   if ( .not. (seap >= 0.99) .and. &
        .not. (icep >= 0.99) .and. &
        .not. (lndp >= 0.99) .and. &
        .not. (snop >= 0.99) ) then
      if (seap > lndp) then
         if (seap > icep) then
            if (seap > snop) then
               isflg = 4
            else
               isflg = 7
            end if
         else
            if (icep > snop) then
               isflg = 5
            else
               isflg = 7
            end if
         end if
      else
         if (lndp > icep) then
            if (lndp > snop) then
               isflg = 6
            else
               isflg = 7
            end if
         else
            if (icep > snop) then
               isflg = 5
            else
               isflg = 7
            end if
         end if
      end if 
   end if
   
   ! 4.0 find nearest grid vegtype and soil type
   !         at obs location
   !-----------------------------------------
   minw=min(w(1),w(2),w(3),w(4))
   if (minw == w(1)) then
      ob_vegtyp  = float(vegtyp (i+1,j+1))
      ob_soiltyp = float(soiltyp(i+1,j+1))
   else if (minw == w(2)) then
      ob_vegtyp  = float(vegtyp (i,j+1))
      ob_soiltyp = float(soiltyp(i,j+1))
   else if (minw == w(3)) then
      ob_vegtyp  = float(vegtyp (i+1,j))
      ob_soiltyp = float(soiltyp(i+1,j))
   else if (minw == w(4)) then
      ob_vegtyp  = float(vegtyp (i,j))
      ob_soiltyp = float(soiltyp(i,j))
   end if

   if (trace_use) call da_trace_exit("da_detsurtyp")    

end subroutine da_detsurtyp