! Description:
!> @file
!!   Subroutines needed by MFASIS-NN for computation of cloud NN input parameters
!!   from model profiles
!
!> @brief
!!   Subroutines needed by MFASIS-NN for computation of cloud NN input parameters
!!   from model profiles
!!
!! @details
!! ...
!
! 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_mfasis_nn_hydro_inpar_ad( &
    err,                &
    nlayers,            &
    prof,               &
    profad,             &
    cloud_columns,      &
    ist,                &
    nch,                &
    cc,                 &
    opdpext,            &
    opdpext_ad,         &
    mfasisnn_coefs,     &
    profiles,           &
    profiles_ad,        &
    aux,                &
    aux_ad,             &
    nn_inpar_idx,       &
    pressure_gradients, &
    sza,                &
    vza,                &
    tauw,               &
    taui,               &
    tauwi_top,          &
    tauwi_bot,          &
    reffw_top,          &
    reffi_top,          &
    reffw_bot,          &
    reffi_bot,          &
    reffwi_top,         &
    reffwi_bot,         &
    psfc,               &
    tauw_top,           & ! TL/AD version, previously optional output in direct subroutine
    tauw_bot,           &
    taui_top,           &
    taui_bot,           &
    tauwi,              &
    logtauw,            &
    tauw_thr,           &
    logtaui,            &
    taui_thr,           &
    dtauw,              &
    dtaui,              &
    dreffw,             &
    dreffi,             &
    fmp,                &
    fbotw,              &
    fboti,              &
!   fbot_X,             &
    tauw_ad,         &  ! TL/AD output
    taui_ad,         &
    tauwi_top_ad,    &
    tauwi_bot_ad,    &
    reffw_top_ad,    &
    reffi_top_ad,    &
    reffw_bot_ad,    &
    reffi_bot_ad,    &
    reffwi_top_ad,   &  
    reffwi_bot_ad,   &
    fpct_ad,         &
!   psfc_in_ad,         &
    psfc_ad,          &
    l_iwvc,        &
    fp_c_ad,            &
    iwvn_c_ad          &
    )

!INTF_OFF
#include "throw.h"
!INTF_ON
  USE rttov_kinds, ONLY : jpim, jprv, jplm

  USE rttov_types, ONLY :          &
      rttov_cloud_columns,         &
      rttov_coef_mfasis_nn,        &
      rttov_profile_internal,      &
      rttov_profile_aux
!INTF_OFF

  USE rttov_const, ONLY :          &
      wcl_opac_deff,               &
      pi,                          &
      nopac_wcl,                   &
      clw_deff_index,              &
      ice_baum_index

  USE rttov_const, ONLY :        &
!X  pi,                          &
!X  deg2rad,                     &
!X  rad2deg,                     &
!X  pi_r,                        &
!X  mfasis_file_type_hydro,      &
!X  mfasis_file_type_aer,        &
!X  qflag_mfasis_nn_limits,      &
    gas_id_watervapour,          &
    gas_mass,                    &
    mair,                        &
    gravity


  USE rttov_mfasis_nn_mod, ONLY :  &
      comp_fbot_ad, iwv2iwvn_ad,   &
      comp_fbot   ,                &
      nn_idx_reffw_top, nn_idx_reffi_top

  USE yomhook, ONLY : lhook, dr_hook, jphook
!INTF_ON
  IMPLICIT NONE

  INTEGER(jpim),                  INTENT(OUT)   :: err
  INTEGER(jpim),                  INTENT(IN)    :: nlayers
  INTEGER(jpim),                  INTENT(IN)    :: prof
  INTEGER(jpim),                  INTENT(IN)    :: profad
  TYPE(rttov_cloud_columns),      INTENT(IN)    :: cloud_columns
  INTEGER(jpim),                  INTENT(IN)    :: ist
  INTEGER(jpim),                  INTENT(IN)    :: nch
  INTEGER(jpim),                  INTENT(IN)    :: cc
  REAL(jprv),                     INTENT(IN)    :: opdpext(:,:,:,:)
  REAL(jprv),                     INTENT(INOUT) :: opdpext_ad(:,:,:,:)
  TYPE(rttov_coef_mfasis_nn),     INTENT(IN)    :: mfasisnn_coefs
  TYPE(rttov_profile_internal),   INTENT(IN)    :: profiles(:)
  TYPE(rttov_profile_internal),   INTENT(INOUT) :: profiles_ad(:)
  TYPE(rttov_profile_aux),        INTENT(IN)    :: aux
  TYPE(rttov_profile_aux),        INTENT(INOUT) :: aux_ad
  INTEGER(jpim),                  INTENT(IN)    :: nn_inpar_idx(:,:)
  LOGICAL(jplm),                  INTENT(IN)    :: pressure_gradients
  REAL(jprv),                     INTENT(IN)    :: sza
  REAL(jprv),                     INTENT(IN)    :: vza
  REAL(jprv),                     INTENT(IN)    :: tauw
  REAL(jprv),                     INTENT(IN)    :: taui
  REAL(jprv),                     INTENT(IN)    :: tauwi_top
  REAL(jprv),                     INTENT(IN)    :: tauwi_bot
  REAL(jprv),                     INTENT(IN)    :: reffw_top
  REAL(jprv),                     INTENT(IN)    :: reffi_top
  REAL(jprv),                     INTENT(IN)    :: reffw_bot
  REAL(jprv),                     INTENT(IN)    :: reffi_bot
  REAL(jprv),                     INTENT(IN)    :: reffwi_top
  REAL(jprv),                     INTENT(IN)    :: reffwi_bot
  REAL(jprv),                     INTENT(IN)    :: psfc
  REAL(jprv),                     INTENT(IN)    :: tauw_top  ! TL/AD version
  REAL(jprv),                     INTENT(IN)    :: tauw_bot
  REAL(jprv),                     INTENT(IN)    :: taui_top
  REAL(jprv),                     INTENT(IN)    :: taui_bot
  REAL(jprv),                     INTENT(IN)    :: tauwi
  REAL(jprv),                     INTENT(IN)    :: tauw_thr
  REAL(jprv),                     INTENT(IN)    :: taui_thr
  REAL(jprv),                     INTENT(IN)    :: logtauw
  REAL(jprv),                     INTENT(IN)    :: logtaui
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: dtauw
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: dtaui
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: dreffw
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: dreffi
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: fmp
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: fbotw
  REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: fboti
! REAL(jprv), DIMENSION(nlayers), INTENT(IN)    :: fbot_X
  REAL(jprv),                     INTENT(INOUT) :: tauw_ad
  REAL(jprv),                     INTENT(INOUT) :: taui_ad
  REAL(jprv),                     INTENT(INOUT) :: tauwi_top_ad
  REAL(jprv),                     INTENT(INOUT) :: tauwi_bot_ad
  REAL(jprv),                     INTENT(INOUT) :: reffw_top_ad
  REAL(jprv),                     INTENT(INOUT) :: reffi_top_ad
  REAL(jprv),                     INTENT(INOUT) :: reffw_bot_ad
  REAL(jprv),                     INTENT(INOUT) :: reffi_bot_ad
  REAL(jprv),                     INTENT(INOUT) :: reffwi_top_ad
  REAL(jprv),                     INTENT(INOUT) :: reffwi_bot_ad
  REAL(jprv),                     INTENT(INOUT) :: fpct_ad
! REAL(jprv),                     INTENT(IN)    :: psfc_in_ad
  REAL(jprv),                     INTENT(INOUT) :: psfc_ad
  INTEGER(jpim),        OPTIONAL, INTENT(IN)  :: l_iwvc    ! >0: compute fp_c/iwvn_c
  REAL(jprv),           OPTIONAL, INTENT(IN)  :: fp_c_ad(1:6) ! pressure ratio between a level and the level below
  REAL(jprv),           OPTIONAL, INTENT(IN)  :: iwvn_c_ad(1:7) ! WV integrated  between a level and the level below

!INTF_END

  REAL(jprv)                     :: fpct
  REAL(jprv)                     :: fac_hi, fac_lo, tau_lo, tau_hi, tau, tau_thr
  REAL(jprv)                     :: tau_thr_ad, tau_thr_v1, tau_ad
  REAL(jprv)                     :: tauw_thr_ad, taui_thr_ad, logtauw_ad, logtaui_ad
  REAL(jprv)                     :: tauw_top_ad, tauw_bot_ad, taui_top_ad, taui_bot_ad, tauwi_ad
  REAL(jprv), DIMENSION(nlayers) :: dtauw_ad, dtaui_ad, fmp_ad, dreffw_ad, dreffi_ad, fbotw_ad, fboti_ad, tauv1
  REAL(jprv), DIMENSION(nlayers) :: ddtaui_ad
  INTEGER(jpim)                  :: j, ii, coli

  REAL(jprv)                     :: min_p_X, max_p_X, range_p_X, qq
  REAL(jprv)                     :: dtau_help_size_X, rel_tauw_thr_X, rel_taui_thr_X
  REAL(jprv)                     :: min_taui_thr_X, min_tauw_thr_X
  REAL(jprv)                     :: ddthr(6), ddccc(6), tau_thr_c(6)
  REAL(jprv)                     :: dd_thr, dd_ccc
  REAL(jprv)                     :: fk_pc(6), pc(0:7)
  REAL(jprv)                     :: iwv_c(0:7) ! WV integrated  between a level and the level below
  REAL(jprv), DIMENSION(nlayers) :: dp_X, diwv_X, dtaubins_X, fbot_X
! REAL(jprv), DIMENSION(nlayers) :: dp_X, diwv_X, dtaubins_X
  REAL(jprv), DIMENSION(nlayers) :: q_mxr

  REAL(jprv)                     :: range_p_ad
  REAL(jprv)                     :: dtau_help_size_ad
  REAL(jprv)                     :: ddthr_ad(6), tau_thr_c_ad(6)
  REAL(jprv)                     :: dd_thr_ad, dd_ccc_ad
  REAL(jprv)                     :: fk_pc_ad(6), pc_ad(0:7), pca_ad, pcb_ad
  REAL(jprv)                     :: iwv_c_ad(0:7) ! WV integrated  between a level and the level below
  REAL(jprv), DIMENSION(nlayers) :: dp_ad, diwv_ad, dtaubins_ad, fbot_ad
  REAL(jprv), DIMENSION(nlayers) :: q_mxr_ad


  REAL(jphook) :: zhook_handle
  !- End of header --------------------------------------------------------

  TRY

  IF (lhook) CALL dr_hook('RTTOV_CALC_MFASIS_NN_HYDRO_INPAR_AD',0_jpim,zhook_handle)

  ! --------------------------------------------------------------------------
  ! Computation of NN input parameters from model profiles:
  ! Optical depths, effective radii, ...
  ! --------------------------------------------------------------------------
  tauwi_ad     = 0

  tau_thr_ad   = 0
  tauw_top_ad  = 0
  taui_top_ad  = 0
  tau_ad       = 0
  tauw_thr_ad  = 0
  taui_thr_ad  = 0
  logtauw_ad   = 0
  logtaui_ad   = 0
  tauw_bot_ad  = 0
  taui_bot_ad  = 0

  dtauw_ad  = 0
  dtaui_ad  = 0
  ddtaui_ad = 0
  fmp_ad    = 0
  dreffw_ad = 0
  dreffi_ad = 0
  fbotw_ad  = 0
  fboti_ad  = 0
  tauv1     = 0



!       tauw_a_ad(ncl) = tauw_ad
!       taui_a_ad(ncl) = taui_ad
!       tauwi_top_a_ad(ncl) = tauwi_top_ad
!       tauwi_bot_a_ad(ncl) = tauwi_bot_ad
!       reffw_top_a_ad(ncl) = reffw_top_ad
!       reffi_top_a_ad(ncl) = reffi_top_ad
!       reffw_bot_a_ad(ncl) = reffw_bot_ad
!       reffi_bot_a_ad(ncl) = reffi_bot_ad
!       reffwi_top_a_ad(ncl) = reffwi_top_ad
!       reffwi_bot_a_ad(ncl) = reffwi_bot_ad
!       psfc_a_ad(ncl) = psfc_ad
!       fpct_a_ad(ncl) = fpct_ad

! psfc_ad       = psfc_ad       + psfc_in_ad

!-------------------------------------------------------------------------------------------------
! threshold optical depth for detecting cloud top
!            tau_thr_v1 = (tauw + taui + tauwi)/2._jprv
!            tau_thr    = tau_thr_v1/(tau_thr_v1 + 1._jprv)
!
!            tau_thr_ad = (tauw_ad + taui_ad + tauwi_ad)/2._jprv
!            tau_thr_ad = tau_thr_ad/(tau_thr_v1 + 1._jprv)  -  &
!                         tau_thr_ad * tau_thr_v1/(tau_thr_v1 + 1._jprv)**2
!
!            tau = 0._jprv
!            tau_ad = 0._jprv
!            DO j = 1, nlay
!               tau = tau + dtauw(j) + dtaui(j)
!               tau_ad = tau_ad + dtauw_ad(j) + dtaui_ad(j)
!               IF ( tau > tau_thr) THEN
!                  fpct    = profiles   (prof)%p_half(j+1) -   (tau    - tau_thr)/(dtauw(j) + dtaui(j))   &
!                                                * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j))
!
!                  fpct_ad = profiles_ad(profad)%p_half(j+1) +( -(tau_ad - tau_thr_ad)/(dtauw(j) + dtaui(j)) &
!                                                * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j)) &
!                                                          +(tau    - tau_thr)/(dtauw(j) + dtaui(j))**2 &
!                                             * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j)) &
!                                                                            * (dtauw_ad(j) + dtaui_ad(j)) &
!                                                          -(tau    - tau_thr)/(dtauw(j) + dtaui(j))    &
!                                 * (profiles_ad(profad)%p_half(j+1) - profiles_ad(profad)%p_half(j)) &
!                                                        )
!                  fpct = fpct / psfc
!                  fpct_ad = fpct_ad / psfc - fpct * psfc0_ad/psfc
!                  EXIT
!               ENDIF
!            ENDDO
!=================================================================================================
  IF(l_iwvc==0) THEN
!=================================================================================================
    tau_thr_v1 = (tauw + taui + tauwi)/2._jprv
    tau_thr    = tau_thr_v1/(tau_thr_v1 + 1._jprv)
    tau = 0._jprv
    tauv1(:) = - 9999999._jprv
    DO j = 1, nlayers
      tau = tau + dtauw(j) + dtaui(j)
      tauv1(j)=tau
      IF(tauv1(j) > tau_thr) THEN
        EXIT
      ENDIF
    ENDDO
    tau_ad = 0
    DO j = nlayers, 1, -1
      IF ( tauv1(j) > tau_thr) THEN
        fpct    = profiles   (prof)%p_half(j+1) -   (tauv1(j)-tau_thr)/(dtauw(j) + dtaui(j))    &
                                       * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j))
        fpct = fpct / psfc

!                  fpct_ad = fpct_ad / psfc - fpct * psfc0_ad/psfc
        psfc_ad = psfc_ad - fpct_ad * fpct/psfc
        fpct_ad  = fpct_ad / psfc

!                  fpct_ad = profiles_ad(profad)%p_half(j+1) +( -(tau_ad - tau_thr_ad)/(dtauw(j) + dtaui(j)) &
!                                                * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j)) &
!                                                          +(tau    - tau_thr)/(dtauw(j) + dtaui(j))**2 &
!                                                * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j)) &
!                                                                            * (dtauw_ad(j) + dtaui_ad(j)) &
!                                                          -(tau    - tau_thr)/(dtauw(j) + dtaui(j))    &
!                                        * (profiles_ad(profad)%p_half(j+1) - profiles_ad(profad)%p_half(j)) &

        tau_ad     =    tau_ad  - fpct_ad/(dtauw(j)+dtaui(j)) &
                                               * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j))
        tau_thr_ad = tau_thr_ad + fpct_ad/(dtauw(j)+dtaui(j)) &
                                               * (profiles(prof)%p_half(j+1) - profiles(prof)%p_half(j))
        dtauw_ad(j)= dtauw_ad(j)+ fpct_ad * (tauv1(j)-tau_thr)/(dtauw(j)+ dtaui(j))**2 &
                                               * (profiles(prof)%p_half(j+1)-profiles(prof)%p_half(j))
        dtaui_ad(j)= dtaui_ad(j)+ fpct_ad * (tauv1(j)-tau_thr)/(dtauw(j)+ dtaui(j))**2 &
                                               * (profiles(prof)%p_half(j+1)-profiles(prof)%p_half(j))
        IF (pressure_gradients) THEN
          profiles_ad(profad)%p_half(j+1) = profiles_ad(profad)%p_half(j+1) - &
              fpct_ad * (tauv1(j)-tau_thr)/(dtauw(j)+dtaui(j)) + fpct_ad
          profiles_ad(profad)%p_half(j  ) = profiles_ad(profad)%p_half(j  ) + &
              fpct_ad * (tauv1(j)-tau_thr)/(dtauw(j)+dtaui(j))
        ENDIF
        fpct_ad = 0
      ENDIF
!              tau_ad = tau_ad + dtauw_ad(j) + dtaui_ad(j)
      dtauw_ad(j) = dtauw_ad(j) + tau_ad
      dtaui_ad(j) = dtaui_ad(j) + tau_ad
    ENDDO
!   tau_thr_ad = (tauw_ad + taui_ad + tauwi_ad)/2._jprv
!   tau_thr_ad = tau_thr_ad/(tau_thr_v1 + 1._jprv)  -  &
!                         tau_thr_ad * tau_thr_v1/(tau_thr_v1 + 1._jprv)**2
    tau_thr_ad = tau_thr_ad/(tau_thr_v1 + 1._jprv)  -  &
                      tau_thr_ad * tau_thr_v1/(tau_thr_v1 + 1._jprv)**2
    tauw_ad    = tauw_ad + tau_thr_ad/2._jprv
    taui_ad    = taui_ad + tau_thr_ad/2._jprv
    tauwi_ad   = tauwi_ad+ tau_thr_ad/2._jprv
    tau_thr_ad = 0
!=================================================================================================
  ELSEIF(l_iwvc>0) THEN
!=================================================================================================
    min_taui_thr_X = 0.2_jprv
    min_tauw_thr_X = 0.2_jprv
    rel_taui_thr_X = 0.2_jprv
    rel_tauw_thr_X = 0.2_jprv


    q_mxr(:) = profiles(prof)%q(:) * gas_mass(gas_id_watervapour) / &
                         (mair * 1.E06_jprv + profiles(prof)%q(:) * gas_mass(gas_id_watervapour))

    dp_X  (:) = profiles(prof)%p_half(2:nlayers+1) - profiles(prof)%p_half(1:nlayers)
    diwv_X(:) = 100._jprv * dp_X *  q_mxr              /gravity
!-------------------------------------------------------------------------------------------------
    min_p_X = 100._jprv
    max_p_X = psfc
    range_p_X = max_p_X - min_p_X

!   dtau_help_size_X = 0.1d0 / prof_pres(i, nz + 1)  ! max value set to 0.01
    dtau_help_size_X = 0.1_jprv / profiles(prof)%p_half(nlayers+1)  ! max value set to 0.01
    dtaubins_X(:) = dtauw(:) + dtaui(:) + dtau_help_size_X * dp_X(:)

    tau_thr_c(6) = 0
    tau_thr_c(5) = taui_top
    tau_thr_c(4) = taui_top + taui_bot

    tau_thr_c(3) = taui_top + taui_bot
    tau_thr_c(2) = taui_top + taui_bot + tauw_top + tauwi_top
    tau_thr_c(1) = taui_top + taui_bot + tauw_top + tauw_bot &
                                + tauwi_top + tauwi_bot
!======================

    ddthr(6) =   rel_taui_thr_X * taui_top
    ddthr(4) =   rel_taui_thr_X * taui_bot
    ddthr(3) =   rel_tauw_thr_X *(tauw_top+tauwi_top)
    ddthr(1) =   rel_tauw_thr_X *(tauw_bot+tauwi_bot)

    ddccc(:) = 0
    ddccc(6) =   min_taui_thr_X  ! ct_i
    ddccc(4) = - min_taui_thr_X  ! cb_i
    ddccc(3) =   min_tauw_thr_X  ! ct_w
    ddccc(1) = - min_tauw_thr_X  ! cb_w
!======================
    do ii=1, 6
      if(ddccc(ii)==0) cycle
      dd_thr = ddthr(ii) * (abs(ddccc(ii)) + ddthr(ii))
      dd_ccc = ddccc(ii)*dd_thr/(ddccc(ii)**2 + dd_thr)
      tau_thr_c(ii) = tau_thr_c(ii) + dd_ccc
    enddo
!=====================               

    do ii=1, 6
      fk_pc(ii)     = (6._jprv -real(ii))*range_p_X / 6._jprv + min_p_X
      tau_thr_c(ii)  = tau_thr_c(ii) + fk_pc(ii) * dtau_help_size_X
    enddo
!=====================

    do ii=1, 6
      call comp_fbot(fbot_X, tau_thr_c(ii), dtaubins_X, nlayers, &
                     mfasisnn_coefs%smoothing_parameter)

      iwv_c(ii) = sum(       fbot_X (:)   * diwv_X(:))        ! WV below level ii
      pc(ii)    = sum((1._jprv-fbot_X(:)) * dp_X  (:))  ! pres.diff btw level ii and toplevel
    enddo

    iwv_c   (7)    = sum(diwv_X(:))
    iwv_c   (0)    = 0
    iwv_c   (1:7)  = iwv_c   (1:7) - iwv_c   (0:6)           ! WV btw level ii and level below

    pc   (7) = profiles   (prof)%p_half(1)
    pc   (0) = psfc

!=====================
    diwv_ad(:)=0
    iwv_c_ad(:)    = 0
    fbot_ad(:)=0 
    dp_ad(:)  =0
    dtau_help_size_ad =0
    fk_pc_ad(:)      =0 
    dd_ccc_ad  =0 
    dd_thr_ad    = 0
    ddthr_ad(:) = 0
    q_mxr_ad(:) = 0
    range_p_ad  = 0
    
    dtaubins_ad = 0
    tau_thr_c_ad=0
    
    pc_ad(:)        =0

!=====================
!=====================
!   pc_tl(7) = profiles_tl(prof)%p_half(1)
!   do ii=1, 7
!     iwvn_c_tl(ii) = iwv2iwvn_tl( pc   (ii-1), pc   (ii), iwv_c   (ii)      &
!                                  ,pc_tl(ii-1), pc_tl(ii), iwv_c_tl(ii))
!   enddo
!=====================
    do ii=1, 7
      call iwv2iwvn_ad (iwvn_c_ad(ii),  & ! in
                         pc   (ii-1), pc   (ii), iwv_c   (ii),      & ! in
                         pca_ad      , pcb_ad  , iwv_c_ad(ii))        ! out
!                        pc_ad(ii-1), pc_ad(ii), iwv_c_ad(ii))        ! out
      pc_ad(ii-1)=pc_ad(ii-1)+pca_ad
      pc_ad(ii)  =pc_ad(ii)  +pcb_ad
    enddo

!   pc_tl(7) = profiles_tl(prof)%p_half(1)
    profiles_ad(profad)%p_half(1) = profiles_ad(profad)%p_half(1) + pc_ad(7)
    pc_ad(7) = 0

!=====================
!   iwv_c_tl(7)    = sum(diwv_tl(:))
!   iwv_c_tl(0)    = 0
!   iwv_c_tl(1:7)  = iwv_c_tl(1:7) - iwv_c_tl(0:6)           ! WV between level ii and level below
!=====================
    iwv_c_ad(0:6)  = iwv_c_ad(0:6) - iwv_c_ad(1:7)
    iwv_c_ad(0)    = 0
    diwv_ad(:)     = diwv_ad(:) + iwv_c_ad(7)
    iwv_c_ad(7) = 0
               

!=====================
!   pc_tl(0) = psfc_tl
!   do ii=1, 6
!     fp_c_tl(ii) = (pc_tl(ii)*pc(ii-1)-pc(ii)*pc_tl(ii-1))/pc(ii-1)**2 ! ratio between level ii and level below
!   enddo
!=====================
    do ii=6, 1, -1
!     fp_c_ad(ii) = (pc_ad(ii)*pc(ii-1)-pc(ii)*pc_ad(ii-1))/pc(ii-1)**2 ! ratio between level ii and level below
      pc_ad(ii)   = pc_ad(ii)   +         fp_c_ad(ii)/pc(ii-1)
      pc_ad(ii-1) = pc_ad(ii-1) - pc(ii)* fp_c_ad(ii)/pc(ii-1)**2
    enddo
    psfc_ad  = psfc_ad + pc_ad(0) 
    pc_ad(0) = 0
!=====================
!   do ii=1, 6
!     call comp_fbot(fbot_X, tau_thr_c(ii), dtaubins_X, nlayers, &
!                    mfasisnn_coefs%smoothing_parameter)
!     call comp_fbot_tl(fbot_tl, tau_thr_c_tl(ii), dtaubins_tl,tau_thr_c(ii), dtaubins_X, nlayers, &
!                       mfasisnn_coefs%smoothing_parameter)

!     iwv_c(ii) = sum(       fbot_X (:)   * diwv_X(:))        ! WV below level ii
!     pc(ii)    = sum((1._jprv-fbot_X(:)) * dp_X  (:))  ! pres.diff between level ii and toplevel

!     iwv_c_tl(ii) = sum(    fbot_tl(:)   * diwv_X(:) +  fbot_X(:)  * diwv_tl(:))        ! WV below level ii
!     pc_tl(ii)    = sum((-fbot_tl(:))    * dp_X(:)   + (1._jprv-fbot_X(:)) * dp_tl(:))  ! pres.diff btw lev ii and toplevel
!   enddo
!=====================
    do ii=1, 6
      call comp_fbot(fbot_X, tau_thr_c(ii), dtaubins_X, nlayers, &
                     mfasisnn_coefs%smoothing_parameter)
      iwv_c(ii) = sum(       fbot_X (:)   * diwv_X(:))        ! WV below level ii
      pc(ii)    = sum((1._jprv-fbot_X(:)) * dp_X  (:))  ! pres.diff between level ii and toplevel

!     pc_tl(ii)    = sum((-fbot_tl(:))    * dp_X(:)   + (1._jprv-fbot_X(:)) * dp_tl(:))  ! pres.diff btw lev ii and toplevel
      fbot_ad(:) = fbot_ad(:) -   dp_X(:)            * pc_ad(ii)
      dp_ad(:)   = dp_ad(:)   +  (1._jprv-fbot_X(:)) * pc_ad(ii)
      pc_ad(ii)  = 0

!     iwv_c_tl(ii) = sum(    fbot_tl(:)   * diwv_X(:) +  fbot_X(:)  * diwv_tl(:))        ! WV below level ii
      fbot_ad(:) = fbot_ad(:) +  diwv_X(:) * iwv_c_ad(ii)
      diwv_ad(:) = diwv_ad(:) +  fbot_X(:) * iwv_c_ad(ii)
      iwv_c_ad(ii)= 0

!     call comp_fbot_tl(fbot_tl, tau_thr_c_tl(ii), dtaubins_tl,tau_thr_c(ii), dtaubins_X, nlayers, &
!                       mfasisnn_coefs%smoothing_parameter)
      call comp_fbot_ad(fbot_ad, tau_thr_c_ad(ii), dtaubins_ad,tau_thr_c(ii), dtaubins_X, nlayers, &
                        mfasisnn_coefs%smoothing_parameter)

    enddo

!=====================
!   do ii=1, 6
!     fk_pc(ii)     = (6._jprv -real(ii))*range_p_X / 6._jprv + min_p_X
!     tau_thr_c(ii)  = tau_thr_c(ii) + fk_pc(ii) * dtau_help_size_X

!     fk_pc_tl(ii)  = (6._jprv -real(ii))*range_p_tl/ 6._jprv
!     tau_thr_c_tl(ii)= tau_thr_c_tl(ii) + fk_pc(ii)   * dtau_help_size_tl  &
!                                                   + fk_pc_tl(ii)* dtau_help_size_X
!   enddo
!=====================
    do ii=1, 6

!     tau_thr_c_tl(ii)= tau_thr_c_tl(ii) + fk_pc(ii)   * dtau_help_size_tl  &
!                                                   + fk_pc_tl(ii)* dtau_help_size_X
      dtau_help_size_ad = dtau_help_size_ad + fk_pc(ii)        * tau_thr_c_ad(ii)
      fk_pc_ad(ii)      = fk_pc_ad(ii)      + dtau_help_size_X * tau_thr_c_ad(ii)

!     fk_pc_tl(ii)  = (6._jprv -real(ii))*range_p_tl/ 6._jprv
      range_p_ad  = range_p_ad  + (6._jprv -real(ii))/ 6._jprv *  fk_pc_ad(ii)
      fk_pc_ad(ii) = 0
    enddo

!=====================
!   do ii=1, 6
!     if(ddccc(ii)==0) cycle
!     dd_thr = ddthr(ii) * (abs(ddccc(ii)) + ddthr(ii))
!     dd_thr_tl = ddthr_tl(ii) * (abs(ddccc(ii)) + 2._jprv* ddthr(ii))
!     qq = ddccc(ii)**2 + dd_thr
!     dd_ccc_tl = ddccc(ii)*dd_thr_tl* (qq - dd_thr)/qq**2
!     tau_thr_c_tl(ii) = tau_thr_c_tl(ii) + dd_ccc_tl
!   enddo
!=====================
    do ii=1, 6
      if(ddccc(ii)==0) cycle
      dd_thr = ddthr(ii) * (abs(ddccc(ii)) + ddthr(ii))
      
!     tau_thr_c_tl(ii) = tau_thr_c_tl(ii) + dd_ccc_tl
      dd_ccc_ad  = dd_ccc_ad + tau_thr_c_ad(ii)

!     dd_ccc_tl = ddccc(ii)*dd_thr_tl* (qq - dd_thr)/qq**2
      qq = ddccc(ii)**2 + dd_thr
      dd_thr_ad = dd_thr_ad + ddccc(ii) * (qq - dd_thr)/qq**2  * dd_ccc_ad
      dd_ccc_ad = 0

!     dd_thr_tl = ddthr_tl(ii) * (abs(ddccc(ii)) + 2._jprv* ddthr(ii))
      ddthr_ad(ii) = ddthr_ad(ii) + (abs(ddccc(ii)) + 2._jprv* ddthr(ii)) *dd_thr_ad
      dd_thr_ad    = 0
    enddo

!=====================
!   ddthr_tl(6) =   rel_taui_thr_X * taui_top_tl
!   ddthr_tl(4) =   rel_taui_thr_X * taui_bot_tl
!   ddthr_tl(3) =   rel_tauw_thr_X *(tauw_top_tl+tauwi_top_tl)
!   ddthr_tl(1) =   rel_tauw_thr_X *(tauw_bot_tl+tauwi_bot_tl)
!=====================
!   ddthr_tl(1) =   rel_tauw_thr_X *(tauw_bot_tl+tauwi_bot_tl)
    tauwi_bot_ad = tauwi_bot_ad + rel_tauw_thr_X * ddthr_ad(1) 
    tauw_bot_ad  = tauw_bot_ad  + rel_tauw_thr_X * ddthr_ad(1)
    ddthr_ad(1)  = 0

!   ddthr_tl(3) =   rel_tauw_thr_X *(tauw_top_tl+tauwi_top_tl)
    tauw_top_ad  = tauw_top_ad  + rel_tauw_thr_X *ddthr_ad(3)
    tauwi_top_ad = tauwi_top_ad + rel_tauw_thr_X *ddthr_ad(3)
    ddthr_ad(3)  = 0

!   ddthr_tl(4) =   rel_taui_thr_X * taui_bot_tl
    taui_bot_ad  = taui_bot_ad  + rel_taui_thr_X *ddthr_ad(4)
    ddthr_ad(4)  = 0

!   ddthr_tl(6) =   rel_taui_thr_X * taui_top_tl
    taui_top_ad  = taui_top_ad  + rel_taui_thr_X *ddthr_ad(6)
    ddthr_ad(6)  = 0

!=====================
!   tau_thr_c_tl(:)= 0
!   tau_thr_c_tl(5)= taui_top_tl
!   tau_thr_c_tl(4)= taui_top_tl + taui_bot_tl
!
!   tau_thr_c_tl(3)= tau_thr_c_tl(4)
!   tau_thr_c_tl(2)= tau_thr_c_tl(3) + tauw_top_tl + tauwi_top_tl
!   tau_thr_c_tl(1)= tau_thr_c_tl(2) + tauw_bot_tl + tauwi_bot_tl
!=====================
    tauwi_bot_ad   = tauwi_bot_ad   + tau_thr_c_ad(1)
    tauw_bot_ad    = tauw_bot_ad    + tau_thr_c_ad(1)
    tau_thr_c_ad(2)= tau_thr_c_ad(2)+ tau_thr_c_ad(1)
    tau_thr_c_ad(1)= 0

    tauwi_top_ad   = tauwi_top_ad   + tau_thr_c_ad(2)
    tauw_top_ad    = tauw_top_ad    + tau_thr_c_ad(2)
    tau_thr_c_ad(3)= tau_thr_c_ad(3)+ tau_thr_c_ad(2)
    tau_thr_c_ad(2)= 0

    tau_thr_c_ad(4)= tau_thr_c_ad(4)+ tau_thr_c_ad(3)
    tau_thr_c_ad(3)= 0

    taui_bot_ad    = taui_bot_ad    + tau_thr_c_ad(4)
    taui_top_ad    = taui_top_ad    + tau_thr_c_ad(4)
    tau_thr_c_ad(4)= 0

    taui_top_ad    = taui_top_ad    + tau_thr_c_ad(5)
    tau_thr_c_ad(5)= 0
    tau_thr_c_ad(:)= 0

!=====================
!   dtaubins_tl(:)= dtauw_tl(:) + dtaui_tl(:) + dtau_help_size_X * dp_tl(:) +dtau_help_size_tl*dp_X(:)
!=====================
    dtauw_ad(:) = dtauw_ad(:) +                    dtaubins_ad(:)
    dtaui_ad(:) = dtaui_ad(:) +                    dtaubins_ad(:)
    dp_ad(:)    = dp_ad(:)    + dtau_help_size_X * dtaubins_ad(:)
    dtau_help_size_ad = dtau_help_size_ad + sum(dp_X(:)*dtaubins_ad(:))

!=====================
!   q_mxr_tl(:) = profiles_tl(prof)%q(:) * gas_mass(gas_id_watervapour) * (1._jprv - q_mxr(:)) / &
!                       (mair * 1.E06_jprv + profiles(prof)%q(:) * gas_mass(gas_id_watervapour))
!
!   dp_tl (:) = 0
!   range_p_tl =0
!   dtau_help_size_tl= 0
!   IF (pressure_gradients) THEN
!     dp_tl (:) = profiles_tl(prof)%p_half(2:nlayers+1) - profiles_tl(prof)%p_half(1:nlayers)
!     range_p_tl= psfc_tl
!     dtau_help_size_tl= -profiles_tl(prof)%p_half(nlayers+1)*dtau_help_size_X /profiles(prof)%p_half(nlayers+1)
!   ENDIF
!
!   diwv_tl(:) = 100._jprv * (dp_tl (:)*q_mxr(:) + dp_X (:)*q_mxr_tl(:))/ gravity
!=====================
    q_mxr_ad(:)  = q_mxr_ad(:) + 100._jprv*dp_X (:)/ gravity  * diwv_ad(:)
    dp_ad   (:)  = dp_ad   (:) + 100._jprv*q_mxr(:)/ gravity  * diwv_ad(:)
    diwv_ad(:)   = 0

    IF (pressure_gradients) THEN
!     dtau_help_size_tl= -profiles_tl(prof)%p_half(nlayers+1)*dtau_help_size_X /profiles(prof)%p_half(nlayers+1)
      profiles_ad(profad)%p_half(nlayers+1) = profiles_ad(profad)%p_half(nlayers+1)     &
                                              - dtau_help_size_X/profiles(prof)%p_half(nlayers+1) * dtau_help_size_ad
      dtau_help_size_ad = 0

!     range_p_tl= psfc_tl
      psfc_ad = psfc_ad + range_p_ad
      range_p_ad = 0

!     DO ii=1, nlayers
!       dp_tl (ii) = profiles_tl(prof)%p_half(ii+1) - profiles_tl(prof)%p_half(ii)
!     ENDDO
      DO ii=1, nlayers
!       dp_ad (ii) =  profiles_ad(prof)%p_half(ii+1) - profiles_ad(prof)%p_half(ii)
        profiles_ad(profad)%p_half(ii+1) = profiles_ad(profad)%p_half(ii+1) + dp_ad (ii)
        profiles_ad(profad)%p_half(ii)   = profiles_ad(profad)%p_half(ii)   - dp_ad (ii)
      ENDDO
    ENDIF! IF (pressure_gradients) THEN

    dtau_help_size_ad= 0
    range_p_ad =0
    dp_ad (:) = 0

!   q_mxr_tl(:) = profiles_tl(prof)%q(:) * gas_mass(gas_id_watervapour) * (1._jprv - q_mxr(:)) / &
!                       (mair * 1.E06_jprv + profiles(prof)%q(:) * gas_mass(gas_id_watervapour))
    profiles_ad(profad)%q(:) = profiles_ad(profad)%q(:) +         &
                                        gas_mass(gas_id_watervapour) * (1._jprv - q_mxr(:)) / &
                                        (mair * 1.E06_jprv + profiles(prof)%q(:) * gas_mass(gas_id_watervapour)) &
                                        * q_mxr_ad(:)
    q_mxr_ad(:) = 0
!=================================================================================================
  ENDIF!(l_iwvc>0) THEN
!=================================================================================================
  IF(l_iwvc<2) THEN
    IF (SIZE(mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffi_top))%auxparams) > 0) THEN
      IF(taui > 0) THEN
        IF ( taui_top > 1E-6 ) THEN
!                    reffi_top_ad = SUM( (    -fboti_ad(:)) * dtaui(:) * (1._jprv-fmp) * dreffi(:)    +       &
!                                        (1._jprv-fboti(:)) * dtaui_ad(:)*(1._jprv-fmp) * dreffi(:)   +       &
!                                        (1._jprv-fboti(:)) * dtaui(:)   *( -fmp_ad) * dreffi(:)   +       &
!                                        (1._jprv-fboti(:)) * dtaui(:) *(1._jprv-fmp) * dreffi_ad(:) ) /taui_top - &
!                                   reffi_top/taui_top * taui_top_ad

          fboti_ad(:)   = fboti_ad(:)    - dtaui(:)           * (1._jprv-fmp) * dreffi(:) * reffi_top_ad/ taui_top
          dtaui_ad(:)   = dtaui_ad(:)    + (1._jprv-fboti(:)) * (1._jprv-fmp) * dreffi(:) * reffi_top_ad/ taui_top
          fmp_ad        = fmp_ad         - (1._jprv-fboti(:)) * dtaui(:)      * dreffi(:) * reffi_top_ad/ taui_top
          dreffi_ad(:)  = dreffi_ad(:)   + (1._jprv-fboti(:)) * (1._jprv-fmp) * dtaui(:)  * reffi_top_ad/ taui_top
  
          taui_top_ad   = taui_top_ad    -                                     reffi_top  * reffi_top_ad/ taui_top
          reffi_top_ad  = 0
        ELSE
          reffi_top_ad = 0._jprv
        ENDIF

        IF( taui_bot > 1E-6 ) THEN
!                    reffi_bot_ad = SUM( fboti_ad(:)     * dtaui(:) *(1._jprv-fmp) * dreffi(:)   +       &
!                                        fboti(:)        * dtaui_ad(:)*(1._jprv-fmp) * dreffi(:)   +       &
!                                        fboti(:)        * dtaui(:)   *( -fmp_ad) * dreffi(:)   +       &
!                                        fboti(:)        * dtaui(:) *(1._jprv-fmp) * dreffi_ad(:)    ) / taui_bot - &
!                                   reffi_bot/taui_bot * taui_bot_ad
          fboti_ad(:)    = fboti_ad(:)    + dtaui(:) *(1._jprv-fmp) * dreffi(:) * reffi_bot_ad/taui_bot
          dtaui_ad(:)    = dtaui_ad(:)    + fboti(:) *(1._jprv-fmp) * dreffi(:) * reffi_bot_ad/taui_bot
          fmp_ad         = fmp_ad         - fboti(:) * dtaui(:)     * dreffi(:) * reffi_bot_ad/taui_bot
          dreffi_ad(:)   = dreffi_ad(:)   + fboti(:) *(1._jprv-fmp) * dtaui(:)  * reffi_bot_ad/taui_bot
          taui_bot_ad    = taui_bot_ad    -                          reffi_bot  * reffi_bot_ad/taui_bot
          reffi_bot_ad = 0
        ELSE
          reffi_bot_ad = 0._jprv 
        ENDIF

!                 taui_bot_ad = SUM( fboti_ad(:) * dtaui   (:) * (1._jprv - fmp (:))    &
!                                  + fboti(:)    * dtaui_ad(:) * (1._jprv - fmp (:))    &
!                                  - fboti(:)    * dtaui   (:) * fmp_ad(:)  )
!                 taui_top_ad = taui_ad - taui_bot_ad
        taui_ad     = taui_ad     + taui_top_ad
        taui_bot_ad = taui_bot_ad - taui_top_ad
        taui_top_ad = 0
        dtaui_ad(:) = dtaui_ad(:) + fboti(:) * (1._jprv - fmp (:))   * taui_bot_ad 
        fboti_ad(:) = fboti_ad(:) + dtaui(:) * (1._jprv - fmp (:))   * taui_bot_ad
        fmp_ad(:)   = fmp_ad(:)   - fboti(:) * dtaui(:)              * taui_bot_ad
        taui_bot_ad = 0

        !---------------------------------------------------------------
        call comp_fbot_ad(fboti_ad, taui_thr_ad, ddtaui_ad, taui_thr, dtaui*(1._jprv - fmp(:)), nlayers, &
                          mfasisnn_coefs%smoothing_parameter)
        !---------------------------------------------------------------
        fboti_ad = 0

!                 taui_top_ad = taui_thr_ad
!                 taui_bot_ad = taui_ad    - taui_thr_ad
!                 ddtaui_ad(:)= dtaui_ad(:)*(1._jprv-fmp) - dtaui(:)*fmp_ad
        fmp_ad      = fmp_ad        - dtaui(:)* ddtaui_ad(:)
        dtaui_ad(:) = dtaui_ad(:)   +  ddtaui_ad(:)*(1._jprv-fmp)                 
        ddtaui_ad(:)= 0

!--------------------------------------------------------
        fac_hi =  mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffi_top))%auxparams(3) &
                  - mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffi_top))%auxparams(4) * MAX( vza, sza ) / 90._jprv
        fac_lo = 0.5_jprv 
        tau_lo = mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffi_top))%auxparams(1)
        tau_hi = mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffi_top))%auxparams(2)
!--------------------------------------------------------
!                 taui_thr_ad = (fac_lo- (fac_lo-fac_hi) * 0.5_jprv*( 1._jprv - COS( logtaui * pi ) ) ) * taui_ad
!                 IF(logtaui > 0.0_jprv .AND. logtaui <  1.0_jprv ) then
!                   logtaui_ad  = taui_ad/(taui + 1E-6) / (LOG(tau_hi) - LOG(tau_lo))
!                   taui_thr_ad = taui_thr_ad - (fac_lo-fac_hi)*0.5_jprv*pi* SIN(logtaui * pi) * taui * logtaui_ad
!                 ENDIF
        ! logtaui = MIN( MAX( (LOG(taui) - LOG(tau_lo)) / (LOG(tau_hi) - LOG(tau_lo)), 0._jprv ), 1._jprv ) 
        IF(logtaui > 0.0_jprv .AND. logtaui <  1.0_jprv ) then
!                   logtaui_ad  = taui_ad/(taui + 1E-6) / (LOG(tau_hi) - LOG(tau_lo))
!                   taui_thr_ad = taui_thr_ad - (fac_lo-fac_hi)*0.5_jprv*pi* SIN(logtaui * pi) * taui * logtaui_ad
          logtaui_ad  = logtaui_ad  - (fac_lo-fac_hi)*0.5_jprv*pi* SIN(logtaui * pi) * taui * taui_thr_ad

          taui_ad     = taui_ad + logtaui_ad / taui / (LOG(tau_hi) - LOG(tau_lo))
          logtaui_ad  = 0
        ENDIF
!                 taui_thr_ad = (fac_lo- (fac_lo-fac_hi) * 0.5_jprv*( 1._jprv - COS( logtaui * pi ) ) ) * taui_ad
        taui_ad = taui_ad + (fac_lo- (fac_lo-fac_hi) * 0.5_jprv*( 1._jprv - COS( logtaui * pi))) * taui_thr_ad
        taui_thr_ad = 0

      ELSE ! taui > 0)
        reffi_bot_ad = 0._jprv
        reffi_top_ad = 0._jprv
      ENDIF
    ELSE
      err = errorstatus_fatal
      THROWM(err.NE.0, "MFASIS-NN has no information on two-layer parameterisation of ice cloud")
    ENDIF

    IF (SIZE(mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffw_top))%auxparams) > 0) THEN
      IF(tauw > 0) THEN
        IF ( tauwi_top > 1E-6 ) THEN
!                    reffwi_top_ad = SUM(-fbotw_ad(:)    * dtaui(:)   * fmp(:) * dreffi(:) +       &
!                                    (1._jprv - fbotw(:))* dtaui_ad(:)* fmp(:) * dreffi(:) +       &
!                                    (1._jprv - fbotw(:))* dtaui(:)   * fmp_ad(:) * dreffi(:) +       &
!                                    (1._jprv - fbotw(:))* dtaui(:)   * fmp(:) * dreffi_ad(:)    ) / tauwi_top -  &
!                                    reffwi_top/tauwi_top * tauwi_top_ad
          fbotw_ad(:)    = fbotw_ad(:)    - dtaui(:)           * fmp(:)  * dreffi(:) * reffwi_top_ad/tauwi_top
          dtaui_ad(:)    = dtaui_ad(:)    +(1._jprv - fbotw(:))* fmp(:)  * dreffi(:) * reffwi_top_ad/tauwi_top 
          fmp_ad(:)      = fmp_ad(:)      +(1._jprv - fbotw(:))* dtaui(:)* dreffi(:) * reffwi_top_ad/tauwi_top
          dreffi_ad(:)   = dreffi_ad(:)   +(1._jprv - fbotw(:))* fmp(:)  * dtaui(:)  * reffwi_top_ad/tauwi_top
          tauwi_top_ad   = tauwi_top_ad   -                              reffwi_top  * reffwi_top_ad/tauwi_top
          reffwi_top_ad = 0
        ELSE
          reffwi_top_ad = 0._jprv !
        ENDIF
!                 tauwi_top_ad = SUM(-fbotw_ad(:) * dtaui(:)   * fmp(:)    + &
!                         (1._jprv - fbotw(:)) * dtaui_ad(:)* fmp(:)    + &
!                         (1._jprv - fbotw(:)) * dtaui(:)   * fmp_ad(:)     )
        fmp_ad(:)   = fmp_ad(:)   + (1._jprv - fbotw(:)) * dtaui(:) * tauwi_top_ad 
        dtaui_ad(:) = dtaui_ad(:) + (1._jprv - fbotw(:)) * fmp(:)   * tauwi_top_ad 
        fbotw_ad(:) = fbotw_ad(:) -            dtaui(:)  * fmp(:)   * tauwi_top_ad 
        tauwi_top_ad = 0


        IF ( tauwi_bot > 1E-6 ) THEN
!                    reffwi_bot_ad = SUM( fbotw_ad(:)    * dtaui(:)   * fmp(:) * dreffi(:) +       &
!                                        fbotw(:)        * dtaui_ad(:)* fmp(:) * dreffi(:) +       &
!                                        fbotw(:)        * dtaui(:)   * fmp_ad(:) * dreffi(:) +       &
!                                        fbotw(:)        * dtaui(:)   * fmp(:) * dreffi_ad(:)    ) / tauwi_bot - &
!                                    reffwi_bot/tauwi_bot * tauwi_bot_ad
          fbotw_ad(:)   = fbotw_ad(:)    +  dtaui(: )* fmp(:)   * dreffi(:) * reffwi_bot_ad/tauwi_bot
          dtaui_ad(:)   = dtaui_ad(:)    +  fbotw(:) * fmp(:)   * dreffi(:) * reffwi_bot_ad/tauwi_bot
          fmp_ad(:)     = fmp_ad(:)      +  fbotw(:) * dtaui(:) * dreffi(:) * reffwi_bot_ad/tauwi_bot
          dreffi_ad(:)  = dreffi_ad(:)   +  fbotw(:) * fmp(:)   * dtaui(:)  * reffwi_bot_ad/tauwi_bot
          tauwi_bot_ad  = tauwi_bot_ad   -                       reffwi_bot * reffwi_bot_ad/tauwi_bot
          reffwi_bot_ad = 0
        ELSE
          reffwi_bot_ad = 0._jprv !
        ENDIF
!                 tauwi_bot_ad = SUM( fbotw_ad(:) * dtaui(:)   * fmp(:)    + &
!                                     fbotw(:)    * dtaui_ad(:)* fmp(:)    + &
!                                     fbotw(:)    * dtaui(:)   * fmp_ad(:)     )
        fbotw_ad(:)  = fbotw_ad(:)  + dtaui(:)   * fmp(:)  * tauwi_bot_ad
        dtaui_ad(:)  = dtaui_ad(:)  + fbotw(:)   * fmp(:)  * tauwi_bot_ad
        fmp_ad(:)    = fmp_ad(:)    + fbotw(:)   * dtaui(:)* tauwi_bot_ad
        tauwi_bot_ad = 0

        IF ( tauw_top > 1E-6 ) THEN
!                    reffw_top_ad = SUM(     -fbotw_ad(:)* dtauw(:)    * dreffw(:) +       &
!                                     (1._jprv-fbotw(:)) * dtauw_ad(:) * dreffw(:) +       &
!                                     (1._jprv-fbotw(:)) * dtauw(:)    * dreffw_ad(:)    ) / tauw_top     -       &
!                                   reffw_top/tauw_top * tauw_top_ad
          fbotw_ad    = fbotw_ad    -                   dtauw(:)  * dreffw(:) * reffw_top_ad/tauw_top
          dtauw_ad(:) = dtauw_ad(:) +(1._jprv-fbotw(:))           * dreffw(:) * reffw_top_ad/tauw_top
          dreffw_ad(:)= dreffw_ad(:)+(1._jprv-fbotw(:))*dtauw(:)              * reffw_top_ad/tauw_top
          tauw_top_ad = tauw_top_ad -                               reffw_top * reffw_top_ad/tauw_top
          reffw_top_ad = 0
        ELSE
          reffw_top_ad = 0._jprv !
        ENDIF

        IF ( tauw_bot > 1E-6 ) THEN
!                    reffw_bot_ad = SUM( fbotw_ad(:)     * dtauw(:)    * dreffw(:) +       &
!                                        fbotw(:)        * dtauw_ad(:) * dreffw(:) +       &
!                                        fbotw(:)        * dtauw(:)    * dreffw_ad(:)    ) / tauw_bot     -       &
!                                   reffw_bot/tauw_bot * tauw_bot_ad
          fbotw_ad(:)  = fbotw_ad(:)  +         dtauw(:) * dreffw(:) * reffw_bot_ad/ tauw_bot
          dtauw_ad(:)  = dtauw_ad(:)  +fbotw(:)          * dreffw(:) * reffw_bot_ad/ tauw_bot
          dreffw_ad(:) = dreffw_ad(:) +fbotw(:)*dtauw(:)             * reffw_bot_ad/ tauw_bot
          tauw_bot_ad  = tauw_bot_ad  -                   reffw_bot  * reffw_bot_ad/ tauw_bot
          reffw_bot_ad = 0
        ELSE
          reffw_bot_ad = 0._jprv !
        ENDIF

!                 tauw_bot_ad = SUM( fbotw_ad(:) * dtauw(:) + fbotw(:) * dtauw_ad(:))
!                 tauw_top_ad = tauw_ad - tauw_bot_ad
        tauw_ad     = tauw_ad     + tauw_top_ad 
        tauw_bot_ad = tauw_bot_ad - tauw_top_ad 
        tauw_top_ad = 0
        fbotw_ad(:) = fbotw_ad(:) + dtauw(:) * tauw_bot_ad
        dtauw_ad(:) = dtauw_ad(:) + fbotw(:) * tauw_bot_ad
        tauw_bot_ad = 0
        !---------------------------------------------------------------
        call comp_fbot_ad(fbotw_ad, tauw_thr_ad, dtauw_ad, tauw_thr, dtauw, nlayers, &
                          mfasisnn_coefs%smoothing_parameter)
        !---------------------------------------------------------------
        fbotw_ad = 0

!                 tauw_top_ad = tauw_thr_ad
!                 tauw_bot_ad = tauw_ad    - tauw_thr_ad

        tauw_thr_ad = tauw_thr_ad + tauw_top_ad
        tauw_top_ad = 0

!-----------------------------------
        fac_hi =  mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffw_top))%auxparams(3) &
                  - mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffw_top))%auxparams(4) * MAX( vza, sza ) / 90._jprv

        fac_lo = 0.5_jprv !
        tau_lo = mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffw_top))%auxparams(1)
        tau_hi = mfasisnn_coefs%nn(nch)%in(nn_inpar_idx(nch, nn_idx_reffw_top))%auxparams(2)
!-----------------------------------
!               tauw_thr_ad = ( fac_lo - (fac_lo-fac_hi) * 0.5_jprv*( 1._jprv - COS( logtauw * pi ) ) ) * tauw_ad
!               IF(logtauw > 0.0_jprv .AND. logtauw <  1.0_jprv ) then
!                 logtauw_ad  = tauw_ad/(tauw + 1E-6) / (LOG(tau_hi) - LOG(tau_lo))
!                 tauw_thr_ad = tauw_thr_ad - (fac_lo-fac_hi)*0.5_jprv*pi* SIN(logtauw * pi) * tauw * logtauw_ad
!               ENDIF
        ! logtauw = MIN( MAX( (LOG(tauw) - LOG(tau_lo)) / (LOG(tau_hi) - LOG(tau_lo)), 0._jprv ), 1._jprv ) 
        IF(logtauw > 0.0_jprv .AND. logtauw <  1.0_jprv ) then
!                 logtauw_ad  = tauw_ad/(tauw + 1E-6) / (LOG(tau_hi) - LOG(tau_lo))
!                 tauw_thr_ad = tauw_thr_ad - (fac_lo-fac_hi)*0.5_jprv*pi* SIN(logtauw * pi) * tauw * logtauw_ad
          logtauw_ad = logtauw_ad   - (fac_lo-fac_hi)*0.5_jprv*pi* SIN(logtauw * pi) * tauw * tauw_thr_ad
          tauw_ad    = tauw_ad      + logtauw_ad / tauw / (LOG(tau_hi) - LOG(tau_lo))
          logtauw_ad = 0
        ENDIF

!           tauw_thr_ad =          (fac_lo- (fac_lo-fac_hi)* 0.5_jprv*( 1._jprv - COS(logtauw*pi))) * tauw_ad
        tauw_ad     = tauw_ad +(fac_lo- (fac_lo-fac_hi)* 0.5_jprv*( 1._jprv - COS(logtauw*pi))) * tauw_thr_ad
        tauw_thr_ad = 0
!---------------------=====================
      ELSE ! tauw > 0)
        reffw_bot_ad = 0._jprv
        reffw_top_ad = 0._jprv
        reffwi_bot_ad = 0._jprv
        reffwi_top_ad = 0._jprv
      ENDIF  ! tauw > 0
    ELSE
      err = errorstatus_fatal
      THROWM(err.NE.0, "MFASIS-NN has no information on two-layer parameterisation of water cloud")
    ENDIF

!            tauwi_ad = SUM(   dtaui_ad(:) * fmp(:)     &
!                           +  dtaui   (:) * fmp_ad(:) ) ! mixed-phase cloud ice in water cloud
!            taui_ad  = SUM(   dtaui_ad(:) * (1._jprv - fmp(:))  &
!                           -  dtaui   (:) * fmp_ad(:) ) ! mixed-phase cloud ice in water cloud
    fmp_ad(:)   = fmp_ad(:)   - dtaui(:)           * taui_ad
    dtaui_ad(:) = dtaui_ad(:) + (1._jprv - fmp(:)) * taui_ad
    taui_ad  = 0
    fmp_ad(:)   = fmp_ad(:)   + dtaui(:)           * tauwi_ad
    dtaui_ad(:) = dtaui_ad(:) + fmp(:)             * tauwi_ad
    tauwi_ad = 0

    tau_thr    =  1._jprv
    !---------------------------------------------------------------
    call comp_fbot_ad(fmp_ad, tau_thr_ad, dtauw_ad, tau_thr, dtauw, nlayers, &
                      mfasisnn_coefs%smoothing_parameter)
    !---------------------------------------------------------------
    fmp_ad = 0
    tau_thr_ad =  0

    DO j = 1, nlayers
      coli = cloud_columns%icldarr(cc,j,prof)
      IF ( coli > 0 ) THEN
!                  tauw_ad = tauw_ad + dtauw_ad(j)
        dtauw_ad(j) = dtauw_ad(j) + tauw_ad
!                  dreffi_ad(j) = aux_ad%hydro_deff(ice_baum_index,j,coli,prof)/2._jprv  !
        aux_ad%hydro_deff(ice_baum_index,j,coli,profad) = &
            aux_ad%hydro_deff(ice_baum_index,j,coli,profad) + dreffi_ad(j)/2._jprv  !
        dreffi_ad(j) = 0

        ! layer eff. radii water
        IF ( dtauw(j) .GT. 0 ) THEN
!                        dreffw_ad(j) = dreffw_ad(j)/dtauw(j) - dreffw(j)/(dtauw(j))    * dtauw_ad(j)
          dtauw_ad(j)  = dtauw_ad(j)    -  dreffw(j)/(dtauw(j))    * dreffw_ad(j)
          dreffw_ad(j) = dreffw_ad(j)/dtauw(j)

          aux_ad%hydro_deff(clw_deff_index,j,coli,profad) = &
              aux_ad%hydro_deff(clw_deff_index,j,coli,profad) + &
              dreffw_ad(j) /2._jprv * opdpext(clw_deff_index,coli,j,ist)
          opdpext_ad(clw_deff_index,coli,j,ist) = &
              opdpext_ad(clw_deff_index,coli,j,ist) + &
              dreffw_ad(j) /2._jprv * aux%hydro_deff(clw_deff_index,j,coli,prof)
          opdpext_ad(1:nopac_wcl,coli,j,ist) = &
              opdpext_ad(1:nopac_wcl,coli,j,ist) + &
              dreffw_ad(j) * wcl_opac_deff(1:nopac_wcl)/2._jprv
        ELSE
          dreffw_ad(j) = 0._jprv !
        ENDIF

        ! layer optical depths
        opdpext_ad(ice_baum_index,coli,j,ist) = &
            opdpext_ad(ice_baum_index,coli,j,ist) + dtaui_ad(j)
        opdpext_ad(1:clw_deff_index,coli,j,ist) = &
            opdpext_ad(1:clw_deff_index,coli,j,ist) + dtauw_ad(j)
      ENDIF

    ENDDO ! DO j = 1, nlay
  ENDIF!(l_iwvc<2) 

  tauw_ad = 0.0_jprv
  taui_ad = 0.0_jprv
  tauwi_top_ad = 0.0_jprv
  tauwi_bot_ad = 0.0_jprv
  reffw_top_ad = 0.0_jprv
  reffi_top_ad = 0.0_jprv
  reffw_bot_ad = 0.0_jprv
  reffi_bot_ad = 0.0_jprv
  reffwi_top_ad = 0.0_jprv
  reffwi_bot_ad = 0.0_jprv
  fpct_ad = 0._jprv

  ! Extract optical depths and effective radii per layer
  dtauw_ad(:) = 0._jprv
  dtaui_ad(:) = 0._jprv
  fmp_ad(:) = 0._jprv
  dreffw_ad(:) = 0._jprv
  dreffi_ad(:) = 0._jprv


  IF (lhook) CALL dr_hook('RTTOV_CALC_MFASIS_NN_HYDRO_INPAR_AD',1_jpim,zhook_handle)

  CATCH

  IF (lhook) CALL dr_hook('RTTOV_CALC_MFASIS_NN_HYDRO_INPAR_AD',1_jpim,zhook_handle)

END SUBROUTINE rttov_calc_mfasis_nn_hydro_inpar_ad
