! Description:
!> @file
!!   Calculate combined optical properties of aerosols and/or clouds.
!
!> @brief
!!   Calculate combined optical properties of aerosols and/or clouds.
!!
!! @details
!!   This subroutine takes the per-particle aerosol and/or hydrometeor
!!   optical properties from rttov_optp_interp and combines them as
!!   required by the selected solver(s).
!!
!!   If Rayleigh multiple-scattering is enabled via the DOM solar solver, 
!!   then it is treated as if it were an additional aerosol. The Rayleigh 
!!   optical properties are calculated at the start of the aerosol section,
!!   and are included with aerosol (if present), and then combined later
!!   with cloud properties (if present).
!!
!!   Not all quantities required for every solver are computed here.
!!
!!   The output layer optical depths in opdpacl and opdpaclsun and the
!!   accumulated optical depths in opdpac and opdpacsun take the local
!!   path angle into account. These are used for Chou-scaling and Rayleigh
!!   single-scattering.
!!
!!   The output layer extinction and scattering optical depths in opdpext_nadir
!!   and opdpsca_nadir are the *nadir* optical depths. These are used for DOM
!!   calculations.
!!
!!   Chou-scaling - thermal only:
!!     outputs required are the total optical depth (with Chou-scaled
!!     scattering op dep) for the local path. The back-scattering parameter
!!     (bpr) is required and this is used to scale the scattering op deps.
!!
!!   Delta-Eddington - thermal only:
!!     outputs required are total extinction coefficient (not optical depth),
!!     summed scattering coefficient, and summed asymmetry parameters weighted
!!     by scattering coefficients. These are computed separately for aerosols
!!     and hydrometeors. These are combined and final SSA and asym are
!!     computed elsewhere.
!!
!!   DOM - thermal:
!!     outputs are the nadir absorption and scattering optical depths and the
!!     first dom_nstreams Legendre coefficients for the combined phase function
!!     for all scattering particles.
!!
!!   Radar - thermal only:
!!     outputs are total extinction coefficient (not optical depth), and total
!!     reflectivity. These are weighted by hydro_frac.
!!
!!   DOM - solar:
!!     outputs are same as DOM thermal, but in addition the explicit phase fn
!!     is required. The phase fn is interpolated for the upward-scattered
!!     beam. The total optical depths are also computed as these are used for
!!     Rayleigh single-scattering calculations where relevant.
!!
!!   MFASIS-NN - solar only:
!!     the only output are the nadir per-particle extinction optical depths.
!!
!!   All layer outputs are computed for the clear (aerosol-only) and (when
!!   relevant) the cloudy(+aerosol) cases. The appropriate values are used
!!   later in each cloud column according to whether the layer is clear or
!!   cloudy in that particular column.
!!
!!
!! @param[in]     nlayers             number of layers in input profile
!! @param[in]     chanprof            specifies channels and profiles to simulate
!! @param[in]     opts                options to configure the simulations
!! @param[in]     cloud_columns       computed cloud column data
!! @param[in]     profiles            internal atmospheric profiles and surface variables
!! @param[in]     do_thermal          flag to indicate if any thermal (emissive) simulations are being performed
!! @param[in]     thermal             per-channel flag to indicate if thermal (emissive) simulations are being performed
!! @param[in]     do_solar            flag to indicate if any solar simulations are being performed
!! @param[in]     solar               per-channel flag to indicate if any solar simulations are being performed
!! @param[in]     coef                optical depth coefficients structure
!! @param[in]     coef_scatt          visible/IR scattering coefficients structure
!! @param[in]     optp_interp         interpolated optical properties
!! @param[in]     angles              geometry structure
!! @param[in]     raytracing          raytracing structure
!! @param[in,out] scatt_aux           structure containing scattering property data
!! @param[in,out] scatt_aux_column    computed optical depths per cloud column and phase functions
!! @param[in,out] transmission        output transmission structure
!
! Copyright:
!    This software was developed within the context of
!    the EUMETSAT Satellite Application Facility on
!    Numerical Weather Prediction (NWP SAF), under the
!    Cooperation Agreement dated 7 September 2021, between
!    EUMETSAT and the Met Office, UK, by one or more partners
!    within the NWP SAF. The partners in the NWP SAF are
!    the Met Office, ECMWF, DWD and MeteoFrance.
!
!    Copyright 2024, EUMETSAT, All Rights Reserved.
!
SUBROUTINE rttov_scatt_optp( &
              nlayers,            &
              chanprof,           &
              opts,               &
              cloud_columns,      &
              profiles,           &
              do_thermal,         &
              thermal,            &
              do_solar,           &
              solar,              &
              coef,               &
              coef_scatt,         &
              optp_interp,        &
              angles,             &
              raytracing,         &
              scatt_aux,          &
              scatt_aux_column,   &
              transmission)

  USE rttov_types, ONLY :  &
       rttov_chanprof,         &
       rttov_options,          &
       rttov_coef,             &
       rttov_profile_internal, &
       rttov_geometry,         &
       rttov_raytracing,       &
       rttov_scatt_aux,        &
       rttov_scatt_aux_column, &
       rttov_coef_scatt,       &
       rttov_optp_interp,      &
       rttov_cloud_columns,    &
       rttov_transmission
  USE rttov_kinds, ONLY : jpim, jplm
!INTF_OFF
  USE yomhook, ONLY : lhook, dr_hook, jphook
  USE rttov_kinds, ONLY : jprv
  USE rttov_const, ONLY :  &
       realtol,                 &
       deg2rad,                 &
       max_optical_depth,       &
       ice_baum_index,          &
       solar_solver_dom,        &
       solar_solver_mfasis_nn,  &
       thermal_solver_dom,      &
       thermal_solver_chou,     &
       thermal_solver_delta_edd
  USE rttov_types, ONLY : &
      rttov_scatt_aerhydro, &
      rttov_optp_interp_aerhydro
!INTF_ON
  IMPLICIT NONE
  INTEGER(KIND=jpim),                INTENT(IN)    :: nlayers
  TYPE(rttov_chanprof),              INTENT(IN)    :: chanprof(:)
  TYPE(rttov_options),               INTENT(IN)    :: opts
  TYPE(rttov_cloud_columns),         INTENT(IN)    :: cloud_columns
  TYPE(rttov_profile_internal),      INTENT(IN)    :: profiles(:)
  LOGICAL(KIND=jplm),                INTENT(IN)    :: do_thermal
  LOGICAL(KIND=jplm),                INTENT(IN)    :: thermal(SIZE(chanprof))
  LOGICAL(KIND=jplm),                INTENT(IN)    :: do_solar
  LOGICAL(KIND=jplm),                INTENT(IN)    :: solar(SIZE(chanprof))
  TYPE(rttov_coef),                  INTENT(IN)    :: coef
  TYPE(rttov_coef_scatt),            INTENT(IN)    :: coef_scatt
  TYPE(rttov_optp_interp), TARGET,   INTENT(IN)    :: optp_interp
  TYPE(rttov_geometry),              INTENT(IN)    :: angles(:)
  TYPE(rttov_raytracing),            INTENT(IN)    :: raytracing
  TYPE(rttov_scatt_aux),             INTENT(INOUT) :: scatt_aux
  TYPE(rttov_scatt_aux_column),      INTENT(INOUT) :: scatt_aux_column
  TYPE(rttov_transmission),          INTENT(INOUT) :: transmission
!INTF_END

  LOGICAL(KIND=jplm) :: do_dom, do_aer_or_ray_dom
  LOGICAL(KIND=jplm) :: do_chou_scaling, do_delta_edd, do_mfasis_nn
  LOGICAL(KIND=jplm) :: do_dom_chan, do_ray_dom_chan, do_mfasis_nn_chan
  LOGICAL(KIND=jplm) :: do_chou_scaling_chan, do_delta_edd_chan, do_chan
  INTEGER(KIND=jpim) :: nchanprof
  INTEGER(KIND=jpim) :: j, prof, chan
  INTEGER(KIND=jpim) :: col, coli, iaer, ihydro
  INTEGER(KIND=jpim) :: lev, lay
  REAL(KIND=jprv)    :: sca, opd, opdsun
  REAL(KIND=jprv)    :: relazi, gamma, ray_lcoef2, ray_phup, cosscata
  REAL(KIND=jprv)    :: lcoef(0:opts%scatt%dom_nstreams)

  TYPE(rttov_scatt_aerhydro),       POINTER :: aer, hydro
  TYPE(rttov_optp_interp_aerhydro), POINTER :: optp_interp_aer
  TYPE(rttov_optp_interp_aerhydro), POINTER :: optp_interp_hydro

  REAL(jphook) :: zhook_handle
!-----End of header------------------------------------------------------------
  IF (LHOOK) CALL DR_HOOK('RTTOV_SCATT_OPTP', 0_jpim, ZHOOK_HANDLE)
  nchanprof = SIZE(chanprof)

  do_chou_scaling = do_thermal .AND. opts%scatt%thermal_solver == thermal_solver_chou
  do_delta_edd = do_thermal .AND. opts%scatt%thermal_solver == thermal_solver_delta_edd
  do_dom = (do_solar .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. &
           (do_thermal .AND. opts%scatt%thermal_solver == thermal_solver_dom)
  do_mfasis_nn = do_solar .AND. opts%scatt%solar_solver == solar_solver_mfasis_nn
  do_aer_or_ray_dom = opts%scatt%aerosols .OR. opts%scatt%rayleigh_multi_scatt

  IF (do_thermal .AND. .NOT. do_delta_edd) scatt_aux%opdpacl = 0._jprv
  IF (do_solar .AND. .NOT. do_mfasis_nn) THEN
    scatt_aux%opdpaclsun = 0._jprv
    scatt_aux%phup = 0._jprv
  ENDIF
  IF (do_dom) THEN
    DO j = 1, nchanprof
      IF (ASSOCIATED(scatt_aux_column%phasefn(j)%lcoef)) &
          scatt_aux_column%phasefn(j)%lcoef = 0._jprv
    ENDDO
  ENDIF

  !----------------------------------------------------------------------------
  ! CALCULATE OPTICAL DEPTHS OF AEROSOLS
  !----------------------------------------------------------------------------

  IF (do_aer_or_ray_dom) THEN

    aer => scatt_aux%aer

    IF (do_thermal .AND. .NOT. do_delta_edd .OR. .NOT. do_mfasis_nn) THEN
      aer%opdpext = 0._jprv
      aer%opdpsca = 0._jprv
    ENDIF
    IF (do_chou_scaling) THEN
      aer%opdpscabpr = 0._jprv
      IF (opts%scatt%chou_tang_mod) aer%opdpsca = 0._jprv
    ENDIF
    IF (do_solar .AND. .NOT. do_mfasis_nn) aer%phtotup = 0._jprv
    IF (do_dom) aer%sca = 0._jprv

    optp_interp_aer => optp_interp%optp_interp_aer

    DO j = 1, nchanprof
      chan = chanprof(j)%chan
      prof = chanprof(j)%prof
      relazi = profiles(prof)%azangle - profiles(prof)%sunazangle

      do_dom_chan = (solar(j) .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. &
                    (thermal(j) .AND. opts%scatt%thermal_solver == thermal_solver_dom)
      do_ray_dom_chan = opts%scatt%rayleigh_multi_scatt .AND. solar(j) .AND. &
                        10000._jprv / coef%ff_cwn(chan) <= opts%rt_all%rayleigh_max_wavelength
      do_chou_scaling_chan = thermal(j) .AND. do_chou_scaling
      do_delta_edd_chan = thermal(j) .AND. do_delta_edd
      do_mfasis_nn_chan = solar(j) .AND. do_mfasis_nn
      do_chan = thermal(j) .OR. solar(j)

      ! Initialisation for Rayleigh phase function:
      !   0.75 * (1 + cos(theta)^2) for scattering angle theta
      IF (do_ray_dom_chan) THEN
        ! Due to strict plane-parallel geometry, ray_phup (used below) is fixed for all
        ! layers so can be calculated outside layer loop, and has zero TL/AD/K
        cosscata = - angles(prof)%coszen * angles(prof)%coszen_sun - &
                     angles(prof)%sinzen * angles(prof)%sinzen_sun * COS(relazi * deg2rad)
        gamma = coef%rayleigh_depol_gamma(chan)
        ray_lcoef2 = 0.5_jprv * (1._jprv - gamma) / (1._jprv + 2._jprv * gamma)
        ray_phup = (1._jprv + 3._jprv * gamma + (1._jprv - gamma) * cosscata**2) * &
                   0.75_jprv / (1._jprv + 2._jprv * gamma)
      ENDIF

      DO lay = 1, nlayers

        !----------------------------------------------------------------------
        ! Calculate combined aerosol parameters for layer
        !----------------------------------------------------------------------
        IF (do_dom_chan) THEN
          lcoef(:) = 0._jprv

          ! If DOM Rayleigh enabled calculate the Rayleigh properties
          IF (do_ray_dom_chan .AND. &
              profiles(prof)%p_half(lay+1) > opts%rt_all%rayleigh_min_pressure) THEN

            aer%opdpsca(lay,j) = scatt_aux%ray_sca(lay,j)
            aer%opdpext(lay,j) = scatt_aux%ray_sca(lay,j)

            ! Ensure Legendre coefficients are weighted consistently with the aerosol ones below
            ! Only zeroth and second Legendre coefs are non-zero
            sca = scatt_aux%ray_sca(lay,j) / raytracing%ltick(lay,prof)
            aer%sca(lay,j) = sca
            lcoef(0) = sca  ! Zeroth lcoef is 1.
            lcoef(2) = sca * ray_lcoef2
            aer%phtotup(lay,j) = ray_phup * sca
          ENDIF
        ENDIF

        IF (opts%scatt%aerosols) THEN
          IF (do_mfasis_nn_chan) THEN
            aer%opdpext_part(:,lay,j) = &
              raytracing%ltick(lay,prof) * optp_interp_aer%ext(:,lay,j)
          ELSEIF (do_chan) THEN
            !------------------------------------------------------------------
            ! Accumulate total aerosol optical parameters (all particles)
            !------------------------------------------------------------------

            IF (do_delta_edd_chan) THEN
              scatt_aux%ext_all(0,lay,j) = SUM(optp_interp_aer%ext(:,lay,j))
              scatt_aux%sca_all(0,lay,j) = SUM(optp_interp_aer%sca(:,lay,j))
              scatt_aux%asm_sum(0,lay,j) = SUM(optp_interp_aer%sca(:,lay,j) * &
                                               optp_interp_aer%asym(:,lay,j))
            ENDIF

            IF (thermal(j) .AND. .NOT. do_delta_edd_chan .OR. solar(j)) THEN
              aer%opdpext(lay,j) = aer%opdpext(lay,j) + &
                SUM(optp_interp_aer%ext(:,lay,j)) * raytracing%ltick(lay,prof)
              aer%opdpsca(lay,j) = aer%opdpsca(lay,j) + &
                SUM(optp_interp_aer%sca(:,lay,j)) * raytracing%ltick(lay,prof)
            ENDIF

            IF (do_chou_scaling_chan) THEN
              aer%opdpscabpr(lay,j) = &
                SUM(optp_interp_aer%sca(:,lay,j) * &
                    optp_interp_aer%bpr(:,lay,j)) * &
                    raytracing%ltick(lay,prof)
            ENDIF

            IF (solar(j)) THEN
              aer%phtotup(lay,j) = aer%phtotup(lay,j) + &
                                   SUM(optp_interp_aer%phfn_int(:,lay,j) * &
                                       optp_interp_aer%sca(:,lay,j))
            ENDIF
          ENDIF
        ENDIF

        !------------------------------------------------------------------
        ! Calculate final phase function Leg. coefs for all aerosol types
        !------------------------------------------------------------------
        IF (do_dom_chan) THEN
          IF (opts%scatt%aerosols) THEN
            aer%sca(lay,j) = aer%sca(lay,j) + SUM(optp_interp_aer%sca(:,lay,j))
            DO iaer = 1, MAX(coef_scatt%optp_aer%ntypes, 1_jpim)
              IF (optp_interp_aer%sca(iaer,lay,j) <= 0) CYCLE
              lcoef(:) = lcoef(:) + optp_interp_aer%lcoef(:,iaer,lay,j) * &
                                    optp_interp_aer%sca(iaer,lay,j)
            ENDDO
          ENDIF
          IF (ABS(aer%sca(lay,j)) > realtol) THEN
            scatt_aux_column%phasefn(j)%lcoef(:,0,lay) = lcoef(:) / aer%sca(lay,j)
          ENDIF
        ENDIF
      ENDDO ! layers

      !------------------------------------------------------------------------
      ! Calculate total aerosol optical depths
      !------------------------------------------------------------------------
      IF (thermal(j)) THEN
        IF (do_chou_scaling_chan) THEN
          ! Chou-scaled optical depth for thermal channels
          aer%opdp(:,j) = aer%opdpext(:,j) - aer%opdpsca(:,j) + aer%opdpscabpr(:,j)
        ELSE
          ! Full optical depth for thermal channels
          aer%opdp(:,j) = aer%opdpext(:,j)
        ENDIF
      ENDIF

      IF (solar(j) .AND. .NOT. do_mfasis_nn_chan) THEN
        WHERE (ABS(aer%sca(:,j)) > realtol)
          scatt_aux%phup(0,:,j) = aer%phtotup(:,j) / aer%sca(:,j)
        ENDWHERE
      ENDIF
    ENDDO ! chanprof
  ENDIF ! do_aer_or_ray_dom


  !----------------------------------------------------------------------------
  ! CALCULATE OPTICAL DEPTHS OF CLOUDS
  !----------------------------------------------------------------------------
  IF (opts%scatt%hydrometeors) THEN

    hydro => scatt_aux%hydro

    IF (do_chou_scaling .AND. opts%scatt%chou_tang_mod) THEN
      hydro%opdpscabpr = 0._jprv
      hydro%opdpsca = 0._jprv
    ENDIF

    optp_interp_hydro => optp_interp%optp_interp_hydro

    DO j = 1, nchanprof
      chan = chanprof(j)%chan
      prof = chanprof(j)%prof
      relazi = profiles(prof)%azangle - profiles(prof)%sunazangle

      do_dom_chan = (solar(j) .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. &
                    (thermal(j) .AND. opts%scatt%thermal_solver == thermal_solver_dom)
      do_chou_scaling_chan = thermal(j) .AND. do_chou_scaling
      do_delta_edd_chan = thermal(j) .AND. do_delta_edd
      do_mfasis_nn_chan = solar(j) .AND. do_mfasis_nn
      do_chan = thermal(j) .OR. solar(j)

      DO lay = 1, nlayers
        IF (do_mfasis_nn_chan) THEN
          hydro%opdpext_part(:,lay,j) = &
            raytracing%ltick(lay,prof) * optp_interp_hydro%ext(:,lay,j)
        ELSEIF (do_chan) THEN
          !------------------------------------------------------------------
          ! Accumulate total cloud optical parameters (all particles)
          !------------------------------------------------------------------

          IF (opts%scatt%radar) THEN
            IF (opts%cloud_overlap%per_hydro_frac) THEN
              scatt_aux%ext_radar(lay,j) = SUM(optp_interp_hydro%ext_radar(:,lay,j) * &
                                               profiles(prof)%hydro_frac(:,lay))
              scatt_aux%zef(lay,j) = SUM(optp_interp_hydro%zef(:,lay,j) * &
                                         profiles(prof)%hydro_frac(:,lay))
            ELSE
              scatt_aux%ext_radar(lay,j) = SUM(optp_interp_hydro%ext_radar(:,lay,j)) * &
                                           profiles(prof)%hydro_frac(1,lay)
              scatt_aux%zef(lay,j) = SUM(optp_interp_hydro%zef(:,lay,j)) * &
                                     profiles(prof)%hydro_frac(1,lay)
            ENDIF
          ENDIF

          IF (do_delta_edd_chan) THEN
            scatt_aux%ext_all(1,lay,j) = SUM(optp_interp_hydro%ext(:,lay,j))
            scatt_aux%sca_all(1,lay,j) = SUM(optp_interp_hydro%sca(:,lay,j))
            scatt_aux%asm_sum(1,lay,j) = SUM(optp_interp_hydro%sca(:,lay,j) * &
                                             optp_interp_hydro%asym(:,lay,j))
          ENDIF

          IF (thermal(j) .AND. .NOT. do_delta_edd_chan .OR. solar(j)) THEN
            hydro%opdpext(lay,j) = &
              SUM(optp_interp_hydro%ext(:,lay,j)) * raytracing%ltick(lay,prof)
            hydro%opdpsca(lay,j) = &
              SUM(optp_interp_hydro%sca(:,lay,j)) * raytracing%ltick(lay,prof)
          ENDIF

          IF (do_chou_scaling_chan) THEN
            hydro%opdpscabpr(lay,j) = &
              SUM(optp_interp_hydro%sca(:,lay,j) * &
                  optp_interp_hydro%bpr(:,lay,j)) * &
                  raytracing%ltick(lay,prof)
          ENDIF

          IF (solar(j)) THEN
            hydro%phtotup(lay,j) = SUM(optp_interp_hydro%phfn_int(:,lay,j) * &
                                     optp_interp_hydro%sca(:,lay,j))
          ENDIF
        ENDIF

        !------------------------------------------------------------------
        ! Calculate final phase function Leg. coefs for all cloud types
        !------------------------------------------------------------------
        IF (do_dom_chan) THEN
          lcoef(:) = 0._jprv
          hydro%sca(lay,j) = SUM(optp_interp_hydro%sca(:,lay,j))
          DO ihydro = 1, MAX(profiles(prof)%nhydro, 1_jpim)
            IF (optp_interp_hydro%sca(ihydro,lay,j) <= 0) CYCLE
            lcoef(:) = lcoef(:) + optp_interp_hydro%lcoef(:,ihydro,lay,j) * &
                                  optp_interp_hydro%sca(ihydro,lay,j)
          ENDDO
          IF (ABS(hydro%sca(lay,j)) > realtol) THEN
            scatt_aux_column%phasefn(j)%lcoef(:,1,lay) = lcoef(:) / hydro%sca(lay,j)
          ENDIF
        ENDIF
      ENDDO ! layers

      !------------------------------------------------------------------------
      ! Calculate total cloud optical depths
      !------------------------------------------------------------------------
      IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN
        IF (do_chou_scaling_chan) THEN
          ! Chou-scaled optical depth for thermal channels
          hydro%opdp(:,j) = hydro%opdpext(:,j) - hydro%opdpsca(:,j) + hydro%opdpscabpr(:,j)
        ELSE
          ! Full optical depth for thermal channels
          hydro%opdp(:,j) = hydro%opdpext(:,j)
        ENDIF
      ENDIF
    ENDDO ! chanprof

    !------------------------------------------------------------------------
    ! Output cloudy transmittances for all channels - direct model only
    !------------------------------------------------------------------------
    ! This is done here (rather than rttov_transmit* for example) because we can
    ! do the calculations for thermal and solar calculations in one place here
    DO j = 1, nchanprof
      IF (.NOT. (thermal(j) .OR. solar(j))) CYCLE

      prof = chanprof(j)%prof
      chan = chanprof(j)%chan

      IF (thermal(j) .AND. do_delta_edd) THEN
        DO lay = 1, nlayers
          transmission%tau_levels_cld(lay+1,j) = transmission%tau_levels_cld(lay,j) + &
            angles(prof)%seczen * scatt_aux%ext_all(1,lay,j) * raytracing%ltick(lay,prof)
        ENDDO
      ELSE IF (solar(j) .AND. do_mfasis_nn) THEN
        DO lay = 1, nlayers
          transmission%tau_levels_cld(lay+1,j) = transmission%tau_levels_cld(lay,j) + &
            raytracing%pathsat(lay,prof) * SUM(hydro%opdpext_part(1:ice_baum_index,lay,j))
        ENDDO
      ELSE
        DO lay = 1, nlayers
          transmission%tau_levels_cld(lay+1,j) = transmission%tau_levels_cld(lay,j) + &
            raytracing%pathsat(lay,prof) * hydro%opdpext(lay,j)
        ENDDO
      ENDIF

      transmission%tau_total_cld(j) = transmission%tau_levels_cld(nlayers+1,j)

      transmission%tau_levels_cld(:,j) = MIN(max_optical_depth, transmission%tau_levels_cld(:,j))
      transmission%tau_levels_cld(:,j) = EXP(-transmission%tau_levels_cld(:,j))
      transmission%tau_total_cld(j)    = MIN(max_optical_depth, transmission%tau_total_cld(j))
      transmission%tau_total_cld(j)    = EXP(-transmission%tau_total_cld(j))
    ENDDO
  ENDIF ! hydrometeors

  !----------------------------------------------------------------------------
  ! CALCULATE TOTAL OPTICAL DEPTHS AND PARAMETERS FOR EACH CLOUD COLUMN
  !----------------------------------------------------------------------------
  DO j = 1, nchanprof
    IF (solar(j) .AND. do_mfasis_nn) CYCLE

    chan = chanprof(j)%chan
    prof = chanprof(j)%prof

    do_dom_chan = (solar(j) .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. &
                  (thermal(j) .AND. opts%scatt%thermal_solver == thermal_solver_dom)
    do_chou_scaling_chan = thermal(j) .AND. do_chou_scaling
    do_delta_edd_chan = thermal(j) .AND. do_delta_edd

    ! For layer-specific quantities store just the non-cloudy and cloudy values
    ! When used later the code looks up the appropriate value for each cloud column

    ! Determine layer total aerosol/cloud optical depths
    IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN
      ! These may be Chou-scaled
      IF (opts%scatt%aerosols) THEN
        scatt_aux%opdpacl(0,:,j) = &
          aer%opdp(:,j) * raytracing%pathsat(:,prof)
      ENDIF
      IF (opts%scatt%hydrometeors) THEN
        scatt_aux%opdpacl(1,:,j) = scatt_aux%opdpacl(0,:,j) + &
          hydro%opdp(:,j) * raytracing%pathsat(:,prof)
      ENDIF
    ENDIF
    IF (solar(j)) THEN
      IF (do_aer_or_ray_dom) THEN
        scatt_aux%opdpaclsun(0,:,j) = &
          aer%opdpext(:,j) * raytracing%patheff(:,prof)
      ENDIF
      IF (opts%scatt%hydrometeors) THEN
        scatt_aux%opdpaclsun(1,:,j) = scatt_aux%opdpaclsun(0,:,j) + &
          hydro%opdpext(:,j) * raytracing%patheff(:,prof)
      ENDIF
    ENDIF

    IF (do_chou_scaling_chan .AND. opts%scatt%chou_tang_mod) THEN
      ! NB opdpacl above is Chou-scaled, here we compute full exinction

      ! Compute layer extinction and scattering optical depths, and bpr values
      IF (opts%scatt%aerosols) THEN
        scatt_aux%opdpext(0,:,j) = aer%opdpext(:,j) * raytracing%pathsat(:,prof)
        scatt_aux%opdpsca(0,:,j) = aer%opdpsca(:,j) * raytracing%pathsat(:,prof)

        WHERE (aer%opdpsca(:,j) > 0._jprv)
          scatt_aux%bpr(0,:,j) = aer%opdpscabpr(:,j) / aer%opdpsca(:,j)
        ELSEWHERE
          scatt_aux%bpr(0,:,j) = 0._jprv
        ENDWHERE
      ELSE
        scatt_aux%opdpext(0,:,j) = 0._jprv
        scatt_aux%opdpsca(0,:,j) = 0._jprv
        scatt_aux%bpr(0,:,j) = 0._jprv
      ENDIF

      IF (opts%scatt%hydrometeors) THEN
        scatt_aux%opdpext(1,:,j) = scatt_aux%opdpext(0,:,j) + &
          hydro%opdpext(:,j) * raytracing%pathsat(:,prof)
        scatt_aux%opdpsca(1,:,j) = scatt_aux%opdpsca(0,:,j) + &
          hydro%opdpsca(:,j) * raytracing%pathsat(:,prof)

        IF (opts%scatt%aerosols) THEN
          WHERE (aer%opdpsca(:,j) + hydro%opdpsca(:,j) > 0._jprv)
            scatt_aux%bpr(1,:,j) = (aer%opdpscabpr(:,j) + hydro%opdpscabpr(:,j)) / &
                                        (aer%opdpsca(:,j) + hydro%opdpsca(:,j))
          ELSEWHERE
            scatt_aux%bpr(1,:,j) = 0._jprv
          ENDWHERE
        ELSE
          WHERE (hydro%opdpsca(:,j) > 0._jprv)
            scatt_aux%bpr(1,:,j) = hydro%opdpscabpr(:,j) / hydro%opdpsca(:,j)
          ELSEWHERE
            scatt_aux%bpr(1,:,j) = 0._jprv
          ENDWHERE
        ENDIF
      ENDIF

    ENDIF

    IF (do_dom_chan) THEN

      ! Determine final layer *nadir* absorption and scattering optical depths
      IF (do_aer_or_ray_dom) THEN
        scatt_aux%opdpext_nadir(0,:,j) = aer%opdpext(:,j)
        scatt_aux%opdpsca_nadir(0,:,j) = aer%opdpsca(:,j)
      ELSE
        scatt_aux%opdpext_nadir(0,:,j) = 0._jprv
        scatt_aux%opdpsca_nadir(0,:,j) = 0._jprv
      ENDIF
      IF (opts%scatt%hydrometeors) THEN
        scatt_aux%opdpext_nadir(1,:,j) = scatt_aux%opdpext_nadir(0,:,j) + &
          hydro%opdpext(:,j)
        scatt_aux%opdpsca_nadir(1,:,j) = scatt_aux%opdpsca_nadir(0,:,j) + &
          hydro%opdpsca(:,j)
      ENDIF

      DO lay = 1, nlayers
        ! Determine combined phase functions for aer+hydro case
        IF (do_aer_or_ray_dom .AND. opts%scatt%hydrometeors) THEN
          IF (aer%sca(lay,j) + hydro%sca(lay,j) > 0._jprv) THEN
            scatt_aux_column%phasefn(j)%lcoef(:,1,lay) = &
                (scatt_aux_column%phasefn(j)%lcoef(:,0,lay) * aer%sca(lay,j) + &
                 scatt_aux_column%phasefn(j)%lcoef(:,1,lay) * hydro%sca(lay,j)) / &
                (aer%sca(lay,j) + hydro%sca(lay,j))
          ENDIF
        ENDIF

        IF (solar(j)) THEN
          IF (opts%scatt%hydrometeors .AND. do_aer_or_ray_dom) THEN
            IF (ABS(hydro%sca(lay,j) + aer%sca(lay,j)) > realtol) THEN
              scatt_aux%phup(1,lay,j) = &
                (aer%phtotup(lay,j) + hydro%phtotup(lay,j)) / &
                (hydro%sca(lay,j) + aer%sca(lay,j))
            ENDIF
          ELSEIF (opts%scatt%hydrometeors) THEN
            IF (ABS(hydro%sca(lay,j)) > realtol) THEN
              scatt_aux%phup(1,lay,j) = hydro%phtotup(lay,j) / hydro%sca(lay,j)
            ENDIF
          ENDIF
        ENDIF

      ENDDO ! layers
    ENDIF ! do dom

    DO col = 0, cloud_columns%ncolumn(prof)
      IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN
        opd = 0._jprv
        scatt_aux_column%opdpac(1,col,j) = 0._jprv
      ENDIF
      IF (solar(j)) THEN
        opdsun = 0._jprv
        scatt_aux_column%opdpacsun(1,col,j) = 0._jprv
      ENDIF
      IF (col == 0) THEN
        IF (do_aer_or_ray_dom) THEN
          DO lay = 1, nlayers
            lev = lay + 1
            IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN
              opd = opd + scatt_aux%opdpacl(0,lay,j)
              scatt_aux_column%opdpac(lev,col,j) = opd
            ENDIF
            IF (solar(j)) THEN
              opdsun = opdsun + scatt_aux%opdpaclsun(0,lay,j)
              scatt_aux_column%opdpacsun(lev,col,j) = opdsun
            ENDIF
          ENDDO ! layers
        ELSE
          IF (thermal(j) .AND. .NOT. do_delta_edd_chan) scatt_aux_column%opdpac(:,col,j) = 0._jprv
          IF (solar(j)) scatt_aux_column%opdpacsun(:,col,j) = 0._jprv
        ENDIF
      ELSE
        DO lay = 1, nlayers
          lev = lay + 1
          coli = cloud_columns%icldarr(col,lay,prof) ! This is 0 or 1 for clear(+aer)/cloud(+aer)
          IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN
            opd = opd + scatt_aux%opdpacl(coli,lay,j)
            scatt_aux_column%opdpac(lev,col,j) = opd
          ENDIF
          IF (solar(j)) THEN
            opdsun = opdsun + scatt_aux%opdpaclsun(coli,lay,j)
            scatt_aux_column%opdpacsun(lev,col,j) = opdsun
          ENDIF
        ENDDO ! layers
      ENDIF ! col == 0
    ENDDO ! col
  ENDDO ! channels

  IF (LHOOK) CALL DR_HOOK('RTTOV_SCATT_OPTP', 1_jpim, ZHOOK_HANDLE)

END SUBROUTINE rttov_scatt_optp
