
  subroutine binned_cross_section_Raman( &
    maxlay,           temp,             jmax,            &
    nsmax,            ns_elas,          wavelength_ns,   &
    wavenr_ns,        wavenrshift_ns,   dnu,             &
    relN2,            relO2,            relAr,           &
    csray_N2,         csray_O2,         csray_Ar,        &
    eps_N2,           eps_O2,           linestr_Sto_N2,  &
    linestr_aSt_N2,   linestr_Sto_O2,   linestr_aSt_O2,  &
    dwn_Sto_N2,       dwn_aSt_N2,       dwn_Sto_O2,      &
    dwn_aSt_O2,       csram_Sto_N2,     csram_aSt_N2,    &
    csram_Sto_O2,     csram_aSt_O2,     csram_N2,        &
    csram_O2,         csram_air,        scat_type,       &
    weightCabRam_air)

!-----------------------------------------------------------
! Subroutine:
!    csrambin
!
! Requires:
!    binner
!    shift_column
!
!
! Purpose:
!
! Count count_ns only for ns2=ns2check and k=kcheck
!------------------------------------------------------------------

  implicit none

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

  integer,intent(in):: &
    jmax,   &  !maximum rotational qunatum number J
    maxlay, &  !nr of layers
    nsmax,  &     !nr of wavelengths
    ns_elas  

  real(kind=dp),intent(in):: &
    eps_N2(0:nsmax),   &
    eps_O2(0:nsmax),   &
    csray_N2(0:nsmax), &
    csray_O2(0:nsmax), &
    csray_Ar(0:nsmax)

! wavenrshifts
  real(kind=dp),intent(in):: &
    dwn_Sto_N2(0:jmax), &
    dwn_aSt_N2(0:jmax), &
    dwn_Sto_O2(0:jmax), &
    dwn_aSt_O2(0:jmax)

  real(kind=dp),dimension(maxlay,0:jmax,0:nsmax),intent(in):: &   
    csram_Sto_N2, &
    csram_aSt_N2, &  
    csram_Sto_O2, &  
    csram_aSt_O2

! wavelengths around incident light (ns2_elas) [nm]
  real(kind=dp),intent(in):: &
    wavelength_ns(0:nsmax),  &
    wavenr_ns(0:nsmax),      & 
    wavenrshift_ns(0:nsmax), &
    dnu

  real(kind=dp),intent(in):: & 
     temp(maxlay), & !temperature [K] for each layer
     relN2, & !N2 volume mixing ratio
     relO2, & !O2 volume mixing ratio
     relAr

  real(kind=dp),dimension(0:jmax,maxlay),intent(in):: &
    linestr_Sto_N2, &
    linestr_aSt_N2, &
    linestr_Sto_O2, &
    linestr_aSt_O2 

  integer::ns1,ns2,k,j,ns

  real(kind=dp):: &
    testfac

  integer:: &
    kcheck,  &
    ns2check, &
    ns1r,     &
    count_ns(0:nsmax)

  real(kind=dp),dimension(maxlay,0:nsmax,0:nsmax),intent(out):: &
    csram_N2, &
    csram_O2, &
    csram_air    !Raman cross section for air
    
  integer,intent(out):: &
    scat_type(0:nsmax,0:nsmax) !0=nothing; 1=Cab; 2=Ram; 3=Cab+ram

  real(kind=dp),intent(out):: &
    weightCabRam_air(maxlay,0:nsmax) !how many inel. light 
    !added top Cabbannes scattering for each ns2 

  logical,parameter:: &
    verbose = .false.,            &
    energy_conservation = .true.

!-----------------------------------------------------------------------
! Reset Raman cross section
!----------------------------------------------------------------------

  do ns2=0,nsmax
     do ns1=0,nsmax
        do k=1,maxlay
           csram_N2(k,ns1,ns2)  = 0.d0
           csram_O2(k,ns1,ns2)  = 0.d0
           csram_air(k,ns1,ns2) = 0.d0
        enddo
        scat_type(ns1,ns2) = 0
     enddo
  enddo

!-----------------------------------------------------------------------
! Bin Raman lines and count (non-zero) Raman lines per bin
!
! incoming wavelength ns_elas
!-----------------------------------------------------------------------

! select layer to count Raman lines;
! count_ns is determined for this layer
  kcheck = maxlay

! count_ns is determined for wavenumber index ns2check
  ns2check = 0

! reset counter
  do ns1=0,nsmax
     count_ns(ns1) = 0
  enddo

! loop over incident wavenr
  do ns2=0,nsmax

!    reset column
     do ns1=0,nsmax
        do k=1,maxlay
           csram_N2(k,ns1,ns2)  = 0.d0
           csram_O2(k,ns1,ns2)  = 0.d0
           csram_air(k,ns1,ns2) = 0.d0
         enddo
     enddo

!-------------------------------------------------------------
!    bin Raman lines per column
!-------------------------------------------------------------

!    Scattering type
     scat_type(0:nsmax,ns2) = 2 !Raman scattering
     scat_type(ns_elas,ns2) = 1 !Cabannes scattering

!    Initial value, pure Cabannes scattering for elements 
!    (ns_elas,ns2)
     weightCabRam_air(1:maxlay,ns2) = 1.d0
     
     call bin_Raman_lines_per_column(  &
          ns2,            ns2check,      kcheck,         &
          jmax,           maxlay,        nsmax,          &
          ns_elas,        wavelength_ns, wavenr_ns,      &
          wavenrshift_ns, dnu,           relN2,          &
          relO2,          dwn_Sto_N2,    dwn_aSt_N2,     &
          dwn_Sto_O2,     dwn_aSt_O2,    csram_Sto_N2,   &
          csram_aSt_N2,   csram_Sto_O2,  csram_aSt_O2,   &
          csram_N2,       csram_O2,      csram_air,      &
          count_ns)

!--------------------------------------------------------------
!    Special test case
!--------------------------------------------------------------

     goto 444

     if(ns2.eq.ns_elas)then
        write(*,*) 'csram = 0.0 0.025 0.05 0.800 0.05 0.025 0.00'
     endif

     do k=1,maxlay

        csram_N2(k,0:nsmax,ns2)  = 0.d0
        csram_O2(k,0:nsmax,ns2)  = 0.d0
        csram_air(k,0:nsmax,ns2) = 0.d0

        do ns=-2,2
           if(abs(ns).eq.0) testfac = 0.800d0
           if(abs(ns).eq.1) testfac = 0.050d0
           if(abs(ns).eq.2) testfac = 0.025d0
           csram_N2(k,ns_elas+ns,ns2) = testfac*csray_N2(ns2)
           csram_O2(k,ns_elas+ns,ns2) = testfac*csray_O2(ns2)
           csram_air(k,ns_elas+ns,ns2) = &
              relN2*csram_N2(k,ns_elas+ns,ns2) + &
                relO2*csram_O2(k,ns_elas+ns,ns2)
        enddo
        
     enddo

 444 continue  

     if(verbose)then
        do ns1=0,nsmax
           write(*,'(F8.2,I5,3D16.8,I4,f8.3)') &
             wavenrshift_ns(ns1),         &
             count_ns(ns1),               &
             csram_N2(kcheck,ns1,ns2),    &
             csram_O2(kcheck,ns1,ns2),    &
             csram_air(kcheck,ns1,ns2),   &
             scat_type(ns1,ns2),          &
             weightCabRam_air(kcheck,ns2)
        enddo
        write(*,*) '---'
      endif

!---------------------------------------------------------------
!    Energy conservation
!    Make sure that sum(csram)+cscab = csray
!    (might be obsolete)
!---------------------------------------------------------------

     if(energy_conservation)then

        call energy_conservation_column( &
          ns2,       ns_elas,   nsmax,        &    
          maxlay,    relN2,     relO2,        &
          eps_N2,    eps_O2,    csray_N2,     &  
          csray_O2,  csram_N2,  csram_O2,     &  
          csram_air, scat_type,  weightCabRam_air)

      else
 
         write(*,*) 'csrambin.f: skip energy conservation'

      endif

!-------------------------------------------------------------
!    shift elements in column so that
!        1 1 1 1 1           3 2 1 0 0
!        2 2 2 2 2           4 3 2 1 0
!        3 3 3 3 3  becomes  5 4 3 2 1
!        4 4 4 4 4           0 5 4 3 2
!        5 5 5 5 5           0 0 5 4 3 
!    and determine cabram weights (correct for border-effects)
!-------------------------------------------------------------

     call shift_elements_for_each_column( &
          ns2,       ns_elas,      nsmax,        &    
          maxlay,    relN2,        relO2,        &
          relAr,     eps_N2,       eps_O2,       &
          csray_N2,  csray_O2,     csray_Ar,     &
          csram_N2,  csram_O2,     csram_air,    &
          scat_type,  weightCabRam_air)
     
     if(verbose)then
        do ns1=0,nsmax
           write(*,'(F8.2,I5,3D16.8,I4,F8.3)') &
             wavenrshift_ns(ns1),         &
             count_ns(ns1),               &
             csram_N2(kcheck,ns1,ns2),    &
             csram_O2(kcheck,ns1,ns2),    &
             csram_air(kcheck,ns1,ns2),   &
             scat_type(ns1,ns2),          &
             weightCabRam_air(kcheck,ns2)
        enddo
        write(*,*) '===================================---'
    endif

  enddo

  return
  end subroutine binned_cross_section_Raman
