#ifdef RTTOV
subroutine da_rttov_k(inst, nchanl, nprofiles, nlevels, &
                      con_vars, aux_vars, &
                      tb, rad_xb, rad_ovc, emissivity, emissivity_out)

   !---------------------------------------------------------------------------
   !  PURPOSE: interface to the forward and k matrix subroutine of RTTOV
   !---------------------------------------------------------------------------

   implicit none

   integer             , intent (in) :: inst, nchanl, nprofiles, nlevels
   type (con_vars_type), intent (inout) :: con_vars(nprofiles)
   type (aux_vars_type), intent (in) :: aux_vars(nprofiles)
   real                , intent (inout) :: tb(nchanl,nprofiles)
   real                , intent (inout) :: rad_xb(nchanl,nprofiles)
   real                , intent (inout) :: rad_ovc(nchanl,nlevels-1,nprofiles)
   real                , intent (in)  :: emissivity(nchanl*nprofiles)
   real                , intent (out) :: emissivity_out(nchanl*nprofiles)

   ! local variables
   integer             :: nchanprof, asw
   integer             :: n, i, j, k
   integer             :: alloc_status(10)
  
   ! RTTOV input parameters
   type (rttov_chanprof), allocatable :: chanprof(:)
   type (profile_type),   allocatable :: profiles(:), profiles_k(:)
   logical, allocatable               :: calcemis(:)

   ! RTTOV out parameters
   integer             :: errorstatus

   ! RTTOV inout parameters
   real,                  allocatable :: emissivity_k(:), emissivity_out_k(:)
   type (radiance_type)               :: radiance, radiance_k
   type (transmission_type)           :: transmission, transmission_k

   call da_trace_entry("da_rttov_k")

   nchanprof = nchanl*nprofiles

   ! generate the chanprof array
   allocate (chanprof(nchanprof))
   do n = 1, nprofiles
      chanprof((n-1)*nchanl+1:n*nchanl)%prof = n
      chanprof((n-1)*nchanl+1:n*nchanl)%chan = (/ (j, j=1,nchanl) /)
   end do

   alloc_status (:) = 0

   ! allocate input profile arrays with the number of levels.
   asw = 1  ! allocate
   allocate (profiles(nprofiles), stat=alloc_status(1))
   call rttov_alloc_prof(        &
      & errorstatus,             &
      & nprofiles,               &
      & profiles,                &
      & nlevels,                 &
      & opts(inst),              &
      & asw,                     &
      & coefs = coefs(inst),     &  ! mandatory if either opts%addclouds or opts%addaerosl is true
      & init = .true.            )  ! additionally initialize profiles structure
   if ( errorstatus /= errorstatus_success .or. alloc_status(1) /= 0 ) then
     call da_error(__FILE__,__LINE__, &
        (/"memory allocation error for profile arrays"/))
   end if

   asw = 1  ! allocate
   allocate (profiles_k(nchanprof), stat=alloc_status(2))
   call rttov_alloc_prof(        &
      & errorstatus,             &
      & nchanprof,               &
      & profiles_k,              &
      & nlevels,                 &
      & opts(inst),              &
      & asw,                     &
      & coefs = coefs(inst),     &  ! mandatory if either opts%addclouds or opts%addaerosl is true
      & init = .true.            )  ! additionally initialize profiles structure
   if ( errorstatus /= errorstatus_success .or. alloc_status(2) /= 0 ) then
     call da_error(__FILE__,__LINE__, &
        (/"memory allocation error for profile_k arrays"/))
   end if

   do n = 1, nprofiles
      profiles(n) % p(:)       = coefs(inst)%coef%ref_prfl_p(:)
      profiles(n) % t(:)       = con_vars(n)%t(:)
      profiles(n) % q(:)       = con_vars(n)%q(:)
      if ( opts(inst) % ozone_data ) profiles(n) % o3(:)      = 0.0 !con_vars(n)%o3(:)
      if ( opts(inst) % co2_data )   profiles(n) % co2(:)     = 0.0 !con_vars(n)%co2(:)
      if ( opts(inst) % clw_data )   profiles(n) % clw(:)     = 0.0 !con_vars(n)%clw(:)
      if ( opts(inst) % n2o_data )   profiles(n) % n2o(:)     = 0.0
      if ( opts(inst) % co_data )    profiles(n) % co(:)      = 0.0
      if ( opts(inst) % co_data )    profiles(n) % ch4(:)     = 0.0
      if ( opts(inst) % addaerosl )  profiles(n) % aerosols(:,:)   = 0.0
      if ( opts(inst) % addclouds )  then
         profiles(n) % cloud(:,:)      = 0.0
         profiles(n) % cfrac(:,:)      = 0.0
         profiles(n) % idg             = 1
         profiles(n) % ish             = 1
      end if

      profiles(n)% skin % surftype   = aux_vars(n) % surftype
      if ( profiles(n)% skin % surftype == 1 ) then
         if ( opts(inst) % addsolar ) then
            ! if refelcted solar radiation is to be included in the SWIR channels, then
            ! specification of fresh or salt water needs to be provided
            profiles(n) % skin % watertype = 1
         end if
      end if

      ! for microwave channels, land/sea-ce emissivity is computed
      ! from coefs in prof%skin%fastem, if calcemis = True
      if ( coefs(inst)%coef%id_sensor == sensor_id_mw ) then
         if ( profiles(n) % skin % surftype == 2 ) then
            profiles(n) % skin % fastem (1) = 2.2
            profiles(n) % skin % fastem (2) = 3.7
            profiles(n) % skin % fastem (3) = 122.0
            profiles(n) % skin % fastem (4) = 0.0
            profiles(n) % skin % fastem (5) = 0.15
         else if ( profiles(n) % skin % surftype == 0 ) then
            profiles(n) % skin % fastem (1) = 3.0
            profiles(n) % skin % fastem (2) = 5.0
            profiles(n) % skin % fastem (3) = 15.0
            profiles(n) % skin % fastem (4) = 0.1
            profiles(n) % skin % fastem (5) = 0.3
         end if
      end if

      profiles(n) % skin % t          = aux_vars (n) % surft    
      !profiles(n) % skin % fastem (:) = 0.0 ! aux_vars (n) % fastem (:)

      profiles(n) % s2m % t   = aux_vars (n) % t2m
      profiles(n) % s2m % q   = aux_vars (n) % q2m
      profiles(n) % s2m % o   = 0.0                 ! o3, never used
      profiles(n) % s2m % p   = con_vars (n) % ps
      profiles(n) % s2m % u   = aux_vars (n) % u10
      profiles(n) % s2m % v   = aux_vars (n) % v10

      profiles(n) % zenangle  = aux_vars (n) % satzen
      profiles(n) % elevation = 0.001* aux_vars(n) % elevation   ! km
      profiles(n) % latitude  = aux_vars(n) % rlat

      if ( opts(inst) % addsolar ) then
         profiles(n) % azangle     = aux_vars (n) % satazi
         profiles(n) % sunzenangle = aux_vars (n) % solzen     !50.0
         profiles(n) % sunazangle  = aux_vars (n) % solazi     !86.0
         profiles(n) % s2m % wfetc = 100000.0  ! m
      end if

      profiles(n) % Be          = 0.35   ! optional, for zeeman effect for ssmis and amsua
      profiles(n) % cosbk       = 0.0    ! optional, for zeeman effect for ssmis and amsua

      profiles(n) % ctp         = 500.0  ! hPa, optional, for simple cloud
      profiles(n) % cfraction   = 0.0    ! 0-1, optional, for simple cloud
   end do

   allocate (emissivity_k(nchanprof), stat=alloc_status(3))
   allocate (emissivity_out_k(nchanprof), stat=alloc_status(4))
   allocate (calcemis(nchanprof), stat=alloc_status(5))
   if ( any( alloc_status /= 0 ) ) then
      call da_error(__FILE__,__LINE__, &
         (/"memory allocation error for emissivity arrays"/))
   end if

   ! allocate radiance results arrays with number of channels
   asw = 1 ! allocate
   call rttov_alloc_rad( &
      & errorstatus,     &
      & nchanprof,       &
      & radiance,        &
      & nlevels-1,       &
      & asw,             &
      & init = .true. )
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
         (/"memory allocation error for radiance arrays"/))
   end if

   asw = 1 ! allocate
   call rttov_alloc_rad( &
      & errorstatus,     &
      & nchanprof,       &
      & radiance_k,      &
      & nlevels-1,       &
      & asw,             &
      & init = .true. )
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
         (/"memory allocation error for radiance_k arrays"/))
   end if

   ! allocate transmittance structure
   asw = 1 ! allocate
   call rttov_alloc_transmission( &
      & errorstatus,              &
      & transmission,             &
      & nlevels-1,                &
      & nchanprof,                &
      & asw,                      &
      & init = .true. )
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
         (/"memory allocation error for transmittance arrays"/))
   end if

   asw = 1 ! allocate
   call rttov_alloc_transmission( &
      & errorstatus,              &
      & transmission_k,           &
      & nlevels-1,                &
      & nchanprof,                &
      & asw,                      &
      & init = .true. )
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
         (/"memory allocation error for transmittance K arrays"/))
   end if

   ! calculate emissivity where the input emissivity value is less than 0.01
   calcemis(:) = emissivity(:) < 0.01

   emissivity_k(:) = 0.0
   radiance_k % bt_clear(:)   = 1.0
   radiance_k % clear(:)      = 1.0
   radiance_k % overcast(:,:) = 1.0

   !-----------------------------------
   !  calling RTTOV forward and K model
   !----------------------------------

   call rttov_k(      &
      & errorstatus,       &   ! out
      & chanprof,          &   ! in     chanprof(nchanprof)
      & opts(inst),        &   ! in
      & profiles,          &   ! in     profiles(nprof)
      & profiles_k,        &   ! inout  profiles_k(nchanprof)
      & coefs(inst),       &   ! in
      & calcemis,          &   ! in     calcemis(nchanprof)
      & emissivity,        &   ! in     emissivity(nchanprof)
      & emissivity_k,      &   ! inout  emissivity_k(nchanprof)
      & emissivity_out,    &   ! out    emissivity_out(nchanprof)
      & emissivity_out_k,  &   ! inout  emissivity_out_k(nchanprof)
      & transmission,      &   ! inout
      & transmission_k,    &   ! inout
      & radiance,          &   ! inout
      & radiance_k )           ! inout
   
   if ( print_detail_rad .or. errorstatus /= errorstatus_success ) then
      WRITE (UNIT=stderr,FMT=*) 'rttov_k error code = ', errorstatus
      WRITE (UNIT=stderr,FMT=*) 'nchanprof    = ', nchanprof
      WRITE (UNIT=stderr,FMT=*) 'nprofiles    = ', nprofiles
      WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%t = ', profiles(1)%s2m%t
      WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%q = ', profiles(1)%s2m%q
      WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%o = ', profiles(1)%s2m%o
      WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%p = ', profiles(1)%s2m%p
      WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%u = ', profiles(1)%s2m%u
      WRITE (UNIT=stderr,FMT=*) 'profiles%s2m%v = ', profiles(1)%s2m%v
      WRITE (UNIT=stderr,FMT=*) 'profiles%skin%surftype = ', profiles(1)%skin%surftype
      WRITE (UNIT=stderr,FMT=*) 'profiles%skin%t        = ', profiles(1)%skin%t
      WRITE (UNIT=stderr,FMT=*) 'profiles%skin%fastem   = ', profiles(1)%skin%fastem
      WRITE (UNIT=stderr,FMT=*) 'profiles%zenangle = ', profiles(1)%zenangle
      WRITE (UNIT=stderr,FMT=*) 'profiles%azangle = ', profiles(1)%azangle
      WRITE (UNIT=stderr,FMT=*) 'profiles%p   = ', profiles(1)%p
      WRITE (UNIT=stderr,FMT=*) 'profiles%t   = ', profiles(1)%t
      WRITE (UNIT=stderr,FMT=*) 'profiles%q   = ', profiles(1)%q
      WRITE (UNIT=stderr,FMT=*) 'calcemis     = ', calcemis
      WRITE (UNIT=stderr,FMT=*) 'emissivity   = ', emissivity
      WRITE (UNIT=stderr,FMT=*) 'emissivity_out = ', emissivity_out
      WRITE (UNIT=stderr,FMT=*) 'radiance = ', radiance%bt_clear
      if ( errorstatus /= errorstatus_success ) call da_error(__FILE__,__LINE__,(/"Problem in rttov_k"/))
   end if

   do n = 1, nprofiles
      tb(1:nchanl,n) = radiance % bt_clear((n-1)*nchanl+1:n*nchanl)
      rad_xb(1:nchanl,n) = radiance % clear((n-1)*nchanl+1:n*nchanl)
      do k = 1, nlevels-1
         rad_ovc(1:nchanl,k,n) = radiance % overcast(k,(n-1)*nchanl+1:n*nchanl)
      end do	
      if ( use_rttov_kmatrix ) then
         do i = 1, nchanl
            do k = 1, nlevels
               con_vars(n) % t_jac(i,k) = profiles_k((n-1)*nchanl+i)%t(k)
               con_vars(n) % q_jac(i,k) = profiles_k((n-1)*nchanl+i)%q(k)
            end do
            con_vars(n) % ps_jac(i) = profiles_k((n-1)*nchanl+i)%s2m%p
         end do
      end if
   end do

   deallocate (calcemis)
   deallocate (chanprof)
   deallocate (emissivity_k)
   deallocate (emissivity_out_k)

   asw = 0 ! deallocation
   ! deallocate radiance arrays
   call rttov_alloc_rad (errorstatus,nchanprof,radiance,nlevels-1,asw)
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
        (/"radiance deallocation error"/))
   end if
   call rttov_alloc_rad (errorstatus,nchanprof,radiance_k,nlevels-1,asw)
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
        (/"radiance K deallocation error"/))
   end if

   ! deallocate transmission arrays
   call rttov_alloc_transmission (errorstatus,transmission,nlevels-1,nchanprof,asw)
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
        (/"transmission deallocation error"/))
   end if
   call rttov_alloc_transmission (errorstatus,transmission_k,nlevels-1,nchanprof,asw)
   if ( errorstatus /= errorstatus_success ) then
      call da_error(__FILE__,__LINE__, &
        (/"transmission K deallocation error"/))
   end if

   ! deallocate profile arrays
   call rttov_alloc_prof (errorstatus,nprofiles,profiles,nlevels,opts(inst),asw)
   deallocate(profiles,stat=alloc_status(6))
   if ( errorstatus /= errorstatus_success .or. alloc_status(6) /= 0 ) then
      call da_error(__FILE__,__LINE__, &
        (/"profile deallocation error"/))
   end if
   call rttov_alloc_prof (errorstatus,nchanprof,profiles_k,nlevels,opts(inst),asw)
   deallocate(profiles_k,stat=alloc_status(7))
   if ( errorstatus /= errorstatus_success .or. alloc_status(7) /= 0 ) then
      call da_error(__FILE__,__LINE__, &
        (/"profile deallocation error"/))
   end if

   call da_trace_exit("da_rttov_k")
   
end subroutine da_rttov_k
#endif