
  subroutine dasray_intf( &
       n,                 nsmax,        ns_elas,          &
       maxintstr,         maxstr,       maxstr_internal,  &
       ngeomax,           mmax,         maxleg,           & 
       maxlay,            nstokes,      accuracy,         &
       dnu,               tau,          omega_Ray,        & 
       omega_Cab,         om_air,       scat_type,        &
       weightCabRam_air,  PCab_N2,      PRay_N2,          & 
       PCab_O2,           PRay_O2,      PCab_Ar,          & 
       PRay_Ar,           PCab_air,     PRay_air,         & 
       PRam_air,          maxcoeffs,    maxcoeffstrunc,   &
       surf_type,         alb_Lambert,  gmu,              &
       gwt,               gmumin,       phi,              &
       imuv,              imu0,         plm_0,            &
       plm_m,             plm_p,        u,                & 
       refl)

!--------------------------------------------------------------
! Subroutine:
!   Doubling-adding solver (elastic Rayleigh scattering)
!
! Purpose:
!   Solves radiative transfer in plane-parallel atmosphere
!   in the scalar approximation and includes inelastic 
!   Raman scattering
!
! Author: 
!   Rutger van Deelen
!   SRON - Netherlands Insititute for Space Research
!
! Reference: 
!   R. van Deelen, J. Landgraf, and I. Aben --
!   "Multiple elastic and inelastic light-scattering in the 
!   Earth's atmosphere: A doubling-adding method to include 
!   rotational Raman scattering by air", 
!   J. Quant. Spectrosc. Radiat. Transfer 95, 399-433, 2005
!------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision

! Input ==========================================================

  integer,intent(in):: &
     n,          & !wavenr index, coarse grid
     nsmax,      & !nr of wavenrs 0:nsmax  
     ns_elas,    & !wavenr index corresponding to elastically scattererd light
     maxintstr,  & !Nr of internal streams (must be even and >= 2)
     maxstr,     &
     maxstr_internal, &
     ngeomax,    & !Nr of viewing geometries
     mmax,       & !maximum Fourier moment (molecular scattering=2)
     maxleg,     & !maximum Legendre moment (molecular scatt.=2)
     maxlay,     & !Nr of atmospheric layers
     nstokes       !NR of stokes parameters

! wavenumber binwdith [cm-1]
  real(kind=dp),intent(in):: &
     dnu

! viewing geometries
  real(kind=dp)::  &
     phi(ngeomax)

! Gaussian quadrature
  real(kind=dp)::gmu(maxstr),gwt(maxstr),gmumin 

! points to stream index for viewing angles and solar angles
  integer,intent(in):: &
     imuv(ngeomax), &
     imu0(ngeomax)

! needed for polarization and single scattering solution
  real(kind=dp),dimension(0:maxleg,0:maxleg,maxintstr+4*ngeomax),intent(in):: &
     plm_0, &
     plm_m, &
     plm_p
  real(kind=dp),dimension(maxintstr+4*ngeomax),intent(in):: &
     u 

! doubling-adding accuracy
  real(kind=dp),intent(in):: &
      accuracy

! Internal ========================================================

  logical,parameter:: &
      single_scattering = .false.

  integer:: &
    m,   & !index of Fourier moment m
    k,   & !layer index
    ik,  & !layer index, reverse order
    ikk, & !layer index (fine layers for doubling)
    i,i1,i2,ns,ns1,ns2,lg,i_st,iv,is

!--------------------------------------------------------------
! Optical depth, single scattering albedo, etc.

  real(kind=dp),dimension(maxlay,0:nsmax),intent(in):: &
     tau,       & !total optical depth
     omega_Ray, & !Rayleigh single scattering albedo 
     omega_Cab    !Cabannes single scattering albedo 

  real(kind=dp),intent(in):: &
     om_air(maxlay,0:nsmax,0:nsmax)

! 0=not good; 1=Cab; 2=Ram; 3=mixCab&Ram; 4=Ray
  integer,intent(in):: &
     scat_type(0:nsmax,0:nsmax) 
  
! weight CabRam for eff. phasematrix
  real(kind=dp),intent(in):: &
     weightCabRam_air(maxlay,0:nsmax)

!--------------------------------------------------------------
! phase matrices, depend on i,j, m, k and n
  real(kind=dp),intent(in), &
   dimension(nstokes,nstokes,0:maxleg,maxlay):: &
    PCab_N2,  & !Cabannes phase fct. moments
    PRay_N2,  & !Rayleigh phase fct. moments
    PCab_O2,  & !Cabannes phase fct. moments
    PRay_O2,  & !Rayleigh phase fct. moments
    PCab_Ar,  & !Cabannes phase fct. moments
    PRay_Ar,  & !Rayleigh phase fct. moments
    PCab_air, & !Cabannes phase fct. moments
    PRay_air, & !Rayleigh phase fct. moments
    PRam_air    !Raman    phase fct. moments

! Nr of Legendre terms for each layer
! Nr of coefficient per layer (clear sky = 0:2)
  integer,intent(in):: &
    maxcoeffs(maxlay), &    !all terms 
    maxcoeffstrunc(maxlay)  !truncate sum for diffuse radiation

  real(kind=dp), &
   dimension(0:maxleg):: &
    wl_Cab_air, & !Cabannes phase fct. moments
    wl_Ray_air, & !Rayleigh phase fct. moments
    wl_Ram_air    !Raman    phase fct. moments

! phase matrices, depend on n, m
  real(kind=dp):: &
     Zm_Ray_min(maxstr,maxstr),   &
     Zm_Cab_min(maxstr,maxstr),   &
     Zm_Ram_min(maxstr,maxstr),   &
     Zm_Ray_plus(maxstr,maxstr),  &
     Zm_Cab_plus(maxstr,maxstr),  &
     Zm_Ram_plus(maxstr,maxstr)

! mixture of Cab and Ram, depends on n, m, k, mix for each ns
  real(kind=dp):: &
     Zm_CabRam_min(0:nsmax,maxstr,maxstr),  &
     Zm_CabRam_plus(0:nsmax,maxstr,maxstr) 

!-----------------------------------------------------------

  integer:: &
     ndl !nr of doublings

  real(kind=dp):: &
     dtau(0:nsmax)    !total (sca+abs) optical depth per layer
  
  real(kind=dp):: &
     dtau_A(0:nsmax),  &
     tau_atm(0:nsmax), &
     tau_AB(0:nsmax),  &
     tau_wl(0:nsmax)

!-----------------------------------------------------------
! Reflectance
 
! whole layer ,  single scattering approximation
  real(kind=dp):: &
     R1m_wl(0:nsmax,0:nsmax,ngeomax),     &
     T1m_wl(0:nsmax,0:nsmax,ngeomax),     &
     R1mStar_wl(0:nsmax,0:nsmax,ngeomax), &
     T1mStar_wl(0:nsmax,0:nsmax,ngeomax)

! atm only
  real(kind=dp):: & 
     R1m_atm(0:nsmax,0:nsmax,ngeomax),     &
     T1m_atm(0:nsmax,0:nsmax,ngeomax),     &
     R1mStar_atm(0:nsmax,0:nsmax,ngeomax), &
     T1mStar_atm(0:nsmax,0:nsmax,ngeomax)

! atm+surf
  real(kind=dp):: &
     R1m_sys(0:nsmax,0:nsmax,ngeomax),     &
     T1m_sys(0:nsmax,0:nsmax,ngeomax),     &
     D1m_sys(0:nsmax,0:nsmax,ngeomax),     &
     R1mStar_sys(0:nsmax,0:nsmax,ngeomax), &
     T1mStar_sys(0:nsmax,0:nsmax,ngeomax)

! initial layer to start doubling
  real(kind=dp):: &
     R1m_bl(0:nsmax,0:nsmax,maxstr,maxstr), &
     T1m_bl(0:nsmax,0:nsmax,maxstr,maxstr)

  real(kind=dp):: &
     Rm_A(0:nsmax,0:nsmax,maxstr,maxstr),       & 
     Tm_A(0:nsmax,0:nsmax,maxstr,maxstr),       &
     Rm_AB(0:nsmax,0:nsmax,maxstr,maxstr),      &
     Tm_AB(0:nsmax,0:nsmax,maxstr,maxstr),      &
     RmStar_AB(0:nsmax,0:nsmax,maxstr,maxstr),  &
     TmStar_AB(0:nsmax,0:nsmax,maxstr,maxstr),  &
     Rm_atm(0:nsmax,0:nsmax,maxstr,maxstr),     &
     Tm_atm(0:nsmax,0:nsmax,maxstr,maxstr),     &
     RmStar_atm(0:nsmax,0:nsmax,maxstr,maxstr), &
     TmStar_atm(0:nsmax,0:nsmax,maxstr,maxstr)

   real(kind=dp):: & 
     Rm_sys(0:nsmax,0:nsmax,maxstr,maxstr),     &
     Tm_sys(0:nsmax,0:nsmax,maxstr,maxstr),     &
     RmStar_sys(0:nsmax,0:nsmax,maxstr,maxstr), &
     TmStar_sys(0:nsmax,0:nsmax,maxstr,maxstr)

! refl single scattering
  real(kind=dp):: &
    reflRay_ss(0:nsmax,nstokes,ngeomax),  &
    reflCab_ss(0:nsmax,nstokes,ngeomax),  &
    reflRam_ss(0:nsmax,nstokes,ngeomax),  &
    reflmRay_ss(0:nsmax,nstokes,ngeomax), &
    reflmCab_ss(0:nsmax,nstokes,ngeomax), &
    reflmRam_ss(0:nsmax,nstokes,ngeomax), &
    Zm_Ray_ss(nstokes,ngeomax), &
    Zm_Cab_ss(nstokes,ngeomax)

! aux
  real(kind=dp):: &
    tau_above(0:nsmax), &
    expuv(0:nsmax,ngeomax), &
    expu0(0:nsmax,ngeomax), &
    taut, &
    explay, &
    factor, &
    one_4pi, &
    uv, &
    u0

 real(kind=dp):: &
    deg2rad, &
    rad2deg, &
    pi

 real(kind=dp):: &
    one_d0m, &
    dfac, & !0.5 * (2-d0m) 
    Bmplus(4,ngeomax), &
    Bmmin(4,ngeomax), &
    EpDelta(4), &
    EmDelta(4)

  real(kind=dp):: &
     R_sys_ns(0:nsmax,ngeomax)

  real(kind=dp),intent(out):: &
     refl(0:nsmax,nstokes,ngeomax)

! Surface reflection ==============================================

! Choose type of surface reflection
!   0 = no surface
!   1 = Lambertian surface reflection
!   2 = BDRF

  integer,intent(in):: &
     surf_type 

! Lambertian surface albedo
  real(kind=dp),intent(in):: &
     alb_Lambert

!===========================================================

  pi = 2.d0*asin(1.d0)
  one_4pi = 1.d0/(4.d0*pi)

  rad2deg = 180.d0/pi
  deg2rad = pi/180.d0

!--------------------------------------------------------------------
! Reset
!--------------------------------------------------------------------

  do ns=0,nsmax
     do iv=1,ngeomax
        R_sys_ns(ns,iv) = 0.d0
     enddo
  enddo

!----------------------------------------------------------------------
! Set these to zero. 
! (Have been filled above 
!  in subroutines ss_elastic and ss_inelastic)     
!---------------------------------------------------------------------

  do iv=1,ngeomax
     do i_st = 1,nstokes
        do ns=0,nsmax
           refl(ns,i_st,iv)       = 0.d0
           reflRay_ss(ns,i_st,iv) = 0.d0
           reflCab_ss(ns,i_st,iv) = 0.d0
           reflRam_ss(ns,i_st,iv) = 0.d0
        enddo
     enddo
  enddo

!********************************************************************
! Loop over Fourier moments
!********************************************************************

  do m=0,mmax

!--------------------------------------------------------------------
!    Setup quantities needed for Fourier expansion 
!--------------------------------------------------------------------

     call Fourier_init( &
          m,        ngeomax,  phi,     &
          one_d0m,  dfac,     EpDelta, &
          EmDelta,  Bmplus,   Bmmin)

!********************************************************************
!    Loop over layers with index ik 
!    from bottom to top (reverse order) 
!
!    Optical properties are stored from top to bottom
!    corresponding layer index k = maxlay+1-ik
!********************************************************************

!    reverse order

     do ik=1,maxlay !reverse order

        k = maxlay+1-ik !corresponding k

!--------------------------------------------------------------------
!       Save phase matrix coefficients
!--------------------------------------------------------------------

        do lg=0,maxleg
           wl_Ray_air(lg) = pRay_air(1,1,lg,k)
           wl_Cab_air(lg) = pCab_air(1,1,lg,k)
           wl_Ram_air(lg) = pRam_air(1,1,lg,k)
        enddo

!-------------------------------------------------------------------
!       Determinine number of doublings, 
!       the number of sublayers and their 
!       optical thicknesses for each layer ik, 
!       and for each Fourier component m,
!       and for each incident wavenumber n
!-------------------------------------------------------------------
      
        call nr_doublings( &
             accuracy,   ns_elas,  nsmax,      &    
             k,          maxlay,   m,          &
             maxleg,     gmumin,   wl_Ray_air, &  
             omega_Ray,  tau,      ndl)

!--------------------------------------------------------------------
!       Setup phasefunction moments Zm
!--------------------------------------------------------------------

        call Zm_setup( &
             nsmax,                        &
             k,             maxlay,            m,             &
             mmax,          maxstr,            maxleg,        &
             gmu,           weightCabram_air,  wl_Ray_air,    &
             wl_Cab_air,    wl_Ram_air,        Zm_Ray_min,    &
             Zm_Ray_plus,   Zm_Cab_min,        Zm_Cab_plus,   &
             Zm_Ram_min,    Zm_Ram_plus,       Zm_CabRam_min, &
             Zm_CabRam_plus)

!--------------------------------------------------------------------
!       Setup optical depth dtau for all wavenumbers around
!       wavenumber of incoming light
!--------------------------------------------------------------------

!       set for every layer ik, every m, and every wavenr ns 
        do ns=0,nsmax
           dtau(ns) = tau(k,ns)/(2.d0**dble(ndl)) 
           tau_atm(ns) = 0.d0
        enddo

!********************************************************************
!       Single scattering solution of layer
!********************************************************************

!--------------------------------------------------------------------
!       ss, update
!--------------------------------------------------------------------

        call ss_Zm_setup( &
             m,         k,          nstokes,         &
             ngeomax,   maxstr,     maxintstr,       &
             maxleg,    maxlay,     maxcoeffstrunc,  &
             gmu,       imuv,       imu0,            &
             plm_0,     plm_m,      plm_p,           &
             u,         pRay_air,   Zm_Ray_ss)

        call ss_Zm_setup( &
             m,         k,          nstokes,         &
             ngeomax,   maxstr,     maxintstr,       &
             maxleg,    maxlay,     maxcoeffstrunc,  &
             gmu,       imuv,       imu0,            &
             plm_0,     plm_m,      plm_p,           &
             u,         pCab_air,   Zm_Cab_ss)

        do ns=0,nsmax
           tau_above(ns) = 0.d0
           if(ik.lt.maxlay)then
              do i=ik+1,maxlay
                 tau_above(ns) = tau_above(ns) + tau(maxlay+1-i,ns)
              enddo
           endif
        enddo

        do iv=1,ngeomax

           u0 = gmu(imu0(iv))
           uv = gmu(imuv(iv))

           do ns=0,nsmax

              expu0(ns,iv)  = exp(-tau_above(ns)/u0)
              expuv(ns,iv)  = exp(-tau_above(ns)/uv)
 
              taut   = tau(k,ns)
              explay = exp(-taut * (1.d0/uv + 1.d0/u0))
              factor = one_4pi *  u0/(u0+uv) * &
                          expu0(ns,iv)*expuv(ns,iv) * &
                           (1.d0-explay)

              do i_st=1,nstokes
                 reflmRay_ss(ns,i_st,iv) = &
                        reflmRay_ss(ns,i_st,iv) + &
                         omega_Ray(k,ns) * &
                          Zm_Ray_ss(i_st,iv) * &
                           factor
              enddo 

           enddo

        enddo

!---------------------------------------------------------------------
!
!---------------------------------------------------------------------

        if(single_scattering)then

           call whole_sublayer( &
                ngeomax,        imuv,            imu0,         &
                maxstr,         maxlay,          nsmax,        &
                k,            &
                om_air,         scat_type,       tau,          &
                Zm_Ray_min,     Zm_Ray_plus,     Zm_Cab_min,   &
                Zm_Cab_plus,    Zm_Ram_min,      Zm_Ram_plus,  &
                Zm_CabRam_min,  Zm_CabRam_plus,  gmu,          &
                R1m_wl,         T1m_wl)

!--------------------------------------------------------------------
!          Adding layer A on top of layers below
!--------------------------------------------------------------------

!          Skip adding for first layer

           if(ik.eq.1)then

              do ns=0,nsmax
                 tau_atm(ns) = tau(k,ns)
              enddo

              do iv=1,ngeomax
                 do ns2=0,nsmax
                    do ns1=0,nsmax
                       R1mStar_wl(ns1,ns2,iv)  = R1m_wl(ns1,ns2,iv)
                       R1mStar_wl(ns1,ns2,iv)  = T1m_wl(ns1,ns2,iv)
                       R1m_atm(ns1,ns2,iv)     = R1m_wl(ns1,ns2,iv)
                       T1m_atm(ns1,ns2,iv)     = T1m_wl(ns1,ns2,iv)
                       R1mStar_atm(ns1,ns2,iv) = R1mStar_wl(ns1,ns2,iv)
                       T1mStar_atm(ns1,ns2,iv) = T1mStar_wl(ns1,ns2,iv)
                    enddo
                 enddo
              enddo

           else

              do ns=0,nsmax
                 tau_wl(ns)  = tau(k,ns)   !tau_A  layer ik
                 tau_atm(ns) = tau(k+1,ns) !tau_B  layer ik-1
              enddo

              do iv=1,ngeomax
                 do ns2=0,nsmax
                    do ns1=0,nsmax
                       R1mStar_wl(ns1,ns2,iv) = R1m_wl(ns1,ns2,iv)
                       T1mStar_wl(ns1,ns2,iv) = T1m_wl(ns1,ns2,iv)
                    enddo
                 enddo
              enddo

              call adder1(  &
                   ngeomax,     maxstr,       imuv,         &
                   imu0,        mmax,         nsmax,        &
                   ns_elas,     m,            gmu,          &
                   R1m_wl,      T1m_wl,       R1mStar_wl,   &
                   T1mStar_wl,  tau_wl,       R1m_atm,      &
                   T1m_atm,     R1mStar_atm,  T1mStar_atm,  &
                   tau_atm)

           endif

        else

!********************************************************************
!          Multiple scattering solution of layer
!********************************************************************

!--------------------------------------------------------------------
!          Setup initial sublayer for doubling
!--------------------------------------------------------------------
 
           call basic_sublayer( &
                maxstr,      maxlay,        nsmax,           &
                ns_elas,     n,               &
                k,           om_air,        scat_type,       &
                dtau,        Zm_Ray_min,    Zm_Ray_plus,     &
                Zm_Cab_min,  Zm_Cab_plus,   Zm_Ram_min,      &
                Zm_Ram_plus, Zm_CabRam_min, Zm_CabRam_plus,  &
                gmu,         R1m_bl,        T1m_bl)

           do i2=1,maxstr
              do i1=1,maxstr
                 do ns2=0,nsmax
                    do ns1=0,nsmax
                       Rm_A(ns1,ns2,i1,i2) = R1m_bl(ns1,ns2,i1,i2)
                       Tm_A(ns1,ns2,i1,i2) = T1m_bl(ns1,ns2,i1,i2)
                    enddo
                 enddo
              enddo
           enddo

           do ns=0,nsmax
              dtau_A(ns) = dtau(ns)
           enddo

!----------------------------------------------------------------
!          Start doubling     
!--------------------------------------------------------------------
 
           do ikk = 1,ndl

                 call doublerray( &
                      maxstr,     maxstr_internal,  mmax,   &
                      nsmax,      ns_elas,          dnu,    &
                      m,          one_d0m,          ik,     &
                      ikk,        gmu,              gwt,    &
                      accuracy,   Rm_A,             Tm_A,   &
                      dtau_A,     Rm_AB,            Tm_AB)   
                   
                 do i2=1,maxstr 
                    do i1=1,maxstr    
                       do ns2=0,nsmax
                          do ns1=0,nsmax
                             Rm_A(ns1,ns2,i1,i2) = Rm_AB(ns1,ns2,i1,i2)
                             Tm_A(ns1,ns2,i1,i2) = Tm_AB(ns1,ns2,i1,i2)
                          enddo
                       enddo
                    enddo
                 enddo

                do ns=0,nsmax
                  dtau_A(ns) = 2.d0 * dtau_A(ns)
                enddo
  
           enddo 

!--------------------------------------------------------------------
!          Adding layer on top of layer(s) below
!--------------------------------------------------------------------

!          Skip adding for first layer

           if(ik.eq.1)then

               do ns=0,nsmax
                  tau_atm(ns) = tau(k,ns)
               enddo

!              Total homogenous atmosphere layer
               do i2=1,maxstr
                  do i1=1,maxstr
                     do ns2=0,nsmax
                        do ns1=0,nsmax
                           Rm_atm(ns1,ns2,i1,i2)     = Rm_AB(ns1,ns2,i1,i2)
                           Tm_atm(ns1,ns2,i1,i2)     = Tm_AB(ns1,ns2,i1,i2)
                           RmStar_atm(ns1,ns2,i1,i2) = Rm_AB(ns1,ns2,i1,i2)
                           TmStar_atm(ns1,ns2,i1,i2) = Tm_AB(ns1,ns2,i1,i2)
                        enddo
                     enddo
                  enddo
               enddo

           else

              do ns=0,nsmax
                 tau_AB(ns)  = tau(k,ns)   !tau_A layer ik
                 tau_atm(ns) = tau(k+1,ns) !tau_B layer ik-1
              enddo

!             star
              do i2=1,maxstr
                 do i1=1,maxstr
                    do ns2=0,nsmax
                       do ns1=0,nsmax
                          RmStar_AB(ns1,ns2,i1,i2) = Rm_AB(ns1,ns2,i1,i2)
                          TmStar_AB(ns1,ns2,i1,i2) = Tm_AB(ns1,ns2,i1,i2)
                       enddo
                    enddo
                 enddo
              enddo

              call adderray( &
                   maxstr,     maxstr_internal, mmax,      &
                   nsmax,      ns_elas,         dnu,       &
                   m,          one_d0m,         gmu,       &
                   gwt,        accuracy,        Rm_AB,     &
                   Tm_AB,      RmStar_AB,       TmStar_AB, &
                   tau_AB,     Rm_atm,          Tm_atm,    &
                   RmStar_atm, TmStar_atm,      tau_atm) 

          endif !first layer or not

        endif !multiple/single scat

     enddo !end loop over layers

!**********************************************************************
!    Adding surface
!**********************************************************************

!    total optical depth atmosphere above surface
     do ns=0,nsmax
        tau_atm(ns) = 0.d0
        do k=1,maxlay
           tau_atm(ns) = tau_atm(ns) + tau(k,ns)
        enddo
     enddo

     call add_surf(    &
          single_scattering,  surf_type,        alb_Lambert,   &
          m,                  mmax,             one_d0m,       &
          maxstr,             maxstr_internal,  nsmax,         &
          ns_elas,            dnu,              ngeomax,       &
          gmu,                gwt,              imuv,          &               
          imu0,               accuracy,         tau_atm,       &
          R1m_atm,            T1m_atm,          R1mStar_atm,   &
          T1mStar_atm,        R1m_sys,          T1m_sys,       &
          D1m_sys,            R1mStar_sys,      T1mStar_sys,   &    
          Rm_atm,             Tm_atm,           RmStar_atm,    &
          TmStar_atm,         Rm_sys,           Tm_sys,        &
          RmStar_sys,         TmStar_sys)

!    ss

     if(m.eq.0)then

        do iv=1,ngeomax

           u0 = gmu(imu0(iv))
           uv = gmu(imuv(iv))

           do i_st=1,nstokes
              do ns=0,nsmax

                 expu0(ns,iv) = exp(-tau_atm(ns)/u0)
                 expuv(ns,iv) = exp(-tau_atm(ns)/uv)

!                add Lambertian
                 if(i_st.eq.1)then
                    reflmRay_ss(ns,i_st,iv) = &
                      reflmRay_ss(ns,i_st,iv) + &
                          alb_Lambert * &
                           expu0(ns,iv) * expuv(ns,iv) * &
                            u0 / pi
                 endif

              enddo
           enddo

        enddo

     endif

!---------------------------------------------------------------------
!    Add Fourier moments
!
!    R_sys_ns = sum { Rm_sys_ns ...  }
!                m
!---------------------------------------------------------------------

     call add_Fourier_moment( &
            m,          nsmax,      ns_elas,           &
            ngeomax,    maxstr,     nstokes,           &
            single_scattering, & 
            imuv,       imu0,       phi,               &
            deg2rad,    pi,         dfac,              &
            Bmplus,     Bmmin,      EpDelta,           &
            EmDelta,    R1m_sys,    Rm_sys,            &
            R_sys_ns)

!    ss
     do iv=1,ngeomax
        do i_st=1,nstokes
              do ns=0,nsmax
                 reflRay_ss(ns,i_st,iv) = reflRay_ss(ns,i_st,iv) + &
                                            EpDelta(i_st) * &
                                             dfac * Bmplus(i_st,iv) * &
                                              reflmRay_ss(ns,i_st,iv) + &
                                            EmDelta(i_st) * &
                                             dfac * Bmmin(i_st,iv) * &
                                              reflmRay_ss(ns,i_st,iv)

                 reflCab_ss(ns,i_st,iv) = reflCab_ss(ns,i_st,iv) + &
                                            EpDelta(i_st) * &
                                             dfac * Bmplus(i_st,iv) * &
                                              reflmCab_ss(ns,i_st,iv) + &
                                            EmDelta(i_st) * &
                                             dfac * Bmmin(i_st,iv) * &
                                              reflmCab_ss(ns,i_st,iv)

                 reflRam_ss(ns,i_st,iv) = reflRam_ss(ns,i_st,iv) + &
                                            EpDelta(i_st) * &
                                             dfac * Bmplus(i_st,iv) * &
                                              reflmRam_ss(ns,i_st,iv) + &
                                            EmDelta(i_st) * &
                                             dfac * Bmmin(i_st,iv) * &
                                              reflmRam_ss(ns,i_st,iv)

!                reflectivity
                 reflRay_ss(ns,i_st,iv) = reflRay_ss(ns,i_st,iv) * pi/u0
                 reflCab_ss(ns,i_st,iv) = reflCab_ss(ns,i_st,iv) * pi/u0
                 reflRam_ss(ns,i_st,iv) = reflRam_ss(ns,i_st,iv) * pi/u0

              enddo
           enddo
     enddo

  enddo !end loop over Fourier moments

!-----------------------------------------------------------------------
! for testing only
!-----------------------------------------------------------------------

  do iv=1,ngeomax
     do i_st=1,nstokes
        if(i_st.eq.1)then        
           do ns=0,nsmax
              refl(ns,i_st,iv) = R_sys_ns(ns,iv)
           enddo
        else
           do ns=0,nsmax
              refl(ns,i_st,iv) = 0.d0
           enddo
        endif
     enddo
  enddo

  return
  end subroutine dasray_intf
