
MODULE MSINSCA

  IMPLICIT NONE

  !! PUBLIC  :: CALL_SINSCA, SINSCA, SECSCA
  
CONTAINS

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

  SUBROUTINE pray(cteta, steta, delta, p11, p21)

    USE MCONSTANTS, only : dp

    implicit none

    !--- Input Arguments
    real(kind=dp), intent(in)  :: cteta
    real(kind=dp), intent(in)  :: steta
    real(kind=dp), intent(in)  :: delta
    real(kind=dp), intent(out) :: p11
    real(kind=dp), intent(out) :: p21

    !--- Local Arguments
    real(kind=dp) :: Gdelta
    real(kind=dp) :: Gdeltap

    !******************************************
    !     Compute Rayleigh phase matrix
    !******************************************

    if (delta .ne. 0D0) then

       GDelta  = (1.0D0 - delta) / (1.0D0 + delta / 2.0D0)
       GDeltap = (1.0D0 - 2.0D0*delta) / (1.0D0 - delta)

       p11 = GDelta * 3.0D0 * (1.0D0 + cteta**2.D0) / 4.0D0 &
            + (1.0D0 - GDelta)                                    ! P11
       P21 = - GDelta * 3.0D0 * steta**2.D0 /4.0D0  ! P21=P12

    else

       p11 = 3.0D0 / 4.0D0 *( 1.0D0 + cteta**2.D0 )   !P11
       p21 = 3.0D0 / 4.0D0 *( cteta**2.0D0 - 1.0D0 )  !P21=P12

    endif

    RETURN 

  END SUBROUTINE pray

  !---------------------------------------------------------------------------------------------------
  
  SUBROUTINE CALL_SINSCA

    USE MCONSTANTS, only : dp, undef_dp
    USE MCOMMON, only : verbose, &
         nmat, &
         mode, &
         nsza, &
         sza, &
         nvza, &
         vza, &
         nvaa, &
         vaa, &
         nlambda, &
         wlambda, &
         depol, &
         nlayers, &
         nalt_atm, &
         Svin, &
         F0, &
         nptcle, & 
         ptcle_nang_phasemat_layers, &
         ptcle_u_phasemat_layers, &
         ptcle_phasemat_layers, &
         ptcle_dtau_layers, &
         ptcle_ssa_layers, &
         surface_albedo, &
         surface_type, &
         surface_family, &
         n_aitaui_total, &
         n_aitaui_mono_layers, &
         nmax_aitaui_mono_layers, &
         ai_mono_layers, &
         gas_dtau_mono_layers, &
         gas_ssa_mono_layers, &
         rho_surf, &
         SvR, &
         SvR_sig, &
         flux_out, &
         rt_cputime, &
         thermal, &
         varsol_fact

    USE MUTILITY, only : REVERSE
    USE MLAYERS, only : GET_MONO_LAYERS

    IMPLICIT NONE

    INTEGER :: sinsca_nmu
    REAL(kind=dp), ALLOCATABLE :: sinsca_mu(:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_p11_ptcle(:,:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_21_ptcle(:,:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_tau_gas(:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_ssa_gas(:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_tau_ptcle(:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_ssa_ptcle(:)
    REAL(kind=dp) :: sinsca_Stoke_in(nmat)

    REAL(kind=dp) :: sinsca_surface_albedo

    REAL(kind=dp), ALLOCATABLE :: sinsca_Stoke(:,:,:,:)
    REAL(kind=dp), ALLOCATABLE :: sinsca_rho_surf(:,:)

    INTEGER :: i, j, k, u, iai

    ! time
    REAL :: tcpu1
    REAL :: tcpu2

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

    IF (thermal) THEN
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       WRITE(*,*) ' (sub. call_sinsca)  ERROR            '
       WRITE(*,*) '  Thermal emission is not implemented '
       WRITE(*,*) '  in call_sinsca                      '
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       STOP
    ENDIF

    IF (verbose) THEN
       WRITE (*,*) ''
       WRITE (*, *)                   '  (sub. call_sinsca) set arguments'
       WRITE (*,*) ''
       WRITE (*,FMT='(A,1x,I6,1x,A)') '                    SINSCA will be called ', n_aitaui_total, ' times'
       if (mode.eq.'kdis') &
            WRITE (*,FMT='(A,1x,I6,A)') '                  (maximum number of call per wavelength is ', nmax_aitaui_mono_layers,')'
    ENDIF

    ALLOCATE(sinsca_tau_gas(nlayers), sinsca_ssa_gas(nlayers),&
         sinsca_tau_ptcle(nlayers), sinsca_ssa_ptcle(nlayers), &
         sinsca_Stoke(nmat, nvza, nvaa, nalt_atm), &
         sinsca_rho_surf(nvza,nvaa))

    ! init in case no lambertian surface
    sinsca_surface_albedo = 0.0_dp

    DO j = 1, nlambda

       CALL GET_MONO_LAYERS(j)

       sinsca_nmu = ptcle_nang_phasemat_layers(j)

       ALLOCATE(sinsca_mu(sinsca_nmu), &
            sinsca_p11_ptcle(nlayers,sinsca_nmu), &
            sinsca_21_ptcle(nlayers,sinsca_nmu))

       if (surface_family .eq. 'lambert')  sinsca_surface_albedo = surface_albedo(j)

       ! we revesre mu for it to be sorted in increasing order
       sinsca_mu(:) = REVERSE(sinsca_nmu, ptcle_u_phasemat_layers(1:sinsca_nmu,j))

       if (nptcle .gt. 0) then
          do i = 1, nlayers
             ! we revert P11 and 21 cause mu was reverted to be 
             ! sorted in increasing order
             sinsca_p11_ptcle(i,:)  = REVERSE(sinsca_nmu, ptcle_phasemat_layers(i,1,1:sinsca_nmu,j)) ! P11
             sinsca_21_ptcle(i,:)   = REVERSE(sinsca_nmu, ptcle_phasemat_layers(i,5,1:sinsca_nmu,j)) ! P12 = P21
             sinsca_tau_ptcle(i)    = ptcle_dtau_layers(i,j) 
             sinsca_ssa_ptcle(i)    = ptcle_ssa_layers(i,j) 
          enddo
       else
          sinsca_p11_ptcle(:,:) = 0.0_dp
          sinsca_21_ptcle(:,:) = 0.0_dp
          sinsca_tau_ptcle(:)   = 0.0_dp
          sinsca_ssa_ptcle(:)   = 0.0_dp
       endif

       SvR(:,:,:,:,j)    = 0.0_dp
       flux_out(:,:,:,j) = 0.0_dp

       IF (verbose) then
          if (mode.eq.'kdis') then
             WRITE (*, FMT='(A,2x,F12.5,1x,A, 1x, I5)')  &
                  &'  (sub. call_sinsca) run SINSCA for wlambda=', wlambda(j),' microns with nai=', n_aitaui_mono_layers
          else
             WRITE (*, FMT='(A45,2x,F12.5,1x,A8)')  &
                  &'  (sub. call_sinsca) run SINSCA for wlambda=', wlambda(j),' microns'  
          endif
       END IF

       CALL CPU_TIME(tcpu1)

       do iai = 1, n_aitaui_mono_layers

          ! init result arrays
          sinsca_Stoke(:,:,:,:) = undef_dp 
          sinsca_rho_surf(:,:) = undef_dp
          sinsca_tau_gas(:)   = gas_dtau_mono_layers(:,iai) 
          sinsca_ssa_gas(:)   = gas_ssa_mono_layers(:,iai)

          DO i = 1, nsza

             sinsca_Stoke_in(1:nmat) = Svin(1:nmat) * F0(j) * varsol_fact(i)

             CALL SINSCA(verbose, &
                  wlambda(j), &
                  nmat, &
                  sza(i), nvza, vza, nvaa, vaa, &
                  nlayers, &
                  sinsca_nmu, sinsca_mu, sinsca_p11_ptcle, sinsca_21_ptcle, &
                  sinsca_tau_ptcle, &
                  sinsca_ssa_ptcle, &
                  sinsca_tau_gas,   &
                  sinsca_ssa_gas,   &
                  depol(j), &
                  sinsca_Stoke_in, &
                  surface_family, &
                  sinsca_surface_albedo, &
                  sinsca_Stoke, &
                  sinsca_rho_surf)

             ! implement the common variable for stokes vector for reflected intensities for each geometry
             DO u = 1, nvza
                DO k = 1, nvaa
                   SvR(i,u,k,1,j)    = SvR(i,u,k,1,j) + sinsca_Stoke(1,u,k,1) * ai_mono_layers(iai)
                   ! rho_surf(i,u,k,j) = rho_surf(i,u,k,j) + sinsca_rho_surf(u,k) * ai_mono_layers(iai)
                   if (nmat .gt. 1)then
                      SvR(i,u,k,2,j) = SvR(i,u,k,2,j) + sinsca_Stoke(2,u,k,1) * ai_mono_layers(iai)
                      SvR(i,u,k,3,j) = SvR(i,u,k,3,j) + sinsca_Stoke(3,u,k,1) * ai_mono_layers(iai)
                   endif
                   if (nmat .gt. 3) then 
                      SvR(i,u,k,4,j) = SvR(i,u,k,4,j) + sinsca_Stoke(4,u,k,1) * ai_mono_layers(iai)
                   endif
                ENDDO
             ENDDO

          ENDDO ! on sza

       end do

       CALL CPU_TIME(tcpu2)
       rt_cputime(j) = tcpu2-tcpu1
       ! IF (verbose) THEN
       !    WRITE (*,FMT='(1x,A35,1PE14.5, A9)') '  (sub. call_sinsca) CPU time     ', tcpu2-tcpu1, ' sec(s) -'
       ! ENDIF

       DEALLOCATE(sinsca_mu, &
            sinsca_p11_ptcle, &
            sinsca_21_ptcle)

    ENDDO ! on lambda

    SvR_sig(:,:,:,:,:)  = 0.0_dp

    DEALLOCATE(sinsca_tau_gas, sinsca_ssa_gas,&
         sinsca_tau_ptcle, sinsca_ssa_ptcle, &
         sinsca_Stoke, sinsca_rho_surf)

  END SUBROUTINE CALL_SINSCA

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

  SUBROUTINE SINSCA(verbose, &
       lambda, &
       nmat, &
       tetas, ntetav, tetav, nphiv, phiv, &
       nlay, &
       nmu, mu, p11_ptcle, p21_ptcle, &
       tau_ptcle, &
       ssa_ptcle, &
       tau_gas,   &
       ssa_gas,   &
       delta_ray, &
       I0, &
       surface, &
       surface_albedo, &
       SvR_levels, &
       rho_surf,&
       brdf_arr)

    !  Compute RT under single scattering assumption

    USE MCONSTANTS, only : dp, xpi, undef_dp, deg2rad, tiniestdp
    USE MUTILITY, only : LININTPOL
    USE MSURFACE_BRDF

    IMPLICIT NONE

    LOGICAL, INTENT(IN)       :: verbose
    REAL(KIND=dp), INTENT(IN) :: lambda    
    !    LOGICAL, INTENT(IN)       :: thermal
    INTEGER, INTENT(IN)       :: nmat
    REAL(KIND=dp), INTENT(IN) :: tetas
    INTEGER, INTENT(IN)       :: ntetav
    REAL(KIND=dp), INTENT(IN) :: tetav(ntetav)
    INTEGER, INTENT(IN)       :: nphiv
    REAL(KIND=dp), INTENT(IN) :: phiv(nphiv)
    INTEGER, INTENT(IN)       :: nlay
    INTEGER, INTENT(IN)       :: nmu
    REAL(KIND=dp), INTENT(IN) :: mu(nmu)
    REAL(KIND=dp), INTENT(IN) :: p11_ptcle(nlay,nmu)
    REAL(KIND=dp), INTENT(IN) :: p21_ptcle(nlay,nmu)
    REAL(KIND=dp), INTENT(IN) :: tau_ptcle(nlay)
    REAL(KIND=dp), INTENT(IN) :: ssa_ptcle(nlay)
    REAL(KIND=dp), INTENT(IN) :: tau_gas(nlay)
    REAL(KIND=dp), INTENT(IN) :: ssa_gas(nlay)
    REAL(KIND=dp), INTENT(IN) :: delta_ray
    REAL(KIND=dp), INTENT(IN) :: I0(nmat)
    CHARACTER(len=*), INTENT(IN) :: surface
    REAL(KIND=dp), INTENT(IN)  :: surface_albedo
    REAL(kind=dp), INTENT(OUT) :: rho_surf(ntetav, nphiv) 
    REAL(KIND=dp), INTENT(OUT) :: SvR_levels(nmat, ntetav, nphiv, 0:nlay) 
    REAL(KIND=dp), OPTIONAL, INTENT(IN)  :: brdf_arr(3, ntetav, nphiv)
    !                               brdf1 = brdf11
    !                               brdf2 = brdf21
    !                               brdf3 = brdf31

    ! local arguments
    INTEGER :: itetav
    INTEGER :: iphiv
    INTEGER :: ilay, ilay1
    INTEGER :: i, j

    REAL(kind=dp) :: umu0
    REAL(kind=dp) :: umu
    REAL(kind=dp) :: ctetad
    REAL(kind=dp) :: stetad
    REAL(kind=dp) :: i2, ci2

    REAL(KIND=dp) :: TAU( 0:nlay )
    REAL(KIND=dp) :: tau_ly(nlay)
    REAL(KIND=dp) :: tausca_ly(nlay)
    REAL(KIND=dp) :: ptcle_tausca_ly(nlay)
    REAL(KIND=dp) :: gas_tausca_ly(nlay)
    REAL(KIND=dp) :: ssa_ly(nlay)
    REAL(KIND=dp) :: p11_ly(nlay)
    REAL(KIND=dp) :: p21_ly(nlay)

    REAL(KIND=DP) :: tau_tot

    REAL(kind=dp) :: p11ray
    REAL(kind=dp) :: p21ray

    REAL(KIND=dp) ::  EXP0, EXP1
    REAL(KIND=dp) ::  I_SINSCA, Q_SINSCA
    REAL(KIND=dp) ::  I_SURF, Q_SURF, U_SURF
    REAL(kind=dp) ::  brdf(nmat,nmat)
    REAL(kind=dp) :: coeff


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

    if ( nmat .gt. 1) then
       if ( (I0(2) .ne. 0.0_dp) .or. (I0(3) .ne. 0.0_dp) .or. (I0(nmat) .ne. 0.0_dp) )  then
          write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          write(*,*) ' (sub. sinsca) :   ERROR   '
          WRITE(*,*) '     The -sinsca- RT model only apply for unpolarized '
          WRITE(*,*) '     incoming radiation at the top of the atmopshere (SvR={1,0,0,0}) '
          write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          STOP
       endif
    endif

    ! do j = 1, nmu
    !   print*, mu(j)
    ! end do

    ! mu must be sorted in increasing order
    do i = 1, nmu-1
       if (mu(i) .ge. mu(i+1)) then
          write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          write(*,*) ' (sub. sinsca) :   ERROR'
          write(*,*) '   mu must be sorted in increasing order'
          write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          STOP
       endif
    enddo

    SvR_levels(:,:,:,:) = 0.0_dp
    rho_surf(:,:) = 0.0_dp
    
    ! get the total optical depth of layers (ptcle+gas)
    tau( 0 ) = 0.0_dp
    tau_tot  = 0.0_dp
    do ilay = 1, nlay
       tau_ly(ilay)    = tau_ptcle(ilay) + tau_gas(ilay) 
       tausca_ly(ilay) = (tau_ptcle(ilay) * ssa_ptcle(ilay)) + &
            &            (tau_gas(ilay)   * ssa_gas(ilay))
       if (tau_ly(ilay).gt.0.0) then  
          ssa_ly(ilay)    = tausca_ly(ilay) / tau_ly(ilay)
       else
          ssa_ly(ilay)    = 1.0
       end if
       ptcle_tausca_ly(ilay) = tau_ptcle(ilay) * ssa_ptcle(ilay)
       gas_tausca_ly(ilay)   = tau_gas(ilay)   * ssa_gas(ilay)
       tau(ilay)             = tau(ilay-1) + tau_ly(ilay) 
       tau_tot               = tau_tot + tau_ly(ilay)  
    end do

    ! start to loop over geometries
    umu0 = COS(tetas* deg2rad)

    DO itetav = 1, ntetav

       UMU = COS(tetav(itetav) * deg2rad)

       coeff = 4.0_dp*XPI * ( 1.0_dp + UMU/UMU0 ) 

       DO iphiv = 1, nphiv

          ! Get the cosine of scattering angle
          ctetad = -UMU0 * UMU + SQRT( 1.0_dp - UMU0**2.0_dp ) * &
               SQRT( 1.0_dp - UMU**2.0_dp )  * COS( phiv(iphiv)*deg2rad )
          if (nmat .gt. 1) then
             STETAD  = SQRT(1.0_dp - CTETAD**2.0_dp)              
             if (STETAD.eq.0.0_dp) then
                CI2  = 1.0_dp
             else
                CI2  = ( - UMU0*SQRT( 1.0_dp - UMU**2.0_dp ) &
                     - UMU*SQRT( 1.0_dp - UMU0**2.0_dp ) * COS( phiv(iphiv)*deg2rad ) ) / &
                     STETAD
                IF (ci2 .GT.  1.0D0) ci2 =  1.0D0
                IF (ci2 .LT. -1.0D0) ci2 = -1.0D0                
             endif
             !--- The sign (-) is set to make a clockwise direction rotation to go from 
             !    the scattering plan to the local meridian plan
             if (phiv(iphiv).le.180.0_dp.and.phiv(iphiv).ge.0.0_dp)   i2 =  ACOS(ci2) 
             if (phiv(iphiv).gt.180.0_dp.and.phiv(iphiv).le.360.0_dp) i2 = -ACOS(ci2)
          else
             i2 = undef_dp
          endif

          ! get the phase matrix terms P11 and P21 for 
          ! each layer of the atmosphere (gas + particles)
          p11_ly(:) = 0.0_dp
          p21_ly(:) = 0.0_dp
          CALL pray(ctetad, stetad, delta_ray, p11ray, p21ray)
          do ilay = 1, nlay
             ! ptcle contribution
             if (ptcle_tausca_ly(ilay).gt.0.0_dp) then
                p11_ly(ilay) = ptcle_tausca_ly(ilay) / tausca_ly(ilay) *    &
                     LININTPOL(p11_ptcle(ilay,:), mu(:), nmu, ctetad)
                if (nmat.gt.1) then
                   p21_ly(ilay) = ptcle_tausca_ly(ilay) / tausca_ly(ilay) * &
                        LININTPOL(p21_ptcle(ilay,:), mu(:), nmu, ctetad)
                endif
             endif
             ! add rayleigh contribution
             if (gas_tausca_ly(ilay).gt.0.0_dp) then
                p11_ly(ilay) =  p11_ly(ilay) + &
                     gas_tausca_ly(ilay) / tausca_ly(ilay) * p11ray
                if (nmat.gt.1) then
                   p21_ly(ilay) =  p21_ly(ilay) + &
                        gas_tausca_ly(ilay) / tausca_ly(ilay) * p21ray
                endif
             end if
          end do

          !IF( ABS( UMU+UMU0 ).LE. tiniestdp * 10.0) THEN
          IF( ABS( UMU+UMU0 ).LE.1d-10) THEN

             ! print*, "here"
             ! print*, ABS( UMU+UMU0 )
             ! print*, tiniestdp

             !  Calculate downward intensity when
             !  UMU=UMU0, Eq. STWL (65e)

             DO ilay1 = 1, nlay
                I_SINSCA = 0.0_dp
                Q_SINSCA = 0.0_dp
                EXP0     = EXP( -TAU( ilay1) /UMU0 )
                DO ilay = 1, ilay1
                   I_SINSCA = I_SINSCA + SSA_LY( ilay )*P11_LY( ilay )*( TAU( ilay ) - TAU( ilay -1) )
                   if (nmat .gt. 1) Q_SINSCA = Q_SINSCA + SSA_LY( ilay )*P21_ly( ilay )*( TAU( ilay ) - TAU( ilay -1) )
                END DO

                SvR_levels(1,itetav,iphiv, ilay1) = I0(1) / ( 4.0_dp * XPI * UMU0  ) * EXP0 * I_SINSCA
                if (nmat .gt. 1) then
                   SvR_levels(2,itetav,iphiv,ilay1) = I0(1) / ( 4.0_dp * XPI * UMU0  )  * EXP0 * Q_SINSCA * COS(2.0_dp * i2 )
                   SvR_levels(3,itetav,iphiv,ilay1) = -I0(1) / ( 4.0_dp * XPI * UMU0  )  * EXP0 * Q_SINSCA * SIN(2.0_dp * i2 )
                endif
                if (nmat .eq. 4) SvR_levels(4,itetav,iphiv,ilay1) = 0.0_dp
                
             END DO

             I_SURF = 0.0
             Q_SURF = 0.0
             U_SURF = 0.0

          ELSE IF(UMU.GT.0.0_dp) THEN

             ! Upward intensity

             ! compute surface contribution
             IF (surface.eq.'lambert') then
                IF(surface_albedo.eq.0.0_dp) then
                   I_SURF = 0.0_dp
                ELSE
                   I_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * surface_albedo
                ENDIF
                IF (nmat.gt.1) then
                   Q_SURF = 0.0_dp
                   U_SURF = 0.0_dp
                END IF
             ELSE IF (surface.eq.'brdf_arr') then
                I_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * brdf_arr(1,itetav,iphiv)
                if (nmat .gt.1) then
                   Q_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * brdf_arr(2,itetav,iphiv)
                   U_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * brdf_arr(3,itetav,iphiv)
                end if
             ELSE IF (surface.eq.'brdf') then
                CALL surface_brdf(nmat, lambda, tetas, tetav(itetav), phiv(iphiv)+180.0, brdf(:,:))
                I_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * brdf(1,1)
                IF (nmat.gt.1) then
                   Q_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * brdf(2,1)
                   U_SURF = I0(1) * UMU0 / XPI * exp(-tau_tot/UMU0) * brdf(3,1)
                ENDIF
             ELSE
                write(*,*) ""
                write(*,*) " (sinsca) Error"
                write(*,*) "  Surface not implemented"
                write(*,*) ""
                STOP
             ENDIF

             ! compute intensity for each atmospheric layer
             DO ilay1 = 0, nlay

                ! atmospheric contribution
                if (ilay1<nlay) then
                   I_SINSCA = 0.0_dp
                   Q_SINSCA = 0.0_dp
                   EXP0     = EXP( -TAU( ilay1) /UMU0 )
                   DO ilay = ilay1+1, nlay
                      ! ** Upward intensity, Eq. STWL (65b)
                      EXP1 = EXP( -( ( TAU( ilay ) - TAU( ilay1 ) )/UMU + TAU( ilay )/UMU0 ) )
                      I_SINSCA = I_SINSCA + SSA_LY( ilay )*P11_LY( ilay )*( EXP0 - EXP1 )
                      if (nmat .gt. 1) Q_SINSCA = Q_SINSCA + SSA_LY( ilay )*P21_ly( ilay )*( EXP0 - EXP1 )
                      EXP0 = EXP1
                   ENDDO
                   SvR_levels(1,itetav,iphiv, ilay1) = I0(1) / coeff * I_SINSCA 
                   if (nmat .gt. 1) then
                      SvR_levels(2,itetav,iphiv,ilay1) = I0(1) / coeff * Q_SINSCA * COS(2.0_dp * i2 )
                      SvR_levels(3,itetav,iphiv,ilay1) = -I0(1) / coeff * Q_SINSCA * SIN(2.0_dp * i2 )
                   endif
                   if (nmat .eq. 4) SvR_levels(4,itetav,iphiv,ilay1) = 0.0_dp
                end if
                ! add surface contribution
                EXP0 = EXP( - ( tau_tot-TAU( ilay1) ) / UMU )
                SvR_levels(1,itetav,iphiv, ilay1) = SvR_levels(1,itetav,iphiv, ilay1) + I_SURF * EXP0
                if (nmat .gt. 1) then
                   SvR_levels(2,itetav,iphiv,ilay1) = SvR_levels(2,itetav,iphiv,ilay1) + Q_SURF * EXP0
                   SvR_levels(3,itetav,iphiv,ilay1) = SvR_levels(3,itetav,iphiv,ilay1) + U_SURF * EXP0
                endif
                if (nmat .eq. 4) SvR_levels(4,itetav,iphiv,ilay1) = 0.0_dp

             END DO

             rho_surf(itetav,iphiv) = I_SURF / I0(1) / UMU0 * XPI *  EXP( - ( tau_tot / UMU ))

          ELSE

             ! downward intensity
             ! compute intensity for each atmospheric layer
             DO ilay1 = 1, nlay

                ! atmospheric contribution
                I_SINSCA = 0.0_dp
                Q_SINSCA = 0.0_dp
                EXP0     = EXP( -TAU( ilay1) /UMU0 )
                DO ilay = ilay1, 1, -1
                   ! ** Upward intensity, Eq. STWL (65b)
                   EXP1 = EXP( -( ( TAU( ilay-1 ) - TAU( ilay1 ) )/UMU + TAU( ilay-1 )/UMU0 ) )
                   I_SINSCA = I_SINSCA + SSA_LY( ilay )*P11_LY( ilay )*( EXP0 - EXP1 )
                   if (nmat .gt. 1) Q_SINSCA = Q_SINSCA + SSA_LY( ilay )*P21_ly( ilay )*( EXP0 - EXP1 )
                   EXP0 = EXP1
                ENDDO
                SvR_levels(1,itetav,iphiv, ilay1) = I0(1) / coeff * I_SINSCA 
                if (nmat .gt. 1) then
                   SvR_levels(2,itetav,iphiv,ilay1) = I0(1) / coeff * Q_SINSCA * COS(2.0_dp * i2 )
                   SvR_levels(3,itetav,iphiv,ilay1) = -I0(1) / coeff * Q_SINSCA * SIN(2.0_dp * i2 )
                endif
                if (nmat .eq. 4) SvR_levels(4,itetav,iphiv,ilay1) = 0.0_dp


             END DO

          END IF


       END DO
    END DO


  END SUBROUTINE SINSCA

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

  SUBROUTINE SECSCA(verbose, &
       lambda, &
       nmat, &
       tetas, ntetav, tetav, nphiv, phiv, &
       nlay, &
       tau_gas, &
       ssa_gas, &
       tau_ptcle,&
       ssa_ptcle,&
       ftrunc_ptcle,&
       nmu, &
       mu, &
       p11_delta, &
       secsca_mu_eq, &
       secsca_neg_phas, &
       secsca_norm_phas, &
       I0, &
       dsecsca)
    
    ! Compute IMS correction for all atm levels
    ! This routine is largely inspirated by the cdisort.c file from libratran
    ! This corresponds to the implemenetation of the method described in
    ! Buras, R.; Dowling, T. & Emde, C.
    ! New secondary-scattering correction in DISORT with increased efficiency for forward scattering
    ! jqsrt, 2011, 112, 2028-2034

       USE MCONSTANTS, only : dp, &
       xpi, &
       undef_dp, &
       deg2rad, &
       tiniestdp, &
       min_ctheta_secsca, & 
       nf_delta_phase_secsca 
    
    USE MUTILITY, only : LININTPOL

    IMPLICIT NONE

    LOGICAL, INTENT(IN)       :: verbose
    REAL(KIND=dp), INTENT(IN) :: lambda    
    INTEGER, INTENT(IN)       :: nmat
    REAL(KIND=dp), INTENT(IN) :: tetas
    INTEGER, INTENT(IN)       :: ntetav
    REAL(KIND=dp), INTENT(IN) :: tetav(ntetav)
    INTEGER, INTENT(IN)       :: nphiv
    REAL(KIND=dp), INTENT(IN) :: phiv(nphiv)
    INTEGER, INTENT(IN)       :: nlay
    REAL(KIND=dp), INTENT(IN) :: tau_gas(nlay)
    REAL(KIND=dp), INTENT(IN) :: ssa_gas(nlay)    
    REAL(KIND=dp), INTENT(IN) :: tau_ptcle(nlay)
    REAL(KIND=dp), INTENT(IN) :: ssa_ptcle(nlay)
    REAL(KIND=dp), INTENT(IN) :: ftrunc_ptcle(nlay)
    INTEGER, INTENT(IN)       :: nmu
    REAL(KIND=dp), INTENT(IN) :: mu(nmu)
    REAL(KIND=dp), INTENT(IN) :: p11_delta(nlay,nmu) 
    REAL(KIND=dp), INTENT(IN) :: secsca_mu_eq(nlay,nf_delta_phase_secsca)
    LOGICAL, INTENT(IN)       :: secsca_neg_phas(nlay,nf_delta_phase_secsca)
    REAL(KIND=dp), INTENT(IN) :: secsca_norm_phas(nlay)
    REAL(KIND=dp), INTENT(IN)  :: I0(nmat)
    REAL(KIND=dp), INTENT(OUT) :: dsecsca(nmat, ntetav, nphiv, 0:nlay)
    !-----------


    ! local arguments
    INTEGER :: itetav
    INTEGER :: iphiv
    INTEGER :: ilay, ilay1
    INTEGER :: i, j

    REAL(kind=dp) :: umu0
    REAL(kind=dp) :: umu
    REAL(kind=dp) :: ctetad
    REAL(kind=dp) :: stetad

    REAL(KIND=dp) :: tausca_ly

    REAL(KIND=dp) :: TAU( 0:nlay )
    REAL(KIND=dp) :: tau_ly(nlay)
    REAL(KIND=dp) :: ssa_ly(nlay)
    
    REAL(KIND=DP) :: tau_tot

    REAL(KIND=DP) :: wbar( nlay )
    REAL(KIND=DP) :: fbar( nlay )
    REAL(KIND=DP) :: stau( nlay )

    REAL(KIND=DP) :: wbar_tmp
    REAL(KIND=DP) :: fbar_tmp
    REAL(KIND=DP) :: stau_tmp
    REAL(KIND=DP) :: umu0p
    REAL(KIND=DP) :: pspike
    REAL(KIND=DP) :: pspike1
    REAL(KIND=DP) :: pspike2
    REAL(KIND=DP) :: tiny = 1.e-4

    LOGICAL :: do_it(nlay)
    
    !-----------

    if (nmat.gt.1) then
       write(*,*) ""
       write(*,*) ""
       write(*,*) " (secsca) ERROR"
       write(*,*) "          Not implemented for nmat > 1"
       write(*,*) ""
       STOP
    end if

    ! mu must be sorted in increasing order
    do i = 1, nmu-1
       if (mu(i) .ge. mu(i+1)) then
          write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          write(*,*) ' (sub. secsca) :   ERROR'
          write(*,*) '         mu must be sorted in increasing order'
          write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          STOP
       endif
    enddo
   
    dsecsca(:,:,:,:) = 0.0_dp

    ! get the total optical depth of layers (ptcle+gas)
    tau( 0 ) = 0.0_dp
    tau_tot  = 0.0_dp
    do ilay = 1, nlay
       tau_ly(ilay)    = tau_ptcle(ilay) + tau_gas(ilay) 
       tausca_ly       = (tau_ptcle(ilay) * ssa_ptcle(ilay)) + &
            &            (tau_gas(ilay)   * ssa_gas(ilay)) 
       ssa_ly(ilay)    = tausca_ly / tau_ly(ilay)
       tau(ilay)       = tau(ilay-1) + tau_ly(ilay) 
       tau_tot         = tau_tot + tau_ly(ilay)  
    end do

    ! Calculate vertically averaged value of single scattering albedo and separated
    ! fraction f, eq. STWL (A.15)
    do_it(:) = .true.
    wbar(:) = 0.0
    fbar(:) = 0.0
    stau(:) = 0.0
    do ilay = 1, nlay
       do i = 1, ilay
          wbar(ilay) = wbar(ilay) + tau_ly(i) * ssa_ly(i) 
          fbar(ilay) = fbar(ilay) +                              &
               (tau_ptcle(i) * ssa_ptcle(i) * ftrunc_ptcle(i)) + &
               (tau_gas(i)   * ssa_gas(i))
          stau(ilay) = stau(ilay) + tau_ly(i)
       end do
       if ( (wbar(ilay) .lt. tiny ) .or. &
            (fbar(ilay) .lt. tiny ) .or. &
            (stau(ilay) .lt. tiny ) ) then
          do_it(ilay) = .false.
       end if
       fbar(ilay) = fbar(ilay) / wbar(ilay)
       wbar(ilay) = wbar(ilay) / stau(ilay)          
    end do
    
    
    
    ! start to loop over geometries
    umu0 = COS(tetas* deg2rad)

    DO itetav = 1, ntetav

       UMU = COS(tetav(itetav) * deg2rad)

       if (UMU.lt.0.0) then

          DO iphiv = 1, nphiv

             ! Get the cosine of scattering angle
             ctetad = -UMU0 * UMU + SQRT( 1.0_dp - UMU0**2.0_dp ) * &
                  SQRT( 1.0_dp - UMU**2.0_dp )  * COS( phiv(iphiv)*deg2rad )

             IF (ctetad.ge.min_ctheta_secsca) THEN

                do ilay = 1, nlay

                   if ( do_it(ilay)) then

                      pspike1 = LININTPOL(p11_delta(ilay,:), mu(:), nmu, ctetad)

                      pspike2 = calc_phase_squared(nmu, &
                           ctetad, &
                           nf_delta_phase_secsca, &
                           mu(:), &
                           p11_delta(ilay,:), &
                           secsca_mu_eq(ilay,:), &
                           secsca_neg_phas(ilay,:), &
                           secsca_norm_phas(ilay))                   

                      pspike = 2.*pspike1 - pspike2

                      umu0p = umu0 / (1.-fbar(ilay)*wbar(ilay))

                      ! Calculate IMS correction term, eq. STWL (A.13)
                      dsecsca(1, itetav, iphiv, ilay) = I0(1) / (4.0 * xpi) *          &
                           (fbar(ilay)*wbar(ilay))**2.0 / (1.-fbar(ilay)*wbar(ilay)) * &
                           pspike * &
                           xi_func(-UMU,umu0p,tau(ilay))
                     
                      !print*, UMU0, UMU, FBAR(ilay), WBAR(ilay), PSPIKE1, pspike2, &
                      !     xi_func(-UMU,umu0p,tau(ilay)), dsecsca(1, itetav, iphiv, ilay)

                   end if

                end do

             end IF

          END DO
       endif

    END DO





  END SUBROUTINE SECSCA

  !==================================================================
  
  function xi_func(umu1,umu2,tau)

    USE MCONSTANTS, only : dp
    
    real(kind=dp) :: xi_func

    real(kind=dp), intent(in) :: umu1
    real(kind=dp), intent(in) :: umu2
    real(kind=dp), intent(in) :: tau

    real(kind=dp) :: exp1,x1

    x1   = (umu2-umu1)/(umu2*umu1)
    exp1 = exp(-tau/umu1)

    if (x1 .ne. 0.) then
       xi_func = ((tau*x1-1.)*exp(-tau/umu2)+exp1)/(x1*x1*umu1*umu2)
    else 
       xi_func =  tau*tau*exp1/(2.*umu1*umu2)
    end if

  end function xi_func
  
 !==================================================================
 
  function calc_phase_squared(nmu, &
       ctheta, &
       nf, &
       mu, &
       p11, &
       mu_eq, &
       neg_phas, &
       norm_phas)

    USE MCONSTANTS, only : dp, xpi
    USE MUTILITY, only : LOCATE

    REAL (KIND=dp) :: calc_phase_squared

    integer, intent(in) :: nmu
    real(kind=dp), intent(in) :: ctheta
    integer, intent(in) :: nf
    real(kind=dp), intent(in) :: mu(nmu)
    real(kind=dp), intent(in) :: p11(nmu)
    real(kind=dp), intent(in) :: mu_eq(nf)
    logical, intent(in)       :: neg_phas(nf)
    real(kind=dp), intent(in) :: norm_phas

    integer:: j=0
    integer:: k=0
    integer:: it=0

    real(kind=dp) :: pspike2=0.0
    real(kind=dp) :: stheta=0.0
    real(kind=dp) :: smueq=0.0
    real(kind=dp) :: phint=0.0
    real(kind=dp) :: mumin=0.0
    real(kind=dp) :: mumax=0.0
    integer :: imin=0
    integer :: imax=0
    real(kind=dp) :: D=0.0
    real(kind=dp) :: C=0.0
    real(kind=dp) :: DDp=0.0
    real(kind=dp) :: Cp=0.0
    logical :: cutting=.false.

    stheta = sqrt( 1.0 - ctheta * ctheta )

    do j = 1, nf

       ! special case: second scattering angle does not depend on
       ! azimuth of first scattering angle
       if ( (ctheta.eq.1.0) .or. (mu_eq(j).eq.1.0)) then

          it = LOCATE( mu , nmu, mu_eq(j) * ctheta )

          phint = xpi * ( p11(it)              &
		       + ( mu_eq(j)*ctheta - MU(it) )  &
		       / ( MU(it+1) - MU(it) )         &
		       * ( P11(it+1) - P11(it) ) )

          if (ctheta.eq.1.0) phint = phint / 2.0

       else

          phint = 0.0

          smueq = sqrt ( 1. - MU_EQ(j)*MU_EQ(j) )

          ! locate integration borders 
          mumin = ctheta *  MU_EQ(j) - stheta * smueq
          mumax = ctheta *  MU_EQ(j) + stheta * smueq

          ! cut where mu_1 = mu_2 
          if (MU_EQ(j) .lt. mumax) then
             mumax = MU_EQ(j)
             cutting=.true.
          else
             cutting=.false.
          end if

          if (mumin.lt.mumax) then

             imin = locate ( mu, nmu, mumin)
             imax = locate ( mu, nmu, mumax)

             k=imin
             ! assuming SPF is linear in mu
             D = ( P11(k+1) - P11(k) ) / ( MU(k+1) - MU(k) )
             C = P11(k) - MU(k) * D

             phint = phint + ( D * ctheta * MU_EQ(j) + C ) * xpi / 2.0

             do k = imin+1, imax

                DDp = ( P11(k+1) - P11(k) ) / ( MU(k+1) - MU(k) )
                Cp = P11(k) - MU(k) * DDp

                phint = phint +                                      & 
                      ( DDp - D ) * sqrt ( 1.0 - ctheta * ctheta     &
                                           - MU_EQ(j) * MU_EQ(j)             &
                                           + 2.0 * ctheta * MU_EQ(j) * MU(k) &
                                           - MU(k) * MU(k) )                 &
                     + ( ( DDp - D )* ctheta * MU_EQ(j) + Cp - C ) * &
                     asin ( ( ctheta * MU_EQ(j) - MU(k) )            & 
                     / ( smueq * stheta ) ) 

                D=DDp
                C=Cp

             end do

             if (cutting.eqv..true.) then
                phint = phint +                              &
                     - D * sqrt ( 1.0 - ctheta * ctheta      &
                                 + 2.0 * MU_EQ(j) * MU_EQ(j) * &
                                 ( ctheta - 1.0 ) )            &
                     - ( D * ctheta * MU_EQ(j) + C ) *       &
                     asin ( ( ctheta - 1.0 ) * MU_EQ(j)      &
                     / ( smueq * stheta ) )
             else
                phint = phint + ( D * ctheta * MU_EQ(j) + C ) * xpi / 2.0
             end if

          endif

       end if

       if ( (j.eq.1) .or. (j.eq.nf) ) then 
          if ( NEG_PHAS(j).eqv..true. ) then
             pspike2 = pspike2 - 0.5 * phint
          else
             pspike2 = pspike2 + 0.5 * phint
          end if
       else
          if ( NEG_PHAS(j).eqv..true. ) then
             pspike2 = pspike2 - phint
          else
             pspike2 = pspike2 + phint
          end if
       end if

    end do

	pspike2 = pspike2 * norm_phas

    calc_phase_squared = pspike2

  end function calc_phase_squared

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



END MODULE MSINSCA
