
  subroutine clearsky_optical_properties( &
       wavelength_ns, nsmax,      ns_elas,    &
       maxlay,        nstokes,    maxleg,     &
       relN2,         relO2,      relAr,      &
       csabs_O3,      csabs_O4,   csabs_NO2,  &    
       csRay_N2,      csRay_O2,   csRay_Ar,   &
       csram_N2,      csram_O2,   csram_air,  &
       eps_N2,        eps_O2,     eps_Ar,     &
       dvO4,          dvO3,       dvNO2,      &
       dvair,         dz,         tautot,     &
       omega_Ray,     omega_Cab,  om_air,     &
       pRay_N2,       pCab_N2,    pRay_O2,    &
       pCab_O2,       pRay_Ar,    pCab_Ar,    &
       pRay_air,      pCab_air,   pRam_air,   &
       maxcoeffs)

!-------------------------------------------------------------------
!
! Subroutine:
!   clearsky_optical_properties
!
! Purpose:
!   Return optical properties of clear sky atmosphere 
!   (optical thickness; single scattering albedo; phase functions; 
!    expansion coefficient matrices)
!!
! Assumption: phase matrix does not change over wavelength range
! wavelength_ns(0 ... nsmax)
!
! omega = (csram+cscab) / csray
!
! For phase matrix and phase function see Table A.1 in 
! appendix A of "Rotational Raman scattering  of polarized light in   
! the Earth's atmosphere..." - Landgraf et al., jqsrt 2004
!-------------------------------------------------------------------
!
! The Rayleigh, Cabannes and Raman phase matrix P has the shape:          
!
!        ( P11  P12   0    0  )
!   P =  ( P21  P22   0    0  )   where P11 = phase function
!        (  0    0   P33  P34 )
!        (  0    0   P43  P44 )
!
!        P11 = [ A [1 + (costheta)**2] + C] / N
!        P12 =  -A [1 - (costheta)**2] + C] / N 
!        P22 =   A [1 + (costheta)**2] / N                                    
!        P33 =  2A [costheta] / N   
!        P44 =   B [costheta] / N  
!        P21 =  P12 
!        P43 =  -P34 = 0   
!
!  Rayleigh:
!
!    A = 3(45 + eps)            
!    B = 30(9 - eps)  
!    C = 36 eps     
!    N = 180 + 40 eps 
!
!  Cabannes:   
!  
!    A = 3(180 + eps) 
!    B = 30(36 - eps) 
!    C = 36 eps      
!    N = 40(18 + eps) 
!
!  Raman:  
!    A =   3 
!    B =  -30 (erratum Landgraf et al. -10 -> -30)       
!    C =  36  
!    N =  40  
!
!  where eps is the molecular anistropy, 
!  related to the king factor (Bates 1984) as Fking = 1 + (2/9)eps 
! 
!  The phase matrix is expanded in generalized spherical function as 
!  follows:
!                 2
!    P11       = sum [alpha1(ns)] * p(l,0,0)
!                l=0  
!
!                 2
!    P22 + P33 = sum [alpha2(ns) + alpha3(ns)] * p(l,2,+2)
!                l=2  
!  
!                 2
!    P22 - P33 = sum [alpha2(ns) - alpha3(ns)] * p(l,2,-2)
!                l=2  
!
!                 2
!    P44       = sum [alpha4(ns)] * p(l,0,0)
!                l=0  
!                             
!                 2
!    P12 = P21 = sum [beta1(ns)] * p(l,0,+2)
!                l=2                    
!
!                 2
!    P34       = sum [beta2(ns)] * p(l,0,+2)
!                l=2   
!                      
!    P43       = -P34
!
!  where generalized spherical functions p(l,m,n) are given by
!  [de Haan et al 1987, A&A 183, 371-391]:
!
!    l    p(l,m,n)
!   ---------------------------------------------------------
!    0    p(0,0,0)  = 1
!    1    p(1,0,0)  = costheta
!    2    p(2,0,0)  = 0.50*[3(costheta)**2 - 1]
!    0    p(0,0,+2) = p(0,0,-2) =  0
!    1    p(1,0,+2) = p(1,0,-2) =  0
!    2    p(2,0,+2) = p(2,0,-2) = -0.25*sqrt(6)*[1 - (costheta)**2]
!    0    p(0,2,+2) = p(0,2,-2) = 0 
!    1    p(1,2,+2) = p(0,2,-2) = 0
!    2    p(2,2,+2) = 0.25*[1 + 2costheta + costheta**2] 
!    2    p(2,2,-2) = 0.25*[1 - 2costheta + costheta**2] 
!
!  see also subroutine "gsf.f90"
!            
!  Expansion coefficient matrix for Rayleigh, Cabannes and Raman 
!  scattering is returned as phaseRay, phaseCab and phaseRam :
!
!    ( alpha1  beta1     0       0     )
!    ( beta1   alpha2    0       0     )
!    (   0       0     alpha3   beta2  )
!    (   0       0     -beta2   alpha4 )
! 
!-------------------------------------------------------------------

  implicit none  

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

  integer,intent(in):: &
    nsmax,    & !nr of wavenumbers (incident light)
    ns_elas,  & !wavelength_ns elastic scattering
    maxlay,   &
    nstokes,  &          
    maxleg    

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

  real(kind=dp),intent(in):: & 
     relN2, & !N2 volume mixing ratio
     relO2, & !O2 volume mixing ratio
     relAr    !Ar volume mixing ratio

  real(kind=dp),intent(in):: &
    dz(maxlay),              & !layer tickness [m]
    csabs_O3(maxlay,0:nsmax),  & !O3 csabs section [cm^2]
    csabs_O4(maxlay,0:nsmax),  & !O4 csabs section [cm^5 molec-2]
    csabs_NO2(maxlay,0:nsmax), & !NO2 csabs section [cm^2]
    csRay_N2(0:nsmax),         & !Rayleigh csabs section [cm2] for N2
    csRay_O2(0:nsmax),         & !Rayleigh csabs section [cm2] for O2
    csRay_Ar(0:nsmax),         & !Rayleigh csabs section [cm2] for Ar
    eps_N2(0:nsmax),           & !epsilon for N2
    eps_O2(0:nsmax),           & !epsilon for O2
    eps_Ar(0:nsmax),           & !epsilon for Ar = 0
    dvO4(maxlay),            & !diff. O4 column density [part^2./cm^5]
    dvO3(maxlay),            & !diff. O3 column density [part./cm^2]
    dvNO2(maxlay),           & !diff. NO2 column density [part./cm^2]
    dvair(maxlay)              !diff. O2 column density [part./cm^2]   

  real(kind=dp),intent(in):: &
    csram_N2(maxlay,0:nsmax,0:nsmax), &
    csram_O2(maxlay,0:nsmax,0:nsmax), &
    csram_air(maxlay,0:nsmax,0:nsmax) !Cab+Ram

  real(kind=dp)::  &
    elasfac_N2,       & !elastic fraction
    elasfac_O2,       &
    elasfac_Ar,       &
    weight_N2,     &
    weight_O2,     &
    weight_Ar,     &
    tausca_Ray_N2, &
    tausca_Ray_O2, &
    tausca_Ray_Ar

! phase matrix for vector radiative transfer, 
! elements are divided by 2m+1
  real(kind=dp):: &
    A_Ray_N2, & !Rayleigh
    A_Ray_O2, &
    A_Ray_Ar, &     
    B_Ray_N2, &
    B_Ray_O2, &
    B_Ray_Ar, &       
    C_Ray_N2, &
    C_Ray_O2, &
    C_Ray_Ar, &    
    N_Ray_N2, &
    N_Ray_O2, &
    N_Ray_Ar, &    
    A_Cab_N2, & !Cabannes     
    A_Cab_O2, &
    A_Cab_Ar, &   
    B_Cab_N2, &
    B_Cab_O2, &
    B_Cab_Ar, &    
    C_Cab_N2, &
    C_Cab_O2, &
    C_Cab_Ar, &        
    N_Cab_N2, &
    N_Cab_O2, &
    N_Cab_Ar, &
    A_Ram,    & !Raman
    B_Ram,    &
    C_Ram,    &
    N_Ram

  real(kind=dp):: &
    phaseRay_air(4,4,0:2), & !Rayleigh phase fct. moments
    phaseRay_N2(4,4,0:2),  & !Rayleigh phase fct. moments
    phaseRay_O2(4,4,0:2),  & !Rayleigh phase fct. moments
    phaseRay_Ar(4,4,0:2),  & !Rayleigh phase fct. moments
    phaseCab_air(4,4,0:2), & !Cabannes phase fct. moments
    phaseCab_N2(4,4,0:2),  & !Cabannes phase fct. moments
    phaseCab_O2(4,4,0:2),  & !Cabannes phase fct. moments
    phaseCab_Ar(4,4,0:2),  & !Cabannes phase fct. moments
    phaseRam(4,4,0:2)        !Raman phase fct. moments

  integer::ns,k,m,lg,i1,i2,i_st,j_st

  real(kind=dp),dimension(maxlay,0:nsmax):: &
    tauabs, & !total abs. optical depth
    tausca,&    !total scat. optical depth
    bscat_Ray, & !beta_scat [m-1]
    bscat_Cab  !beta_scat [m-1]

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

  real(kind=dp),intent(out):: &
    om_air(maxlay,0:nsmax,0:nsmax) !watch out, different for each layer k

 real(kind=dp), &
  dimension(nstokes,nstokes,0:maxleg,maxlay),intent(out):: &
    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
 
 integer::maxcoeffs(maxlay) !number of phase fct. coeff. 

 logical,parameter::scalar =.true.

 logical,parameter:: &
   verbose_input  = .false., & !show setup expansion coefficients
   verbose_output = .false.    !show selected expans.coefficients

 do ns=0,nsmax

!-------------------------------------------------------------------
!    Optical depths and single scattering albedos
!-------------------------------------------------------------------

     elasfac_N2 = (18.d0+eps_N2(ns))/(18.d0+4.d0*eps_N2(ns))
     elasfac_O2 = (18.d0+eps_O2(ns))/(18.d0+4.d0*eps_O2(ns))
     elasfac_Ar =  1.d0
         
!    Total optical depth (Rayleigh + absorption) and
!    single scattering albedo for Cabannes and Rayleigh scattering

     do k = 1,maxlay
            
        tausca_Ray_N2 = csRay_N2(ns) * relN2 * dvair(k)
        tausca_Ray_O2 = csRay_O2(ns) * relO2 * dvair(k) 
        tausca_Ray_Ar = csRay_Ar(ns) * relAr * dvair(k)          


        tausca(k,ns) = tausca_Ray_N2 + &
                       tausca_Ray_O2 + &
                       tausca_Ray_Ar

        tauabs(k,ns) = csabs_O3(k,ns)  * dvO3(k)  + &
                       csabs_O4(k,ns)  * dvO4(k)  + &
                       csabs_NO2(k,ns) * dvNO2(k)
                   
        tautot(k,ns) = tauabs(k,ns) + tausca(k,ns)
        
!        if(ns.eq.ns_elas) &
!           write(*,'(i4,f10.4,4d12.4)') & 
!             k,wavelength_ns(ns),tausca(k,ns_elas),tauabs(k,ns_elas),tautot(k,ns_elas)

        omega_Ray(k,ns)    = 0.d0
        omega_Cab(k,ns)    = 0.d0
        om_air(k,0:nsmax,ns) = 0.d0

        if(tautot(k,ns).ne.0.d0)then

           omega_Ray(k,ns) = tausca(k,ns)/tautot(k,ns)
           omega_Cab(k,ns) =(elasfac_N2*tausca_Ray_N2 + &
                             elasfac_O2*tausca_Ray_O2 + &
                             elasfac_Ar*tausca_Ray_Ar)/tautot(k,ns)

           om_air(k,0:nsmax,ns) = &
             csram_air(k,0:nsmax,ns)*dvair(k)/tautot(k,ns)
        endif
               
        bscat_Ray(k,ns) = omega_Ray(k,ns)*tautot(k,ns)/dz(k) ![m-1]
        bscat_Cab(k,ns) = omega_Cab(k,ns)*tautot(k,ns)/dz(k) ![m-1]

     enddo

!--------------------------------------------------------------------
!    Setup Rayleigh phase matrix, 
!    assume independent of ns (calculate for ns_elas)
!-------------------------------------------------------------------

     if(ns.eq.ns_elas)then

!       Save omega, tautot and pRay for selected wavelength       
!       weight = omega_Ray_species(k)/omega_Ray_air(k,ns), 
!       which does not depend on altitude
      
        weight_N2 = csRay_N2(ns)*relN2/(csRay_N2(ns)*relN2 + &
                                        csRay_O2(ns)*relO2 + &
                                        csRay_Ar(ns)*relAr)
        weight_O2 = csRay_O2(ns)*relO2/(csRay_N2(ns)*relN2 + &
                                        csRay_O2(ns)*relO2 + &
                                        csRay_Ar(ns)*relAr)
        weight_Ar = csRay_Ar(ns)*relAr/(csRay_N2(ns)*relN2 + &
                                        csRay_O2(ns)*relO2 + &
                                        csRay_Ar(ns)*relAr)
     
!       parameters A,B,C and normalization constant N in 
!       Table A.1 landgraf et al.
               
        A_Ray_N2 = 3.d0 * (45.d0 + eps_N2(ns))
        A_Ray_O2 = 3.d0 * (45.d0 + eps_O2(ns))
        A_Ray_Ar = 3.d0 *  45.d0 
            
        B_Ray_N2 = 30.d0 * (9.d0 - eps_N2(ns))
        B_Ray_O2 = 30.d0 * (9.d0 - eps_O2(ns))
        B_Ray_Ar = 30.d0 *  9.d0
            
        C_Ray_N2 = 36.d0 * eps_N2(ns)
        C_Ray_O2 = 36.d0 * eps_O2(ns)
        C_Ray_Ar =  0.d0
            
        N_Ray_N2 = 180.d0 + 40.d0 * eps_N2(ns)
        N_Ray_O2 = 180.d0 + 40.d0 * eps_O2(ns)
        N_Ray_Ar = 180.d0 + 40.d0 * eps_Ar(ns)
            
!       Rayleigh phase matrix elements 
!       in generalized spherical function expansion
      
        call RayCaB_Ram_EXPANSION_coeff( &
             A_Ray_N2, B_Ray_N2, C_Ray_N2, N_Ray_N2, phaseRay_N2)
            
        call RayCaB_Ram_EXPANSION_coeff( &
             A_Ray_O2, B_Ray_O2, C_Ray_O2, N_Ray_O2, phaseRay_O2)
            
        call RayCaB_Ram_EXPANSION_coeff( &
             A_Ray_Ar, B_Ray_Ar, C_Ray_Ar, N_Ray_Ar, phaseRay_Ar)
            
        do m=0,2
           do i2=1,4
              do i1=1,4
                 phaseRay_air(i1,i2,m) = &
                      weight_N2 * phaseRay_N2(i1,i2,m) + &
                      weight_O2 * phaseRay_O2(i1,i2,m) + &
                      weight_Ar * phaseRay_Ar(i1,i2,m)
              enddo
           enddo
        enddo

     endif !ns_elas

!--------------------------------------------------------------------
!    Setup Cabannes phase matrix, assume independent of ns
!    (Calculate for ns_elas)
!--------------------------------------------------------------------
      
     if(ns.eq.ns_elas)then

!       weight = omega_Cab_species(k)/omega_Cab(k,ns) 
!       which does not depend on altitude
      
        weight_N2 = elasfac_N2 * csRay_N2(ns) * relN2 / &
                   (elasfac_N2 * csRay_N2(ns) * relN2 + &
                    elasfac_O2 * csRay_O2(ns) * relO2 + &
                    elasfac_Ar * csRay_Ar(ns) * relAr)
        weight_O2 = elasfac_O2 * csRay_O2(ns) * relO2 / &
                   (elasfac_N2 * csRay_N2(ns) * relN2 + &
                    elasfac_O2 * csRay_O2(ns) * relO2 + &
                    elasfac_Ar * csRay_Ar(ns) * relAr)
        weight_Ar = elasfac_Ar * csRay_Ar(ns) * relAr / & 
                   (elasfac_N2 * csRay_N2(ns) * relN2 + &
                    elasfac_O2 * csRay_O2(ns) * relO2 + &
                    elasfac_Ar * csRay_Ar(ns) * relAr)
            
!       parameters A,B,C and normalization constant N in 
!       Table A.1 landgraf et al.

        A_Cab_N2 = 3.d0 * (180.d0 + eps_N2(ns))
        A_Cab_O2 = 3.d0 * (180.d0 + eps_O2(ns))
        A_Cab_Ar = 3.d0 *  180.d0 
            
        B_Cab_N2 = 30.d0 * (36.d0 - eps_N2(ns))
        B_Cab_O2 = 30.d0 * (36.d0 - eps_O2(ns))
        B_Cab_Ar = 30.d0 *  36.d0
            
        C_Cab_N2 = 36.d0 * eps_N2(ns)
        C_Cab_O2 = 36.d0 * eps_O2(ns)
        C_Cab_Ar =  0.d0
            
        N_Cab_N2 = 40.d0 * (18.d0 + eps_N2(ns))
        N_Cab_O2 = 40.d0 * (18.d0 + eps_O2(ns))
        N_Cab_Ar = 40.d0 *  18.d0
            
!       Cabannes phase matrix elements 
!       in generalized spherical function expansion
      
        call RayCaB_Ram_EXPANSION_coeff( &
             A_Cab_N2, B_Cab_N2, C_Cab_N2, N_Cab_N2, phaseCab_N2)
            
        call RayCaB_Ram_EXPANSION_coeff( &
             A_Cab_O2, B_Cab_O2, C_Cab_O2, N_Cab_O2, phaseCab_O2)
            
        call RayCaB_Ram_EXPANSION_coeff( &
             A_Cab_Ar, B_Cab_Ar, C_Cab_Ar, N_Cab_Ar, phaseCab_Ar)
            
        do m=0,2
           do i2=1,4
              do i1=1,4
                 phaseCab_air(i1,i2,m) = &
                      weight_N2 * phaseCab_N2(i1,i2,m) + &
                      weight_O2 * phaseCab_O2(i1,i2,m) + &
                      weight_Ar * phaseCab_Ar(i1,i2,m)
              enddo
           enddo
        enddo

     endif !ns_elas only

  enddo !ns
            
!--------------------------------------------------------------------
! Raman phase matrix, phase function 
! independent of wavelength and species
!--------------------------------------------------------------------

! parameters A,B,C and normalization constant N in 
! Table A.1 landgraf et al.

  A_Ram =   3.d0 
  B_Ram = -30.d0 !correction, used to be -10.d0 in paper Landgraf 2004
  C_Ram =  36.d0 
  N_Ram =  40.d0

  call RayCaB_Ram_EXPANSION_coeff( &
       A_Ram, B_Ram, C_Ram, N_Ram, phaseRam)

  if(verbose_input)then
     write(*,*) 'Setup Expansion coefficients for Rayleigh scattering'
     write(*,'(A4,6A9)') &
       'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
     do lg=0,2
        write(*,'(I4,6F9.5)')   &
          lg,                   &
          phaseRay_air(1,1,lg), &
          phaseRay_air(2,2,lg), &
          phaseRay_air(3,3,lg), &
          phaseRay_air(4,4,lg), &
          phaseRay_air(1,2,lg), &
          phaseRay_air(3,4,lg)
      enddo  
      write(*,*) &
        'Setup Expansion coefficients for Cabannes scattering'
      write(*,'(A4,6A9)') &
        'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
      do lg=0,2
         write(*,'(I4,6F9.5)')   &
           lg,                   &
           phaseCab_air(1,1,lg), &
           phaseCab_air(2,2,lg), &
           phaseCab_air(3,3,lg), &
           phaseCab_air(4,4,lg), &
           phaseCab_air(1,2,lg), &
           phaseCab_air(3,4,lg)
       enddo        
       write(*,*) & 
         'Setup Expansion coefficients for Raman scattering'
       write(*,'(A4,6A9)') & 
         'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
       do lg=0,2 
          write(*,'(I4,6F9.5)')  &
            lg,                  &
            phaseRam(1,1,lg),    &
            phaseRam(2,2,lg),    &
            phaseRam(3,3,lg),    &
            phaseRam(4,4,lg),    &
            phaseRam(1,2,lg),    &
            phaseRam(3,4,lg)
       enddo
  endif

!--------------------------------------------------------------------
! Select phase matrix elements 
! for (4,4,lgmax) -> (nstokes,nstokes,maxleg) 
!--------------------------------------------------------------------

  do k = 1,maxlay
     maxcoeffs(k) = 2
  enddo
      
  do k=1,maxlay
   do lg = 0,maxleg   
    do i_st = 1,nstokes
     do j_st = 1,nstokes
      if(lg.lE.2)then
         pRay_N2( j_st,i_st,lg,k) = phaseRay_N2( j_st,i_st,lg)
         pCab_N2( j_st,i_st,lg,k) = phaseCab_N2( j_st,i_st,lg)
         pRay_O2( j_st,i_st,lg,k) = phaseRay_O2( j_st,i_st,lg)
         pCab_O2( j_st,i_st,lg,k) = phaseCab_O2( j_st,i_st,lg)
         pRay_Ar( j_st,i_st,lg,k) = phaseRay_Ar( j_st,i_st,lg)
         pCab_Ar( j_st,i_st,lg,k) = phaseCab_Ar( j_st,i_st,lg)
         pRay_air(j_st,i_st,lg,k) = phaseRay_air(j_st,i_st,lg)
         pCab_air(j_st,i_st,lg,k) = phaseCab_air(j_st,i_st,lg)
         pRam_air(j_st,i_st,lg,k) = phaseRam(j_st,i_st,lg)
         if(scalar)then
            pRay_N2( j_st,i_st,lg,k) = 0.d0
            pCab_N2( j_st,i_st,lg,k) = 0.d0
            pRay_O2( j_st,i_st,lg,k) = 0.d0
            pCab_O2( j_st,i_st,lg,k) = 0.d0
            pRay_Ar( j_st,i_st,lg,k) = 0.d0
            pCab_Ar( j_st,i_st,lg,k) = 0.d0
            pRay_air(j_st,i_st,lg,k) = 0.d0
            pCab_air(j_st,i_st,lg,k) = 0.d0
            pRam_air(j_st,i_st,lg,k) = 0.d0
            if((j_st.eq.1).AND.(i_st.eq.1))then
               pRay_N2( j_st,i_st,lg,k) = phaseRay_N2( j_st,i_st,lg)
               pCab_N2( j_st,i_st,lg,k) = phaseCab_N2( j_st,i_st,lg)
               pRay_O2( j_st,i_st,lg,k) = phaseRay_O2( j_st,i_st,lg)
               pCab_O2( j_st,i_st,lg,k) = phaseCab_O2( j_st,i_st,lg)
               pRay_Ar( j_st,i_st,lg,k) = phaseRay_Ar( j_st,i_st,lg)
               pCab_Ar( j_st,i_st,lg,k) = phaseCab_Ar( j_st,i_st,lg)
               pRay_air(j_st,i_st,lg,k) = phaseRay_air(j_st,i_st,lg)
               pCab_air(j_st,i_st,lg,k) = phaseCab_air(j_st,i_st,lg)
               pRam_air(j_st,i_st,lg,k) = phaseRam(j_st,i_st,lg)
            endif
         endif
      else
          pRay_N2( j_st,i_st,lg,k) = 0.d0 
          pCab_N2( j_st,i_st,lg,k) = 0.d0
          pRay_O2( j_st,i_st,lg,k) = 0.d0
          pCab_O2( j_st,i_st,lg,k) = 0.d0
          pRay_Ar( j_st,i_st,lg,k) = 0.d0
          pCab_Ar( j_st,i_st,lg,k) = 0.d0
          pRay_air(j_st,i_st,lg,k) = 0.d0
          pCab_air(j_st,i_st,lg,k) = 0.d0
          pRam_air(j_st,i_st,lg,k) = 0.d0
      endif
     enddo
    enddo
   enddo
  enddo

 goto 44
  write(*,*)'Isotropic phasefunctions'
  do lg = 0,maxleg   
   do i_st = 1,nstokes
    do j_st = 1,nstokes
     do k = 1,maxlay  
        if((i_st.eq.1).AND.(j_st.eq.1).AND.(lg.eq.0))then
           pRay_N2( j_st,i_st,lg,k) = 1.d0
           pCab_N2( j_st,i_st,lg,k) = 1.d0
           pRay_O2( j_st,i_st,lg,k) = 1.d0
           pCab_O2( j_st,i_st,lg,k) = 1.d0
           pRay_Ar( j_st,i_st,lg,k) = 1.d0
           pCab_Ar( j_st,i_st,lg,k) = 1.d0
           pRay_air(j_st,i_st,lg,k) = 1.d0
           pCab_air(j_st,i_st,lg,k) = 1.d0
        else
           pRay_N2( j_st,i_st,lg,k) = 0.d0
           pCab_N2( j_st,i_st,lg,k) = 0.d0
           pRay_O2( j_st,i_st,lg,k) = 0.d0
           pCab_O2( j_st,i_st,lg,k) = 0.d0
           pRay_Ar( j_st,i_st,lg,k) = 0.d0
           pCab_Ar( j_st,i_st,lg,k) = 0.d0
           pRay_air(j_st,i_st,lg,k) = 0.d0
           pCab_air(j_st,i_st,lg,k) = 0.d0
        endif
        if((i_st.eq.1).AND.(j_st.eq.1).AND.(lg.eq.0))then
           pRam_air(j_st,i_st,lg,k) = 1.d0
        else
          pRam_air(j_st,i_st,lg,k) = 0.d0
        endif
      enddo
     enddo
    enddo
   enddo
44 continue

  if(verbose_output)then
     write(*,*) &
       'Output Expansion coefficients for Rayleigh scattering'
     write(*,*) 'lambda [nm] = ',wavelength_ns(ns_elas)
     write(*,'(A4,6A9)') &
       'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
     k = maxlay
     do lg=0,maxleg
        if(nstokes.eq.1)then
           write(*,'(I4,6F9.5)')  &
             lg,                  &
             pRay_air(1,1,lg,k),  &
             0.d0,                &
             0.d0,                &
             0.d0,                &
             0.d0,                &
             0.d0
        endif
        if(nstokes.eq.3)then
           write(*,'(I4,6F9.5)')  &
             lg,                  &
             pRay_air(1,1,lg,k),  &
             pRay_air(2,2,lg,k),  &
             pRay_air(3,3,lg,k),  &
             0.d0,                &
             pRay_air(1,2,lg,k),  &
             0.d0       
        endif
        if(nstokes.eq.4)then
           write(*,'(I4,6F9.5)')  &
             lg,                  &
             pRay_air(1,1,lg,k),  &
             pRay_air(2,2,lg,k),  &
             pRay_air(3,3,lg,k),  &
             pRay_air(4,4,lg,k),  &
             pRay_air(1,2,lg,k),  &
             pRay_air(3,4,lg,k)   
        endif
     enddo            
     write(*,*) &
       'Output Expansion coefficients for Cabannes scattering'
     write(*,'(A4,6A9)') & 
       'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
     k = maxlay
     do lg=0,maxleg
        if(nstokes.eq.1)then
           write(*,'(I4,6F9.5)') &
             lg,                 &
             pCab_air(1,1,lg,k), &
             0.d0,               &
             0.d0,               &
             0.d0,               &
             0.d0,               &
             0.d0                
         endif
         if(nstokes.eq.3)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pCab_air(1,1,lg,k), &
              pCab_air(2,2,lg,k), &
              pCab_air(3,3,lg,k), &
              0.d0,               &
              pCab_air(1,2,lg,k), &
              0.d0
         endif
         if(nstokes.eq.4)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pCab_air(1,1,lg,k), &
              pCab_air(2,2,lg,k), &
              pCab_air(3,3,lg,k), &
              pCab_air(4,4,lg,k), &
              pCab_air(1,2,lg,k), &
              pCab_air(3,4,lg,k)
         endif
      enddo         
      write(*,*) &
        'Output Expansion coefficients for Raman scattering'
      write(*,*) 'lambda [nm] = ',wavelength_ns(ns_elas)
      write(*,'(A4,6A9)') &
        'lg','alpha1','alpha2','alpha3','alpha4','beta1','beta2'
      k = maxlay
      do lg=0,maxleg
         if(nstokes.eq.1)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pRam_air(1,1,lg,k), &
              0.d0,               &
              0.d0,               &
              0.d0,               &
              0.d0,               &
              0.d0 
         endif
         if(nstokes.eq.3)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pRam_air(1,1,lg,k), &
              pRam_air(2,2,lg,k), &
              pRam_air(3,3,lg,k), &
              0.d0,               &
              pRam_air(1,2,lg,k), &
              0.d0
         endif
         if(nstokes.eq.4)then
            write(*,'(I4,6F9.5)') &
              lg,                 &
              pRam_air(1,1,lg,k), &
              pRam_air(2,2,lg,k), &
              pRam_air(3,3,lg,k), &
              pRam_air(4,4,lg,k), &
              pRam_air(1,2,lg,k), &
              pRam_air(3,4,lg,k)
         endif
      enddo
  endif

  return
end subroutine clearsky_optical_properties

!********************************************************************

subroutine RayCaB_Ram_expansion_coeff(A, B, C, N, coeff)

!--------------------------------------------------------------------
! Purpose:
!   Returns expansion coefficient matrix for Rayleigh, Cabannes 
!   scattering in the form of weights, 
!
!             ( alpha1  beta1     0       0     )
!     coeff = ( beta1   alpha2    0       0     )
!             (   0       0     alpha3   beta2  )
!             (   0       0     -beta2   alpha4 )
!
!   where expansion coeffcicients depend on Fourier moment lg=0,..,2
!   and wavelength [see Stam et al. JGR vol 107, D20, 2002]
!
!    lg   alpha1  alpha2  alpha3  alpha4      beta1      beta2
!   ------------------------------------------------------------
!     0      1       0       0       0          0           0
!     1      0       0       0      B/N         0           0
!     2    2A/3N    4A/N     0       0    4A/(N*(6**0.5))   0
!
!   the rest is zero. coeff=coeff(i1,i2,lg) where i1,i2 are index
!   for Stokes component, lg is Fourier moment index 
!--------------------------------------------------------------------

  implicit none
  integer,parameter::dp=selected_real_kind(8) !double precision
  integer::lg,i1,i2
  logical,parameter::verbose = .false.
    real(kind=dp)::A,B,C,N,coeff(4,4,0:2)   

  if(verbose) write(*,*) A,B,C,N

! Set all to zero
  do lg=0,2
     do i2=1,4
        do i1=1,4
           coeff(i1,i2,lg) = 0.d0 
        enddo
     enddo
  enddo

! alpha1
  coeff(1,1,0) = 1.d0
  coeff(1,1,1) = 0.d0
  coeff(1,1,2) = 2.d0*A/(3.d0*N)

! alpha2
  coeff(2,2,0) = 0.d0
  coeff(2,2,1) = 0.d0
  coeff(2,2,2) = 4.d0*A/N

! alpha3
  coeff(3,3,0) = 0.d0
  coeff(3,3,1) = 0.d0
  coeff(3,3,2) = 0.d0

! alpha4
  coeff(4,4,0) = 0.d0
  coeff(4,4,1) = B/N
  coeff(4,4,2) = 0.d0
    
! beta1
  coeff(1,2,0) = 0.d0
  coeff(1,2,1) = 0.d0
  coeff(1,2,2) = 4.d0*A/(dsqrt(6.d0)*N)

  coeff(2,1,0) = coeff(1,2,0)
  coeff(2,1,1) = coeff(1,2,1)
  coeff(2,1,2) = coeff(1,2,2)

! beta2
  coeff(3,4,0) = 0.d0
  coeff(3,4,1) = 0.d0
  coeff(3,4,2) = 0.d0

  coeff(4,3,0) = -coeff(3,4,0)
  coeff(4,3,1) = -coeff(3,4,1)
  coeff(4,3,2) = -coeff(3,4,2)

  return
end subroutine RayCaB_Ram_expansion_coeff
      



