!
! ADA_Module
!
! Module containing the Adding Doubling Adding (ADA) radiative
! transfer solution procedures used in the CRTM.
!
!
! CREATION HISTORY:
!       Written by:     Quanhua Liu, QSS at JCSDA; quanhua.liu@noaa.gov
!                       Yong Han,    NOAA/NESDIS;  yong.han@noaa.gov
!                       Paul van Delst; CIMMS/SSEC; paul.vandelst@noaa.gov
!                       08-Jun-2004

MODULE ADA_Module

  ! ------------------
  ! Environemnt set up
  ! ------------------
  ! Module use statements
  USE RTV_Define
  USE CRTM_Parameters
  USE Type_Kinds
  USE Message_Handler
  USE CRTM_Utility
    
  IMPLICIT NONE
  
  ! --------------------
  ! Default visibilities
  ! --------------------
  ! Everything private by default
  PRIVATE
  
  PUBLIC CRTM_ADA
  PUBLIC CRTM_ADA_TL
  PUBLIC CRTM_ADA_AD
  
  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module
  CHARACTER(*),  PARAMETER :: MODULE_VERSION_ID = &
  '$Id: $'
  
CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                         ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                            ##
!################################################################################
!################################################################################

  
  SUBROUTINE CRTM_ADA(n_Layers, & ! Input  number of atmospheric layers
                             w, & ! Input  layer scattering albedo
                          T_OD, & ! Input  layer optical depth
             cosmic_background, & ! Input  cosmic background radiance
                    emissivity, & ! Input  surface emissivity
                  reflectivity, & ! Input  surface reflectivity matrix 
           direct_reflectivity, & ! Input  surface direct reflectivity 
                           RTV)   ! IN/Output upward radiance and others
! ------------------------------------------------------------------------- !
! FUNCTION:                                                                 !
!   This subroutine calculates IR/MW radiance at the top of the atmosphere  !
!   including atmospheric scattering. The scheme will include solar part.   !
!   The ADA algorithm computes layer reflectance and transmittance as well  !
!   as source function by the subroutine CRTM_Doubling_layer, then uses     !
!   an adding method to integrate the layer and surface components.         !
!                                                                           ! 
!    Quanhua Liu    Quanhua.Liu@noaa.gov                                    !
! ------------------------------------------------------------------------- !
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n_Layers
    INTEGER nZ
    TYPE(RTV_type), INTENT( INOUT ) :: RTV
    REAL (fp), INTENT(IN), DIMENSION( : ) ::  w,T_OD
    REAL (fp), INTENT(IN), DIMENSION( : ) ::  emissivity, direct_reflectivity
    REAL (fp), INTENT(IN), DIMENSION( :,: ) :: reflectivity 
    REAL (fp), INTENT(IN) ::  cosmic_background

    ! -------------- internal variables --------------------------------- !
    !  Abbreviations:                                                     !
    !      s: scattering, rad: radiance, trans: transmission,             !
    !         refl: reflection, up: upward, down: downward                !
    ! --------------------------------------------------------------------!
    REAL (fp), DIMENSION(RTV%n_Angles, RTV%n_Angles) :: temporal_matrix
    REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down 
    REAL (fp), DIMENSION(0:n_Layers) :: total_opt
    INTEGER :: i, j, k, Error_Status
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_ADA'
    CHARACTER(256) :: Message
 
    total_opt(0) = ZERO
    DO k = 1, n_Layers
      total_opt(k) = total_opt(k-1) + T_OD(k)
    END DO

    nZ = RTV%n_Angles
    RTV%s_Layer_Trans = ZERO
    RTV%s_Layer_Refl = ZERO
    RTV%s_Level_Refl_UP = ZERO
    RTV%s_Level_Rad_UP = ZERO
    RTV%s_Layer_Source_UP = ZERO
    RTV%s_Layer_Source_DOWN = ZERO

    RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,n_Layers)=reflectivity(1:RTV%n_Angles,1:RTV%n_Angles)

    IF( RTV%mth_Azi == 0 ) THEN
      RTV%s_Level_Rad_UP(1:RTV%n_Angles,n_Layers ) = emissivity(1:RTV%n_Angles)*RTV%Planck_Surface
    END IF

    IF( RTV%Solar_Flag_true ) THEN
      RTV%s_Level_Rad_UP(1:nZ,n_Layers ) = RTV%s_Level_Rad_UP(1:nZ,n_Layers )+direct_reflectivity(1:nZ)* &
        RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_opt(n_Layers)/RTV%COS_SUN)       
    END IF

    ! UPWARD ADDING LOOP STARTS FROM BOTTOM LAYER TO ATMOSPHERIC TOP LAYER.
    DO 10 k = n_Layers, 1, -1

    ! Compute tranmission and reflection matrices for a layer
    IF(w(k) > SCATTERING_ALBEDO_THRESHOLD) THEN 

    !  ----------------------------------------------------------- !
    !    CALL  multiple-stream algorithm for computing layer       !
    !    transmission, reflection, and source functions.           !
    !  ----------------------------------------------------------- !

    CALL CRTM_AMOM_layer( &
           RTV%n_Streams,            &
           RTV%n_Angles,k,w(k),      &
           T_OD(k),                  &
           total_opt(k-1),           &
           RTV%COS_Angle,            & ! Input
           RTV%COS_Weight,           &
           RTV%Pff(:,:,k),           &
           RTV%Pbb(:,:,k),           & ! Input
           RTV%Planck_Atmosphere(k), & ! Input
           RTV                       ) ! Internal variable  
    !  ----------------------------------------------------------- !
    !    Adding method to add the layer to the present level       !
    !    to compute upward radiances and reflection matrix         !
    !    at new level.                                             !
    !  ----------------------------------------------------------- !

    temporal_matrix = -matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),  &
                       RTV%s_Layer_Refl(1:RTV%n_Angles,1:RTV%n_Angles,k))
    DO i = 1, RTV%n_Angles 
      temporal_matrix(i,i) = ONE + temporal_matrix(i,i)
    END DO

    RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k) = matinv(temporal_matrix, Error_Status)
    IF( Error_Status /= SUCCESS  ) THEN
      WRITE( Message,'("Error in matrix inversion matinv(temporal_matrix, Error_Status) ")' ) 
      CALL Display_Message( ROUTINE_NAME,  &                                                    
                            TRIM(Message), &                                                   
                            Error_Status   )                                          
      RETURN                                                                                    
    END IF
         
    RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k) =   &
     matmul(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k), RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k))
    refl_down(1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),  &
                                  RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))

    RTV%s_Level_Rad_UP(1:RTV%n_Angles,k-1 )=RTV%s_Layer_Source_UP(1:RTV%n_Angles,k)+ &
    matmul(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k),refl_down(1:RTV%n_Angles,k) &
          +RTV%s_Level_Rad_UP(1:RTV%n_Angles,k ))
    RTV%Refl_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k), &
          RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k))
    RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k-1)=RTV%s_Layer_Refl(1:RTV%n_Angles,1:RTV%n_Angles,k) + &
    matmul(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k),RTV%Refl_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)) 

    ELSE
      DO i = 1, RTV%n_Angles 
        RTV%s_Layer_Trans(i,i,k) = exp(-T_OD(k)/RTV%COS_Angle(i))
        RTV%s_Layer_Source_UP(i,k) = RTV%Planck_Atmosphere(k) * (ONE - RTV%s_Layer_Trans(i,i,k) )
        RTV%s_Layer_Source_DOWN(i,k) = RTV%s_Layer_Source_UP(i,k)

      END DO

      ! Adding method
      DO i = 1, RTV%n_Angles 
        RTV%s_Level_Rad_UP(i,k-1 )=RTV%s_Layer_Source_UP(i,k)+ &
        RTV%s_Layer_Trans(i,i,k)*(sum(RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,k)*RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))  &
         +RTV%s_Level_Rad_UP(i,k ))
      ENDDO
        DO i = 1, RTV%n_Angles 
          DO j = 1, RTV%n_Angles 
            RTV%s_Level_Refl_UP(i,j,k-1)=RTV%s_Layer_Trans(i,i,k)*RTV%s_Level_Refl_UP(i,j,k)*RTV%s_Layer_Trans(j,j,k)
          ENDDO
        ENDDO
    ENDIF
    
    10     CONTINUE

    !  Adding reflected cosmic background radiation
    IF( RTV%mth_Azi == 0 ) THEN
       DO i = 1, RTV%n_Angles 
         RTV%s_Level_Rad_UP(i,0)=RTV%s_Level_Rad_UP(i,0)+sum(RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,0))*cosmic_background
       ENDDO
    END IF

    RETURN
    
  END SUBROUTINE CRTM_ADA 
      
  SUBROUTINE CRTM_AMOM_layer(    n_streams, & ! Input, number of streams
                                  nZ, & ! Input, number of angles
                                  KL, & ! Input, KL-th layer 
                       single_albedo, & ! Input, single scattering albedo
                       optical_depth, & ! Input, layer optical depth
                           total_opt, & ! Input, accumulated optical depth from the top to current layer top
                           COS_Angle, & ! Input, COSINE of ANGLES
                          COS_Weight, & ! Input, GAUSSIAN Weights
                                  ff, & ! Input, Phase matrix (forward part)
                                  bb, & ! Input, Phase matrix (backward part)
                         Planck_Func, & ! Input, Planck for layer temperature
                                 RTV)   ! Output, layer transmittance, reflectance, and source 
! ---------------------------------------------------------------------------------------
!   FUNCTION
!    Compute layer transmission, reflection matrices and source function 
!    at the top and bottom of the layer.
!
!   Method and References
!    The transmittance and reflectance matrices is further derived from 
!    matrix operator method. The matrix operator method is referred to the paper by
!
!    Weng, F., and Q. Liu, 2003: Satellite Data Assimilation in Numerical Weather Prediction
!    Model: Part 1: Forward Radiative Transfer and Jacobian Modeling in Cloudy Atmospheres,
!    J. Atmos. Sci., 60, 2633-2646.
!
!   see also ADA method.
!   Quanhua Liu
!   Quanhua.Liu@noaa.gov
! ----------------------------------------------------------------------------------------
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: n_streams,nZ,KL
   TYPE(RTV_type), INTENT( INOUT ) :: RTV
   REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
   REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight 
   REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt

   ! internal variables
   REAL(fp), DIMENSION(nZ,nZ) :: trans, refl, tempo
   REAL(fp) :: s, c, xx
   INTEGER :: i,j,N2,N2_1
   INTEGER :: Error_Status
   REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
   REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)    
   CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer'
   CHARACTER(256) :: Message

   ! for small layer optical depth, single scattering is applied.  
   IF( optical_depth < DELTA_OPTICAL_DEPTH ) THEN
     s = optical_depth * single_albedo
     DO i = 1, nZ
       RTV%Thermal_C(i,KL) = ZERO
       c = s/COS_Angle(i)
       DO j = 1, nZ
         RTV%s_Layer_Refl(i,j,KL) = c * bb(i,j) * COS_Weight(j)
         RTV%s_Layer_Trans(i,j,KL) = c * ff(i,j) * COS_Weight(j)
         IF( i == j ) THEN
           RTV%s_Layer_Trans(i,i,KL) = RTV%s_Layer_Trans(i,i,KL) + &
             ONE - optical_depth/COS_Angle(i)
         END IF
         IF( RTV%mth_Azi == 0 ) THEN
           RTV%Thermal_C(i,KL) = RTV%Thermal_C(i,KL) + &
           ( RTV%s_Layer_Refl(i,j,KL) + RTV%s_Layer_Trans(i,j,KL) )
         END IF
       ENDDO

       IF( RTV%mth_Azi == 0 ) THEN
         RTV%s_Layer_Source_UP(i,KL) = ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func
         RTV%s_Layer_Source_DOWN(i,KL) = RTV%s_Layer_Source_UP(i,KL)
       END IF
     ENDDO

     RETURN

   END IF
   !
   ! for numerical stability, 
   IF( single_albedo < max_albedo ) THEN
     s = single_albedo
   ELSE
     s = max_albedo
   END IF
   !
   ! building phase matrices
   DO i = 1, nZ
     c = s/COS_Angle(i)
     DO j = 1, nZ
       RTV%PM(i,j,KL) = c * bb(i,j) * COS_Weight(j)
       RTV%PP(i,j,KL) = c * ff(i,j) * COS_Weight(j)
     ENDDO
       RTV%PP(i,i,KL) = RTV%PP(i,i,KL) - ONE/COS_Angle(i)
   ENDDO
   RTV%PPM(1:nZ,1:nZ,KL) = RTV%PP(1:nZ,1:nZ,KL) - RTV%PM(1:nZ,1:nZ,KL)
   RTV%i_PPM(1:nZ,1:nZ,KL) = matinv( RTV%PPM(1:nZ,1:nZ,KL), Error_Status )
   IF( Error_Status /= SUCCESS  ) THEN
     WRITE( Message,'("Error in matrix inversion matinv( RTV%PPM(1:nZ,1:nZ,KL), Error_Status ) ")' ) 
     CALL Display_Message( ROUTINE_NAME, &                                                    
                           TRIM(Message), &                                                   
                           Error_Status )                                          
     RETURN                                                                                    
   END IF

   RTV%PPP(1:nZ,1:nZ,KL) = RTV%PP(1:nZ,1:nZ,KL) + RTV%PM(1:nZ,1:nZ,KL)
   RTV%HH(1:nZ,1:nZ,KL) = matmul( RTV%PPM(1:nZ,1:nZ,KL), RTV%PPP(1:nZ,1:nZ,KL) )   
   !
   ! save phase element RTV%HH, call ASYMTX for calculating eigenvalue and vectors.
   tempo = RTV%HH(1:nZ,1:nZ,KL)
   CALL ASYMTX(tempo,nZ,nZ,nZ,RTV%EigVe(1:nZ,1:nZ,KL),RTV%EigVa(1:nZ,KL),Error_Status)
   DO i = 1, nZ
     IF( RTV%EigVa(i,KL) > ZERO ) THEN         
       RTV%EigValue(i,KL) = sqrt( RTV%EigVa(i,KL) )
     ELSE
       RTV%EigValue(i,KL) = ZERO
     END IF
   END DO

   DO i = 1, nZ
     DO j = 1, nZ
       RTV%EigVeVa(i,j,KL) = RTV%EigVe(i,j,KL) * RTV%EigValue(j,KL)
     END DO
   END DO       
   RTV%EigVeF(1:nZ,1:nZ,KL) = matmul( RTV%i_PPM(1:nZ,1:nZ,KL), RTV%EigVeVa(1:nZ,1:nZ,KL) )

   ! compute layer reflection, transmission and source function
   RTV%Gp(1:nZ,1:nZ,KL) = ( RTV%EigVe(1:nZ,1:nZ,KL) + RTV%EigVeF(1:nZ,1:nZ,KL) )/2.0_fp
   RTV%Gm(1:nZ,1:nZ,KL) = ( RTV%EigVe(1:nZ,1:nZ,KL) - RTV%EigVeF(1:nZ,1:nZ,KL) )/2.0_fp
   RTV%i_Gm(1:nZ,1:nZ,KL) = matinv( RTV%Gm(1:nZ,1:nZ,KL), Error_Status)

   IF( Error_Status /= SUCCESS  ) THEN
     WRITE( Message,'("Error in matrix inversion matinv( RTV%Gm(1:nZ,1:nZ,KL), Error_Status) ")' ) 
     CALL Display_Message( ROUTINE_NAME, &                                                    
                           TRIM(Message), &                                                   
                           Error_Status )                                          
     RETURN                                                                                    
   END IF             

   DO i = 1, nZ
     xx = RTV%EigValue(i,KL)*optical_depth
     RTV%Exp_x(i,KL) = exp(-xx)
   END DO

   DO i = 1, nZ
     DO j = 1, nZ
       RTV%A1(i,j,KL) = RTV%Gp(i,j,KL) * RTV%Exp_x(j,KL)
       RTV%A4(i,j,KL) = RTV%Gm(i,j,KL) * RTV%Exp_x(j,KL)
     END DO
   END DO

   RTV%A2(1:nZ,1:nZ,KL) = matmul( RTV%i_Gm(1:nZ,1:nZ,KL), RTV%A1(1:nZ,1:nZ,KL) )
   RTV%A3(1:nZ,1:nZ,KL) = matmul( RTV%Gp(1:nZ,1:nZ,KL), RTV%A2(1:nZ,1:nZ,KL) )
   RTV%A5(1:nZ,1:nZ,KL) = matmul( RTV%A1(1:nZ,1:nZ,KL), RTV%A2(1:nZ,1:nZ,KL) )
   RTV%A6(1:nZ,1:nZ,KL) = matmul( RTV%A4(1:nZ,1:nZ,KL), RTV%A2(1:nZ,1:nZ,KL) )
   RTV%Gm_A5(1:nZ,1:nZ,KL) = RTV%Gm(1:nZ,1:nZ,KL) - RTV%A5(1:nZ,1:nZ,KL)     
   RTV%i_Gm_A5(1:nZ,1:nZ,KL) = matinv(RTV%Gm_A5(1:nZ,1:nZ,KL), Error_Status)
   IF( Error_Status /= SUCCESS  ) THEN
     WRITE( Message,'("Error in matrix inversion matinv(RTV%Gm_A5(1:nZ,1:nZ,KL), Error_Status) ")' ) 
     CALL Display_Message( ROUTINE_NAME, &                                                    
                           TRIM(Message), &                                                   
                           Error_Status )                                          
     RETURN                                                                                    
   END IF
   trans = matmul( RTV%A4(1:nZ,1:nZ,KL) - RTV%A3(1:nZ,1:nZ,KL), RTV%i_Gm_A5(1:nZ,1:nZ,KL) ) 
   refl = matmul( RTV%Gp(1:nZ,1:nZ,KL) - RTV%A6(1:nZ,1:nZ,KL), RTV%i_Gm_A5(1:nZ,1:nZ,KL) )

   ! post processing  
   RTV%s_Layer_Trans(1:nZ,1:nZ,KL) = trans(:,:)
   RTV%s_Layer_Refl(1:nZ,1:nZ,KL) = refl(:,:)
   RTV%s_Layer_Source_UP(:,KL) = ZERO
   IF( RTV%mth_Azi == 0 ) THEN
     DO i = 1, nZ
       RTV%Thermal_C(i,KL) = ZERO
       DO j = 1, n_Streams
         RTV%Thermal_C(i,KL) = RTV%Thermal_C(i,KL) + (trans(i,j) + refl(i,j) )
       END DO
       IF ( i == nZ .AND. nZ == (n_Streams+1) ) THEN
         RTV%Thermal_C(i,KL) = RTV%Thermal_C(i,KL) + trans(nZ,nZ)
       END IF
       RTV%s_Layer_Source_UP(i,KL) = ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func
       RTV%s_Layer_Source_DOWN(i,KL) = RTV%s_Layer_Source_UP(i,KL)
     END DO
   END IF

   !  compute visible part for visible channels during daytime
   IF( RTV%Visible_Flag_true ) THEN
     N2 = 2 * nZ
     N2_1 = N2 - 1
     source_up = ZERO
     source_down = ZERO
     !
     ! Solar source  
     Sfactor = single_albedo*RTV%Solar_irradiance/PI
     IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO
       EXPfactor = exp(-optical_depth/RTV%COS_SUN)
       s_transmittance = exp(-total_opt/RTV%COS_SUN)

       DO i = 1, nZ     
         Solar(i) = -bb(i,nZ+1)*Sfactor
         Solar(i+nZ) = -ff(i,nZ+1)*Sfactor

         DO j = 1, nZ
           V0(i,j) = single_albedo * ff(i,j) * COS_Weight(j)
           V0(i+nZ,j) = single_albedo * bb(i,j) * COS_Weight(j)
           V0(i,j+nZ) = V0(i+nZ,j)
           V0(nZ+i,j+nZ) = V0(i,j)
         ENDDO
         V0(i,i) = V0(i,i) - ONE - COS_Angle(i)/RTV%COS_SUN
         V0(i+nZ,i+nZ) = V0(i+nZ,i+nZ) - ONE + COS_Angle(i)/RTV%COS_SUN
       ENDDO

       V1(1:N2_1,1:N2_1) = matinv(V0(1:N2_1,1:N2_1), Error_Status)
       IF( Error_Status /= SUCCESS  ) THEN
         WRITE( Message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' ) 
         CALL Display_Message( ROUTINE_NAME, &                                                    
                               TRIM(Message), &                                                   
                               Error_Status )                                          
         RETURN                                                                                    
       END IF         

       Solar1(1:N2_1) = matmul( V1(1:N2_1,1:N2_1), Solar(1:N2_1) )
       Solar1(N2) = ZERO
       Sfac2 = Solar(N2) - sum( V0(N2,1:N2_1)*Solar1(1:N2_1) )


       DO i = 1, nZ
         source_up(i) = Solar1(i)
         source_down(i) = EXPfactor*Solar1(i+nZ)
         DO j = 1, nZ
          source_up(i) =source_up(i)-refl(i,j)*Solar1(j+nZ)-trans(i,j)*EXPfactor*Solar1(j)
          source_down(i) =source_down(i) -trans(i,j)*Solar1(j+nZ) -refl(i,j)*EXPfactor*Solar1(j)
         END DO
       END DO
       ! specific treatment for downeward source function
       IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN
         source_down(nZ) =source_down(nZ) +(EXPfactor-trans(nZ,nZ))*Sfac2/V0(N2,N2)
       ELSE
         source_down(nZ) =source_down(nZ) -EXPfactor*Sfac2*optical_depth/COS_Angle(nZ)
       END IF

       source_up(1:nZ) = source_up(1:nZ)*s_transmittance
       source_down(1:nZ) = source_down(1:nZ)*s_transmittance

       RTV%s_Layer_Source_UP(1:nZ,KL) = RTV%s_Layer_Source_UP(1:nZ,KL)+source_up(1:nZ)
       RTV%s_Layer_Source_DOWN(1:nZ,KL) = RTV%s_Layer_Source_DOWN(1:nZ,KL)+source_down(1:nZ)
   END IF

   RETURN

  END SUBROUTINE CRTM_AMOM_layer
  
   SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input  number of atmospheric layers
                                 w, & ! Input  layer scattering albedo
                              T_OD, & ! Input  layer optical depth
                 cosmic_background, & ! Input  cosmic background radiance
                        emissivity, & ! Input  surface emissivity
               direct_reflectivity, & ! Input  direct reflectivity
                               RTV, & ! Input  structure containing forward part results 
              Planck_Atmosphere_TL, & ! Input  tangent-linear atmospheric layer Planck radiance 
                 Planck_Surface_TL, & ! Input  TL surface Planck radiance
                              w_TL, & ! Input  TL layer scattering albedo
                           T_OD_TL, & ! Input  TL layer optical depth
                     emissivity_TL, & ! Input  TL surface emissivity
                   reflectivity_TL, & ! Input  TL  reflectivity
            direct_reflectivity_TL, & ! Input  TL  direct reflectivity
                            Pff_TL, & ! Input  TL forward phase matrix
                            Pbb_TL, & ! Input  TL backward phase matrix
                         s_rad_up_TL) ! Output TL upward radiance 
! ------------------------------------------------------------------------- !
! FUNCTION:                                                                 !
!   This subroutine calculates IR/MW tangent-linear radiance at the top of  !
!   the atmosphere including atmospheric scattering. The structure RTV      !
!   carried in forward part results.                                        !
!   The CRTM_ADA_TL algorithm computes layer tangent-linear reflectance and !
!   transmittance as well as source function by the subroutine              !
!   CRTM_Doubling_layer as source function by the subroutine                !
!   CRTM_Doubling_layer, then uses                                          !
!   an adding method to integrate the layer and surface components.         !
!                                                                           ! 
!    Quanhua Liu    Quanhua.Liu@noaa.gov                                    !
! ------------------------------------------------------------------------- !
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n_Layers
      TYPE(RTV_type), INTENT(IN) :: RTV
      REAL (fp), INTENT(IN), DIMENSION( : ) ::  w,T_OD
      REAL (fp), INTENT(IN), DIMENSION( : ) ::  emissivity,direct_reflectivity
      REAL (fp), INTENT(IN) ::  cosmic_background

      REAL (fp),INTENT(IN),DIMENSION( :,:,: ) ::  Pff_TL, Pbb_TL
      REAL (fp),INTENT(IN),DIMENSION( : ) ::  w_TL,T_OD_TL
      REAL (fp),INTENT(IN),DIMENSION( 0: ) ::  Planck_Atmosphere_TL
      REAL (fp),INTENT(IN) ::  Planck_Surface_TL
      REAL (fp),INTENT(IN),DIMENSION( : ) ::  emissivity_TL
      REAL (fp),INTENT(IN),DIMENSION( :,: ) :: reflectivity_TL 
      REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_TL 
      REAL (fp),INTENT(INOUT),DIMENSION( : ) :: direct_reflectivity_TL
   ! -------------- internal variables --------------------------------- !
   !  Abbreviations:                                                     !
   !      s: scattering, rad: radiance, trans: transmission,             !
   !         refl: reflection, up: upward, down: downward                !
   ! --------------------------------------------------------------------!
      REAL (fp), DIMENSION(RTV%n_Angles,RTV%n_Angles) :: temporal_matrix_TL

      REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down 
      REAL (fp), DIMENSION( RTV%n_Angles ) :: s_source_up_TL,s_source_down_TL,refl_down_TL
 
      REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles ) :: s_trans_TL,s_refl_TL,Refl_Trans_TL 
      REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles ) :: s_refl_up_TL,Inv_Gamma_TL,Inv_GammaT_TL
      REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_TL
      INTEGER :: i, j, k, nZ
!
       nZ = RTV%n_Angles
       
       total_opt(0) = ZERO
       total_opt_TL(0) = ZERO
       DO k = 1, n_Layers
         total_opt(k) = total_opt(k-1) + T_OD(k)
         total_opt_TL(k) = total_opt_TL(k-1) + T_OD_TL(k)
       END DO
       
       Refl_Trans_TL = ZERO
       s_rad_up_TL = ZERO
       s_refl_up_TL = reflectivity_TL
       IF( RTV%mth_Azi == 0 ) THEN
       s_rad_up_TL = emissivity_TL * RTV%Planck_Surface + emissivity * Planck_Surface_TL
       END IF
 
       IF( RTV%Solar_Flag_true ) THEN
         s_rad_up_TL = s_rad_up_TL+direct_reflectivity_TL*RTV%COS_SUN*RTV%Solar_irradiance/PI  &
                     * exp(-total_opt(n_Layers)/RTV%COS_SUN) &
                     - direct_reflectivity * RTV%Solar_irradiance/PI  &
                     * total_opt_TL(n_Layers) * exp(-total_opt(n_Layers)/RTV%COS_SUN)
       END IF 
     
       DO 10 k = n_Layers, 1, -1
         s_source_up_TL = ZERO
         s_source_down_TL = ZERO
         s_trans_TL = ZERO
         s_refl_TL = ZERO
         Inv_GammaT_TL = ZERO
         Inv_Gamma_TL = ZERO
         refl_down_TL = ZERO
!
!      Compute tranmission and reflection matrices for a layer
      IF(w(k) > SCATTERING_ALBEDO_THRESHOLD) THEN 
       !  ----------------------------------------------------------- !
       !    CALL Doubling algorithm to computing forward and tagent   !
       !    layer transmission, reflection, and source functions.     !
       !  ----------------------------------------------------------- !

      call CRTM_AMOM_layer_TL(RTV%n_Streams,RTV%n_Angles,k,w(k),T_OD(k),total_opt(k-1), & !Input
         RTV%COS_Angle(1:RTV%n_Angles),RTV%COS_Weight(1:RTV%n_Angles)                   , & !Input
         RTV%Pff(:,:,k), RTV%Pbb(:,:,k),RTV%Planck_Atmosphere(k)                        , & !Input
         w_TL(k),T_OD_TL(k),total_opt_TL(k-1),Pff_TL(:,:,k)                             , & !Input
         Pbb_TL(:,:,k),Planck_Atmosphere_TL(k),RTV                                      , & !Input
         s_trans_TL,s_refl_TL,s_source_up_TL,s_source_down_TL)                              !Output
   
!         Adding method
         temporal_matrix_TL = -matmul(s_refl_up_TL,RTV%s_Layer_Refl(1:RTV%n_Angles,1:RTV%n_Angles,k))  &
                            - matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),s_refl_TL)

         temporal_matrix_TL = matmul(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k),temporal_matrix_TL)
         Inv_Gamma_TL = -matmul(temporal_matrix_TL,RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k))

         Inv_GammaT_TL = matmul(s_trans_TL, RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k))  &
                       + matmul(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k), Inv_Gamma_TL)

         refl_down(1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),  &
                               RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))
         refl_down_TL(:) = matmul(s_refl_up_TL,RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k)) &
                           + matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),s_source_down_TL(:))
         s_rad_up_TL(1:RTV%n_Angles)=s_source_up_TL(1:RTV%n_Angles)+ &
         matmul(Inv_GammaT_TL,refl_down(:,k)+RTV%s_Level_Rad_UP(1:RTV%n_Angles,k))  &
         +matmul(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k),refl_down_TL(1:RTV%n_Angles)+s_rad_up_TL(1:RTV%n_Angles))

         Refl_Trans_TL = matmul(s_refl_up_TL,RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k))  &
                       + matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),s_trans_TL)

         s_refl_up_TL=s_refl_TL+matmul(Inv_GammaT_TL,RTV%Refl_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k))  &
                     +matmul(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k),Refl_Trans_TL)

         Refl_Trans_TL = ZERO
         
      ELSE

         DO i = 1, RTV%n_Angles
           s_trans_TL(i,i) = -T_OD_TL(k)/RTV%COS_Angle(i) * RTV%s_Layer_Trans(i,i,k)
           s_source_up_TL(i) = Planck_Atmosphere_TL(k) * (ONE - RTV%s_Layer_Trans(i,i,k) ) &
                             - RTV%Planck_Atmosphere(k) * s_trans_TL(i,i)
           s_source_down_TL(i) = s_source_up_TL(i)
         ENDDO

!         Adding method
        DO i = 1, RTV%n_Angles 
        s_rad_up_TL(i)=s_source_up_TL(i) &
        +s_trans_TL(i,i)*(sum(RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,k)  &
        *RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))+RTV%s_Level_Rad_UP(i,k)) &
        +RTV%s_Layer_Trans(i,i,k)  &
        *(sum(s_refl_up_TL(i,1:RTV%n_Angles)*RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k)  &
        +RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,k)*s_source_down_TL(1:RTV%n_Angles))+s_rad_up_TL(i))

        ENDDO
                                       
        DO i = 1, RTV%n_Angles 
        DO j = 1, RTV%n_Angles 
        s_refl_up_TL(i,j)=s_trans_TL(i,i)*RTV%s_Level_Refl_UP(i,j,k)  &
        *RTV%s_Layer_Trans(j,j,k) &
        +RTV%s_Layer_Trans(i,i,k)*s_refl_up_TL(i,j)*RTV%s_Layer_Trans(j,j,k)  &
        +RTV%s_Layer_Trans(i,i,k)*RTV%s_Level_Refl_UP(i,j,k)*s_trans_TL(j,j)
        ENDDO
        ENDDO

      ENDIF
   10     CONTINUE
!
!  Adding reflected cosmic background radiation
    IF( RTV%mth_Azi == 0 ) THEN
      DO i = 1, RTV%n_Angles 
      s_rad_up_TL(i)=s_rad_up_TL(i)+sum(s_refl_up_TL(i,:))*cosmic_background
      ENDDO
    END IF
!
      RETURN
      END SUBROUTINE CRTM_ADA_TL
      
      SUBROUTINE CRTM_AMOM_layer_TL( n_streams, & ! Input, number of streams
                                            nZ, & ! Input, number of angles
                                            KL, & ! Input, KL-th layer 
                                 single_albedo, & ! Input, single scattering albedo
                                 optical_depth, & ! Input, layer optical depth
                                     total_opt, & ! Input, accumulated optical depth from the top to current layer top
                                     COS_Angle, & ! Input, COSINE of ANGLES
                                    COS_Weight, & ! Input, GAUSSIAN Weights
                                            ff, & ! Input, Phase matrix (forward part)
                                            bb, & ! Input, Phase matrix (backward part)
                                   Planck_Func, & ! Input, Planck for layer temperature
                              single_albedo_TL, & ! Input, tangent-linear single albedo
                              optical_depth_TL, & ! Input, TL layer optical depth
                                  total_opt_TL, & ! Input, accumulated TL optical depth from the top to current layer top
                                         ff_TL, & ! Input, TL forward Phase matrix
                                         bb_TL, & ! Input, TL backward Phase matrix
                                Planck_Func_TL, & ! Input, TL Planck for layer temperature
                                           RTV, & ! Input, structure containing forward results 
                                      trans_TL, & ! Output, layer tangent-linear trans 
                                       refl_TL, & ! Output, layer tangent-linear refl 
                                  source_up_TL, & ! Output, layer tangent-linear source_up 
                                source_down_TL)   ! Output, layer tangent-linear source_down 
                                                                
! ---------------------------------------------------------------------------------------
!   FUNCTION
!    Compute TL layer transmission, reflection matrices and source function 
!    at the top and bottom of the layer.
!
!   see also ADA method.
!   Quanhua Liu
!   Quanhua.Liu@noaa.gov
! ----------------------------------------------------------------------------------------
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: n_streams,nZ,KL
     TYPE(RTV_type), INTENT( IN ) :: RTV
     REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
     REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight 
     REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt
     !
     ! internal variables
     REAL(fp) :: s, c, s_TL,c_TL,xx_TL
     INTEGER :: i,j
     INTEGER :: Error_Status
     !
     ! Tangent-Linear Part
     REAL(fp), INTENT(OUT), DIMENSION( :,: ) :: trans_TL,refl_TL
     REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL,source_down_TL
     
     REAL(fp), INTENT(IN) :: single_albedo_TL
     REAL(fp), INTENT(IN) :: optical_depth_TL,Planck_Func_TL,total_opt_TL
     REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff_TL,bb_TL

     REAL(fp), DIMENSION(nZ) :: Exp_x_TL,EigVa_TL,EigValue_TL
     REAL(fp), DIMENSION(nZ,nZ) :: i_Gm_TL, Gm_TL, Gp_TL, EigVe_TL, EigVeF_TL, EigVeVa_TL
     REAL(fp), DIMENSION(nZ,nZ) :: A1_TL,A2_TL,A3_TL,A4_TL,A5_TL,A6_TL,Gm_A5_TL,i_Gm_A5_TL
     REAL(fp), DIMENSION(nZ,nZ) :: HH_TL,PM_TL,PP_TL,PPM_TL,i_PPM_TL,PPP_TL

     REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
     REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) 
     REAL(fp) :: EXPfactor_TL,Sfactor_TL,s_transmittance_TL,Solar_TL(2*nZ),V0_TL(2*nZ,2*nZ),Solar1_TL(2*nZ)
     REAL(fp) :: Sfac2_TL
     REAL(fp), DIMENSION( nZ ) :: thermal_up_TL,thermal_down_TL
     REAL(fp), DIMENSION(nZ,nZ) :: trans, refl
     INTEGER :: N2, N2_1
     REAL(fp) :: Thermal_C_TL
     CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer_TL'
     CHARACTER(256) :: Message
     !
     ! for small layer optical depth, single scattering is applied.  
     IF( optical_depth < DELTA_OPTICAL_DEPTH ) THEN
       s = optical_depth * single_albedo
       s_TL = optical_depth_TL * single_albedo + optical_depth * single_albedo_TL
       DO i = 1, nZ
         Thermal_C_TL = ZERO
         c = s/COS_Angle(i)
         c_TL = s_TL/COS_Angle(i)
         DO j = 1, nZ
           refl_TL(i,j) = (c_TL * bb(i,j) + c*bb_TL(i,j) )* COS_Weight(j)
           trans_TL(i,j) = (c_TL * ff(i,j) + c*ff_TL(i,j) )* COS_Weight(j)
           IF( i == j ) THEN
             trans_TL(i,j) = trans_TL(i,j) - optical_depth_TL/COS_Angle(i)
           END IF
         
           IF( RTV%mth_Azi == 0 ) THEN
             Thermal_C_TL = Thermal_C_TL + refl_TL(i,j) + trans_TL(i,j)
           END IF
         ENDDO
         IF( RTV%mth_Azi == 0 ) THEN
           source_up_TL(i) = -Thermal_C_TL * Planck_Func + &
             ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func_TL
           source_down_TL(i) = source_up_TL(i)
         END IF
       ENDDO
       RETURN
     END IF

     !
     ! for numerical stability, 
     IF( single_albedo < max_albedo ) THEN
       s = single_albedo
       s_TL = single_albedo_TL
     ELSE
       s = max_albedo
       s_TL = 0.0_fp
     END IF
     !
     ! building TL phase matrices
     DO i = 1, nZ
       c = s/COS_Angle(i)
       c_TL = s_TL/COS_Angle(i)
       DO j = 1, nZ
         PM_TL(i,j) = (c_TL * bb(i,j) + c*bb_TL(i,j) ) * COS_Weight(j)
         PP_TL(i,j) = (c_TL * ff(i,j) + c*ff_TL(i,j) ) * COS_Weight(j)
       END DO
     ENDDO

     PPM_TL(:,:) = PP_TL(:,:) - PM_TL(:,:)
     i_PPM_TL(:,:) = - matmul( RTV%i_PPM(1:nZ,1:nZ,KL), matmul(PPM_TL(:,:),RTV%i_PPM(1:nZ,1:nZ,KL)) )     
     PPP_TL(:,:) = PP_TL(:,:) + PM_TL(:,:)
     HH_TL(:,:) = matmul( PPM_TL(:,:), RTV%PPP(1:nZ,1:nZ,KL) )+matmul( RTV%PPM(1:nZ,1:nZ,KL), PPP_TL(:,:) )
     !
     ! compute TL eigenvectors EigVe, and eigenvalues EigVa
     CALL ASYMTX_TL(nZ,RTV%EigVe(1:nZ,1:nZ,KL),RTV%EigVa(1:nZ,KL),HH_TL, &
          EigVe_TL,EigVa_TL,Error_Status)

     DO i = 1, nZ
       IF( RTV%EigVa(i,KL) > ZERO ) THEN
         EigValue_TL(i) = 0.5_fp*EigVa_TL(i)/RTV%EigValue(i,KL)
       ELSE
         EigValue_TL(i) = ZERO 
       END IF
     END DO
     EigVeVa_TL = ZERO
     
     DO i = 1, nZ
       DO j = 1, nZ
         EigVeVa_TL(i,j) = EigVe_TL(i,j) * RTV%EigValue(j,KL)+RTV%EigVe(i,j,KL) * EigValue_TL(j)
       END DO
     END DO     
     EigVeF_TL(:,:) = matmul( i_PPM_TL(:,:), RTV%EigVeVa(1:nZ,1:nZ,KL) )  &
                    + matmul( RTV%i_PPM(1:nZ,1:nZ,KL), EigVeVa_TL(:,:) )
     !  
     ! compute TL reflection and transmission matrices, TL source function     
     Gp_TL(:,:) = ( EigVe_TL(:,:) + EigVeF_TL(:,:) )/2.0_fp
     Gm_TL(:,:) = ( EigVe_TL(:,:) - EigVeF_TL(:,:) )/2.0_fp
     i_Gm_TL = -matmul( RTV%i_Gm(1:nZ,1:nZ,KL), matmul(Gm_TL,RTV%i_Gm(1:nZ,1:nZ,KL)) )
     DO i = 1, nZ
       xx_TL = EigValue_TL(i)*optical_depth+RTV%EigValue(i,KL)*optical_depth_TL
       Exp_x_TL(i) = -xx_TL*RTV%Exp_x(i,KL)
     END DO

     DO i = 1, nZ
       DO j = 1, nZ
         A1_TL(i,j) = Gp_TL(i,j)* RTV%Exp_x(j,KL)+ RTV%Gp(i,j,KL)* Exp_x_TL(j)
         A4_TL(i,j) = Gm_TL(i,j)* RTV%Exp_x(j,KL)+ RTV%Gm(i,j,KL)* Exp_x_TL(j)
       END DO
     END DO        
     A2_TL(:,:) = matmul(i_Gm_TL(:,:),RTV%A1(1:nZ,1:nZ,KL))+matmul(RTV%i_Gm(1:nZ,1:nZ,KL),A1_TL(:,:))          
     A3_TL(:,:) = matmul(Gp_TL(:,:),RTV%A2(1:nZ,1:nZ,KL))+matmul(RTV%Gp(1:nZ,1:nZ,KL),A2_TL(:,:))
     A5_TL(:,:) = matmul(A1_TL(:,:),RTV%A2(1:nZ,1:nZ,KL))+matmul(RTV%A1(1:nZ,1:nZ,KL),A2_TL(:,:))
     A6_TL(:,:) = matmul(A4_TL(:,:),RTV%A2(1:nZ,1:nZ,KL))+matmul(RTV%A4(1:nZ,1:nZ,KL),A2_TL(:,:))
  
     Gm_A5_TL(:,:) = Gm_TL(:,:) - A5_TL(:,:)
     i_Gm_A5_TL(:,:) = -matmul( RTV%i_Gm_A5(1:nZ,1:nZ,KL),matmul(Gm_A5_TL,RTV%i_Gm_A5(1:nZ,1:nZ,KL)))
     !
     ! T = matmul( RTV%A4(:,:,KL) - RTV%A3(:,:,KL), RTV%i_Gm_A5(:,:,KL) )
     trans_TL = matmul( A4_TL(:,:) - A3_TL(:,:), RTV%i_Gm_A5(1:nZ,1:nZ,KL) )  &
          + matmul( RTV%A4(1:nZ,1:nZ,KL) - RTV%A3(1:nZ,1:nZ,KL), i_Gm_A5_TL(:,:) )
     refl_TL = matmul( Gp_TL(:,:) - A6_TL(:,:), RTV%i_Gm_A5(1:nZ,1:nZ,KL) )  &
          + matmul( RTV%Gp(1:nZ,1:nZ,KL) - RTV%A6(1:nZ,1:nZ,KL), i_Gm_A5_TL(:,:) )

     trans(1:nZ,1:nZ) = RTV%s_Layer_Trans(1:nZ,1:nZ,KL)
     refl(1:nZ,1:nZ) = RTV%s_Layer_Refl(1:nZ,1:nZ,KL)
     Source_UP_TL = ZERO
     Source_DOWN_TL = ZERO    
     !
     ! Thermal part
     IF( RTV%mth_Azi == 0 ) THEN  
       DO i = 1, nZ
         Thermal_C_TL = ZERO
         DO j = 1, n_Streams 
           Thermal_C_TL = Thermal_C_TL + (trans_TL(i,j) + refl_TL(i,j))
         ENDDO
         IF(i == nZ .AND. nZ == (n_Streams+1)) THEN
           Thermal_C_TL = Thermal_C_TL + trans_TL(nZ,nZ)
         ENDIF
         thermal_up_TL(i) = -Thermal_C_TL * Planck_Func  &
           + ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func_TL
         thermal_down_TL(i) = thermal_up_TL(i)
       ENDDO
     END IF
     
     !
     ! for visible channels at daytime
     IF( RTV%Visible_Flag_true ) THEN 
       N2 = 2 * nZ
       N2_1 = N2 - 1
       V0 = ZERO
       V1 = ZERO
       Solar = ZERO
       Solar1 = ZERO
       Sfac2 = ZERO
       V0_TL = ZERO
       Solar_TL = ZERO
       Solar1_TL = ZERO
       Sfac2_TL = ZERO
       !
       ! Solar source  
       Sfactor = single_albedo*RTV%Solar_irradiance/PI
       Sfactor_TL = single_albedo_TL*RTV%Solar_irradiance/PI
               
       IF( RTV%mth_Azi == 0 ) THEN
         Sfactor = Sfactor/TWO
         Sfactor_TL = Sfactor_TL/TWO
       END IF
       EXPfactor = exp(-optical_depth/RTV%COS_SUN)
       EXPfactor_TL = -optical_depth_TL/RTV%COS_SUN*EXPfactor
      
       s_transmittance = exp(-total_opt/RTV%COS_SUN)
       s_transmittance_TL = -total_opt_TL/RTV%COS_SUN*s_transmittance
       
       DO i = 1, nZ     
         Solar(i) = -bb(i,nZ+1)*Sfactor
         Solar_TL(i) = -bb_TL(i,nZ+1)*Sfactor-bb(i,nZ+1)*Sfactor_TL
         Solar(i+nZ) = -ff(i,nZ+1)*Sfactor
         Solar_TL(i+nZ) = -ff_TL(i,nZ+1)*Sfactor-ff(i,nZ+1)*Sfactor_TL
         DO j = 1, nZ
           V0(i,j) = single_albedo * ff(i,j) * COS_Weight(j)
           V0_TL(i,j) = single_albedo_TL*ff(i,j)*COS_Weight(j)+single_albedo*ff_TL(i,j)*COS_Weight(j)
           V0(i+nZ,j) = single_albedo * bb(i,j) * COS_Weight(j)
           V0_TL(i+nZ,j) = single_albedo_TL*bb(i,j)*COS_Weight(j)+single_albedo*bb_TL(i,j)*COS_Weight(j)
           V0(i,j+nZ) = V0(i+nZ,j)
           V0_TL(i,j+nZ) = V0_TL(i+nZ,j)
           V0(nZ+i,j+nZ) = V0(i,j)
           V0_TL(nZ+i,j+nZ) = V0_TL(i,j)
         ENDDO
         V0(i,i) = V0(i,i) - ONE - COS_Angle(i)/RTV%COS_SUN
         V0(i+nZ,i+nZ) = V0(i+nZ,i+nZ) - ONE + COS_Angle(i)/RTV%COS_SUN
       ENDDO

       V1(1:N2_1,1:N2_1) = matinv(V0(1:N2_1,1:N2_1), Error_Status)
       IF( Error_Status /= SUCCESS  ) THEN
         WRITE( Message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' ) 
         CALL Display_Message( ROUTINE_NAME, &                                                    
                               TRIM(Message), &                                                   
                               Error_Status )                                          
         RETURN                                                                                    
       END IF           
       
       Solar1(1:N2_1) = matmul( V1(1:N2_1,1:N2_1), Solar(1:N2_1) )
       
       Solar(1:N2_1) =  matmul( V1(1:N2_1,1:N2_1),Solar(1:N2_1) )
       Solar1_TL(1:N2_1) = matmul( V0_TL(1:N2_1,1:N2_1),Solar(1:N2_1) )
       Solar1_TL(1:N2_1) = -matmul(  V1(1:N2_1,1:N2_1),Solar1_TL(1:N2_1) )  &
                         + matmul( V1(1:N2_1,1:N2_1), Solar_TL(1:N2_1) )
       
       Solar1(N2) = ZERO
       Solar1_TL(N2) = ZERO
       Sfac2 = Solar(N2) - sum( V0(N2,1:N2_1)*Solar1(1:N2_1) )
       Sfac2_TL = Solar_TL(N2) - sum( V0_TL(N2,1:N2_1)*Solar1(1:N2_1) )  &
                - sum( V0(N2,1:N2_1)*Solar1_TL(1:N2_1) )
               
       DO i = 1, nZ
         source_up(i) = Solar1(i)
         source_up_TL(i) = Solar1_TL(i)
         source_down(i) = EXPfactor*Solar1(i+nZ)
         source_down_TL(i) = EXPfactor_TL*Solar1(i+nZ)+EXPfactor*Solar1_TL(i+nZ)
         DO j = 1, nZ
           source_up(i) =source_up(i)-refl(i,j)*Solar1(j+nZ)-trans(i,j)*EXPfactor*Solar1(j)
           source_up_TL(i) =source_up_TL(i)-refl_TL(i,j)*Solar1(j+nZ) -refl(i,j)*Solar1_TL(j+nZ) &
           - trans_TL(i,j)*EXPfactor*Solar1(j) - trans(i,j)*EXPfactor_TL*Solar1(j) -trans(i,j)*EXPfactor*Solar1_TL(j)
           source_down(i) =source_down(i) -trans(i,j)*Solar1(j+nZ) -refl(i,j)*EXPfactor*Solar1(j)
           source_down_TL(i) =source_down_TL(i) -trans_TL(i,j)*Solar1(j+nZ) -trans(i,j)*Solar1_TL(j+nZ) &
           -refl_TL(i,j)*EXPfactor*Solar1(j) -refl(i,j)*EXPfactor_TL*Solar1(j) -refl(i,j)*EXPfactor*Solar1_TL(j)
         END DO
       END DO
       !
       ! specific treatment for downeward source function
       IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN
         source_down(nZ) =source_down(nZ) +(EXPfactor-trans(nZ,nZ))*Sfac2/V0(N2,N2)
         source_down_TL(nZ) =source_down_TL(nZ) +(EXPfactor_TL-trans_TL(nZ,nZ))*Sfac2/V0(N2,N2) &
          +(EXPfactor-trans(nZ,nZ))*Sfac2_TL/V0(N2,N2)-(EXPfactor-trans(nZ,nZ))*Sfac2*V0_TL(N2,N2)/V0(N2,N2)/V0(N2,N2)
       ELSE
         source_down(nZ) =source_down(nZ) -EXPfactor*Sfac2*optical_depth/COS_Angle(nZ)
         source_down_TL(nZ) =source_down_TL(nZ) -EXPfactor_TL*Sfac2*optical_depth/COS_Angle(nZ)  &
         -EXPfactor*Sfac2_TL*optical_depth/COS_Angle(nZ)-EXPfactor*Sfac2*optical_depth_TL/COS_Angle(nZ)
       END IF
        
       ! source_up(1:nZ) = source_up(1:nZ)*s_transmittance
        source_up_TL(1:nZ) = source_up_TL(1:nZ)*s_transmittance+source_up(1:nZ)*s_transmittance_TL
       ! source_down(1:nZ) = source_down(1:nZ)*s_transmittance
        source_down_TL(1:nZ) = source_down_TL(1:nZ)*s_transmittance + source_down(1:nZ)*s_transmittance_TL
     END IF

     source_up_TL(:) = source_up_TL(:) + thermal_up_TL(:)
     source_down_TL(:) = source_down_TL(:) + thermal_down_TL(:)
    
     RETURN

     END SUBROUTINE CRTM_AMOM_layer_TL      
  
   SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input  number of atmospheric layers
                                 w, & ! Input  layer scattering albedo
                              T_OD, & ! Input  layer optical depth
                 cosmic_background, & ! Input  cosmic background radiance
                        emissivity, & ! Input  surface emissivity
               direct_reflectivity, & ! surface direct reflectivity
                               RTV, & ! Input  structure containing forward results 
                       s_rad_up_AD, & ! Input  adjoint upward radiance 
              Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance
                 Planck_Surface_AD, & ! Output AD surface Planck radiance
                              w_AD, & ! Output AD layer scattering albedo
                           T_OD_AD, & ! Output AD layer optical depth
                     emissivity_AD, & ! Output AD surface emissivity
                   reflectivity_AD, & ! Output AD surface reflectivity
            direct_reflectivity_AD, & ! Output AD surface direct reflectivity
                            Pff_AD, & ! Output AD forward phase matrix
                            Pbb_AD)   ! Output AD backward phase matrix
! ------------------------------------------------------------------------- !
! FUNCTION:                                                                 !
!   This subroutine calculates IR/MW adjoint radiance at the top of         !
!   the atmosphere including atmospheric scattering. The structure RTV      !
!   carried in forward part results.                                        !
!   The CRTM_ADA_AD algorithm computes layer tangent-linear reflectance and !
!   transmittance as well as source function by the subroutine              !
!   CRTM_Doubling_layer as source function by the subroutine                !
!   CRTM_Doubling_layer, then uses                                          !
!   an adding method to integrate the layer and surface components.         !
!                                                                           ! 
!    Quanhua Liu    Quanhua.Liu@noaa.gov                                    !
! ------------------------------------------------------------------------- !
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n_Layers
      TYPE(RTV_type), INTENT(IN) :: RTV
      REAL (fp), INTENT(IN), DIMENSION( : ) ::  w,T_OD
      REAL (fp), INTENT(IN), DIMENSION( : ) ::  emissivity,direct_reflectivity
      REAL (fp), INTENT(IN) ::  cosmic_background

      REAL (fp),INTENT(INOUT),DIMENSION( :,:,: ) ::  Pff_AD, Pbb_AD
      REAL (fp),INTENT(INOUT),DIMENSION( : ) ::  w_AD,T_OD_AD
      REAL (fp),INTENT(INOUT),DIMENSION( 0: ) ::  Planck_Atmosphere_AD
      REAL (fp),INTENT(INOUT) ::  Planck_Surface_AD
      REAL (fp),INTENT(INOUT),DIMENSION( : ) ::  emissivity_AD,direct_reflectivity_AD
      REAL (fp),INTENT(INOUT),DIMENSION( :,: ) :: reflectivity_AD 
      REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_AD 
   ! -------------- internal variables --------------------------------- !
   !  Abbreviations:                                                     !
   !      s: scattering, rad: radiance, trans: transmission,             !
   !         refl: reflection, up: upward, down: downward                !
   ! --------------------------------------------------------------------!
      REAL (fp), DIMENSION(RTV%n_Angles,RTV%n_Angles) :: temporal_matrix_AD

      REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down 
      REAL (fp), DIMENSION( RTV%n_Angles ) :: s_source_up_AD,s_source_down_AD,refl_down_AD
 
      REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles) :: s_trans_AD,s_refl_AD,Refl_Trans_AD
      REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles) :: s_refl_up_AD,Inv_Gamma_AD,Inv_GammaT_AD
      REAL (fp) :: sum_s_AD, sums_AD, xx
      REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_AD
      INTEGER :: i, j, k,nZ
!
      nZ = RTV%n_Angles
      
       total_opt_AD = ZERO
       total_opt(0) = ZERO
       DO k = 1, n_Layers
         total_opt(k) = total_opt(k-1) + T_OD(k)
       END DO
!
       s_trans_AD = ZERO
       Planck_Atmosphere_AD = ZERO
       Planck_Surface_AD = ZERO
       s_refl_up_AD = ZERO

      Pff_AD = ZERO
      Pbb_AD = ZERO
!      T_OD_AD = ZERO
!  Adding reflected cosmic background radiation

      DO i = 1, RTV%n_Angles 
      sum_s_AD = s_rad_up_AD(i)*cosmic_background
      DO j = 1, RTV%n_Angles 
      s_refl_up_AD(i,j) = sum_s_AD
      ENDDO
      ENDDO

!
       DO 10 k = 1, n_Layers 
       s_source_up_AD = ZERO
       s_source_down_AD = ZERO
       s_trans_AD = ZERO
!
!      Compute tranmission and reflection matrices for a layer
      IF(w(k) > SCATTERING_ALBEDO_THRESHOLD) THEN 

        refl_down(1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),  &
                               RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))

        s_refl_AD = s_refl_up_AD
        Inv_GammaT_AD = matmul(s_refl_up_AD,transpose(RTV%Refl_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)))
        Refl_Trans_AD = matmul(transpose(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k)),s_refl_up_AD)

        s_refl_up_AD=matmul(Refl_Trans_AD,transpose(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)))
        s_trans_AD=matmul(transpose(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k)),Refl_Trans_AD)
 
        s_source_up_AD(1:RTV%n_Angles) = s_rad_up_AD(1:RTV%n_Angles)

        DO i = 1, RTV%n_Angles 
        sums_AD = s_rad_up_AD(i)
        DO j = 1, RTV%n_Angles 
        Inv_GammaT_AD(i,j)=Inv_GammaT_AD(i,j)+sums_AD*(refl_down(j,k)+RTV%s_Level_Rad_UP(j,k))
        ENDDO 
        ENDDO

        refl_down_AD(1:RTV%n_Angles)=matmul(transpose(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k)),s_rad_up_AD(1:RTV%n_Angles))
        s_rad_up_AD(1:RTV%n_Angles)=matmul(transpose(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k)),s_rad_up_AD(1:RTV%n_Angles))

        DO i = 1, RTV%n_Angles 
        sums_AD = refl_down_AD(i)
        DO j = 1, RTV%n_Angles 
        s_refl_up_AD(i,j)=s_refl_up_AD(i,j)+sums_AD*RTV%s_Layer_Source_DOWN(j,k)
        ENDDO 
        ENDDO

        s_source_down_AD=matmul(transpose(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k)),refl_down_AD(:)) 

        s_trans_AD=s_trans_AD+matmul(Inv_GammaT_AD,transpose(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k)))
        Inv_Gamma_AD= matmul(transpose(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)),Inv_GammaT_AD)

        temporal_matrix_AD= -matmul(Inv_Gamma_AD,transpose(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k)))
        temporal_matrix_AD=matmul(transpose(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k)),temporal_matrix_AD)

        s_refl_up_AD=s_refl_up_AD-matmul(temporal_matrix_AD,transpose(RTV%s_Layer_Refl(1:RTV%n_Angles,1:RTV%n_Angles,k)))
        s_refl_AD=s_refl_AD-matmul(transpose(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k)),temporal_matrix_AD)
       !  ----------------------------------------------------------- !
       !    CALL Doubling algorithm to computing forward and tagent   !
       !    layer transmission, reflection, and source functions.     !
       !  ----------------------------------------------------------- !

      call CRTM_AMOM_layer_AD(RTV%n_Streams,RTV%n_Angles,k,w(k),T_OD(k),total_opt(k-1)      , & !Input
         RTV%COS_Angle,RTV%COS_Weight,RTV%Pff(:,:,k),RTV%Pbb(:,:,k),RTV%Planck_Atmosphere(k), & !Input
         s_trans_AD,s_refl_AD,s_source_up_AD,s_source_down_AD,RTV,w_AD(k),T_OD_AD(k)        , &
         total_opt_AD(k-1),Pff_AD(:,:,k),Pbb_AD(:,:,k),Planck_Atmosphere_AD(k))  !Output

      ELSE
        DO i = 1, RTV%n_Angles 
        DO j = 1, RTV%n_Angles 
        s_trans_AD(j,j)=s_trans_AD(j,j)+RTV%s_Layer_Trans(i,i,k)*RTV%s_Level_Refl_UP(i,j,k)*s_refl_up_AD(i,j)
        s_trans_AD(i,i)=s_trans_AD(i,i)+s_refl_up_AD(i,j)*RTV%s_Level_Refl_UP(i,j,k)*RTV%s_Layer_Trans(j,j,k)
        s_refl_up_AD(i,j)=RTV%s_Layer_Trans(i,i,k)*s_refl_up_AD(i,j)*RTV%s_Layer_Trans(j,j,k)
        ENDDO
        ENDDO
!         Adding method
       DO i = 1, RTV%n_Angles 

         s_source_up_AD(i)=s_rad_up_AD(i)
         s_trans_AD(i,i)=s_trans_AD(i,i)+s_rad_up_AD(i)*(sum(RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,k)  &
         *RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))+RTV%s_Level_Rad_UP(i,k))

         sum_s_AD=RTV%s_Layer_Trans(i,i,k)*s_rad_up_AD(i)
         DO j = 1, RTV%n_Angles 
          s_refl_up_AD(i,j)=s_refl_up_AD(i,j)+sum_s_AD*RTV%s_Layer_Source_DOWN(j,k)
         ENDDO
                                            
         sum_s_AD=RTV%s_Layer_Trans(i,i,k)*s_rad_up_AD(i)
         DO j = 1, RTV%n_Angles 
          s_source_down_AD(j)=s_source_down_AD(j)+sum_s_AD*RTV%s_Level_Refl_UP(i,j,k)
         ENDDO
         s_rad_up_AD(i)=RTV%s_Layer_Trans(i,i,k)*s_rad_up_AD(i)

       ENDDO

         DO i = RTV%n_Angles, 1, -1
           s_source_up_AD(i) = s_source_up_AD(i) +  s_source_down_AD(i)
           s_trans_AD(i,i) = s_trans_AD(i,i) - RTV%Planck_Atmosphere(k) * s_source_up_AD(i)
           Planck_Atmosphere_AD(k) = Planck_Atmosphere_AD(k) + s_source_up_AD(i) * (ONE - RTV%s_Layer_Trans(i,i,k) )
           s_source_up_AD(i) = ZERO

           T_OD_AD(k)=T_OD_AD(k)-s_trans_AD(i,i)/RTV%COS_Angle(i)*RTV%s_Layer_Trans(i,i,k)
         ENDDO

      ENDIF
   10     CONTINUE

! 

       IF( RTV%Solar_Flag_true ) THEN         
         xx = RTV%Solar_irradiance/PI * exp(-total_opt(n_Layers)/RTV%COS_SUN)
         total_opt_AD(n_Layers) = total_opt_AD(n_Layers)  &
            - xx*sum(direct_reflectivity(1:RTV%n_Angles)*s_rad_up_AD(1:RTV%n_Angles))
         direct_reflectivity_AD = direct_reflectivity_AD + s_rad_up_AD * RTV%COS_SUN * xx
       END IF 
 
        emissivity_AD = s_rad_up_AD * RTV%Planck_Surface
        Planck_Surface_AD = sum(emissivity(:) * s_rad_up_AD(:) )

     
      reflectivity_AD = s_refl_up_AD
!
       s_rad_up_AD = ZERO
       s_refl_up_AD = ZERO
 
       
       DO k = n_Layers, 1, -1
         T_OD_AD(k) = T_OD_AD(k) + total_opt_AD(k)
         total_opt_AD(k-1) = total_opt_AD(k-1) + total_opt_AD(k)
       END DO
       
      RETURN
      END SUBROUTINE CRTM_ADA_AD

!
!
     SUBROUTINE CRTM_AMOM_layer_AD(  n_streams, & ! Input, number of streams
                                            nZ, & ! Input, number of angles
                                            KL, & ! Input, KL-th layer 
                                 single_albedo, & ! Input, single scattering albedo
                                 optical_depth, & ! Input, layer optical depth
                                     total_opt, & ! Input, 
                                     COS_Angle, & ! Input, COSINE of ANGLES
                                    COS_Weight, & ! Input, GAUSSIAN Weights
                                            ff, & ! Input, Phase matrix (forward part)
                                            bb, & ! Input, Phase matrix (backward part)
                                   Planck_Func, & ! Input, Planck for layer temperature
                                      trans_AD, & ! Input, layer tangent-linear trans 
                                       refl_AD, & ! Input, layer tangent-linear refl 
                                  source_up_AD, & ! Input, layer tangent-linear source_up 
                                source_down_AD, & ! Input, layer tangent-linear source_down 
                                           RTV, & ! Input, structure containing forward results 
                              single_albedo_AD, & ! Output adjoint single scattering albedo
                              optical_depth_AD, & ! Output AD layer optical depth
                                  total_opt_AD, & ! Output AD accumulated optical depth ftom TOA to current layer top
                                         ff_AD, & ! Output AD forward Phase matrix
                                         bb_AD, & ! Output AD backward Phase matrix
                                Planck_Func_AD)   ! Output AD Planck for layer temperature
! ---------------------------------------------------------------------------------------
!   FUNCTION
!    Compute AD layer transmission, reflection matrices and source function 
!    at the top and bottom of the layer.
!
!   see also ADA method.
!   Quanhua Liu
!   Quanhua.Liu@noaa.gov
! ----------------------------------------------------------------------------------------
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: n_streams,nZ,KL
     TYPE(RTV_type), INTENT( IN ) :: RTV
     REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
     REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight 
     REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt

     ! internal variables
     REAL(fp) :: s, c, s_AD,c_AD,xx_AD
     INTEGER :: i,j
     INTEGER :: Error_Status

     ! Tangent-Linear Part
     REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD
     REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD
     REAL(fp), INTENT( INOUT ) :: single_albedo_AD
     REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD,total_opt_AD
     REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD

     REAL(fp), DIMENSION(nZ) :: Exp_x_AD,EigVa_AD,EigValue_AD
     REAL(fp), DIMENSION(nZ,nZ) :: i_Gm_AD, Gm_AD, Gp_AD, EigVe_AD, EigVeF_AD, EigVeVa_AD
     REAL(fp), DIMENSION(nZ,nZ) :: A1_AD,A2_AD,A3_AD,A4_AD,A5_AD,A6_AD,Gm_A5_AD,i_Gm_A5_AD
     REAL(fp), DIMENSION(nZ,nZ) :: HH_AD,PM_AD,PP_AD,PPM_AD,i_PPM_AD,PPP_AD
    
     REAL(fp), DIMENSION(nZ) :: thermal_up_AD, thermal_down_AD
     REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
     REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) 
     REAL(fp) :: EXPfactor_AD,Sfactor_AD,s_transmittance_AD,Solar_AD(2*nZ),V0_AD(2*nZ,2*nZ),Solar1_AD(2*nZ)
     REAL(fp) :: V1_AD(2*nZ,2*nZ),Sfac2_AD
     REAL(fp), DIMENSION(nZ,nZ) :: trans, refl
     INTEGER :: N2, N2_1
     REAL(fp) :: Thermal_C_AD
     CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer_AD'
     CHARACTER(256) :: Message
          
     s_AD = ZERO
     c_AD = ZERO
     !
     ! for small layer optical depth, single scattering is applied.  
     IF( optical_depth < DELTA_OPTICAL_DEPTH ) THEN
       s = optical_depth * single_albedo
       DO i = 1, nZ         
         c = s/COS_Angle(i)
         IF( RTV%mth_Azi == 0 ) THEN
           source_up_AD(i) = source_up_AD(i) + source_down_AD(i)
           source_down_AD(i) = ZERO
           Planck_Func_AD = Planck_Func_AD + (ONE - RTV%Thermal_C(i,KL))*source_up_AD(i)
           Thermal_C_AD = -source_up_AD(i) * Planck_Func
         END IF

         DO j = 1, nZ   
           IF( RTV%mth_Azi == 0 ) THEN
             refl_AD(i,j) = refl_AD(i,j) + Thermal_C_AD
             trans_AD(i,j) = trans_AD(i,j) + Thermal_C_AD
           END IF
           IF( i == j ) THEN
             optical_depth_AD = optical_depth_AD - trans_AD(i,j)/COS_Angle(i)
           END IF
            
           c_AD = c_AD + trans_AD(i,j) * ff(i,j) * COS_Weight(j)
           ff_AD(i,j) = ff_AD(i,j) + c * trans_AD(i,j) * COS_Weight(j)
           c_AD = c_AD + refl_AD(i,j) * bb(i,j) * COS_Weight(j)
           bb_AD(i,j) = bb_AD(i,j) + c * refl_AD(i,j) * COS_Weight(j)
         ENDDO

         source_up_AD(i) = ZERO
         s_AD = s_AD + c_AD/COS_Angle(i)
         c_AD = ZERO
       ENDDO
       optical_depth_AD = optical_depth_AD + s_AD * single_albedo
       single_albedo_AD = single_albedo_AD + optical_depth * s_AD
       RETURN
     END IF

     trans(1:nZ,1:nZ) = RTV%s_Layer_Trans(1:nZ,1:nZ,KL)
     refl(1:nZ,1:nZ) = RTV%s_Layer_Refl(1:nZ,1:nZ,KL)

     thermal_up_AD(:) = source_up_AD(:)
     thermal_down_AD(:) = source_down_AD(:)

     IF( RTV%Visible_Flag_true ) THEN 
       N2 = 2 * nZ
       N2_1 = N2 - 1

       ! forward part  start   ********
       source_up = ZERO
       source_down = ZERO
       Solar_AD = ZERO
       Solar1_AD = ZERO
       Sfactor_AD = ZERO
       Sfac2_AD = ZERO
       EXPfactor_AD = ZERO
       s_transmittance_AD = ZERO
       V0_AD = ZERO
       V1_AD = ZERO
       !
       ! Solar source  
       Sfactor = single_albedo*RTV%Solar_irradiance/PI
       IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO
       EXPfactor = exp(-optical_depth/RTV%COS_SUN)
       s_transmittance = exp(-total_opt/RTV%COS_SUN)
       
       DO i = 1, nZ     
         Solar(i) = -bb(i,nZ+1)*Sfactor
         Solar(i+nZ) = -ff(i,nZ+1)*Sfactor

         DO j = 1, nZ
           V0(i,j) = single_albedo * ff(i,j) * COS_Weight(j)
           V0(i+nZ,j) = single_albedo * bb(i,j) * COS_Weight(j)
           V0(i,j+nZ) = V0(i+nZ,j)
           V0(nZ+i,j+nZ) = V0(i,j)
         ENDDO
       V0(i,i) = V0(i,i) - ONE - COS_Angle(i)/RTV%COS_SUN
       V0(i+nZ,i+nZ) = V0(i+nZ,i+nZ) - ONE + COS_Angle(i)/RTV%COS_SUN
       ENDDO
       
       V1(1:N2_1,1:N2_1) = matinv(V0(1:N2_1,1:N2_1), Error_Status)
       IF( Error_Status /= SUCCESS  ) THEN
         WRITE( Message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' ) 
         CALL Display_Message( ROUTINE_NAME, &                                                    
                               TRIM(Message), &                                                   
                               Error_Status )                                          
         RETURN                                                                                    
       END IF               
       
       Solar1(1:N2_1) = matmul( V1(1:N2_1,1:N2_1), Solar(1:N2_1) )
       Solar1(N2) = ZERO
       Sfac2 = Solar(N2) - sum( V0(N2,1:N2_1)*Solar1(1:N2_1) )

       DO i = 1, nZ
         source_up(i) = Solar1(i)
         source_down(i) = EXPfactor*Solar1(i+nZ)
         DO j = 1, nZ
           source_up(i) =source_up(i)-refl(i,j)*Solar1(j+nZ)-trans(i,j)*EXPfactor*Solar1(j)
           source_down(i) =source_down(i) -trans(i,j)*Solar1(j+nZ) -refl(i,j)*EXPfactor*Solar1(j)
         END DO
       END DO
       ! specific treatment for downeward source function
       IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN
         source_down(nZ) =source_down(nZ) +(EXPfactor-trans(nZ,nZ))*Sfac2/V0(N2,N2)
       ELSE
         source_down(nZ) =source_down(nZ) -EXPfactor*Sfac2*optical_depth/COS_Angle(nZ)
       END IF
        
       ! forward part end  ********      
       !
       s_transmittance_AD = s_transmittance_AD+sum (source_down(1:nZ)*source_down_AD(1:nZ) )
       source_down_AD(1:nZ) = source_down_AD(1:nZ)*s_transmittance
       s_transmittance_AD = s_transmittance_AD + sum (source_up(1:nZ)*source_up_AD(1:nZ) )
       source_up_AD(1:nZ) = source_up_AD(1:nZ)*s_transmittance
       !
       ! specific treatment for downeward source function
       IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN
         V0_AD(N2,N2)=V0_AD(N2,N2)-(EXPfactor-trans(nZ,nZ))*Sfac2*source_down_AD(nZ)/V0(N2,N2)/V0(N2,N2)
         Sfac2_AD = Sfac2_AD+(EXPfactor-trans(nZ,nZ))*source_down_AD(nZ)/V0(N2,N2)
         EXPfactor_AD = EXPfactor_AD+source_down_AD(nZ)*Sfac2/V0(N2,N2)
         trans_AD(nZ,nZ) = trans_AD(nZ,nZ)-source_down_AD(nZ)*Sfac2/V0(N2,N2)
       ELSE
         optical_depth_AD = optical_depth_AD -EXPfactor*Sfac2*source_down_AD(nZ)/COS_Angle(nZ)
         Sfac2_AD = Sfac2_AD-EXPfactor*source_down_AD(nZ)*optical_depth/COS_Angle(nZ)
         EXPfactor_AD = EXPfactor_AD-source_down_AD(nZ)*Sfac2*optical_depth/COS_Angle(nZ)
       END IF

       DO i = nZ, 1, -1
         DO j = nZ, 1, -1
           Solar1_AD(j)=Solar1_AD(j)-refl(i,j)*EXPfactor*source_down_AD(i)
           EXPfactor_AD = EXPfactor_AD -refl(i,j)*source_down_AD(i)*Solar1(j)
           refl_AD(i,j) = refl_AD(i,j) -source_down_AD(i)*EXPfactor*Solar1(j)
           Solar1_AD(j+nZ) = Solar1_AD(j+nZ) -trans(i,j)*source_down_AD(i)
           trans_AD(i,j) = trans_AD(i,j) -source_down_AD(i)*Solar1(j+nZ)
           
           Solar1_AD(j)=Solar1_AD(j)-trans(i,j)*EXPfactor*source_up_AD(i)
           EXPfactor_AD = EXPfactor_AD - trans(i,j)*source_up_AD(i)*Solar1(j)
           trans_AD(i,j)=trans_AD(i,j) - source_up_AD(i)*EXPfactor*Solar1(j)
           Solar1_AD(j+nZ) = Solar1_AD(j+nZ) -refl(i,j)*source_up_AD(i)
           refl_AD(i,j) = refl_AD(i,j) -source_up_AD(i)*Solar1(j+nZ)
         END DO
         
         Solar1_AD(i+nZ) = Solar1_AD(i+nZ) + EXPfactor * source_down_AD(i)
         EXPfactor_AD = EXPfactor_AD + source_down_AD(i)*Solar1(i+nZ)
         Solar1_AD(i) = Solar1_AD(i) + source_up_AD(i)          
       END DO

       Solar1_AD(1:N2_1)=Solar1_AD(1:N2_1) -Sfac2_AD*V0(N2,1:N2_1)
       V0_AD(N2,1:N2_1)=V0_AD(N2,1:N2_1) -Sfac2_AD*Solar1(1:N2_1)
       Solar_AD(N2) = Solar_AD(N2) + Sfac2_AD

       Solar1_AD(N2) = ZERO
       Solar_AD(1:N2_1)=Solar_AD(1:N2_1)+matmul( transpose(V1(1:N2_1,1:N2_1)),Solar1_AD(1:N2_1) )
       Solar(1:N2_1) =  matmul( V1(1:N2_1,1:N2_1),Solar(1:N2_1) )                  
       Solar1_AD(1:N2_1) = -matmul( transpose(V1(1:N2_1,1:N2_1)),Solar1_AD(1:N2_1) ) 
       DO i = 1, N2_1
       DO j = 1, N2_1
         V0_AD(i,j)=V0_AD(i,j)+Solar1_AD(i)*Solar(j)
       END DO
       END DO        
        
       ! Solar source            
       DO i = nZ, 1, -1                
         DO j = nZ, 1, -1    
           V0_AD(i,j)=V0_AD(i,j) + V0_AD(nZ+i,j+nZ)
           V0_AD(i+nZ,j)=V0_AD(i+nZ,j) + V0_AD(i,j+nZ)
           bb_AD(i,j)=bb_AD(i,j) + single_albedo*V0_AD(i+nZ,j)*COS_Weight(j)
           single_albedo_AD=single_albedo_AD + V0_AD(i+nZ,j)*bb(i,j)*COS_Weight(j)
           ff_AD(i,j)=ff_AD(i,j) + single_albedo*V0_AD(i,j)*COS_Weight(j)
           single_albedo_AD=single_albedo_AD +V0_AD(i,j)*ff(i,j)*COS_Weight(j)
         ENDDO
       
         Sfactor_AD = Sfactor_AD -ff(i,nZ+1)*Solar_AD(i+nZ)
         ff_AD(i,nZ+1) = ff_AD(i,nZ+1) -Solar_AD(i+nZ)*Sfactor
         Sfactor_AD = Sfactor_AD -bb(i,nZ+1)*Solar_AD(i)
         bb_AD(i,nZ+1)=bb_AD(i,nZ+1) - Solar_AD(i)*Sfactor
       ENDDO
        
       total_opt_AD = total_opt_AD -s_transmittance_AD/RTV%COS_SUN*s_transmittance
       optical_depth_AD = optical_depth_AD -EXPfactor_AD/RTV%COS_SUN*EXPfactor
       
       IF( RTV%mth_Azi == 0 ) THEN
         Sfactor_AD = Sfactor_AD/TWO
       END IF
             
       single_albedo_AD = single_albedo_AD + Sfactor_AD*RTV%Solar_irradiance/PI
       
     END IF

    ! Thermal part
     IF( RTV%mth_Azi == 0 ) THEN
       DO i = nZ, 1, -1
         thermal_up_AD(i) = thermal_up_AD(i) + thermal_down_AD(i)
         thermal_down_AD(i) = ZERO
         Planck_Func_AD = Planck_Func_AD + ( ONE - RTV%Thermal_C(i,KL) ) * thermal_up_AD(i)
         Thermal_C_AD = -thermal_up_AD(i) * Planck_Func

         IF ( i == nZ .AND. nZ == (n_Streams+1) ) THEN
           trans_AD(nZ,nZ) = trans_AD(nZ,nZ) + Thermal_C_AD
         END IF
       
         DO j = n_Streams, 1, -1
           trans_AD(i,j) = trans_AD(i,j) + Thermal_C_AD
           refl_AD(i,j) = refl_AD(i,j) + Thermal_C_AD
         ENDDO
         thermal_up_AD(i) = ZERO
       ENDDO
     END IF
!     
     i_Gm_A5_AD = matmul( transpose(RTV%Gp(1:nZ,1:nZ,KL)-RTV%A6(1:nZ,1:nZ,KL)),refl_AD )
     Gp_AD(:,:) = matmul( refl_AD, transpose(RTV%i_Gm_A5(1:nZ,1:nZ,KL)) )       
     
     A6_AD = - GP_AD  
     i_Gm_A5_AD = i_Gm_A5_AD + matmul( transpose(RTV%A4(1:nZ,1:nZ,KL)-RTV%A3(1:nZ,1:nZ,KL)),trans_AD )
     A4_AD(:,:) = matmul( trans_AD, transpose(RTV%i_Gm_A5(1:nZ,1:nZ,KL)) )
     A3_AD = - A4_AD
     Gm_A5_AD = -matmul( transpose(RTV%i_Gm_A5(1:nZ,1:nZ,KL)) ,matmul( i_Gm_A5_AD,transpose(RTV%i_Gm_A5(1:nZ,1:nZ,KL)) ) )       
     Gm_AD = Gm_A5_AD
     A5_AD = - Gm_A5_AD

     A4_AD = A4_AD + matmul( A6_AD(:,:), transpose(RTV%A2(1:nZ,1:nZ,KL)) )
     A2_AD = matmul( transpose(RTV%A4(1:nZ,1:nZ,KL)),A6_AD(:,:) )

     A1_AD = matmul( A5_AD(:,:), transpose(RTV%A2(1:nZ,1:nZ,KL)) )
     A2_AD = A2_AD + matmul( transpose(RTV%A1(1:nZ,1:nZ,KL)), A5_AD(:,:) )

     Gp_AD = Gp_AD + matmul( A3_AD(:,:), transpose(RTV%A2(1:nZ,1:nZ,KL)) )
     A2_AD = A2_AD + matmul( transpose(RTV%Gp(1:nZ,1:nZ,KL)),A3_AD(:,:) )

     i_Gm_AD = matmul( A2_AD(:,:), transpose(RTV%A1(1:nZ,1:nZ,KL)) )
     A1_AD = A1_AD + matmul( transpose(RTV%i_Gm(1:nZ,1:nZ,KL)), A2_AD(:,:) )

     Exp_x_AD = ZERO
     
     DO i = nZ, 1, -1
       DO j = nZ, 1, -1
         Gm_AD(i,j) = Gm_AD(i,j) + A4_AD(i,j)* RTV%Exp_x(j,KL)
         Exp_x_AD(j) = Exp_x_AD(j) + RTV%Gm(i,j,KL)*A4_AD(i,j)
         Gp_AD(i,j) = Gp_AD(i,j) + A1_AD(i,j)* RTV%Exp_x(j,KL)
         Exp_x_AD(j) = Exp_x_AD(j) + RTV%Gp(i,j,KL)*A1_AD(i,j)
       END DO
     END DO
    
     DO i = nZ, 1, -1
       xx_AD = -Exp_x_AD(i)*RTV%Exp_x(i,KL)
       Exp_x_AD(i) = ZERO
       EigValue_AD(i) = xx_AD*optical_depth
       optical_depth_AD = optical_depth_AD + RTV%EigValue(i,KL)*xx_AD
     END DO

     Gm_AD = Gm_AD -matmul( transpose(RTV%i_Gm(1:nZ,1:nZ,KL)), matmul( i_Gm_AD, transpose(RTV%i_Gm(1:nZ,1:nZ,KL)) ) )
 
     EigVe_AD(:,:) = Gm_AD(:,:)/2.0_fp           
     EigVeF_AD(:,:) = - Gm_AD(:,:)/2.0_fp         

     EigVe_AD = EigVe_AD + Gp_AD(:,:)/2.0_fp
     EigVeF_AD = EigVeF_AD + Gp_AD(:,:)/2.0_fp

     i_PPM_AD(:,:) = matmul( EigVeF_AD(:,:), transpose(RTV%EigVeVa(1:nZ,1:nZ,KL)) )
     EigVeVa_AD(:,:) = matmul( transpose(RTV%i_PPM(1:nZ,1:nZ,KL)), EigVeF_AD(:,:) )           

     DO i = nZ, 1, -1
       DO j = nZ, 1, -1              
         EigVe_AD(i,j)=EigVe_AD(i,j)+EigVeVa_AD(i,j)* RTV%EigValue(j,KL)
         EigValue_AD(j) = EigValue_AD(j)+RTV%EigVe(i,j,KL)*EigVeVa_AD(i,j)
       END DO
     END DO
    
     DO i = nZ, 1, -1
       IF( RTV%EigVa(i,KL) > ZERO ) THEN
         EigVa_AD(i) = 0.5_fp*EigValue_AD(i)/RTV%EigValue(i,KL)
       ELSE
         EigValue_AD(i) = ZERO
         EigVa_AD(i) = ZERO 
       END IF
     END DO

     ! compute eigenvectors EigVe, and eigenvalues EigVa
     CALL ASYMTX_AD(nZ,RTV%EigVe(1:nZ,1:nZ,KL),RTV%EigVa(1:nZ,KL), &
          EigVe_AD,EigVa_AD,HH_AD,Error_Status) 

     PPM_AD(:,:) = matmul( HH_AD(:,:), transpose(RTV%PPP(1:nZ,1:nZ,KL)) )
     PPP_AD(:,:) = matmul( transpose(RTV%PPM(1:nZ,1:nZ,KL)), HH_AD(:,:) )

     PP_AD = PPP_AD
     PM_AD = PPP_AD
   
     PPM_AD(:,:) = PPM_AD(:,:)-matmul( transpose(RTV%i_PPM(1:nZ,1:nZ,KL)),matmul(i_PPM_AD(:,:),transpose(RTV%i_PPM(1:nZ,1:nZ,KL))) )

     PP_AD = PP_AD + PPM_AD
     PM_AD = PM_AD - PPM_AD

     IF( single_albedo < max_albedo ) THEN
       s = single_albedo
     ELSE
       s = max_albedo
     END IF 
     
       c_AD = ZERO
       s_AD = ZERO
     DO i = nZ, 1, -1
       c = s/COS_Angle(i) 
       DO j = nZ, 1, -1
       c_AD = c_AD + PP_AD(i,j) * ff(i,j) * COS_Weight(j)
       ff_AD(i,j) = ff_AD(i,j) + c * PP_AD(i,j) * COS_Weight(j)
       c_AD = c_AD + PM_AD(i,j) * bb(i,j) * COS_Weight(j)
       bb_AD(i,j) = bb_AD(i,j) + c * PM_AD(i,j) * COS_Weight(j)
       END DO
       s_AD = s_AD + c_AD/COS_Angle(i)
       c_AD = ZERO
     ENDDO
!
     IF( single_albedo < max_albedo ) THEN
       s = single_albedo
       single_albedo_AD = s_AD + single_albedo_AD       
     ELSE
       s = max_albedo
       s_AD = 0.0_fp
     END IF
!       
     RETURN

     END SUBROUTINE CRTM_AMOM_layer_AD
  
  
END MODULE ADA_Module
