subroutine da_intpsfc_tem (val, ho, po, to, height, pressure, temp, kts, kte)

   !-----------------------------------------------------------------------
   ! Purpose: Correct temperature between two levels.
   !
   ! Reference: 
   ! ---------
   ! The use of surface observations in four dimensional data assmiilation
   !  Using a mesoscale model.
   !  Ruggiero et al., 1996, Monthly Weather Review, Volume 124, 1018-1033
   !
   !----------------------------------------------------------------------------

   implicit none

   real,    intent (out) :: val
   integer, intent (in)  :: kts, kte
   real,    intent (in)  :: ho, po, to
   real,    intent (in)  :: height(kts:kte)
   real,    intent (in)  :: pressure(kts:kte)
   real,    intent (in)  :: temp(kts:kte)

   real    :: prs_mb(kts:kte)
   ! calculated but never used
   !      real    :: dth_12, dth_21, dth_sfc, dth_obs
   !     real    :: dhe_12, dhe_21, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
   real    :: dth_21, dth_sfc, dth_obs
   real    :: dhe_12, dhe_sfc1, dhe_obs1, dhe_sfc2, dhe_obs2
   real    :: th_100mb, th_200mb, th_obs, th_sfc
   real    :: th_obs_int, th_sfc_int
   real    :: pdif, rcp
   integer :: k_100mb, k_200mb, jk
   integer :: inc_100mb, inc_200mb

   if (trace_use_dull) call da_trace_entry("da_intpsfc_tem")

   rcp = gas_constant/cp

   ! 1.Find levels: model surface + 100hpa and model surface + 200hpa ar obs loc
   ! ===========================================================================

   ! 1.1 Convert model pressure profile from Pa to hPa  

   prs_mb = pressure / 100.0

   ! 1.2  Find levels surface + 100hPA

   inc_100mb = 100.0
   k_100mb   = kts

   do jk =  kts+1, kte
      pdif = prs_mb (kts) - prs_mb (jk)
      if (pdif .GE. inc_100mb) then
         k_100mb = jk
         exit
      end if
   end do

   ! 1.2  Find levels surface + 200hPA

   inc_200mb = 200.0
   k_200mb   = kts

   do jk =  kts+1, kte
      pdif = prs_mb (kts) - prs_mb (jk)
      if (pdif .GE. inc_200mb) then
         k_200mb = jk
         exit
      end if
   end do

   ! 1.3  Check consistency 

   if ((k_100mb .LE. kts) .OR. (k_200mb .LE. kts) .OR. &
        (k_200mb .LT. k_100mb)) then
      write (unit=message(1),fmt='(A)') ' Cannot find sfc + 200hPa and sfc + 100hPa' 
      write (unit=message(2),fmt='(A,I2,A,F10.3)') ' P (',k_200mb,') = ',prs_mb (k_200mb) 
      write (unit=message(3),fmt='(A,I2,A,F10.3)') ' P (',k_100mb,') = ',prs_mb (k_100mb) 
      write (unit=message(4),fmt='(A,F10.3)')      ' P_SFC  = ',         prs_mb (kts) 
      call da_warning(__FILE__,__LINE__,message(1:4))
      val = missing_r
      return
   end if

   ! 2.  potential temperature 
   ! =========================

   ! 2.1 Potential temperature at 100hPa above model surface

   th_100mb = temp (k_100mb) * (1000.0 / prs_mb (k_100mb))**rcp

   ! 2.2 Potential temperature at 200hPa above model surface

   th_200mb = temp (K_200mb) * (1000.0 / prs_mb (k_200mb))**rcp

   ! 2.3 Potential temperature at observation location

   th_obs   = to * (1000.0 / (po/100.0)) ** rcp


   ! 3.  lapse rate between surface+200hpa and surface+100hpa
   ! =========================================================

   ! 3.1 Potential temperature increment
    
   dth_21 = th_100mb - th_200mb
   ! never used
   ! dth_12 = th_200mb - th_100mb

   ! 3.1 Height increments
    
   ! never used
   ! dhe_21   = height (k_100mb)- height (k_200mb)
   dhe_sfc1 = height (k_100mb)- height (kts)
   dhe_obs1 = height (k_100mb)- ho

   dhe_12   = height (k_200mb)- height (k_100mb)
   dhe_sfc2 = height (k_200mb)- height (kts)
   dhe_obs2 = height (k_200mb)- ho

   ! 3.2 Extrapolated potential temperature at model surface and observation loc

   th_sfc_int = th_100mb + (dth_21/dhe_12) * dhe_sfc1 
   th_obs_int = th_100mb + (dth_21/dhe_12) * dhe_obs1 


   ! 4.  Bring temperature onto model surface
   ! ========================================

   ! 4.1 Difference at observation locations

   dth_obs = th_obs_int - th_obs

   ! 4.2 Difference at model surface

   dth_sfc = (dhe_sfc2/dhe_obs2) * dth_obs

   ! 4.3 Potentiel temperature brought to model surface

   th_sfc = th_sfc_int - dth_sfc

   ! 4.3 Corresponding Temperature

   val  = th_sfc * (prs_mb (kts) / 1000.0)**rcp

   if (trace_use_dull) call da_trace_exit("da_intpsfc_tem")

end subroutine da_intpsfc_tem