! Description:
!> @file
!!   K of calculation of the fraction of solar radiance that is
!!   reflected by a wind roughened water surface.
!
!> @brief
!!   K of calculation of the fraction of solar radiance that is
!!   reflected by a wind roughened water surface.
!!
!! @param[in]     opts           options to configure the simulations
!! @param[in]     chanprof       specifies channels and profiles to simulate
!! @param[in]     solar          flag to indicate channels with solar contribution
!! @param[in]     profiles       input atmospheric profiles and surface variables
!! @param[in,out] profiles_k     Jacobian wrt atmospheric profile and surface variables
!! @param[in]     emis_refl      input/output emissivities/reflectances
!! @param[in]     sunglint       internal structure for sea surface BRDF model variables
!! @param[in,out] sunglint_k     Jacobian wrt sea surface BRDF model variables
!! @param[in]     raytracing     raytracing structure
!! @param[in,out] raytracing_k   Jacobian wrt raytracing variables
!
! 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_calc_sunglint_k( &
              opts,         &
              chanprof,     &
              solar,        &
              profiles,     &
              profiles_k,   &
              emis_refl,    &
              sunglint,     &
              sunglint_k,   &
              raytracing,   &
              raytracing_k)

  USE rttov_types, ONLY : &
      rttov_options,          &
      rttov_chanprof,         &
      rttov_profile_internal, &
      rttov_emis_refl,        &
      rttov_raytracing,       &
      rttov_sunglint
  USE rttov_kinds, ONLY : jplm
!INTF_OFF
  USE rttov_kinds, ONLY : jprv, jpim
  USE rttov_const, ONLY : &
      realtol,                  &
      deg2rad, pi,              &
      gravity,                  &
      surftype_sea,             &
      min_windsp,               &
      max_sol_zen,              &
      max_exp_exponent,         &
      min_exponent,             &
      nk_elf, cm_elf,           &
      km_elf, x0_elf,           &
      a0_elf, ap_elf,           &
      solar_refl_model_elfouhaily
  USE yomhook, ONLY : lhook, dr_hook, jphook
!INTF_ON
  IMPLICIT NONE

  TYPE(rttov_options),          INTENT(IN)    :: opts
  TYPE(rttov_chanprof),         INTENT(IN)    :: chanprof(:)
  LOGICAL(jplm),                INTENT(IN)    :: solar(SIZE(chanprof))
  TYPE(rttov_profile_internal), INTENT(IN)    :: profiles(:)
  TYPE(rttov_profile_internal), INTENT(INOUT) :: profiles_k(SIZE(chanprof))
  TYPE(rttov_emis_refl),        INTENT(IN)    :: emis_refl(:)
  TYPE(rttov_raytracing),       INTENT(IN)    :: raytracing
  TYPE(rttov_raytracing),       INTENT(INOUT) :: raytracing_k
  TYPE(rttov_sunglint),         INTENT(IN)    :: sunglint
  TYPE(rttov_sunglint),         INTENT(INOUT) :: sunglint_k
!INTF_END

  INTEGER(jpim) :: j, i, isurf, prof, nchanprof, nlayers
  REAL(jprv)    :: csi_k                                            ! Angle between the zenith and the proJection of the
                                                                         ! incident ray on the x-z plane
  REAL(jprv)    :: alfa_k                                           ! Angle between the incident ray and the x-z plane
  REAL(jprv)    :: c_shad_k                                         ! The average magnitude of tan(alfa)
  REAL(jprv)    :: p_prime_k                                        ! The probability density of tan(alfa)
  REAL(jprv)    :: pxy_gammaxy_k                                    ! The joint probability density of the along-view and
                                                                         ! cross view slope
  REAL(jprv)    :: gamma_o_k                                        ! The mean square of the along-view (X axis) slope
  REAL(jprv)    :: gamma_p_k                                        ! The mean square of the cross-view (Y axis) slope
  REAL(jprv)    :: g_shad_k                                         ! Normalization function
  REAL(jprv)    :: gammax_k                                         ! The x-slope of the facet
  REAL(jprv)    :: theta_fi                                         ! First order shadowing factor for the reflected ray
  REAL(jprv)    :: theta_csi                                        ! First ordesr shadowing factor for the incident ray
  REAL(jprv)    :: q_shad_k                                         ! Second order shadowing factor
  REAL(jprv)    :: q_shad_a_k                                       ! Second order shadowing factor
  REAL(jprv)    :: q_shad_b_k                                       ! Second order shadowing factor
  REAL(jprv)    :: zensat_k                                         ! Zenith angle of satellite viewing angle at surface
  REAL(jprv)    :: zensun_k                                         ! Zenith angle of sun at surface
  REAL(jprv)    :: fac1_k
  REAL(jprv)    :: a_shad_k
  REAL(jprv)    :: b_shad_k
  REAL(jprv)    :: lambda_a_k
  REAL(jprv)    :: lambda_b_k
  REAL(jprv)    :: k
  REAL(jprv)    :: windsp_k
  REAL(jprv)    :: wangl_k
  REAL(jprv)    :: sqrt2pi
  REAL(jprv)    :: exp_max_exp_exponent_r

  REAL(jprv)    :: karr(0:nk_elf), sk2_k, sk2o_k(0:nk_elf), sk2p_k(0:nk_elf)
  REAL(jprv)    :: sqrt_k_kp, cm_c_pow25, c_cp_pow25
  REAL(jprv)    :: cos2wangl, cos2wangl_k
  REAL(jprv)    :: windsp, k0, x, kp, gamma, sigma
  REAL(jprv)    :: k0_k, x_k, kp_k, gamma_k, sigma_k
  REAL(jprv)    :: cp, omega, alpha_p, ustar, am
  REAL(jprv)    :: cp_k, omega_k, alpha_p_k, ustar_k, alpha_m_k, am_k
  REAL(jprv)    :: fp, fp_k, bl_k, bh_k
  REAL(jprv)    :: omega_c_k, lpm_k, jp_k, gamma_exp_k, fpexp_k, dk_k
  REAL(jprv)    :: cosh_tmp, cp_r, c_r, kp_r, sigma_r2, sigma_r3

  REAL(jphook) :: zhook_handle
!-----End of header-------------------------------------------------------------
  IF (LHOOK) CALL DR_HOOK('RTTOV_CALC_SUNGLINT_K', 0_jpim, ZHOOK_HANDLE)

  nchanprof = SIZE(chanprof)
  nlayers = profiles(1)%nlayers
  sqrt2pi = SQRT(2 * pi)
  k = 3.3_jprv
  exp_max_exp_exponent_r = EXP(-max_exp_exponent)

  IF (opts%surface%solar_sea_refl_model == solar_refl_model_elfouhaily) THEN
    ! nk_elf = 318 : 0-1 @ 0.1 [10], 1-20 @ 0.5 [38], 20-100 @ 1 [80], 100-2000 @ 10 [190]
    DO i = 0, 10
      karr(i) = 0.1_jprv * i
    ENDDO
    DO i = 1, 38
      karr(i+10) = 1 + 0.5_jprv * i
    ENDDO
    DO i = 1, 80
      karr(i+48) = 20 + 1._jprv * i
    ENDDO
    DO i = 1, 190
      karr(i+128) = 100 + 10._jprv * i
    ENDDO
  ENDIF

  DO j = 1, nchanprof
    prof = chanprof(j)%prof
    IF (profiles(prof)%sunzenangle < 0.0 .OR. &
        profiles(prof)%sunzenangle > max_sol_zen) CYCLE

    DO isurf = 1, profiles(prof)%nsurfaces

      IF (profiles(prof)%skin(isurf)%surftype /= surftype_sea) CYCLE

      IF (.NOT. (solar(j) .AND. emis_refl(isurf)%calc_brdf(j))) CYCLE

      windsp_k      = 0._jprv
      wangl_k       = 0._jprv
      q_shad_k      = 0._jprv
      q_shad_a_k    = 0._jprv
      q_shad_b_k    = 0._jprv
      g_shad_k      = 0._jprv
      p_prime_k     = 0._jprv
      pxy_gammaxy_k = 0._jprv
      fac1_k        = 0._jprv
      gamma_o_k     = 0._jprv
      gammax_k      = 0._jprv
      alfa_k        = 0._jprv
      csi_k         = 0._jprv
      zensat_k      = 0._jprv
      zensun_k      = 0._jprv
      c_shad_k      = 0._jprv
      gamma_p_k     = 0._jprv
      lambda_a_k    = 0._jprv
      lambda_b_k    = 0._jprv
      b_shad_k      = 0._jprv
      a_shad_k      = 0._jprv

      IF (sunglint%s(prof,isurf)%windsp <= min_windsp) THEN

        sunglint_k%s(j,isurf)%glint = 0._jprv

      ELSE

        ! ------------------------------------------------------------------
        ! Yoshimori wave facet reflectance model
        ! ------------------------------------------------------------------

        ! Compute effective distribution function

        q_shad_k      = q_shad_k + &
            sunglint_k%s(j,isurf)%glint * sunglint%s(prof,isurf)%g_shad * sunglint%s(prof,isurf)%p_prime * &
            sunglint%s(prof,isurf)%pxy_gammaxy / sunglint%s(prof,isurf)%fac1
        g_shad_k      = g_shad_k + &
            sunglint_k%s(j,isurf)%glint * sunglint%s(prof,isurf)%q_shad * sunglint%s(prof,isurf)%p_prime * &
            sunglint%s(prof,isurf)%pxy_gammaxy / sunglint%s(prof,isurf)%fac1
        p_prime_k     = p_prime_k + &
            sunglint_k%s(j,isurf)%glint * sunglint%s(prof,isurf)%q_shad * sunglint%s(prof,isurf)%g_shad * &
            sunglint%s(prof,isurf)%pxy_gammaxy / sunglint%s(prof,isurf)%fac1
        pxy_gammaxy_k = pxy_gammaxy_k + &
            sunglint_k%s(j,isurf)%glint * sunglint%s(prof,isurf)%q_shad * sunglint%s(prof,isurf)%g_shad * &
            sunglint%s(prof,isurf)%p_prime / sunglint%s(prof,isurf)%fac1
        fac1_k        = fac1_k - &
            sunglint_k%s(j,isurf)%glint * sunglint%s(prof,isurf)%q_shad * sunglint%s(prof,isurf)%g_shad * &
            sunglint%s(prof,isurf)%p_prime * sunglint%s(prof,isurf)%pxy_gammaxy / sunglint%s(prof,isurf)%fac1**2

        ! Compute the probability density that the slope of the facet at
        ! a certain point is gammax when the incident ray and the ray from
        ! the reflection of the incident ray on the local surface do not
        ! intersect with any other surface.

        gamma_o_k     = gamma_o_k + pxy_gammaxy_k * sunglint%s(prof,isurf)%pxy_gammaxy * &
            (sunglint%s(prof,isurf)%gammax**2 / sunglint%s(prof,isurf)%gamma_o**3 - 1._jprv / sunglint%s(prof,isurf)%gamma_o)
        gammax_k      = gammax_k - &
            pxy_gammaxy_k * sunglint%s(prof,isurf)%pxy_gammaxy * sunglint%s(prof,isurf)%gammax / sunglint%s(prof,isurf)%gamma_o**2
        pxy_gammaxy_k = 0._jprv
        alfa_k        = alfa_k - fac1_k * 2._jprv * sunglint%s(prof,isurf)%fac1 * TAN(sunglint%s(prof,isurf)%alfa)
        csi_k         = csi_k - &
            fac1_k * sunglint%s(prof,isurf)%fac1 * TAN((sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat) / 2._jprv)
        zensat_k      = zensat_k + &
            fac1_k * sunglint%s(prof,isurf)%fac1 * TAN((sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat) / 2._jprv)
        fac1_k        = 0._jprv

        csi_k         = csi_k - g_shad_k * 0.5_jprv * TAN(sunglint%s(prof,isurf)%zensat) / &
            COS(0.5_jprv * (sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat))**2
        zensat_k      = zensat_k + 0.5_jprv * g_shad_k * TAN(sunglint%s(prof,isurf)%zensat) / &
            COS(0.5_jprv * (sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat))**2
        zensat_k      = zensat_k - g_shad_k * TAN(0.5_jprv * (sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat)) / &
            COS(sunglint%s(prof,isurf)%zensat)**2
        g_shad_k      = 0._jprv

        alfa_k        = alfa_k - p_prime_k * sunglint%s(prof,isurf)%p_prime * TAN(sunglint%s(prof,isurf)%alfa) / &
            (COS(sunglint%s(prof,isurf)%alfa) * sunglint%s(prof,isurf)%c_shad)**2
        c_shad_k      = c_shad_k + p_prime_k * sunglint%s(prof,isurf)%p_prime * &
            ((TAN(sunglint%s(prof,isurf)%alfa))**2 / sunglint%s(prof,isurf)%c_shad**3 - 1._jprv / sunglint%s(prof,isurf)%c_shad)
        p_prime_k     = 0._jprv

        IF (COS(sunglint%s(prof,isurf)%csi) + COS(sunglint%s(prof,isurf)%zensat) >= 0._jprv) THEN
          zensat_k  = zensat_k - c_shad_k * sunglint%s(prof,isurf)%gamma_p * SIN(sunglint%s(prof,isurf)%zensat)
          csi_k     = csi_k - c_shad_k * sunglint%s(prof,isurf)%gamma_p * SIN(sunglint%s(prof,isurf)%csi)
          gamma_p_k = gamma_p_k + c_shad_k * (COS(sunglint%s(prof,isurf)%csi) + COS(sunglint%s(prof,isurf)%zensat))
          c_shad_k  = 0._jprv
        ELSE
          zensat_k  = zensat_k + c_shad_k * sunglint%s(prof,isurf)%gamma_p * SIN(sunglint%s(prof,isurf)%zensat)
          csi_k     = csi_k + c_shad_k * sunglint%s(prof,isurf)%gamma_p * SIN(sunglint%s(prof,isurf)%csi)
          gamma_p_k = gamma_p_k - c_shad_k * (COS(sunglint%s(prof,isurf)%csi) + COS(sunglint%s(prof,isurf)%zensat))
          c_shad_k  = 0._jprv
        ENDIF

        ! Direct model - first order shadowing

        IF (ABS(sunglint%s(prof,isurf)%zensat) > 0._jprv) THEN
          IF (1._jprv / TAN(ABS(sunglint%s(prof,isurf)%zensat)) - &
              sunglint%s(prof,isurf)%gammax * sunglint%s(prof,isurf)%zensat / ABS(sunglint%s(prof,isurf)%zensat) >= 0._jprv) THEN
            theta_fi = 1._jprv
          ELSE
            theta_fi = 0._jprv
          ENDIF
        ELSE
          theta_fi = 1._jprv
        ENDIF
        IF (ABS(sunglint%s(prof,isurf)%csi) > 0._jprv) THEN
          IF (1._jprv / TAN(ABS(sunglint%s(prof,isurf)%csi)) + &
              sunglint%s(prof,isurf)%gammax * sunglint%s(prof,isurf)%csi / ABS(sunglint%s(prof,isurf)%csi) >= 0._jprv) THEN
            theta_csi = 1._jprv
          ELSE
            theta_csi = 0._jprv
          ENDIF
        ELSE
          theta_csi = 1._jprv
        ENDIF
        IF (ABS(sunglint%s(prof,isurf)%csi) > pi / 2._jprv) theta_csi = 0._jprv

        ! Second order shadowing (the facet cannot be seen)

        IF ((sunglint%s(prof,isurf)%zensat * sunglint%s(prof,isurf)%csi) <= 0._jprv) THEN
          IF (ABS(sunglint%s(prof,isurf)%zensat) < realtol .AND. ABS(sunglint%s(prof,isurf)%csi) < realtol) THEN
            q_shad_k = 0._jprv
          ELSE IF (sunglint%s(prof,isurf)%zensat >= ABS(sunglint%s(prof,isurf)%csi)) THEN
            q_shad_a_k = q_shad_a_k + q_shad_k * theta_fi
            q_shad_k   = 0._jprv
            lambda_a_k = lambda_a_k - q_shad_a_k / (1._jprv + sunglint%s(prof,isurf)%lambda_a)**2
            q_shad_a_k = 0._jprv
            a_shad_k   = a_shad_k - lambda_a_k * EXP(-sunglint%s(prof,isurf)%a_shad**2 / 2._jprv) / &
                (sqrt2pi * sunglint%s(prof,isurf)%a_shad**2)
            lambda_a_k = 0._jprv
            zensat_k   = zensat_k - a_shad_k / (SIN(sunglint%s(prof,isurf)%zensat)**2 * sunglint%s(prof,isurf)%gamma_o)
            gamma_o_k  = gamma_o_k - a_shad_k / (TAN(sunglint%s(prof,isurf)%zensat) * sunglint%s(prof,isurf)%gamma_o**2)
            a_shad_k   = 0._jprv
          ELSE IF (ABS(sunglint%s(prof,isurf)%csi) > sunglint%s(prof,isurf)%zensat) THEN
            q_shad_b_k = q_shad_b_k + q_shad_k * theta_csi
            q_shad_k   = 0._jprv
            lambda_b_k = lambda_b_k - q_shad_b_k / (1._jprv + sunglint%s(prof,isurf)%lambda_b)**2
            q_shad_b_k = 0._jprv
            b_shad_k   = b_shad_k - lambda_b_k * EXP(-sunglint%s(prof,isurf)%b_shad**2 / 2._jprv) / &
                (sqrt2pi * sunglint%s(prof,isurf)%b_shad**2)
            lambda_b_k = 0._jprv
            IF (sunglint%s(prof,isurf)%csi >= 0._jprv) THEN
              csi_k     = csi_k - b_shad_k / (SIN(sunglint%s(prof,isurf)%csi)**2 * sunglint%s(prof,isurf)%gamma_o)
              gamma_o_k = gamma_o_k - b_shad_k / (TAN(sunglint%s(prof,isurf)%csi) * sunglint%s(prof,isurf)%gamma_o**2)
              b_shad_k  = 0._jprv
            ELSE
              csi_k     = csi_k + b_shad_k  / (SIN(sunglint%s(prof,isurf)%csi)**2 * sunglint%s(prof,isurf)%gamma_o)
              gamma_o_k = gamma_o_k + b_shad_k / (TAN(sunglint%s(prof,isurf)%csi) * sunglint%s(prof,isurf)%gamma_o**2)
              b_shad_k  = 0._jprv
            ENDIF
          ENDIF
        ELSE
          q_shad_k   = q_shad_k * theta_fi * theta_csi
          lambda_a_k = lambda_a_k - q_shad_k / (1._jprv + sunglint%s(prof,isurf)%lambda_a + sunglint%s(prof,isurf)%lambda_b)**2
          lambda_b_k = lambda_b_k - q_shad_k / (1._jprv + sunglint%s(prof,isurf)%lambda_a + sunglint%s(prof,isurf)%lambda_b)**2
          q_shad_k   = 0._jprv
          b_shad_k   = b_shad_k - lambda_b_k * EXP(-sunglint%s(prof,isurf)%b_shad**2 / 2._jprv) / &
              (sqrt2pi * sunglint%s(prof,isurf)%b_shad**2)
          lambda_b_k = 0._jprv
          a_shad_k   = a_shad_k - lambda_a_k * EXP(-sunglint%s(prof,isurf)%a_shad**2 / 2._jprv) / &
              (sqrt2pi * sunglint%s(prof,isurf)%a_shad**2)
          lambda_a_k = 0._jprv
          IF (sunglint%s(prof,isurf)%csi >= 0._jprv) THEN
            csi_k     = csi_k - b_shad_k / (SIN(sunglint%s(prof,isurf)%csi)**2 * sunglint%s(prof,isurf)%gamma_o)
            gamma_o_k = gamma_o_k - b_shad_k / (TAN(sunglint%s(prof,isurf)%csi) * sunglint%s(prof,isurf)%gamma_o**2)
            b_shad_k  = 0._jprv
          ELSE IF (sunglint%s(prof,isurf)%csi < 0._jprv) THEN
            csi_k     = csi_k + b_shad_k / (SIN(sunglint%s(prof,isurf)%csi)**2 * sunglint%s(prof,isurf)%gamma_o)
            gamma_o_k = gamma_o_k + b_shad_k / (TAN(sunglint%s(prof,isurf)%csi) * sunglint%s(prof,isurf)%gamma_o**2)
            b_shad_k  = 0._jprv
          ENDIF
          zensat_k  = zensat_k - a_shad_k / (SIN(sunglint%s(prof,isurf)%zensat)**2 * sunglint%s(prof,isurf)%gamma_o)
          gamma_o_k = gamma_o_k - a_shad_k / (TAN(sunglint%s(prof,isurf)%zensat) * sunglint%s(prof,isurf)%gamma_o**2)
          a_shad_k  = 0._jprv
        ENDIF

        ! First order shadowing (the slope of the facet is negative)

        csi_k    = csi_k + &
            gammax_k / (2._jprv * (COS((sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat) / 2._jprv))**2)
        zensat_k = zensat_k - &
            gammax_k / (2._jprv * (COS((sunglint%s(prof,isurf)%csi - sunglint%s(prof,isurf)%zensat) / 2._jprv))**2)
        gammax_k = 0._jprv

        ! Obtain angles csi and alfa

        IF (sunglint%s(prof,isurf)%csi + sunglint%s(prof,isurf)%zensat >= 0._jprv) THEN
          csi_k    = csi_k + sunglint_k%s(j,isurf)%omega / 2._jprv
          zensat_k = zensat_k + sunglint_k%s(j,isurf)%omega / 2._jprv
        ELSE
          csi_k    = csi_k - sunglint_k%s(j,isurf)%omega / 2._jprv
          zensat_k = zensat_k - sunglint_k%s(j,isurf)%omega / 2._jprv
        ENDIF
        IF (ABS((SIN(sunglint%s(prof,isurf)%zensun) * SIN(pi - sunglint%s(prof,isurf)%dazng * deg2rad))**2 - 1._jprv) > realtol) THEN
          zensun_k = zensun_k + alfa_k * (COS(sunglint%s(prof,isurf)%zensun) * SIN(pi - sunglint%s(prof,isurf)%dazng * deg2rad) / &
              SQRT(1._jprv - (SIN(sunglint%s(prof,isurf)%zensun) * SIN(pi - sunglint%s(prof,isurf)%dazng * deg2rad))**2))
        ENDIF
        alfa_k = 0._jprv
        zensun_k   = zensun_k + &
            csi_k * COS(pi - sunglint%s(prof,isurf)%dazng * deg2rad) / (COS(sunglint%s(prof,isurf)%zensun)**2 * &
            (1._jprv + (TAN(sunglint%s(prof,isurf)%zensun) * COS(pi - sunglint%s(prof,isurf)%dazng * deg2rad))**2))
        csi_k = 0._jprv


        ! ------------------------------------------------------------------
        ! Wave spectrum parameterisation
        ! ------------------------------------------------------------------

        cos2wangl = COS(2 * sunglint%s(prof,isurf)%wangl)
        cos2wangl_k = 0._jprv

        IF (opts%surface%solar_sea_refl_model == solar_refl_model_elfouhaily) THEN

          ! Elfouhaily et al wave spectrum parameterisation

          ! AD init

          am_k = 0._jprv
          cp_k = 0._jprv
          alpha_p_k = 0._jprv
          alpha_m_k = 0._jprv
          omega_k = 0._jprv
          gamma_k = 0._jprv
          sigma_k = 0._jprv
          kp_k = 0._jprv

          ! Direct model

          windsp = MAX(sunglint%s(prof,isurf)%windsp, 0.6_jprv)

          ! Active values independent of wavenumber (k)
          k0 = gravity / windsp**2
          x = k0 * profiles(prof)%near_surface(isurf)%wind_fetch
          kp = k0 * sunglint%s(prof,isurf)%omega_c**2
          kp_r = 1._jprv / kp

          IF (sunglint%s(prof,isurf)%omega_c < 1) THEN
            gamma = 1.7_jprv
          ELSE
            gamma = 1.7_jprv + 6 * LOG10(sunglint%s(prof,isurf)%omega_c)
          ENDIF
          sigma = 0.08_jprv * (1 + 4 / sunglint%s(prof,isurf)%omega_c**3)
          sigma_r2 = 1._jprv / sigma**2
          sigma_r3 = 1._jprv / sigma**3

          cp = SQRT(gravity * (1 + (kp / km_elf)**2) * kp_r)
          cp_r = 1._jprv / cp
          omega = windsp * cp_r
          alpha_p = 6E-3_jprv * SQRT(omega)

          ! Air friction velocity
          ustar = windsp / 25

          ! Angular spread
          am = 0.13_jprv * ustar / cm_elf

          ! AD model

          ! Trapezium integration
          sk2o_k(1:nk_elf-1) = 0.25_jprv * gamma_o_k * (karr(2:nk_elf) - karr(0:nk_elf-2)) / sunglint%s(prof,isurf)%gamma_o
          sk2o_k(0)          = 0.25_jprv * gamma_o_k * (karr(1) - karr(0)) / sunglint%s(prof,isurf)%gamma_o
          sk2o_k(nk_elf)     = 0.25_jprv * gamma_o_k * (karr(nk_elf) - karr(nk_elf-1)) / sunglint%s(prof,isurf)%gamma_o

          sk2p_k(1:nk_elf-1) = 0.25_jprv * gamma_p_k * (karr(2:nk_elf) - karr(0:nk_elf-2)) / sunglint%s(prof,isurf)%gamma_p
          sk2p_k(0)          = 0.25_jprv * gamma_p_k * (karr(1) - karr(0)) / sunglint%s(prof,isurf)%gamma_p
          sk2p_k(nk_elf)     = 0.25_jprv * gamma_p_k * (karr(nk_elf) - karr(nk_elf-1)) / sunglint%s(prof,isurf)%gamma_p

          DO i = 1, nk_elf
            k = karr(i)
            sqrt_k_kp = SQRT(k * kp_r)
            c_r = 1._jprv / sunglint%c(i,prof,isurf)

            sk2_k = 0._jprv
            dk_k = 0._jprv

            IF (sunglint%sk2(i,prof,isurf) > 0._jprv) THEN
              ! Angular spreading factor
              sk2_k = sk2_k + sk2p_k(i) * 0.5_jprv * (1 - 0.5_jprv * sunglint%dk(i,prof,isurf) * cos2wangl)
              dk_k = dk_k - sk2p_k(i) * 0.25_jprv * sunglint%sk2(i,prof,isurf) * cos2wangl
              cos2wangl_k = cos2wangl_k - sk2p_k(i) * 0.25_jprv * sunglint%sk2(i,prof,isurf) * sunglint%dk(i,prof,isurf)

              sk2_k = sk2_k + sk2o_k(i) * 0.5_jprv * (1 + 0.5_jprv * sunglint%dk(i,prof,isurf) * cos2wangl)
              dk_k = dk_k + sk2o_k(i) * 0.25_jprv * sunglint%sk2(i,prof,isurf) * cos2wangl
              cos2wangl_k = cos2wangl_k + sk2o_k(i) * 0.25_jprv * sunglint%sk2(i,prof,isurf) * sunglint%dk(i,prof,isurf)

              cm_c_pow25 = (cm_elf * c_r)**2.5
              c_cp_pow25 = (sunglint%c(i,prof,isurf) * cp_r)**2.5
              cosh_tmp = 1._jprv / COSH(a0_elf + ap_elf * c_cp_pow25 + am * cm_c_pow25)**2
              am_k = am_k + dk_k * cm_c_pow25 * cosh_tmp
              cp_k = cp_k - dk_k *  2.5_jprv * ap_elf * c_cp_pow25 * cp_r * cosh_tmp
            ELSE
              sk2o_k(i) = 0._jprv
              sk2p_k(i) = 0._jprv
            ENDIF

            ! Combined spectrum
            bl_k = sk2_k / k
            bh_k = sk2_k / k

            ! Short-wave spectrum
            alpha_m_k = alpha_m_k + 0.5_jprv * bh_k * cm_elf * sunglint%fm(i,prof,isurf) * c_r

            ! Long-wave spectrum
            fp = sunglint%lpm(i,prof,isurf) * sunglint%jp(i,prof,isurf) * sunglint%fpexp(i,prof,isurf)

            alpha_p_k = alpha_p_k + 0.5_jprv * bl_k * cp * fp * c_r
            cp_k      = cp_k + 0.5_jprv * bl_k * alpha_p * fp * c_r
            fp_k      = 0.5_jprv * bl_k * alpha_p * cp * c_r

            lpm_k   = fp_k * sunglint%jp(i,prof,isurf) * sunglint%fpexp(i,prof,isurf)
            jp_k    = fp_k * sunglint%lpm(i,prof,isurf) * sunglint%fpexp(i,prof,isurf)
            fpexp_k = fp_k * sunglint%lpm(i,prof,isurf) * sunglint%jp(i,prof,isurf)

            omega_k = omega_k - (sunglint%fpexp(i,prof,isurf) / SQRT(10._jprv)) * &
                      (sqrt_k_kp - 1) * fpexp_k
            kp_k    = kp_k + (sunglint%fpexp(i,prof,isurf) / SQRT(10._jprv)) * &
                      0.5_jprv * omega * sqrt_k_kp * fpexp_k * kp_r

            IF (sunglint%gamma_exp(i,prof,isurf) >= exp_max_exp_exponent_r) THEN
              IF (sunglint%gamma_exp(i,prof,isurf) > min_exponent) THEN
                gamma_exp_k = jp_k * sunglint%jp(i,prof,isurf) * LOG(gamma)
                gamma_k = gamma_k + jp_k * sunglint%jp(i,prof,isurf) * &
                          sunglint%gamma_exp(i,prof,isurf) / gamma
              ELSE
                jp_k = 0._jprv
              ENDIF

              kp_k    = kp_k + gamma_exp_k * sunglint%gamma_exp(i,prof,isurf) * &
                          (sqrt_k_kp - 1) * 0.5_jprv * sqrt_k_kp * kp_r * sigma_r2
              sigma_k = sigma_k + gamma_exp_k * sunglint%gamma_exp(i,prof,isurf) * &
                          (sqrt_k_kp - 1)**2 * sigma_r3
            ELSE
              jp_k = 0._jprv
            ENDIF

            kp_k = kp_k - 5 * 2 * lpm_k * kp * sunglint%lpm(i,prof,isurf) / (4 * k**2)
          ENDDO

          ! Angular spread
          ustar_k = 0.13_jprv * am_k / cm_elf

          ! Air friction velocity
          IF (ustar < cm_elf) THEN
            ustar_k = ustar_k + 1E-2_jprv * alpha_m_k / ustar
          ELSE
            ustar_k = ustar_k + 1E-2_jprv * 3 * alpha_m_k / ustar
          ENDIF
          windsp_k = ustar_k / 25

          omega_k = omega_k + 6E-3_jprv * 0.5_jprv * alpha_p_k / SQRT(omega)

          windsp_k = windsp_k + omega_k * cp_r
          cp_k = cp_k - windsp * omega_k * cp_r**2

          kp_k = kp_k + 0.5_jprv  * gravity * cp_k * (1 / km_elf**2 - kp_r**2) * cp_r

          omega_c_k = -0.08_jprv * 4 * 3 * sigma_k / sunglint%s(prof,isurf)%omega_c**4
          IF (sunglint%s(prof,isurf)%omega_c < 1) THEN
            gamma_k = 0._jprv
          ELSE
            ! log10(x) = log(x) / log(10)
            omega_c_k = omega_c_k + 6 * gamma_k / (sunglint%s(prof,isurf)%omega_c * LOG(10._jprv))
          ENDIF

          k0_k = kp_k * sunglint%s(prof,isurf)%omega_c**2
          omega_c_k = omega_c_k + 2 * k0 * kp_k * sunglint%s(prof,isurf)%omega_c

          x_k = 0.84_jprv * (-0.75_jprv) * (TANH((x / x0_elf)**0.4))**(-1.75) * &
                (1._jprv / COSH((x / x0_elf)**0.4)**2) * &
                (0.4_jprv * omega_c_k / (x**0.6_jprv * x0_elf**0.4_jprv))

          k0_k = k0_k + x_k * profiles(prof)%near_surface(isurf)%wind_fetch
          profiles_k(j)%near_surface(isurf)%wind_fetch = profiles_k(j)%near_surface(isurf)%wind_fetch + x_k * k0

          windsp_k = windsp_k - 2 * gravity * k0_k / windsp**3

          IF (sunglint%s(prof,isurf)%windsp < 0.6_jprv) windsp_k = 0._jprv

        ENDIF ! wave spectrum option

        wangl_k = wangl_k - 2 * cos2wangl_k * SIN(2 * sunglint%s(prof,isurf)%wangl)


        ! --------------------------------------------------------------------
        ! Calculate profile variables
        ! --------------------------------------------------------------------

        ! Satellite and solar zenith angles at surface
        IF (ABS(raytracing%pathsun(nlayers,prof) - 1._jprv) > realtol) THEN
          raytracing_k%pathsun(nlayers,j) = raytracing_k%pathsun(nlayers,j) + zensun_k / &
              (raytracing%pathsun(nlayers,prof) * SQRT(raytracing%pathsun(nlayers,prof)**2 - 1._jprv))
        ENDIF
        zensun_k = 0._jprv
        IF (ABS(raytracing%pathsat(nlayers,prof) - 1._jprv) > realtol) THEN
          raytracing_k%pathsat(nlayers,j) = raytracing_k%pathsat(nlayers,j) + zensat_k / &
              (raytracing%pathsat(nlayers,prof) * SQRT(raytracing%pathsat(nlayers,prof)**2 - 1._jprv))
        ENDIF
        zensat_k = 0._jprv

        ! Compute the angle between the wind direction and the U axis

        IF (sunglint%s(prof,isurf)%windsp > min_windsp) THEN
          profiles_k(j)%near_surface(isurf)%wind_v10m = profiles_k(j)%near_surface(isurf)%wind_v10m +      &
              profiles(prof)%near_surface(isurf)%wind_u10m * wangl_k / &
              (profiles(prof)%near_surface(isurf)%wind_u10m**2 + profiles(prof)%near_surface(isurf)%wind_v10m**2)
          profiles_k(j)%near_surface(isurf)%wind_u10m = profiles_k(j)%near_surface(isurf)%wind_u10m -      &
              profiles(prof)%near_surface(isurf)%wind_v10m * wangl_k / &
              (profiles(prof)%near_surface(isurf)%wind_u10m**2 + profiles(prof)%near_surface(isurf)%wind_v10m**2)
        ELSE
          wangl_k = 0._jprv
        ENDIF

        ! Compute the wind speed

        IF (sunglint%s(prof,isurf)%windsp > min_windsp) THEN
          profiles_k(j)%near_surface(isurf)%wind_u10m = profiles_k(j)%near_surface(isurf)%wind_u10m + &
              windsp_k * profiles(prof)%near_surface(isurf)%wind_u10m / sunglint%s(prof,isurf)%windsp
          profiles_k(j)%near_surface(isurf)%wind_v10m = profiles_k(j)%near_surface(isurf)%wind_v10m + &
              windsp_k * profiles(prof)%near_surface(isurf)%wind_v10m / sunglint%s(prof,isurf)%windsp
          windsp_k = 0._jprv
        ELSE
          profiles_k(j)%near_surface(isurf)%wind_u10m = profiles_k(j)%near_surface(isurf)%wind_u10m + windsp_k / SQRT(2._jprv)
          profiles_k(j)%near_surface(isurf)%wind_v10m = profiles_k(j)%near_surface(isurf)%wind_v10m + windsp_k / SQRT(2._jprv)
          windsp_k = 0._jprv
        ENDIF

      ENDIF ! wind speed > minimum
    ENDDO ! nsurfaces
  ENDDO ! nchanprof
  IF (LHOOK) CALL DR_HOOK('RTTOV_CALC_SUNGLINT_K', 1_jpim, ZHOOK_HANDLE)
END SUBROUTINE rttov_calc_sunglint_k
