
MODULE MINOUT

  IMPLICIT NONE

  ! suppress those attribute for f2py
  ! PRIVATE :: SKIP_COM, &
  !      READ_COM, &
  !      UPCASE, &
  !      LOWCASE, &
  !      COUNTCOL, &
  !      COUNTLINE
  ! PUBLIC :: READ_DATA, WRITE_DATA

CONTAINS

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

  SUBROUTINE READ_DATA 

    USE MCOMMON 
    USE MCONSTANTS
    USE MUTILITY, only : TRIMCAT, REVERSE, LININTPOL
    USE nr, ONLY : SORT

    IMPLICIT NONE

    ! local variables
    CHARACTER (LEN=max_len) :: filename_tmp
    CHARACTER (LEN=max_len) :: filename_kdis
    CHARACTER (LEN=max_len) :: the_char
    CHARACTER (LEN=max_len) :: the_char2
    LOGICAL                 :: test_file
    LOGICAL                 :: read_prop 

    INTEGER :: i, k, j, l, n

    INTEGER :: isp, ip, it, ic

    TYPE (FICH) :: finput        = FICH (''        , 11)
    TYPE (FICH) :: fatm          = FICH ('atm_'                  , 12)
    TYPE (FICH) :: fopt          = FICH ('opt_'                  , 13)
    TYPE (FICH) :: fad_spec      = FICH ('ad_spec.dat'           , 14)
    TYPE (FICH) :: fod_spec      = FICH ('od_spec.dat'           , 15)
    TYPE (FICH) :: fmcrad1d_spec = FICH ('mcrad1d_spec.dat'      , 16)
    TYPE (FICH) :: fpotter_spec  = FICH ('potter_spec.dat'       , 17)
    TYPE (FICH) :: fdfit_spec    = FICH ('dfit_spec.dat'         , 18)
    TYPE (FICH) :: focean_spec   = FICH ('ocean_spec.dat'        , 19)
    TYPE (FICH) :: fdoad_spec    = FICH ('doad_spec.dat'         , 20)
    TYPE (FICH) :: fprop         = FICH ('prop_'                 , 21)
    TYPE (FICH) :: frefind       = FICH ('refind_'               , 22)
    TYPE (FICH) :: fsurfalb      = FICH ('surfalb_'              , 23)
    TYPE (FICH) :: fgasconc      = FICH ('uniform_gases_conc.dat', 24)
    TYPE (FICH) :: fkdis         = FICH ('kdis_'                 , 25)
    TYPE (FICH) :: fsolrad       = FICH ('solrad_'               , 26)
    TYPE (FICH) :: fvdist        = FICH ('vdist_'                , 27)
    TYPE (FICH) :: ffilter       = FICH ('filter_'               , 28)
    TYPE (FICH) :: fsurfbrdf     = FICH ('surfbrdf_'             , 29)
    TYPE (FICH) :: fbaum_de      = FICH ('_baum_diameters.dat'   , 30)
    type (FICH) :: fbaum_mu      = FICH ('_baum_mu.dat'          , 31)
    TYPE (FICH) :: fbaum_wvl     = FICH ('_baum_wvl.dat'         , 32)
    TYPE (FICH) :: fbaum         = FICH ('baum_cirrus_'          , 33)

    INTEGER :: nlamb
    INTEGER :: nelem
    REAL(kind=dp), allocatable :: lamb(:)
    INTEGER, allocatable :: nteta(:)
    REAL(kind=dp), allocatable :: salb(:)

    REAL(kind=dp), allocatable :: klamb(:)
    REAL(kind=dp), allocatable :: k2(:)
    REAL(kind=dp), allocatable :: k1(:)
    REAL(kind=dp), allocatable :: cv(:)
    REAL(kind=dp), allocatable :: refind_r(:)
    REAL(kind=dp), allocatable :: refind_i(:)

    REAL(kind=dp) :: cste_alb

    LOGICAL, ALLOCATABLE:: flag_kdiswvl(:)

    LOGICAL :: flag
    INTEGER :: imaterial

    INTEGER :: ntmp                                 
    CHARACTER(LEN=max_len), ALLOCATABLE :: char_arr_tmp(:) 
    REAL(kind=dp), ALLOCATABLE :: dbl_arr_tmp(:) 
    INTEGER, ALLOCATABLE :: logic_arr_tmp(:) 

    logical :: flag_user_vdist

    CHARACTER (LEN=max_len) :: command
    LOGICAL :: dir_exist
    INTEGER :: st

    character(len=8)  :: char8
    character(len=10) :: char10
    character(len=5)  :: char5

    INTEGER :: iwvl
    INTEGER :: ide

    CHARACTER(LEN=max_len) :: baum_type
    INTEGER :: baum_nlamb_file
    INTEGER :: baum_nde_file
    REAL(kind=dp), allocatable :: baum_lamb_file(:)
    REAL(kind=dp), allocatable :: baum_de_file(:)
    INTEGER :: baum_nide
    INTEGER :: baum_nilamb
    integer, allocatable :: baum_ide(:)
    integer, allocatable :: baum_ilamb(:)

    REAL(kind=dp) :: tmp_dp

    LOGICAL :: flag_chnstrtmp
    LOGICAL :: flag_chnstr

    CHARACTER (LEN=max_len) :: artdeco_path_local

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

    CALL CPU_TIME(start_tcpu)
    CALL DATE_AND_TIME(char8, char10, char5, start_date)

    ! init keywords 
    verbose        = .FALSE.
    opt_only       = .FALSE.
    betal_only     = .FALSE.
    thermal        = .FALSE.
    do_rt          = .TRUE.
    sinsca_corint      = .FALSE.
    sinsca_corint_brdf = .FALSE.
    warning        = .TRUE.
    no_rayleigh    = .FALSE.
    thermal_only   = .FALSE. 
    flux_only      = .FALSE.
    print_betal    = .FALSE.
    print_recomp   = .FALSE.
    print_aitaui   = .FALSE.
    od_no_check    = .FALSE.
    is_secsca      = .FALSE.
    bottom_f0      = .FALSE.

    ! init some values
    flag_user_vdist = .FALSE.
    mode         = undef_c
    trunc_method = undef_c
    rt_model     = undef_c
    nmat         = undef_i
    hg_flag      = .FALSE.
    flag_mie     = .FALSE.
    flag_hm      = .FALSE.
    flag_baum    = .FALSE.

    !------------
    ! start reading main input file


    finput%name  = artdeco_infile_name
    filename_tmp = artdeco_infile_name
    ! filename_tmp = TRIMCAT (artdeco_path, dir_input)
    ! filename_tmp = TRIMCAT (filename_tmp, finput%name)

    OPEN (UNIT = finput%unit, FILE = filename_tmp, STATUS = 'old')

    the_char = READ_COM (finput%unit)
    the_char = LOWCASE(the_char, len(the_char)) 
    ! get global keywords 
    IF (INDEX(the_char, 'verbose') .ne. 0) verbose = .true.
    IF (INDEX(the_char, 'opt_only') .ne. 0) opt_only = .true.
    IF (INDEX(the_char, 'betal_only') .ne. 0) betal_only = .true.
    IF (INDEX(the_char, 'nowarn') .ne. 0) warning = .false.
    IF (INDEX(the_char, 'no_rayleigh') .ne. 0) no_rayleigh = .true.
    IF (INDEX(the_char, 'thermal_only') .ne. 0) thermal_only = .true.
    IF (INDEX(the_char, 'flux_only') .ne. 0) flux_only = .true.
    IF (INDEX(the_char, 'print_betal') .ne. 0) print_betal = .true.
    IF (INDEX(the_char, 'print_recomp') .ne. 0) print_recomp = .true.
    IF (INDEX(the_char, 'print_aitaui') .ne. 0) print_aitaui = .true.
    IF (INDEX(the_char, 'od_no_check') .ne. 0) od_no_check = .true.

    IF (betal_only .or. opt_only) do_rt = .FALSE.

    IF (verbose) THEN
       WRITE (*,*) '==========================================='
       WRITE (*,*) '           Running ARTDECO ', version
       WRITE (*,*) ''
       WRITE (*,*)'  (sub. read_data) read artdeco_in.dat'
    ENDIF

    IF (verbose) THEN
       WRITE(*,*) ''
       WRITE(*,*)' >>>>>>               keywords : ', TRIM(ADJUSTL(the_char))
    ENDIF

    ! read artdeco_path 
    artdeco_path = READ_COM(finput%unit)
    artdeco_path = TRIM (ADJUSTL(artdeco_path))
    if (verbose) then
       WRITE(*,*) ' '
       WRITE(*,*) ' >>>>>>               ARTDECO PATH : ', artdeco_path
    ENDIF

    ! read root name for output file
    out_root = READ_COM(finput%unit)
    out_root = TRIM (ADJUSTL(out_root))
    if (verbose) then
       WRITE(*,*) ' '
       WRITE(*,*) ' >>>>>>               Output root name  : ', TRIM (ADJUSTL(out_root))
    ENDIF

    ! get the working mode 
    CALL SKIP_COM(finput%unit)
    READ (UNIT=finput%unit, FMT=*) mode
    mode = LOWCASE ( TRIM(ADJUSTL(mode)), len((TRIM(ADJUSTL(mode)))) )

    if ((mode.ne.'mono') .and. (mode.ne.'kdis')) then
       write(*,*) ''
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       WRITE(*,*) ' (in sub. read_data)  : ERROR '
       write(*,*) ' mode must be "mono" or "kdis"'
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       write(*,*) ''
       STOP
    endif
    IF (verbose) THEN
       if (mode.eq.'mono') then
          WRITE(*,*) ''
          WRITE(*,*)' >>>>>>               mode : ', TRIM(ADJUSTL(mode))
          if (do_rt.and.warning) then
             write(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : WARNING                  '
             write(*,*) '       In monochromatic mode, ARTDECO does not   '
             write(*,*) '       account for gas absorption (and emission) '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          endif
       endif
    ENDIF
    if (mode.eq.'kdis')then
       BACKSPACE(finput%unit)
       READ (UNIT=finput%unit, FMT=*) mode, kdis_model
       mode       = LOWCASE ( TRIM(ADJUSTL(mode)), len(TRIM(ADJUSTL(mode))) )
       kdis_model = LOWCASE ( TRIM(ADJUSTL(kdis_model)), len(TRIM(ADJUSTL(kdis_model))) )
       IF (verbose) THEN
          WRITE(*,*) ''
          WRITE(*,FMT='(1x,A29,1x,A6,1x, A30)')' >>>>>>               mode : ', &
               TRIM(ADJUSTL(mode)), TRIM(ADJUSTL(kdis_model))
       ENDIF
    endif

    filter = .TRUE.
    ! --- read filter description
    call SKIP_COM(finput%unit)
    ntmp = COUNTLINE(finput%unit,finput%name)
    ALLOCATE(filter_name(ntmp))
    DO i = 1, ntmp
       READ (finput%unit, FMT=*) filter_name(i)
       filter_name(i) = LOWCASE (TRIM(ADJUSTL(filter_name(i))), len(TRIM(ADJUSTL(filter_name(i)))))
       if (filter_name(i) .eq. 'none') then
          nfilter = 0
          filter = .FALSE.
       endif
    END DO

    IF (filter) then
       write(*,*) ''
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       WRITE(*,*) ' (in sub. read_data)  : ERROR '
       write(*,*) '                        No filter use anymore '
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       write(*,*) ''
       STOP
    end IF
                 
    if (do_rt.eqv..FALSE.) then
       ! no need for filters
       ! in that case
       nfilter = 0
       filter  = .FALSE.
       DEALLOCATE(filter_name)
    else
       IF (filter) then
          if (mode.ne.'kdis') then
             write(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             write(*,*) '  You must use kdis mode if '
             write(*,*) '  you use filter transmission'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             write(*,*) ''
             STOP
          endif
          nfilter = ntmp
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A,1x,A)')      ' >>>>>>               filter(s) : ', TRIM(ADJUSTL(filter_name(1)))
             if (nfilter.ge.2) then
                DO i = 2, nfilter
                   WRITE(*,FMT='(1x,A,1x,A)')' >>>>>>                           ', TRIM(ADJUSTL(filter_name(i)))
                ENDDO
             endif
          ENDIF
       ELSE
          DEALLOCATE(filter_name)
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A)')' >>>>>>               no filter '
          ENDIF
       END IF
    endif

    ! get working wavelength (mono)
    ! or the bondary wavelengths (kdis)
    lambda_bound(:) = undef_dp 
    CALL SKIP_COM(finput%unit)
    if (filter) then
       ! we skip the wavelength description because
       ! it will be given by the filter
       READ (UNIT=finput%unit, FMT=*) the_char
    else
       if (mode.eq.'mono') then
          nlambda = COUNTCOL(finput%unit)
          ALLOCATE(wlambda(nlambda))
          wlambda(:)     = undef_dp
          READ (UNIT=finput%unit, FMT=*) (wlambda(i), i=1, nlambda)
          CALL SORT(wlambda)
          DO I = 1, nlambda
             IF(wlambda(i).lt.0.0_dp) then
                write(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR                     '
                write(*,*) '    Wavelengths set in artdeco_in.dat            '
                write(*,*) '    must be > 0.0                                 '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) ''
                STOP
             END IF
          ENDDO
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A35,1x,f10.5)')' >>>>>>               wavelengths : ', wlambda(1)
             if (nlambda.ge.2) then
                DO I = 2, nlambda
                   WRITE(*,FMT='(37x,f10.5)') wlambda(i)
                ENDDO
             endif
          ENDIF
       else
          if (COUNTCOL(finput%unit).ne.2) then
             write(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR                     '
             write(*,*) '    Two wavelengths must be provided  in          '
             write(*,*) '    artdeco_in.dat in kdis mode                   '
             write(*,*) '    --> limit wavelengths to work with            '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             write(*,*) ''
             STOP
          endif
          READ (UNIT=finput%unit, FMT=*) lambda_bound(1), lambda_bound(2)
          if ((lambda_bound(1).ne.-1.0_dp) .and. (lambda_bound(2).ne.-1.0_dp)) CALL SORT(lambda_bound)
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A44,1x,f10.5,3x,f10.5)')' >>>>>>               bondary wavelengths : ', &
                  & lambda_bound(1), lambda_bound(2)
          ENDIF
       end if
    end if

    ! --- read particle description
    call SKIP_COM(finput%unit)
    ntmp = COUNTLINE(finput%unit,finput%name)
    nptcle = ntmp
    ALLOCATE(char_arr_tmp(ntmp))
    char_arr_tmp(:) = undef_c
    DO i = 1, ntmp
       READ (finput%unit, FMT=*) char_arr_tmp(i)    
       char_arr_tmp(i) = LOWCASE (TRIM(ADJUSTL(char_arr_tmp(i))), len(TRIM(ADJUSTL(char_arr_tmp(i)))))
       if (i.gt.1) then
          DO j = 1, i-1
             if(char_arr_tmp(i).eq.char_arr_tmp(j)) then
                write(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR                     '
                write(*,*) '   Problem with the particle list in artdeco_in.dat '
                write(*,*) '   The particle ', LOWCASE (TRIM(ADJUSTL(char_arr_tmp(i))), len(TRIM(ADJUSTL(char_arr_tmp(i)))))        
                write(*,*) '   appears more than once '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) ''
                STOP
             end if
          ENDDO
       endif
       if (char_arr_tmp(i).eq.'none') then
          nptcle = 0
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A34)')' >>>>>>               no particle '
          ENDIF
       endif
    ENDDO
    IF (nptcle .gt. 0) THEN
       ALLOCATE(ptcle_type(nptcle),       &
            ptcle_hg(nptcle),             &
            ptcle_opt_interp(nptcle),     &
            ptcle_vdist_type(nptcle),     &
            ptcle_h_vdist(nptcle),        &
            ptcle_sd_gauss_vdist(nptcle), &
            ptcle_tau_reflamb(nptcle), &
            ptcle_user_vdist_nalt(nptcle),& 
            ptcle_user_vdist_alt(nptcle,max_user_vdist_nalt), &
            ptcle_user_vdist(nptcle,max_user_vdist_nalt) )
       ! init arrays
       ptcle_type(:)                 = undef_c
       ptcle_hg(:)                   = .FALSE.
       ptcle_opt_interp(:)           = .TRUE.    
       ptcle_vdist_type(:)           = undef_c
       ptcle_h_vdist(:)              = undef_dp
       ptcle_sd_gauss_vdist(:)       = undef_dp
       ptcle_tau_reflamb(:)          = undef_dp
       ptcle_user_vdist_nalt(:)      = undef_i
       ptcle_user_vdist_alt(:,:)     = undef_dp
       ptcle_user_vdist(:,:)         = undef_dp
       IF (verbose) THEN
          WRITE(*,*) ''
          WRITE(*,FMT='(1x,A35,1x,A30)')      ' >>>>>>               particle(s) : ', TRIM(ADJUSTL(char_arr_tmp(1)))
          if (nptcle.ge.2) then
             DO i=2, nptcle
                WRITE(*,FMT='(1x,A35,1x,A30)')' >>>>>>                             ', TRIM(ADJUSTL(char_arr_tmp(i)))
             ENDDO
          endif
       ENDIF

       IF (OPT_ONLY) THEN
          ! we only need the particle type
          DO i=1, nptcle
             ptcle_type(i) = char_arr_tmp(i)
             ptcle_opt_interp(i) = .FALSE. 
          ENDDO
       ENDIF

       IF (BETAL_ONLY) then
          ! back to the first ptcle line to read
          DO i=1, nptcle
             BACKSPACE(finput%unit)
          END DO
          DO i=1, nptcle
             READ (UNIT=finput%unit,FMT=*) ptcle_type(i), the_char

             ptcle_type(i)       = LOWCASE(TRIM(ADJUSTL(ptcle_type(i))), len(TRIM(ADJUSTL(ptcle_type(i)))))             

             the_char  = LOWCASE(TRIM(ADJUSTL(the_char)), len(TRIM(ADJUSTL(the_char))))

             if (the_char.eq.'yes')  then
                ptcle_opt_interp(i) = .TRUE.
             elseif (the_char.eq.'no') then
                ptcle_opt_interp(i) = .FALSE.
             else
                write(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  :  ERROR                      '
                write(*,*) '   In artdeco_in.dat for the particle definition :  '
                write(*,*) '   Whether we interpolate or not optical properties '
                write(*,*) '   is undefined for particle', &
                     LOWCASE(TRIM(ADJUSTL(ptcle_type(i))), len(TRIM(ADJUSTL(ptcle_type(i)))))        
                write(*,*) '   "interp." should be "yes" or "no"                '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) ''
                STOP
             endif
          END DO
       ENDIF

       IF (DO_RT) THEN
          ! back to the first ptcle line to read
          DO i=1, nptcle
             BACKSPACE(finput%unit)
          END DO
          DO i=1, nptcle
             READ (UNIT=finput%unit,FMT=*) ptcle_type(i), the_char, the_char2, ptcle_tau_reflamb(i), ptcle_vdist_type(i)

             ptcle_type(i)       = LOWCASE(TRIM(ADJUSTL(ptcle_type(i))), len(TRIM(ADJUSTL(ptcle_type(i)))) )             
             ptcle_vdist_type(i) = LOWCASE(TRIM(ADJUSTL(ptcle_vdist_type(i))), len(TRIM(ADJUSTL(ptcle_vdist_type(i)))))

             the_char  = LOWCASE(TRIM(ADJUSTL(the_char)), len(TRIM(ADJUSTL(the_char))))
             the_char2 = LOWCASE(TRIM(ADJUSTL(the_char2)), len(TRIM(ADJUSTL(the_char2))))

             if (the_char.eq.'yes')  then
                ptcle_opt_interp(i) = .TRUE.
             elseif (the_char.eq.'no') then
                ptcle_opt_interp(i) = .FALSE.
             else
                write(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  :  ERROR                      '
                write(*,*) '   In artdeco_in.dat for the particle definition :  '
                write(*,*) '   Whether we interpolate or not optical properties '
                write(*,*) '   is undefined for particle', &
                     LOWCASE(TRIM(ADJUSTL(ptcle_type(i))), len(TRIM(ADJUSTL(ptcle_type(i))))) 
                write(*,*) '   "interp." should be "yes" or "no"                '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) ''
                STOP
             endif

             if (the_char2.eq.'yes')  then
                ptcle_hg(i) = .TRUE.
             elseif (the_char2.eq.'no') then
                ptcle_hg(i) = .FALSE.
             else
                write(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  :  ERROR                      '
                write(*,*) '   In artdeco_in.dat for the particle definition :  '
                write(*,*) '   Whether we approximate the phase function with   '
                write(*,*) '   a H-G one is undefined for particle', &
                     LOWCASE(TRIM(ADJUSTL(ptcle_type(i))), len(TRIM(ADJUSTL(ptcle_type(i))))) 
                write(*,*) '   "H-G" should be "yes" or "no"                    '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) ''
                STOP
             endif

             if (ptcle_hg(i)) hg_flag = .TRUE.

             BACKSPACE(finput%unit)
             SELECT CASE (ptcle_vdist_type(i))
             CASE('gauss')
                READ (UNIT=finput%unit,FMT=*) ptcle_type(i), the_char, the_char2, ptcle_tau_reflamb(i), &
                     ptcle_vdist_type(i), &
                     ptcle_h_vdist(i), ptcle_sd_gauss_vdist(i)
             CASE('sh')
                ptcle_sd_gauss_vdist(i) = undef_dp
                READ (UNIT=finput%unit,FMT=*) ptcle_type(i), the_char, the_char2, ptcle_tau_reflamb(i), &
                     ptcle_vdist_type(i), &
                     ptcle_h_vdist(i)
             CASE('homogeneous')
                ptcle_sd_gauss_vdist(i) = undef_dp
                READ (UNIT=finput%unit,FMT=*) ptcle_type(i), the_char, the_char2, ptcle_tau_reflamb(i), &
                     ptcle_vdist_type(i), &
                     ptcle_h_vdist(i), ptcle_sd_gauss_vdist(i)
             CASE('layer')
                ptcle_sd_gauss_vdist(i) = undef_dp
                READ (UNIT=finput%unit,FMT=*) ptcle_type(i), the_char, the_char2, ptcle_tau_reflamb(i), &
                     ptcle_vdist_type(i), &
                     ptcle_h_vdist(i)
             CASE('user')
                flag_user_vdist = .TRUE.
                ptcle_sd_gauss_vdist(i) = undef_dp
                ptcle_h_vdist(i)        = undef_dp
                READ (UNIT=finput%unit,FMT=*) the_char            
             END SELECT

             ptcle_type(i)       = LOWCASE(TRIM(ADJUSTL(ptcle_type(i))), len(TRIM(ADJUSTL(ptcle_type(i))))) 
             ptcle_vdist_type(i) = LOWCASE(TRIM(ADJUSTL(ptcle_vdist_type(i))), len(TRIM(ADJUSTL(ptcle_vdist_type(i)))))

          ENDDO
       ENDIF
    ENDIF ! nptcle > 0

    DEALLOCATE(char_arr_tmp)

    IF (opt_only .eqv. .FALSE.) THEN
       !-----------------------------
       ! truncation method to be used
       the_char = READ_COM (finput%unit)
       trunc_method = LOWCASE (TRIM(ADJUSTL(the_char)), len(TRIM(ADJUSTL(the_char))))
       IF ( trunc_method .ne. 'none' .and. &
            trunc_method .ne. 'dfit' .and. &
            trunc_method .ne. 'potter' .and. &
            trunc_method .ne. 'dm') THEN
          WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          WRITE(*,*) ' (in sub. read_data)  : ERROR '
          WRITE(*,*) ' The truncation method must be "none", "dfit" "dm" or "potter"'
          WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          WRITE(*,*) ''
          STOP 
       ENDIF
       IF (verbose) THEN
          WRITE(*,*) ''
          WRITE(*,*)' >>>>>>               truncation method : ', TRIM(ADJUSTL(trunc_method))
       ENDIF

       ! read the number of betal to decompose the phase matrix on
       CALL SKIP_COM(finput%unit)
       READ(finput%unit, FMT=*) nbetal_in
       IF (verbose) THEN
          WRITE(*,*) ''
          WRITE(*,FMT='(1x, A30, 3x, i5)')' >>>>>>               nbetal : ', nbetal_in
       ENDIF

       IF  (do_rt) THEN

          ! radiative transfer code to be used
          the_char = READ_COM (finput%unit)
          rt_model = LOWCASE (TRIM(ADJUSTL(the_char)), len(TRIM(ADJUSTL(the_char))))
          IF ( rt_model .ne. 'sinsca' .and.  & 
               rt_model .ne. 'disort' .and.  & 
               rt_model .ne. 'doad'   .and.  & 
               rt_model .ne. 'addoub' .and.  & 
               rt_model .ne. 'mcrad1d') THEN
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' The radiative transfer model must be '
             WRITE(*,*) ' "sinsca", "disort", "mcrad1d", "addoub" or "doad"'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ELSE
             IF (verbose) THEN
                WRITE(*,*) ''
                WRITE(*,*)' >>>>>>               RTE solver : ', TRIM(ADJUSTL(rt_model))
             ENDIF
          ENDIF
          if((rt_model.eq.'mcrad1d').and.(mode.eq.'kdis'))then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' K-distribution treatment is not implemented with the mcrad1d RTE solver'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif
          if((rt_model.ne.'disort').and.(flux_only))then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' "flux_only" option only available with "disort" RTE solver'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif

          ! TMS (single scattering correction) ?
          the_char = READ_COM (finput%unit)
          the_char = LOWCASE (the_char, len(the_char)) 
          the_char = TRIM(ADJUSTL(the_char))
          IF ( the_char .ne. 'yes' .and.   & 
               the_char .ne. 'no' ) THEN
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' Correct value for -TMS correction- must be "yes" or "no" '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          IF (INDEX(the_char, 'yes') .ne. 0) THEN
             sinsca_corint      =.true.
             sinsca_corint_brdf =.true.
          else
             sinsca_corint     =.false.
             sinsca_corint_brdf=.false.
          END IF
          if ( (sinsca_corint)       .and. &
               (rt_model.eq.'disort').and. &
               (trunc_method.eq.'potter')) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR                          '
             WRITE(*,*) '    You can not apply TMS after computation with       '
             WRITE(*,*) '    DISORT RTE solver and Potter truncation.           '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          endif
          if ((sinsca_corint).and.(nbetal_in.ne.-1)) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  Intensity correction can not be correctly performed  '
             WRITE(*,*) '  if the nbetal is not set to default (= -1 )'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          endif
          ! if ((sinsca_corint) .and. (nptcle.eq.0)) then
          !    if (verbose.and.warning) then
          !       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          !       WRITE(*,*) ' (in sub. read_data)  : WARNING '
          !       WRITE(*,*) '  No intensity correction could be performed when '
          !       WRITE(*,*) '  no particle are present.'
          !       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          !       WRITE(*,*) ''
          !    endif
          !    sinsca_corint = .false.
          ! end if
          ! if ( (trunc_method .eq. 'none') .and. (sinsca_corint .eqv. .true.)) then
          !    if (verbose.and.warning) then
          !       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          !       WRITE(*,*) ' (in sub. read_data)  : WARNING '
          !       WRITE(*,*) '  No intensity correction could be performed when '
          !       WRITE(*,*) '  no truncation of the phase matrix is done.'
          !       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          !       WRITE(*,*) ''
          !    endif
          !    sinsca_corint = .false.
          ! endif
          if ( (rt_model .eq. 'sinsca') .and. (sinsca_corint .eqv. .true.)) then
             if (verbose.and.warning) then
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : WARNING '
                WRITE(*,*) '  No intensity correction could be performed when '
                WRITE(*,*) '  using "sinsca" to solve RTE.'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
             endif
             sinsca_corint = .false.
          endif
          if ( (rt_model .eq. 'mcrad1d') .and. (sinsca_corint .eqv. .true.)) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  No intensity correction is performed when '
             WRITE(*,*) '  using "mcrad1d" to solve RTE.'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          endif
          if ( (rt_model .eq. 'addoub') .and. (sinsca_corint .eqv. .true.)) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  No intensity correction is performed when '
             WRITE(*,*) '  using "addoub" to solve RTE.'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          endif

          ! do we account for thermal IR emission in-situ
          the_char = READ_COM (finput%unit)
          the_char = LOWCASE (the_char, len(the_char)) 
          the_char = TRIM(ADJUSTL(the_char))
          IF ( the_char .ne. 'yes' .and.   & 
               the_char .ne. 'no' ) THEN
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' Correct value for -thermal- must be "yes" or "no" '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          IF (INDEX(the_char, 'yes') .ne. 0) THEN
             thermal = .true.
             IF (rt_model .eq. 'mcrad1d')THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  Thermal IR in-situ emission can '
                WRITE(*,*) '  not be accounted for when calling MCRAD1D'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             ENDIF
             IF (sinsca_corint) THEN
                IF (warning) THEN
                   WRITE(*,*) ''
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : WARNING '
                   WRITE(*,*) '  The TMS correction is only applied to '
                   WRITE(*,*) '  solar photons, not those from '
                   WRITE(*,*) '  in-situ thermal emission'
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                ENDIF
             END IF
             IF (rt_model .eq. 'sinsca')THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  Thermal IR in-situ emission is '
                WRITE(*,*) '  not be accounted for when calling SINSCA'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             ENDIF
             IF (rt_model .eq. 'addoub') THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  Thermal IR in-situ emission can '
                WRITE(*,*) '  not be accounted for when calling ADDOUB'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             ENDIF
             IF (rt_model .eq. 'doad') THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  Thermal IR in-situ emission can '
                WRITE(*,*) '  not be accounted for when calling DOAD'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             ENDIF
             ! IF (mode .eq. 'mono') THEN
             !    WRITE(*,*) ''
             !    WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             !    WRITE(*,*) ' (in sub. read_data)  : ERROR '
             !    WRITE(*,*) '  Thermal IR in-situ emission can '
             !    WRITE(*,*) '  only be accounted for in "kdis" mode'
             !    WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             !    WRITE(*,*) ''
             !    STOP
             ! ENDIF
             ! IF ((mode .eq. 'mono').and.(thermal_only.eqv..FALSE.)) THEN
             !    WRITE(*,*) ''
             !    WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             !    WRITE(*,*) ' (in sub. read_data)  : ERROR '
             !    WRITE(*,*) " In 'mono' mode, you must use the option 'thermal_only' "
             !    WRITE(*,*) " if thermal IR in-situ emission is on "
             !    WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             !    WRITE(*,*) ''
             !    STOP
             ! ENDIF
          ENDIF
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,*)' >>>>>>               Thermal emission : ', TRIM(ADJUSTL(the_char))
          ENDIF

          ! polarization included (4) or scalar (1) or (3) for I, Q, U
          CALL SKIP_COM(finput%unit)
          READ(finput%unit, FMT=*) nmat
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,FMT='(1x, A29, 3x, i1)')' >>>>>>               nmat : ', nmat
          ENDIF
          IF ( nmat .ne. 1 .and.  & 
               nmat .ne. 3 .and.  & 
               nmat .ne. 4) THEN
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' nmat must be 1, 3 or 4'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          IF ((nmat .ne. 1) .and. (rt_model.eq.'disort'))THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  When calling DISORT nmat should be 1'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          ENDIF
          IF ((nmat.gt.1).and. (hg_flag)) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  Henyey-Greenstein approx. is not allowed '
             WRITE(*,*) '  for polarized RT (nmat>1)    '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          ENDIF

          ! --- read the atmosphere model parameters
          the_char = READ_COM (finput%unit)
          atm_type = LOWCASE ( TRIM(ADJUSTL(the_char)), len( TRIM(ADJUSTL(the_char))))
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,*)' >>>>>>               Atmosphere profile : ', TRIM(ADJUSTL(atm_type))
          ENDIF

          ! --- read surface description
          CALL SKIP_COM(finput%unit)
          READ(finput%unit, FMT=*) surface_family, surface_type
          surface_family = LOWCASE ( TRIM (ADJUSTL(surface_family)) , len( TRIM (ADJUSTL(surface_family)) ))
          surface_type = LOWCASE ( TRIM (ADJUSTL(surface_type)) , len( TRIM (ADJUSTL(surface_type)) ))
          IF ((surface_family .ne. 'lambert') .and. (surface_family .ne. 'brdf')) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' The surface family must be "lambert" or "brdf"'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          if ((surface_family.eq.'lambert') .and. (surface_type.eq.'cste')) then
             BACKSPACE(finput%unit)
             READ(finput%unit, FMT=*) surface_family, surface_type, cste_alb
             surface_family = LOWCASE ( TRIM (ADJUSTL(surface_family)) , len( TRIM (ADJUSTL(surface_family)) ))
             surface_type = LOWCASE ( TRIM (ADJUSTL(surface_type)) , len( TRIM (ADJUSTL(surface_type)) ))
          endif

          ! get the surface temperature
          CALL SKIP_COM(finput%unit)
          READ(finput%unit, FMT=*) surface_temp

          if (verbose) then
             WRITE(*,*) ' '
             if (surface_type.eq.'cste') then
                WRITE(*,'(1x, A32, 1x, A7, 3x, A7, 1x,A1, 1x, f10.5)') ' >>>>>>               surface : ', &
                     TRIM (ADJUSTL(surface_family)),TRIM (ADJUSTL(surface_type)), '=', cste_alb
             else
                WRITE(*,*) ' >>>>>>               surface: ', &
                     TRIM (ADJUSTL(surface_family)), ' ',TRIM (ADJUSTL(surface_type))
             endif
             WRITE(*,'(1x, A43, 3x, f10.5)') ' >>>>>>               surface temperature = ', surface_temp
          endif

          ! read incident radiation spec
          CALL SKIP_COM(finput%unit)
          READ (UNIT=finput%unit,FMT=*) solrad
          if (mode.ne.'mono') then
             if (verbose) then
                WRITE(*,*) ' '
                WRITE(*,*) ' >>>>>>               Incident radiation field : ', TRIM (ADJUSTL(solrad))
             ENDIF
          endif

          get_possol = .false.
          !--- read geometry
          CALL SKIP_COM(finput%unit)
          nsza = COUNTCOL(finput%unit)
          ALLOCATE(sza(nsza))
          sza(:) = undef_dp
          READ (UNIT=finput%unit, FMT=*) (sza(i), i=1,nsza)
          CALL SORT(sza)
          do i=1, nsza
             ! if no SZA is provided: we will read lon, lat and time
             if (sza(i).eq.-1) then
                get_possol = .true.
                EXIT
             end if
             if ((sza(i).gt.70.0).and.warning) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : WARNING '
                WRITE(*,*) '  Solar zenithal angle must be <= 70 deg'
                WRITE(*,*) '  for the plan-parallel geometry approximation'
                WRITE(*,*) '  to remain accurate '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
             endif
          enddo
          
          CALL SKIP_COM(finput%unit)
          ntmp = COUNTLINE(finput%unit,finput%name)
          if (ntmp.gt.1) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  You can only provide one set of lon, lat '
             WRITE(*,*) '   --> modify ', TRIM(finput%name)
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          end if
          READ (UNIT=finput%unit, FMT=*) lon, lat
          if (lat.eq.-1) then
             lat = lat_ray_default
          else
             if ((lat.gt.90.0) .or.  (lat.lt.0.0) ) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  Latitude must be 0 <= lat <= 90 '
                WRITE(*,*) '   --> modify ', TRIM(finput%name)
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             end if
          end if
          if ( (lat.ne.-1) .and. ((lon.gt.180.0) .or. (lon.lt.-180.0) )) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  Longitude must be -180 <= lat <= 180 '
             WRITE(*,*) '   --> modify ', TRIM(finput%name)
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          end if

          ! get variables needed to compute solar position

          if (get_possol) then 
             CALL SKIP_COM(finput%unit)
             nsza = COUNTLINE(finput%unit,finput%name)
             DEALLOCATE(sza)
             ALLOCATE(sza(nsza))
             ALLOCATE(h_tu(nsza), day(nsza))
             !print*, nsza
             do i = 1, nsza
                READ (UNIT=finput%unit, FMT=*) day(i), h_tu(i) 
                !print*, day(i), h_tu(i) 
             end do
          else
             CALL SKIP_COM(finput%unit)
             ntmp = COUNTLINE(finput%unit,finput%name)
             do i = 1, ntmp
                READ (UNIT=finput%unit, FMT=*) the_char
             end do
          endif

          CALL SKIP_COM(finput%unit)
          nvza = COUNTCOL(finput%unit)
          ALLOCATE(vza(nvza))
          vza(:) = undef_dp
          READ (UNIT=finput%unit, FMT=*) (vza(i), i=1,nvza)
          CALL SORT(vza)
          ! view zenith angles must be in decreasing order in DISORT
          vza(:) = REVERSE(nvza, vza(:)) 
          do i=1, nvza
             if ((vza(i).gt.70.0).and.warning) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : WARNING '
                WRITE(*,*) '  View zenithal angle must be <= 70 deg'
                WRITE(*,*) '  for the plan-parallel geometry approximation'
                WRITE(*,*) '  to remain accurate '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
             endif
          enddo
          CALL SKIP_COM(finput%unit)
          nvaa = COUNTCOL(finput%unit)
          ALLOCATE(vaa(nvaa))
          vaa(:) = undef_dp
          READ (UNIT=finput%unit, FMT=*) (vaa(i), i=1,nvaa)
          CALL SORT(vaa)

          IF (verbose) THEN
             if (get_possol.eqv..false.) then
                WRITE(*,*) ''
                WRITE(*,FMT='(1x,A35,1x,f10.5)')' >>>>>>                       sza : ', sza(1)
                DO I = 2, nsza
                   WRITE(*,FMT='(37x,f10.5)') sza(i)
                ENDDO
             else
                WRITE(*,*) ''
                WRITE(*,FMT='(1x,A35,2(1x,f10.5,1x),I5,1x,f10.5)') &
                     ' >>>>>>    lon, lat, day, time    : ', lon,lat,day(1), h_tu(1) 
                DO I = 2, nsza
                WRITE(*,FMT='(A37,2(1x,f10.5,1x),I5,1x,f10.5)')&
                     ' >>>>>>    lon, lat, day, time    : ', lon,lat,day(i), h_tu(i) 
                ENDDO                
             end if
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A35,1x,f10.5)')' >>>>>>                       vza : ', vza(1)
             DO I = 2, nvza
                WRITE(*,FMT='(37x,f10.5)') vza(i)
             ENDDO
             WRITE(*,*) ''
             WRITE(*,FMT='(1x,A35,1x,f10.5)')' >>>>>>                       vaa : ', vaa(1)
             DO I = 2, nvaa
                WRITE(*,FMT='(37x,f10.5)') vaa(i)
             ENDDO
          ENDIF

       ENDIF !(do_rt)
    ENDIF  ! (opt_only .eqv. .FALSE.)
    CLOSE(UNIT=finput%unit)
    ! END READING artdeco_in.dat
    if (verbose) then
       WRITE(*,*) ' '
       WRITE(*,*) '  (sub. read_data) finish read artdeco_in.dat'
       WRITE(*,*) ' '
    ENDIF

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! Create the output directory
    if ( (opt_only.eqv..false.) .and. f2py_mode.eqv..false.) then
       dir_root_out = TRIMCAT (artdeco_path, dir_out)
       dir_root_out = TRIMCAT (dir_root_out, out_root)
       if (betal_only.eqv..false.) then
          dir_root_out = TRIMCAT (dir_root_out, '_')
          dir_root_out = TRIMCAT (dir_root_out, rt_model)
       end if
       dir_root_out = TRIMCAT (dir_root_out, '/')
       if (verbose) then
          WRITE(*,*) ' '
          WRITE(*,*) '  (sub. read_data) create the output directory ', TRIM(ADJUSTL(dir_root_out))
          WRITE(*,*) ' '
       ENDIF
       INQUIRE(FILE=dir_root_out, EXIST=dir_exist)
       if (dir_exist) then
          command = 'rm -rf'
          command = TRIM(ADJUSTL(command))//" "//dir_root_out       
          CALL SYSTEM(command, status=st)
       endif
       command = 'mkdir'
       command = TRIM(ADJUSTL(command))//" "//dir_root_out
       CALL SYSTEM(command, status=st)
    endif

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! read the filter(s) description
    if (filter) then
       if (verbose) then
          WRITE(*,*) ' '
          WRITE(*,*) '  (sub. read_data) Read filter file(s)'
          WRITE(*,*) ' '
       ENDIF
       ! first look for nmaxlamb_filter
       nmaxlamb_filter = 0
       DO i = 1, nfilter
          filename_tmp = TRIMCAT(artdeco_path, dir_filter)
          filename_tmp = TRIMCAT(filename_tmp, ffilter%name)
          filename_tmp = TRIMCAT(filename_tmp, filter_name(i))
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          INQUIRE(FILE=filename_tmp, EXIST=test_file)
          IF ( test_file .eqv. .FALSE. ) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' You need to provide a filter file ', TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) ' for filter ', TRIM(ADJUSTL(filter_name(i)))
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          OPEN (UNIT = ffilter%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(ffilter%unit) 
          READ (UNIT=ffilter%unit,FMT=*) ntmp
          if (ntmp.gt.nmaxlamb_filter) nmaxlamb_filter = ntmp
          CLOSE(UNIT = ffilter%unit)
       END DO
       ALLOCATE(rt_cputime_filter(nfilter), &
            nlamb_filter(nfilter), &
            lamb_filter(nfilter, nmaxlamb_filter), &
            trans_filter(nfilter, nmaxlamb_filter), &
            bound_lamb_filter(2,nfilter))
       rt_cputime_filter(:) = undef_dp
       nlamb_filter(:)   = undef_i
       lamb_filter(:, :) = undef_dp
       trans_filter(:,:) = undef_dp
       bound_lamb_filter(:,:) = undef_dp
       DO i = 1, nfilter
          filename_tmp = TRIMCAT(artdeco_path, dir_filter)
          filename_tmp = TRIMCAT(filename_tmp, ffilter%name)
          filename_tmp = TRIMCAT(filename_tmp, filter_name(i))
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          IF (verbose) THEN
             WRITE(*,*)  '                   --> read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          OPEN (UNIT = ffilter%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(ffilter%unit) 
          READ (UNIT=ffilter%unit,FMT=*) nlamb_filter(i)
          CALL SKIP_COM(ffilter%unit) 
          DO k = 1, nlamb_filter(i)
             READ (UNIT=ffilter%unit,FMT=*) lamb_filter(i, k), trans_filter(i, k)
             if (k .gt.1) then
                if (lamb_filter(i, k) .le. lamb_filter(i, k-1)) then
                   WRITE(*,*) ''
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) ' In filter file ', TRIM(ADJUSTL(filename_tmp))
                   WRITE(*,*) ' Wavelengths must be sorted in increasing order '
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                   STOP                                           
                endif
             endif
          END DO
          CLOSE(UNIT = ffilter%unit)
          bound_lamb_filter(1,i) = lamb_filter(i, 1)
          bound_lamb_filter(2,i) = lamb_filter(i, nlamb_filter(i))
       END DO
    endif

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! NOTE: if k-distribution will be used, it will constrain the 
    ! working wavelength grid, so we need to read k-dis files first
    ! READ K-DISTRIBUTION FILES
    IF (mode.eq.'kdis') then
       if (verbose) then
          WRITE(*,*) ' '
          WRITE(*,*) '  (sub. read_data) Start to read k-distribution input files'
          WRITE(*,*) ' '
       ENDIF

       filename_tmp = TRIMCAT(artdeco_path, dir_kdis)
       filename_tmp = TRIMCAT(filename_tmp, fkdis%name)
       filename_tmp = TRIMCAT(filename_tmp, kdis_model)
       filename_tmp = TRIMCAT(filename_tmp, '_def.dat')
       filename_kdis = filename_tmp
          
       IF (verbose) THEN
          WRITE(*,*)  '    --> read  ', TRIM(ADJUSTL(filename_tmp))
       ENDIF
       INQUIRE(FILE=filename_tmp, EXIST=test_file)
       IF ( test_file .eqv. .FALSE. ) THEN
          WRITE(*,*) ''
          WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          WRITE(*,*) ' (in sub. read_data)  : ERROR '
          WRITE(*,*) ' You need to provide k-distribution definition a file '
          WRITE(*,*) ' Required file is :'
          WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
          WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          WRITE(*,*) ''
          STOP 
       ENDIF
       OPEN (UNIT = fkdis%unit, FILE = filename_tmp, STATUS = 'old')
       CALL SKIP_COM(fkdis%unit) 
       READ (UNIT=fkdis%unit,FMT=*) kdis_nmaxai
       !print*, kdis_nmaxai
       CALL SKIP_COM(fkdis%unit) 
       READ (UNIT=fkdis%unit,FMT=*) ntmp
       ALLOCATE(char_arr_tmp(ntmp), logic_arr_tmp(ntmp), dbl_arr_tmp(ntmp))
       char_arr_tmp(:)  = undef_c
       logic_arr_tmp(:) = undef_i
       dbl_arr_tmp(:)   = undef_dp
       kdis_nsp   = 0
       kdis_nsp_c = 0
       CALL SKIP_COM(fkdis%unit) 
       DO i = 1, ntmp
          READ (UNIT=fkdis%unit,FMT=*) char_arr_tmp(i), logic_arr_tmp(i), dbl_arr_tmp(i)
          if (logic_arr_tmp(i).eq.1) then
             kdis_nsp_c = kdis_nsp_c + 1
          else
             kdis_nsp = kdis_nsp + 1 
          endif
       ENDDO
       if (kdis_nsp.gt.0) then
          ALLOCATE(kdis_species(kdis_nsp), kdis_fcont(kdis_nsp))
          kdis_species(:) = undef_c
          kdis_fcont(:)   = undef_dp
       endif
       if (kdis_nsp_c.gt.0) then
          ALLOCATE(kdis_species_c(kdis_nsp_c), kdis_fcont_c(kdis_nsp_c))
          kdis_species_c(:) = undef_c
          kdis_fcont_c(:)   = undef_dp
       endif
       k=0
       j=0
       DO i = 1, ntmp
          if (logic_arr_tmp(i).eq.1) then
             k=k+1
             kdis_species_c(k) = TRIM(ADJUSTL(char_arr_tmp(i)))
             kdis_fcont_c(k)   = dbl_arr_tmp(i)
          else
             j=j+1
             kdis_species(j) = TRIM(ADJUSTL(char_arr_tmp(i)))
             kdis_fcont(j)   = dbl_arr_tmp(i)
          endif
       ENDDO
       DEALLOCATE(char_arr_tmp, logic_arr_tmp, dbl_arr_tmp)
       CALL SKIP_COM(fkdis%unit) 
       READ (UNIT=fkdis%unit,FMT=*) kdis_nwvl
       !print*, kdis_nwvl
       ALLOCATE(kdis_wvlband(3,kdis_nwvl))
       kdis_wvlband(:,:) = undef_dp
       CALL SKIP_COM(fkdis%unit) 
       DO i = 1, kdis_nwvl
          READ (UNIT=fkdis%unit,FMT=*) j, kdis_wvlband(1,i), kdis_wvlband(2,i), kdis_wvlband(3,i)
          !print*, j, kdis_wvlband(1,i), kdis_wvlband(2,i), kdis_wvlband(3,i)
          if (i.gt.1) then
             if (kdis_wvlband(1,i).le.kdis_wvlband(1,i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '  Wavelength must be sorted in increasing order '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             endif
             if ((kdis_wvlband(2,i).gt.kdis_wvlband(3,i-1)).and.warning) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : WARNING '
                WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,FMT='(1x,A52,1x,i3,1x,A3,1x,i3)') &
                     '  There seems to be a wavelength gap between bands ', i-1, 'and', i
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''                
             endif
             if ((kdis_wvlband(2,i).lt.kdis_wvlband(3,i-1)).and.warning) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : WARNING '
                WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,FMT='(1x,A58,1x,i3,1x,A3,1x,i3)') &
                     '  There seems to be a wavelength covering between bands ', i-1, 'and', i
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''                
             endif
          endif ! (i.gt.1)
       ENDDO
       ! read PRESSURE
       CALL SKIP_COM(fkdis%unit) 
       READ (UNIT=fkdis%unit,FMT=*) kdis_np
       ALLOCATE(kdis_p(kdis_np))
       kdis_p(:) = undef_dp
       CALL SKIP_COM(fkdis%unit) 
       DO i = 1, kdis_np
          READ (UNIT=fkdis%unit,FMT=*) kdis_p(i)
          !print*, i, kdis_np, kdis_p(i)
          if (i.gt.1) then
             if (kdis_p(i).le.kdis_p(i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '  Pressures must be sorted in increasing order '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             endif
          ENDIF
       ENDDO
       ! read Temperature
       CALL SKIP_COM(fkdis%unit) 
       READ (UNIT=fkdis%unit,FMT=*) kdis_nt
       ALLOCATE(kdis_t(kdis_nt))
       kdis_t(:) = undef_dp
       CALL SKIP_COM(fkdis%unit) 
       DO i = 1, kdis_nt
          READ (UNIT=fkdis%unit,FMT=*) kdis_t(i)
          !print*, i, kdis_nt, kdis_t(i)
          if (i.gt.1) then
             if (kdis_t(i).le.kdis_t(i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '  Temperatures must be sorted in increasing order '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             endif
          endif
       ENDDO
       ! read concentration
       if (kdis_nsp_c.gt.0) then
          CALL SKIP_COM(fkdis%unit) 
          READ (UNIT=fkdis%unit,FMT=*) kdis_nc
          if ((kdis_nc.eq.0).and.(kdis_nsp_c.gt.0)) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) '  At least one gas was declared with concentration dependence '
             WRITE(*,*) '  There should be a concentration grid provided. '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
          endif
          if (kdis_nc.gt.0) THEN
             ALLOCATE(kdis_c(kdis_nc))
             kdis_c(:) = undef_dp
             CALL SKIP_COM(fkdis%unit) 
             DO i = 1, kdis_nc
                READ (UNIT=fkdis%unit,FMT=*) kdis_c(i)
                !print*, i, kdis_nc, kdis_c(i)
                if (i.gt.1) then
                   if (kdis_c(i).le.kdis_c(i-1)) then
                      WRITE(*,*) ''
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (in sub. read_data)  : ERROR '
                      WRITE(*,*) '  In ', TRIM(ADJUSTL(filename_tmp))
                      WRITE(*,*) '  Concentrations must be sorted in increasing order '
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ''
                      STOP 
                   endif
                endif
             ENDDO
          ENDIF
       endif
       CLOSE(unit=fkdis%unit) 
       if(kdis_nsp.gt.0)then 
          ALLOCATE(kdis_nai(kdis_nsp,kdis_nwvl), &
               kdis_ki(kdis_nsp,kdis_nwvl,kdis_nmaxai,kdis_np,kdis_nt),&                  
               kdis_ai(kdis_nsp,kdis_nwvl,kdis_nmaxai), &
               kdis_xsect(kdis_nsp,kdis_nwvl))
          kdis_nai(:,:)         = undef_i
          kdis_ki(:,:,:,:,:)    = undef_dp
          kdis_ai(:,:,:)        = undef_dp
          kdis_xsect(:,:)       = undef_dp
       endif
       if(kdis_nsp_c.gt.0)then 
          ALLOCATE(kdis_nai_c(kdis_nsp_c,kdis_nwvl), &
               kdis_ki_c(kdis_nsp_c,kdis_nwvl,kdis_nmaxai,kdis_np,kdis_nt,kdis_nc),&                  
               kdis_ai_c(kdis_nsp_c,kdis_nwvl,kdis_nmaxai), &
               kdis_xsect_c(kdis_nsp_c,kdis_nwvl))
          kdis_nai_c(:,:)         = undef_i
          kdis_ki_c(:,:,:,:,:,:)  = undef_dp
          kdis_ai_c(:,:,:)        = undef_dp
          kdis_xsect_c(:,:)       = undef_dp
       endif
       ! read ai,ki
       gas_species : do isp = 1, kdis_nsp
          filename_tmp = TRIMCAT(artdeco_path, dir_kdis)
          filename_tmp = TRIMCAT(filename_tmp, fkdis%name)
          filename_tmp = TRIMCAT(filename_tmp, kdis_model)
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, kdis_species(isp))
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          
          IF (verbose) THEN
             WRITE(*,*)  '    --> read  ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          INQUIRE(FILE=filename_tmp, EXIST=test_file)
          IF ( test_file .eqv. .FALSE. ) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' You need to provide k-distribution  definition a file '
             WRITE(*,*) ' Required file is :'
             WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          OPEN (UNIT = fkdis%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fkdis%unit) 
          DO iwvl = 1, kdis_nwvl
             READ (UNIT=fkdis%unit,FMT=*) j, kdis_nai(isp,iwvl), kdis_xsect(isp,iwvl)
          ENDDO
          DO iwvl = 1, kdis_nwvl
             if (kdis_nai(isp,iwvl).gt.1) then
                CALL SKIP_COM(fkdis%unit) 
                ! read ai
                READ (UNIT=fkdis%unit,FMT=*) (kdis_ai(isp,iwvl,i), i=1, kdis_nai(isp,iwvl))
                ! read ki
                DO it = 1, kdis_nt
                   DO ip = 1, kdis_np
                      READ (UNIT=fkdis%unit,FMT=*) (kdis_ki(isp,iwvl,i,ip,it), i=1, kdis_nai(isp,iwvl))
                   END DO
                END DO
             endif
          ENDDO
          CLOSE(UNIT = fkdis%unit)
       enddo gas_species
       gas_species_c : do isp = 1, kdis_nsp_c
          filename_tmp = TRIMCAT(artdeco_path, dir_kdis)
          filename_tmp = TRIMCAT(filename_tmp, fkdis%name)
          filename_tmp = TRIMCAT(filename_tmp, kdis_model)
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, kdis_species_c(isp))
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          IF (verbose) THEN
             WRITE(*,*)  '    --> read  ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          INQUIRE(FILE=filename_tmp, EXIST=test_file)
          IF ( test_file .eqv. .FALSE. ) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' You need to provide k-distribution  definition a file '
             WRITE(*,*) ' Required file is :'
             WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          OPEN (UNIT = fkdis%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fkdis%unit) 
          DO iwvl = 1, kdis_nwvl
             READ (UNIT=fkdis%unit,FMT=*) j, kdis_nai_c(isp,iwvl), kdis_xsect_c(isp,iwvl)
          ENDDO
          DO iwvl = 1, kdis_nwvl
             if (kdis_nai_c(isp,iwvl).gt.1) then
                CALL SKIP_COM(fkdis%unit) 
                ! read ai
                READ (UNIT=fkdis%unit,FMT=*) (kdis_ai_c(isp,iwvl,i), i=1, kdis_nai_c(isp,iwvl))
                ! read ki
                DO ic = 1, kdis_nc
                   DO it = 1, kdis_nt
                      DO ip = 1, kdis_np
                         READ (UNIT=fkdis%unit,FMT=*) (kdis_ki_c(isp,iwvl,i,ip,it,ic), i=1, kdis_nai_c(isp,iwvl))
                      END DO
                   END DO
                END DO
             endif
          ENDDO
          CLOSE(UNIT = fkdis%unit)
       enddo gas_species_c
       if (verbose) then
          WRITE(*,*) ' '
          WRITE(*,*) '  (sub. read_data) finish read k-distribution input files'
          WRITE(*,*) ' '
       ENDIF

       ! set wlambda(:) and lambda_ikdis(:)
       ! used for all scattering purpose
       ALLOCATE(flag_kdiswvl(kdis_nwvl))
       flag_kdiswvl(:) = .FALSE.
       nlambda = 0

       if (filter) then

          ! test whether the filter range is covered by the k-dis wavelengths
          DO k = 1, nfilter
             IF ( (bound_lamb_filter(1,k).lt.kdis_wvlband(2,1)) .or. &
                  (bound_lamb_filter(2,k).gt.kdis_wvlband(3,kdis_nwvl)) ) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' Wavelengths range of filter ', ADJUSTL(TRIM(filter_name(k)))
                WRITE(*,*) ' is out of the k-distribution covered range'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
          END DO

          ! get working wavelengths   
          DO i= 1, kdis_nwvl
             DO k = 1, nfilter
                if ( ((kdis_wvlband(2,i).ge.bound_lamb_filter(1,k)) .and. (kdis_wvlband(2,i).le.bound_lamb_filter(2,k))) .or. &
                     ((kdis_wvlband(3,i).ge.bound_lamb_filter(1,k)) .and. (kdis_wvlband(3,i).le.bound_lamb_filter(2,k))) ) then
                   nlambda = nlambda + 1
                   flag_kdiswvl(i) = .true.
                   exit
                endif
             END DO
          END DO

       else

          if (lambda_bound(1).eq.-1.0_dp) lambda_bound(1) = kdis_wvlband(2,1)
          if (lambda_bound(2).eq.-1.0_dp) lambda_bound(2) = kdis_wvlband(3,kdis_nwvl)
          IF (lambda_bound(1).lt.kdis_wvlband(2,1)) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' The minimum boundary wavelength is out of the k-distribution '
             WRITE(*,*) ' wavelength range.'
             WRITE(*,FMT='(1x,A20,1x,f10.5)') ' Minimum boundary = ',lambda_bound(1)
             WRITE(*,FMT='(1x,A20,1x,f10.5)') ' kdis minimum     = ',kdis_wvlband(2,1)
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''                
             STOP   
          ENDIF
          IF (lambda_bound(2).gt.kdis_wvlband(3,kdis_nwvl)) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' The maximum boundary wavelength is out of the k-distribution '
             WRITE(*,*) ' wavelength range.'
             WRITE(*,FMT='(1x,A20,1x,f10.5)') ' maximum boundary = ',lambda_bound(2)
             WRITE(*,FMT='(1x,A20,1x,f10.5)') ' kdis maximum     = ',kdis_wvlband(3,kdis_nwvl)
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''                
             STOP   
          ENDIF

          DO i = 1, kdis_nwvl
             if ((kdis_wvlband(3,i).gt.lambda_bound(1)) .and. (kdis_wvlband(2,i).lt.lambda_bound(2))) then
                nlambda = nlambda + 1
                flag_kdiswvl(i) = .true.
             endif
          END DO

       endif ! filter

       ALLOCATE(wlambda(nlambda), lambda_ikdis(nlambda))
       k = 0
       DO i= 1, kdis_nwvl
          if (flag_kdiswvl(i)) then
             k = k+1
             wlambda(k)       = kdis_wvlband(1,i)
             lambda_ikdis(k) = i             
          endif
       END DO

       DEALLOCATE(flag_kdiswvl)

    ENDIF !(kdis mode)

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! check if reflamb is equal to any working wavelength
    ireflamb = 0
    if (do_rt) then
       DO j = 1, nlambda
          IF ( (wlambda(j).ge.(reflamb-EPS_LAMB)) .and. &
               (wlambda(j).le.(reflamb+EPS_LAMB)) ) then
             ireflamb = j
             exit
          endif
       ENDDO
    endif

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! READ RTE SOLVER SPEC FILES
    if(do_rt) then
       IF (rt_model .eq. 'mcrad1d')THEN
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fmcrad1d_spec%name)
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          OPEN (UNIT = fmcrad1d_spec%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_nphot
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_nmax_interact
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_vrm
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_epsilon_ddis
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_vrm_n_firstcp, mc_vrm_n_sccp, mc_vrm_n_lecp
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_wspl_crit
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_nsplmax
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_wrr_crit
          CALL SKIP_COM(fmcrad1d_spec%unit)
          READ(fmcrad1d_spec%unit, FMT=*) mc_wrr_min
          IF ((mc_vrm .ne. 0) .and. (mc_vrm .ne. 1)) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' mc_vrm must be 0 or 1 '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif
          IF ((mc_epsilon_ddis .gt. 1) .or. (mc_epsilon_ddis .lt. 0)) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' mc_epsilon_ddis must be >0 and <1'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif
          CLOSE(UNIT = fmcrad1d_spec%unit)
       ENDIF
       IF(rt_model .eq. 'doad') THEN
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fdoad_spec%name)
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          OPEN (UNIT = fdoad_spec%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fdoad_spec%unit)
          READ(fdoad_spec%unit, FMT=*) doad_eps
          CALL SKIP_COM(fdoad_spec%unit)
          READ(fdoad_spec%unit, FMT=*) doad_nmug
          CALL SKIP_COM(fdoad_spec%unit)
          READ(fdoad_spec%unit, FMT=*) doad_nfoumx
          CALL SKIP_COM(fdoad_spec%unit)
          READ(fdoad_spec%unit, FMT=*) doad_ncoef_min_brdf
          CLOSE(UNIT =  fdoad_spec%unit)
       ENDIF
       IF(rt_model .eq. 'addoub') THEN
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fad_spec%name)          
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          OPEN (UNIT = fad_spec%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fad_spec%unit)
          READ(fad_spec%unit, FMT=*) ad_epsilon
          CALL SKIP_COM(fad_spec%unit)
          READ(fad_spec%unit, FMT=*) ad_nmug
          CALL SKIP_COM(fad_spec%unit)
          READ(fad_spec%unit, FMT=*) ad_mstop
          CALL SKIP_COM(fad_spec%unit)
          READ(fad_spec%unit, FMT=*) ad_nfoumax
          CLOSE(UNIT =  fad_spec%unit)
       ENDIF
       IF(rt_model .eq. 'disort') THEN
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fod_spec%name)
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          OPEN (UNIT = fod_spec%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fod_spec%unit)
          READ(fod_spec%unit, FMT=*) od_nstr
          CALL SKIP_COM(fod_spec%unit)
          READ(fod_spec%unit, FMT=*) od_accur
          CLOSE(fod_spec%unit)
          
          ! check nstr 
          flag_chnstrtmp = .TRUE.
          flag_chnstr    = .FALSE.
          DO WHILE (flag_chnstrtmp)
             do i = 1, nsza
                CALL CHECK_NSTR(od_nstr, cos(sza(i) / 180.0 * xpi), 2000, flag_chnstrtmp)  
                if (flag_chnstrtmp) then
                   flag_chnstr = .TRUE.
                   EXIT
                end if
             end do
          END DO
          IF (flag_chnstr) THEN
             if (verbose.and.warning) then
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : WARNING '
                WRITE(*,*) '  Number of streams for DISORT was changed automatically'
                WRITE(*,*) '  NSTR = ',od_nstr 
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
             endif
          ENDIF
          
       ENDIF
    endif !do_rt

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if (opt_only.eqv..false.) then
       IF (nptcle .gt. 0) then

          ! set nmaxbetal which is a dimensioning variable

          IF (nbetal_in .eq. -1) then
             IF (do_rt.eqv..FALSE.) then
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' You need to set nbetal in artdeco_in.dat'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             ! No default value for mcrad1d & sinsca
             IF ( ( (rt_model .eq. 'mcrad1d') .or. (rt_model .eq. 'sinsca') ) &
                  .and.(trunc_method.ne.'none')) THEN
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' You need to set nbetal in artdeco_in.dat'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             ! If no nbetal specification was given in artdeco_in.dat 
             ! we use the number of computationnal gauss angles 
             ! specified for the RT models
             IF(rt_model .eq. 'disort')  nbetal_in = od_nstr       ! 1 time nstr
             IF(rt_model .eq. 'doad')    nbetal_in = 2 * doad_nmug ! optimum in adding with 2 time nstr
             IF(rt_model .eq. 'addoub')  nbetal_in = 2 * ad_nmug   ! optimum in adding with 2 time nstr

          ELSE

             IF(nbetal_in .lt. 2) then
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' Required numbers of Bl for Legendre polynomial expension'
                WRITE(*,*) ' should be greater or equal to 2'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
          ENDIF

          ! set the dimensioning variable
          if (nbetal_in.ge.0) then
             nmaxbetal = nbetal_in
          else
             nmaxbetal = 0
          endif

       ELSE

          nbetal_in = 0
          IF(rt_model .eq. 'sinsca')  nmaxbetal = 0
          IF(rt_model .eq. 'mcrad1d') nmaxbetal = 0
          IF(rt_model .eq. 'disort')  nmaxbetal = od_nstr       ! 1 time nstr
          IF(rt_model .eq. 'doad')    nmaxbetal = 2 * doad_nmug ! optimum in adding with 2 time nstr
          IF(rt_model .eq. 'addoub')  nmaxbetal = 2 * ad_nmug   ! optimum in adding with 2 time nstr

       ENDIF ! nptcle gt 0
    endif

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! START TO RETREIVE PARTICLE OPTICAL PROPERTIES
    ! AND PARTICLES PROPERTIES
    nmaterial = 0
    IF (nptcle.gt.0) THEN

       flag_store_opt = .FALSE.

       ALLOCATE(ptcle_nelem(nptcle),&
            ptcle_flag_store_opt(nptcle),&
            ptcle_opt_exist_reflamb(nptcle),&
            ptcle_opt_exist(nptcle,nlambda), &
            ptcle_ilambda(nptcle, nlambda), &
            ptcle_ireflamb(nptcle))
       ptcle_nelem(:)             = undef_i
       ptcle_flag_store_opt(:)    = .FALSE.
       ptcle_opt_exist_reflamb(:) = .FALSE. 
       ptcle_opt_exist(:,:)       = .FALSE.
       ptcle_ilambda(:,:)         = undef_i
       ptcle_ireflamb(:)          = undef_i

       ! get nlamb_max_store 
       !     nang_max_store
       !     ptcle_opt_exist_reflamb(:)
       !     ptcle_opt_exist(:,:)      
       nlamb_max_store = 0  
       nang_max_store  = 0

       IF (verbose) WRITE(*,*) ''
       DO i = 1, nptcle
          nlamb = 0 ! need to init that value
          filename_tmp = TRIMCAT (artdeco_path, dir_opt)
          filename_tmp = TRIMCAT(filename_tmp, fopt%name)
          filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
          filename_tmp = TRIMCAT(filename_tmp, '.dat')          
          INQUIRE(FILE=filename_tmp, EXIST=test_file) 
          if (test_file.eqv..FALSE.) then
             ! no optical properties at all for that particle
             IF (verbose) THEN
                WRITE(*,*) '  (sub. read_data) no file ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             IF (ptcle_opt_interp(i)) THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' You need to provide an optical properties file   '
                WRITE(*,*) ' for particle ', TRIM(ADJUSTL(ptcle_type(i)))
                WRITE(*,*) ' if you want interpolated optical properties (interp. = yes)'
                WRITE(*,*) ' Required file is :'
                WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
          else
             ptcle_flag_store_opt(i)=.TRUE.
             flag_store_opt = .TRUE.
             IF (verbose) THEN
                WRITE(*,*) '  (sub. read_data) find file ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             ! check wavelengths of optical properties stored in the file
             OPEN (UNIT = fopt%unit, FILE = filename_tmp, STATUS = 'OLD')
             the_char = READ_COM (fopt%unit) ! skip model def.
             the_char = READ_COM (fopt%unit) ! skip material name
             CALL SKIP_COM(fopt%unit)
             READ(fopt%unit, FMT=*) nelem ! read number phase matrix elements stored
             CALL SKIP_COM(fopt%unit)
             READ(fopt%unit, FMT=*) nlamb ! read number of wavelengths in the file
             ALLOCATE(lamb(nlamb), nteta(nlamb))
             CALL SKIP_COM(fopt%unit)
             DO l = 1, nlamb
                READ(fopt%unit, FMT=*) lamb(l), nteta(l)
                if (l.gt.1) then
                   if (lamb(l).le.lamb(l-1)) then
                      WRITE(*,*) ''
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (in sub. read_data)  : ERROR '
                      WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                      WRITE(*,*) ' Wavelengths must be sorted in increasing order'
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ''
                      STOP 
                   endif
                endif
             ENDDO
             CLOSE(fopt%unit)
             ! check for existing optical properties for the reference wavelength
             do l = 1, nlamb
                IF ( (reflamb.ge.(lamb(l)-EPS_LAMB)) .and. &
                     (reflamb.le.(lamb(l)+EPS_LAMB)) ) THEN
                   ptcle_ireflamb(i)           = l         
                   ptcle_opt_exist_reflamb(i)  = .TRUE.
                   EXIT 
                END IF
             end do
             ! see if the stored wavelengths match the required one
             DO j = 1, nlambda
                DO l = 1, nlamb
                   IF ( (wlambda(j).ge.(lamb(l)-EPS_LAMB)) .and. &
                        (wlambda(j).le.(lamb(l)+EPS_LAMB)) ) THEN
                      ptcle_ilambda(i,j)   = l  
                      ptcle_opt_exist(i,j) = .TRUE.  
                      EXIT
                   ENDIF
                END DO
             END DO
             if (nlamb_max_store.lt.nlamb) &
                  nlamb_max_store = nlamb
             if (nang_max_store.lt.MAXVAL(nteta(:))) &
                  nang_max_store = MAXVAL(nteta(:))
             DEALLOCATE(lamb, nteta)          
          endif
       END DO ! loop on ptcles

       IF (flag_store_opt) THEN

          ALLOCATE(ptcle_opt_store(nptcle, ptcle_nopt, nlamb_max_store), &
               ptcle_nlamb_phasemat_store(nptcle), & 
               ptcle_nang_phasemat_store(nptcle,nlamb_max_store), &
               ptcle_lamb_phasemat_store(nptcle, nlamb_max_store), &
               ptcle_u_phasemat_store(nptcle, nang_max_store, nlamb_max_store), &
               ptcle_phasemat_store(nptcle, 6, nang_max_store, nlamb_max_store))
          ptcle_opt_store(:,:,:)          = undef_dp
          ptcle_nlamb_phasemat_store(:)   = undef_i 
          ptcle_nang_phasemat_store(:,:)  = undef_i
          ptcle_lamb_phasemat_store(:, :) = undef_dp
          ptcle_u_phasemat_store(:,:,:)   = undef_dp
          ptcle_phasemat_store(:,:,:,:)   = undef_dp

          ! read the optical properties file entirely
          IF (verbose) WRITE(*,*) ''
          DO i = 1, nptcle

             ptcle_nlamb_phasemat_store(i)  = 0
             ptcle_nang_phasemat_store(i,:) = 0

             nlamb = 0 ! need to init that value
             filename_tmp = TRIMCAT (artdeco_path, dir_opt)
             filename_tmp = TRIMCAT(filename_tmp, fopt%name)
             filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
             filename_tmp = TRIMCAT(filename_tmp, '.dat')          
             INQUIRE(FILE=filename_tmp, EXIST=test_file) 
             if (test_file) then
                IF (verbose) THEN
                   WRITE (*,*)'  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp)) 
                ENDIF
                ! check optical properties stored in the file
                OPEN (UNIT = fopt%unit, FILE = filename_tmp, STATUS = 'OLD')
                the_char = READ_COM (fopt%unit) ! skip model def.
                the_char = READ_COM (fopt%unit) ! skip material name
                CALL SKIP_COM(fopt%unit)
                READ(fopt%unit, FMT=*) ptcle_nelem(i) ! read number phase matrix elements stored
                CALL SKIP_COM(fopt%unit)
                READ(fopt%unit, FMT=*) ptcle_nlamb_phasemat_store(i) ! read number of wavelengths in the file
                CALL SKIP_COM(fopt%unit)
                DO l = 1, ptcle_nlamb_phasemat_store(i)
                   READ(fopt%unit, FMT=*) ptcle_lamb_phasemat_store(i, l), &
                        ptcle_nang_phasemat_store(i,l), &
                        ptcle_opt_store(i, 1, l), &
                        ptcle_opt_store(i, 2, l), &
                        ptcle_opt_store(i, 3, l)
                ENDDO

                if (ptcle_hg(i) .eqv. .FALSE.) then
                   ! read phase matrices
                   CALL SKIP_COM(fopt%unit)
                   DO l = 1, ptcle_nlamb_phasemat_store(i)
                      CALL SKIP_COM(fopt%unit)
                      CALL SKIP_COM(fopt%unit)
                      if (ptcle_nelem(i).eq.4) then
                         DO n = 1, ptcle_nang_phasemat_store(i,l)
                            READ(fopt%unit, FMT=*)  ptcle_u_phasemat_store(i, n, l), &
                                 ptcle_phasemat_store(i, 1, n, l),ptcle_phasemat_store(i, 4, n, l), &
                                 ptcle_phasemat_store(i, 5, n, l),ptcle_phasemat_store(i, 6, n, l)
                            ptcle_phasemat_store(i, 2, n, l) = ptcle_phasemat_store(i, 1, n, l)
                            ptcle_phasemat_store(i, 3, n, l) = ptcle_phasemat_store(i, 4, n, l)
                         ENDDO
                      else
                         DO n = 1, ptcle_nang_phasemat_store(i,l)
                            READ(fopt%unit, FMT=*)  ptcle_u_phasemat_store(i, n, l), &
                                 ptcle_phasemat_store(i, 1, n, l),ptcle_phasemat_store(i, 2, n, l),&
                                 ptcle_phasemat_store(i, 3, n, l),ptcle_phasemat_store(i, 4, n, l),&
                                 ptcle_phasemat_store(i, 5, n, l),ptcle_phasemat_store(i, 6, n, l)
                         ENDDO
                      end if
                   END DO
                end if

                CLOSE(UNIT=fopt%unit)
             endif ! opt file exist
          END DO ! loop on nptcle
       END IF ! flag_store_opt

       ALLOCATE(ptcle_material(nptcle), &
            ptcle_opt_model(nptcle),    &
            ptcle_param_mie3(nptcle, mie3_npar), &
            ptcle_nphot_hm(nptcle), &  
            ptcle_param_hm(nptcle,hm_nmaxpar),&
            ptcle_param_baum(nptcle))
       ptcle_material(:)      = undef_c
       ptcle_opt_model(:)     = undef_c
       ptcle_param_mie3(:, :) = undef_dp
       ptcle_nphot_hm(:)      = undef_i
       ptcle_param_hm(:,:)    = undef_dp
       ptcle_param_baum(:)    = undef_dp

       ! eventually get particle properties to compute
       ! new optical properties
       IF (verbose) WRITE(*,*) ''
       DO i = 1, nptcle

          ! Do we need to read prop_[...].dat ? 
          read_prop = .false.
          IF (ptcle_opt_interp(i).eqv..false.) THEN
             DO j = 1, nlambda
                IF (ptcle_opt_exist(i,j) .eqv. .false.) then
                   read_prop = .true.
                   EXIT
                endif
             ENDDO
             if (do_rt) then
                if (ptcle_opt_exist_reflamb(i).eqv..false.) read_prop = .true.
             endif
          END IF

          ! READ prop_[...].dat files if necessary
          IF (read_prop .eqv. .true.) THEN
             filename_tmp = TRIMCAT(artdeco_path, dir_prop)
             filename_tmp = TRIMCAT(filename_tmp, fprop%name)
             filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             IF (verbose) THEN
                WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             INQUIRE(FILE=filename_tmp, EXIST=test_file)
             IF ( test_file .eqv. .FALSE. ) THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' You need to provide a properties definition file '
                WRITE(*,*) ' for particle ', TRIM(ADJUSTL(ptcle_type(i)))
                WRITE(*,*) ' Required file is :'
                WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             OPEN (UNIT = fprop%unit, FILE = filename_tmp, STATUS = 'old')
             CALL SKIP_COM(fprop%unit)
             READ(fprop%unit, FMT=*) ptcle_material(i)
             ptcle_material(i) = LOWCASE(TRIM(ADJUSTL(ptcle_material(i))), len(TRIM(ADJUSTL(ptcle_material(i)))))
             CALL SKIP_COM(fprop%unit)
             READ(fprop%unit, FMT=*) ptcle_opt_model(i)
             ptcle_opt_model(i) = LOWCASE(TRIM(ADJUSTL(ptcle_opt_model(i))), len(TRIM(ADJUSTL(ptcle_opt_model(i)))))
             IF ( ptcle_opt_model(i) .ne. 'mie3' .and. & 
                  ptcle_opt_model(i) .ne. 'phm' .and. &
                  ptcle_opt_model(i) .ne. 'ihm' .and. &
                  ptcle_opt_model(i) .ne. 'rhm' .and. &
                  ptcle_opt_model(i) .ne. 'baum') THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' Undefined optical model for particle ', TRIM(ADJUSTL(ptcle_type(i)))
                WRITE(*,*) ' Check file: ', TRIM(ADJUSTL(filename_tmp)) 
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             select case (ptcle_opt_model(i))
                ! read Baum parameter
             case('baum')
                ptcle_material(i) = undef_c
                CALL SKIP_COM(fprop%unit)
                READ(fprop%unit, FMT=*) the_char
                if (flag_baum) then
                   if (baum_type.ne. LOWCASE( TRIM(ADJUSTL(the_char)), len(TRIM(ADJUSTL(the_char)))) )  then
                      WRITE(*,*) ' '
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (sub. read_data) :  ERROR              '
                      WRITE(*,*) '   Only one Baum type can be used  '
                      WRITE(*,*) '   within a single ARTDECO run            '
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' '
                      STOP
                   endif
                else
                   baum_type = LOWCASE(TRIM(ADJUSTL(the_char)), len(TRIM(ADJUSTL(the_char))))
                   if ( baum_type.ne.'ghm' .and. &
                        baum_type.ne.'asc' .and. &
                        baum_type.ne.'sc' ) then
                      WRITE(*,*) ' '
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (sub. read_data) :  ERROR               '
                      WRITE(*,*) ' Baum type must be one of the following :'
                      WRITE(*,*) '   asc : Aggregate Solid Columns         '
                      WRITE(*,*) '   ghm : General Habit Mixture           '
                      WRITE(*,*) '   sc  : Solid Columns                   '
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' '
                      STOP                      
                   endif
                endif
                CALL SKIP_COM(fprop%unit)
                READ(fprop%unit, FMT=*) ptcle_param_baum(i)
                flag_baum = .true. 
                ! read parameters for mie3
             case('mie3')
                do k = 1, mie3_npar
                   CALL SKIP_COM(fprop%unit)
                   READ(fprop%unit, FMT=*) ptcle_param_mie3(i,k)
                enddo
                flag_mie = .true. 
                ! read parameters for phm
             case ('phm')
                CALL SKIP_COM(fprop%unit)
                READ(fprop%unit, FMT=*) ptcle_nphot_hm(i)
                do k = 1, phm_npar
                   CALL SKIP_COM(fprop%unit)
                   READ(fprop%unit, FMT=*) ptcle_param_hm(i,k)
                enddo
                flag_hm = .true. 
                ! read parameters for phm
             case('rhm') 
                CALL SKIP_COM(fprop%unit)
                READ(fprop%unit, FMT=*) ptcle_nphot_hm(i)
                do k = 1, rhm_npar
                   CALL SKIP_COM(fprop%unit)
                   READ(fprop%unit, FMT=*) ptcle_param_hm(i,k)
                enddo
                flag_hm = .true. 
             case('ihm') 
                CALL SKIP_COM(fprop%unit)
                READ(fprop%unit, FMT=*) ptcle_nphot_hm(i)
                do k = 1, ihm_npar
                   CALL SKIP_COM(fprop%unit)
                   READ(fprop%unit, FMT=*) ptcle_param_hm(i,k)
                enddo
                flag_hm = .true. 
             end select
             CLOSE(UNIT = fprop%unit)
          ENDIF
          ! End read prop_[...].dat file
       END DO

       ! We may need to read material refractive index in case optical properties has to be computed
       ! How many different material do we need ?
       ! NOTE : if ptcle_material .ne. undef_c, that means that some
       !        optical properties will be computed
       DO i = 1, nptcle 
          if (ptcle_material(i) .ne. undef_c) THEN
             if (i .eq. 1) then
                nmaterial = nmaterial + 1 
             else
                flag=.true.
                DO j = 1, i-1 
                   if (ptcle_material(i) .eq. ptcle_material(j)) flag = .false.
                ENDDO
                if (flag) nmaterial = nmaterial + 1 
             endif
          ENDIF
       ENDDO
       ! then read those refractive index if nmaterial>0
       if (nmaterial.gt.0) then
          if (verbose) then
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read refractive index files '
          endif
          ! get material names
          ALLOCATE(material_name(nmaterial))
          material_name(:)         = undef_c
          imaterial = 0
          DO i = 1, nptcle 
             if (ptcle_material(i) .ne. undef_c) THEN
                if (i .eq. 1) then
                   imaterial = imaterial +1
                   material_name(imaterial) = ptcle_material(i)
                else
                   flag=.true.
                   DO j = 1, i-1 
                      if (ptcle_material(i) .eq. ptcle_material(j)) flag = .false.
                   ENDDO
                   if (flag) then
                      imaterial = imaterial +1
                      material_name(imaterial) = ptcle_material(i)
                   endif
                endif
             ENDIF
          ENDDO
          ! get nlambmax_refind
          nlambmax_refind = 0
          DO i = 1, nmaterial
             filename_tmp = TRIMCAT(artdeco_path, dir_refind)
             filename_tmp = TRIMCAT(filename_tmp, frefind%name)
             filename_tmp = TRIMCAT(filename_tmp, material_name(i))
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             INQUIRE(FILE=filename_tmp, EXIST=test_file)
             IF ( test_file .eqv. .FALSE. ) THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' You need to provide a refractive index file '
                WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' for material ', TRIM(ADJUSTL(material_name(i)))
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             OPEN (UNIT = frefind%unit, FILE = filename_tmp, STATUS = 'OLD')
             CALL SKIP_COM(frefind%unit)
             READ(frefind%unit, FMT=*) ntmp
             CLOSE(frefind%unit)
             if (ntmp.gt.nlambmax_refind) nlambmax_refind = ntmp
          END DO
          ALLOCATE(material_refind(nmaterial, 3, nlambmax_refind), &
               material_nlamb_refind(nmaterial))
          material_refind(:,:,:)   = undef_dp
          material_nlamb_refind(:) = undef_i
          ! read the corresponding file
          DO i = 1, nmaterial
             filename_tmp = TRIMCAT(artdeco_path, dir_refind)
             filename_tmp = TRIMCAT(filename_tmp, frefind%name)
             filename_tmp = TRIMCAT(filename_tmp, material_name(i))
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             IF (verbose) THEN
                WRITE(*,*)  '                   --> read ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             OPEN (UNIT = frefind%unit, FILE = filename_tmp, STATUS = 'OLD')
             CALL SKIP_COM(frefind%unit)
             READ(frefind%unit, FMT=*) material_nlamb_refind(i)
             CALL SKIP_COM(frefind%unit)
             do k = 1,  material_nlamb_refind(i)
                READ(UNIT=frefind%unit, FMT=*) material_refind(i,1,k),material_refind(i,2,k),material_refind(i,3,k)
                if (k .gt.1) then
                   if (material_refind(i,1,k) .le. material_refind(i,1,k-1)) then
                      WRITE(*,*) ''
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (in sub. read_data)  : ERROR '
                      WRITE(*,*) ' In refractive index file ', TRIM(ADJUSTL(filename_tmp))
                      WRITE(*,*) ' Wavelengths must be sorted in increasing order '
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ''
                      STOP                                           
                   endif
                endif
             enddo
             CLOSE(frefind%unit)
          ENDDO
       endif ! (nmaterial.gt.0)

    END IF! nptcle>0

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! if new Baum cirrus optical properties are required, we
    ! need to read the Baum LUT
    if (flag_baum) then
       if (verbose) then
          WRITE(*,*) ' '
          WRITE(*,*) '  (sub. read_data) read Baum cirrus particles LUT  '
       endif

       ! read the effective diameter grid
       filename_tmp = TRIMCAT(artdeco_path, dir_baum)
       filename_tmp = TRIMCAT(filename_tmp, baum_type)
       filename_tmp = TRIMCAT(filename_tmp, '/')
       filename_tmp = TRIMCAT(filename_tmp, fbaum_de%name)
       OPEN (UNIT = fbaum_de%unit, FILE = filename_tmp, STATUS = 'OLD')
       CALL SKIP_COM(fbaum_de%unit)
       READ(fbaum_de%unit, FMT=*)  baum_nde_file
       ALLOCATE(baum_de_file(baum_nde_file))
       CALL SKIP_COM(fbaum_de%unit)
       do i = 1, baum_nde_file
          read(fbaum_de%unit, FMT=*) baum_de_file(i)
          if (i.gt.1) then
             if (baum_de_file(i) .le. baum_de_file(i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' Effective diameters must be sorted in increasing order '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP                                           
             endif
          endif
       end do
       CLOSE(fbaum_de%unit)
       ! read the wavelength grid
       filename_tmp = TRIMCAT(artdeco_path, dir_baum)
       filename_tmp = TRIMCAT(filename_tmp, baum_type)
       filename_tmp = TRIMCAT(filename_tmp, '/')
       filename_tmp = TRIMCAT(filename_tmp, fbaum_wvl%name)
       OPEN (UNIT = fbaum_wvl%unit, FILE = filename_tmp, STATUS = 'OLD')
       CALL SKIP_COM(fbaum_wvl%unit)
       READ(fbaum_wvl%unit, FMT=*)  baum_nlamb_file
       ALLOCATE(baum_lamb_file(baum_nlamb_file))
       CALL SKIP_COM(fbaum_wvl%unit)
       do i = 1, baum_nlamb_file
          read(fbaum_wvl%unit, FMT=*) baum_lamb_file(i)
          if (i.gt.1) then
             if (baum_lamb_file(i) .le. baum_lamb_file(i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' Wavelengths must be sorted in increasing order '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP                                           
             endif
          endif
       end do
       CLOSE(fbaum_wvl%unit)

       ! we look for the files whose LUT
       ! will be needed for further interpolation
       ! We will only read the pair of files
       ! whose values surround the required values 
       ! (we restrict the De and lambda sample to read)
       baum_nide   = 2*nptcle
       if (ireflamb.eq.0) then
          baum_nilamb = 2*nlambda + 2 ! two more wavlength to surround the ref. lamb.
       else
          baum_nilamb = 2*nlambda
       endif
       ALLOCATE(baum_ide(baum_nide), baum_ilamb(baum_nilamb))
       baum_ide(:)   = -1
       baum_ilamb(:) = -1
       do i = 1, nptcle
          if (ptcle_opt_model(i).eq.'baum') then
             if ( (ptcle_param_baum(i).gt.baum_de_file(baum_nde_file)) .or. &
                  (ptcle_param_baum(i).lt.baum_de_file(1))           ) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' For particle ', TRIM(ADJUSTL(ptcle_type(i)))
                WRITE(*,*) ' The required effective diameter is out of'
                WRITE(*,*) ' the LUT range '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP                                           
             endif
             do j = 1, baum_nde_file-1
                if ( (ptcle_param_baum(i).ge.baum_de_file(j)) .and. &
                     (ptcle_param_baum(i).lt.baum_de_file(j+1)) ) then
                   baum_ide(2*i-1)  = j 
                   baum_ide(2*i)    = j+1
                endif
             end do
             if (ptcle_param_baum(i).eq.baum_de_file(baum_nde_file)) then
                baum_ide(2*i-1) = baum_nde_file - 1
                baum_ide(2*i)   = baum_nde_file 
             endif
          endif
          ! print*, TRIM(ADJUSTL(ptcle_type(i))), baum_ide(i)
       end do
       ! check the uniqueness
       !print*, baum_ide
       do i = 1, baum_nide
          do j = i+1, baum_nide
             if ( baum_ide(i).eq.(baum_ide(j)) ) baum_ide(j) = -1
          end do
       end do
       !print*, baum_ide
       nde_baum = 0
       do i = 1, baum_nide
          if (baum_ide(i).ne.-1) nde_baum = nde_baum + 1 
       end do
       CALL SORT(baum_ide)
       !print*, nde_baum
       !print*, baum_ide

       ! look for good wavelength bins
       do i = 1, nlambda
          if ( (wlambda(i).gt.baum_lamb_file(baum_nlamb_file)) .or. &
               (wlambda(i).lt.baum_lamb_file(1))           ) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' For Baum particle properties interpolation '
             WRITE(*,*) ' The required wavelength is out of'
             WRITE(*,*) ' the LUT range '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP                                           
          endif
          do j = 1, baum_nlamb_file-1
             if ( (wlambda(i).ge.baum_lamb_file(j)) .and. &
                  (wlambda(i).lt.baum_lamb_file(j+1)) ) then
                baum_ilamb(2*i-1) = j 
                baum_ilamb(2*i)   = j+1 
             endif
          end do
          if (wlambda(i).eq.baum_lamb_file(baum_nlamb_file)) then
             baum_ilamb(2*i-1) = baum_nlamb_file - 1
             baum_ilamb(2*i)   = baum_nlamb_file
          endif
          !print*, wlambda(i), baum_ilamb(i)             
       end do
       ! print*, '' 
       ! print*, baum_ilamb
       ! Same for reflamb 
       if (ireflamb.eq.0) then
          if ( (reflamb.gt.baum_lamb_file(baum_nlamb_file)) .or. &
               (reflamb.lt.baum_lamb_file(1))           ) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' For Baum particle properties interpolation '
             WRITE(*,*) ' The required wavelength (reflamb) is out of'
             WRITE(*,*) ' the LUT range '
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP                                           
          endif
          do j = 1, baum_nlamb_file-1
             if ( (reflamb.ge.baum_lamb_file(j)) .and. &
                  (reflamb.lt.baum_lamb_file(j+1)) ) then
                baum_ilamb(2*i-1) = j 
                baum_ilamb(2*i)   = j+1 
             endif
          end do
          if (reflamb.eq.baum_lamb_file(baum_nlamb_file)) then
             baum_ilamb(2*i-1) = baum_nlamb_file - 1
             baum_ilamb(2*i)   = baum_nlamb_file
          endif
       end if
       ! check the uniqueness
       do i = 1, baum_nilamb
          do j = i+1, baum_nilamb
             if ( baum_ilamb(i).eq.(baum_ilamb(j)) ) baum_ilamb(j) = -1
          end do
       end do
       !print*, baum_ilamb
       nlamb_baum = 0
       do i = 1, baum_nilamb
          if (baum_ilamb(i).ne.-1) nlamb_baum = nlamb_baum + 1 
       end do
       CALL SORT(baum_ilamb)
       ! print*, nlamb_baum
       ! print*, baum_ilamb
       ! STOP
       ! Read the LUT good bins
       ALLOCATE(lamb_baum(nlamb_baum), &
            de_baum(nde_baum), &
            Cext_baum(nlamb_baum, nde_baum), &
            ssa_baum(nlamb_baum, nde_baum), &
            g_baum(nlamb_baum, nde_baum) )

       ! read teta grid
       filename_tmp = TRIMCAT(artdeco_path, dir_baum)
       filename_tmp = TRIMCAT(filename_tmp, baum_type)
       filename_tmp = TRIMCAT(filename_tmp, '/')
       filename_tmp = TRIMCAT(filename_tmp, fbaum_mu%name)
       OPEN (UNIT = fbaum_mu%unit, FILE = filename_tmp, STATUS = 'OLD')
       CALL SKIP_COM(fbaum_mu%unit)
       READ(fbaum_mu%unit, FMT=*) nmu_baum
       ! print*, nmu_baum
       ALLOCATE(mu_baum(nmu_baum), &
            ptcle_phasemat_baum(nlamb_baum, nde_baum, 6, nmu_baum))
       CALL SKIP_COM(fbaum_mu%unit)
       do i = 1, nmu_baum
          read(fbaum_mu%unit, FMT=*) mu_baum(i)
          if (i.gt.1) then
             if (mu_baum(i) .le. mu_baum(i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' Scattering angles must be sorted in increasing order '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP                                           
             endif
          endif
       end do
       CLOSE(fbaum_mu%unit)

       ! print*, nde_baum
       ! print*, baum_ide
       ! print*, ' '
       ! print*, nlamb_baum
       ! print*, baum_ilamb
       ! print*, ' '

       ! read the optical properties LUT
       do i = 1, nlamb_baum
          iwvl = baum_ilamb(baum_nilamb - nlamb_baum + i)
          lamb_baum(i) = baum_lamb_file(iwvl)

          do j = 1, nde_baum
             ide = baum_ide(baum_nide - nde_baum + j)
             de_baum(j) = baum_de_file(ide)

             ! print*, i, j, iwvl, ide, lamb_baum(i), de_baum(j)
             filename_tmp = TRIMCAT(artdeco_path, dir_baum)
             filename_tmp = TRIMCAT(filename_tmp, baum_type)
             filename_tmp = TRIMCAT(filename_tmp, '/')
             filename_tmp = TRIMCAT(filename_tmp, fbaum%name)
             filename_tmp = TRIMCAT(filename_tmp, baum_type)
             filename_tmp = TRIMCAT(filename_tmp, '_')
             filename_tmp = TRIMCAT(filename_tmp, 'wvl')
             WRITE(the_char,FMT=*) iwvl
             filename_tmp = TRIMCAT(filename_tmp,TRIM(ADJUSTL(the_char)))
             filename_tmp = TRIMCAT(filename_tmp, '_de')
             WRITE(the_char,FMT=*) ide
             filename_tmp = TRIMCAT(filename_tmp,TRIM(ADJUSTL(the_char)))
             filename_tmp = TRIMCAT(filename_tmp,'.dat')

             OPEN (UNIT = fbaum%unit, FILE = filename_tmp, STATUS = 'OLD')
             CALL SKIP_COM(fbaum%unit)
             READ(fbaum%unit, FMT=*) tmp_dp, Cext_baum(i,j), ssa_baum(i,j), g_baum(i,j)
             ! print*, tmp_dp, Cext_baum(i,j), ssa_baum(i,j), g_baum(i,j)
             CALL SKIP_COM(fbaum%unit)
             do k = 1, nmu_baum
                READ(fbaum%unit, FMT=*) ptcle_phasemat_baum(i, j, 1, k), &
                     ptcle_phasemat_baum(i, j, 2, k), &
                     ptcle_phasemat_baum(i, j, 3, k), &
                     ptcle_phasemat_baum(i, j, 4, k), &
                     ptcle_phasemat_baum(i, j, 5, k), &
                     ptcle_phasemat_baum(i, j, 6, k)
             end do
             CLOSE(UNIT=fbaum%unit)
          end do
       end do

       DEALLOCATE(baum_de_file)
       DEALLOCATE(baum_lamb_file)
       DEALLOCATE(baum_ide, baum_ilamb)

    endif

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! READ TRUNCATION MODEL SPEC FILES
    if (opt_only.eqv..false.) then

       ! read the Potter truncation specifications
       if ( trunc_method .eq. 'potter')  then      
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE(*,*) '  (sub. read_data)  read potter_spec.dat'
          ENDIF
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fpotter_spec%name)
          OPEN (UNIT = fpotter_spec%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fpotter_spec%unit)
          READ(fpotter_spec%unit, FMT=*) potter_theta_min, potter_theta_max
          IF (potter_theta_min .eq. potter_theta_max) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' potter_theta_min must be different than potter_theta_max'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif
          IF (potter_theta_min .gt. potter_theta_max) then
             tmp_dp = potter_theta_max
             potter_theta_max = potter_theta_min
             potter_theta_min = tmp_dp
          endif
          CLOSE(UNIT = fpotter_spec%unit)
       endif
       if ( trunc_method .eq. 'dfit')  then      
          IF (verbose) THEN
             WRITE(*,*) ''
             WRITE (*,*) '  (sub. read_data) read dfit_spec.dat'
          ENDIF
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fdfit_spec%name)
          OPEN (UNIT = fdfit_spec%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fdfit_spec%unit)
          READ(fdfit_spec%unit, FMT=*) dfit_thetac
          IF ((dfit_thetac .gt. 180.0_dp ) .or. (dfit_thetac .lt. 0.0_dp )) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' In d-fit truncation thetac must be > 0.0 and < 180.0'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif
          CALL SKIP_COM(fdfit_spec%unit)
          READ(fdfit_spec%unit, FMT=*) i
          IF ((i .ne. 1) .and. (i .ne. 0 )) then
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' In d-fit -fitall- must be 0 or 1'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          endif
          if (i .eq. 0) then
             dfit_fitall = .false.
          else
             dfit_fitall = .true.
          endif
          CLOSE(UNIT = fdfit_spec%unit)
       endif
    end if

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! GET ATMSOPHERIC GAS CHARACTERISTICS
    if (do_rt) then

       ALLOCATE(depol(nlambda))
       depol(:) = undef_dp
              
       ntmp = 0
       if (mode.eq.'kdis') then

          ! Read eventual uniformly mixed gas concentration:
          filename_tmp = TRIMCAT (artdeco_path, dir_input)
          filename_tmp = TRIMCAT (filename_tmp, fgasconc%name)      
          the_char     = filename_tmp 
          INQUIRE(FILE=filename_tmp, EXIST=test_file)
          IF ( test_file .eqv. .FALSE. ) THEN
             IF (verbose) THEN
                WRITE (*,*) '' 
                WRITE (*,*)'  (sub. read_data) no file ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             ntmp = 0
          ELSE
             IF (verbose) THEN
                WRITE (*,*) '' 
                WRITE (*,*)'  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             OPEN (UNIT = fgasconc%unit, FILE = filename_tmp, STATUS = 'old')
             CALL SKIP_COM(fgasconc%unit)
             READ (UNIT=fgasconc%unit, FMT=*) ntmp
             CALL SKIP_COM(fgasconc%unit)
             if (ntmp.gt.0) then
                ALLOCATE(char_arr_tmp(ntmp),dbl_arr_tmp(ntmp))
                DO i=1, ntmp
                   READ (UNIT=fgasconc%unit, FMT=*) char_arr_tmp(i),dbl_arr_tmp(i)
                   char_arr_tmp(i)=LOWCASE(TRIM(ADJUSTL(char_arr_tmp(i))), len(TRIM(ADJUSTL(char_arr_tmp(i)))))
                ENDDO
             endif
             CLOSE(fgasconc%unit)
          endif
       endif! kdis

       ! read atmospheric profil
       filename_tmp = TRIMCAT (artdeco_path, dir_atm)
       filename_tmp = TRIMCAT (filename_tmp, fatm%name)
       filename_tmp = TRIMCAT (filename_tmp, atm_type)
       filename_tmp = TRIMCAT(filename_tmp, '.dat')
       INQUIRE(FILE=filename_tmp, EXIST=test_file)
       IF ( test_file .eqv. .FALSE. ) THEN
          WRITE(*,*) ''
          WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          WRITE(*,*) ' (in sub. read_data)  : ERROR '
          WRITE(*,*) ' You need to provide an atmosphere definition file '
          WRITE(*,*) ' for  ', TRIM(ADJUSTL(atm_type))
          WRITE(*,*) ' Required file is :'
          WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
          WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          WRITE(*,*) ''
          STOP 
       ENDIF
       IF (verbose) THEN
          WRITE (*,*) '' 
          WRITE (*,*)'  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
       ENDIF
       OPEN (UNIT = fatm%unit, FILE = filename_tmp, STATUS = 'old')
       CALL SKIP_COM(fatm%unit)
       READ(fatm%unit, FMT=*) nalt_atm
       CALL SKIP_COM(fatm%unit)
       READ(fatm%unit, FMT=*) ngas_u_atm
       ALLOCATE(alt_atm(nalt_atm), &
            p_atm(nalt_atm), &
            t_atm(nalt_atm))
       alt_atm(:) = undef_dp
       p_atm(:)   = undef_dp
       t_atm(:)   = undef_dp
       if (mode.eq.'kdis') then
          ALLOCATE(u_air_atm(nalt_atm), &
               gas_u_atm(ngas_u_atm+ntmp), &
               u_atm(ngas_u_atm+ntmp, nalt_atm))
          u_air_atm(:) = undef_dp
          gas_u_atm(:) = undef_c
          u_atm(:,:)   = undef_dp
       endif
       CALL SKIP_COM(fatm%unit)
       
       if (ngas_u_atm .gt. 0) then
          DO i = 1, ngas_u_atm
             if (mode.eq.'kdis') then
                READ(UNIT=fatm%unit, FMT=*) gas_u_atm(i)
                gas_u_atm(i)=LOWCASE(TRIM(ADJUSTL(gas_u_atm(i))), len(TRIM(ADJUSTL(gas_u_atm(i)))))
                if (ntmp.gt.0) then
                   do j=1,ntmp
                      IF(gas_u_atm(i).eq.char_arr_tmp(j))THEN
                         WRITE(*,*) ''
                         WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                         WRITE(*,*) ' (in sub. read_data)  : ERROR '
                         WRITE(*,*) ' Conflicting set-up for the gas specy  ', gas_u_atm(i)
                         WRITE(*,*) ' It appears twice (in two files) :'
                         WRITE(*,*) TRIM(ADJUSTL(filename_tmp))
                         WRITE(*,*) ' and'
                         WRITE(*,*) TRIM(ADJUSTL(the_char))
                         WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                         WRITE(*,*) ''
                         STOP 
                      END IF
                   enddo
                endif
             else
                READ(UNIT=fatm%unit, FMT=*) the_char
             endif
          ENDDO
       else
          READ(UNIT=fatm%unit, FMT=*) the_char
       endif
       CALL SKIP_COM(fatm%unit)

       DO i=1, nalt_atm
          if (mode.eq.'kdis') then
             READ (UNIT=fatm%unit, FMT=*) alt_atm(i), p_atm(i), t_atm(i), u_air_atm(i), (u_atm(j,i), j=1,ngas_u_atm)
          else
             READ (UNIT=fatm%unit, FMT=*) alt_atm(i), p_atm(i), t_atm(i)
          endif
          IF (i .gt. 1) THEN
             ! check atmospheric profil sorting
             IF ( alt_atm(i) .ge. alt_atm(i-1) ) THEN
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file :', TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' Altitudes must be unique and in decreasing order'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
          ENDIF
       ENDDO
       CLOSE(UNIT=fatm%unit)

       ! Add concentration profiles for uniformly mixed gas
       if (mode.eq.'kdis') then
          do i = 1, ntmp
             gas_u_atm(i+ngas_u_atm) = char_arr_tmp(i)
             u_atm(i+ngas_u_atm,:)   = u_air_atm(:) * (dbl_arr_tmp(i)*1.0e-6)
          enddo
          ngas_u_atm = ngas_u_atm + ntmp
          if (ntmp.gt.0) DEALLOCATE(char_arr_tmp,dbl_arr_tmp)
       endif

       ! set the number of layers for RTE solvers
       nlayers     = nalt_atm - 1

       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       ! eventually read particles vertical distribution
       if (flag_user_vdist) then
          if (verbose) then
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read particle vertical distribution files'
             WRITE(*,*) ' '
          ENDIF
          do i = 1, nptcle
             if (ptcle_vdist_type(i).eq.'user') then
                filename_tmp = TRIMCAT(artdeco_path, dir_vdist)
                filename_tmp = TRIMCAT(filename_tmp, fvdist%name)
                filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
                filename_tmp = TRIMCAT(filename_tmp, '.dat')
                INQUIRE(FILE=filename_tmp, EXIST=test_file)
                IF ( test_file .eqv. .FALSE. ) THEN
                   WRITE(*,*) ''
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) ' You need to provide a file for vertical distribution  '
                   WRITE(*,*) ' of  ', TRIM(ADJUSTL(ptcle_type(i)))
                   WRITE(*,*) ' Required file is :'
                   WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                   STOP 
                ENDIF
                IF (verbose) THEN
                   WRITE(*,*)  '    --> read  ', TRIM(ADJUSTL(filename_tmp))
                ENDIF
                OPEN (UNIT = fvdist%unit, FILE = filename_tmp, STATUS = 'old')
                CALL SKIP_COM(fvdist%unit) 
                READ (UNIT=fvdist%unit,FMT=*) ptcle_user_vdist_nalt(i)
                IF(ptcle_user_vdist_nalt(i).gt.max_user_vdist_nalt) then 
                   WRITE(*,*) '  '
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) '    The number of sampled altitude in file ', TRIM(filename_tmp)
                   WRITE(*,*) '    must be <= max_user_vdist_nalt defined in artdecot_constant.90'
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) '  '
                   STOP
                ENDIF
                CALL SKIP_COM(fvdist%unit) 
                READ (UNIT=fvdist%unit,FMT=*) ptcle_user_vdist_alt(i,1), ptcle_user_vdist(i,1)   
                DO k = 2, ptcle_user_vdist_nalt(i)
                   READ (UNIT=fvdist%unit,FMT=*) ptcle_user_vdist_alt(i,k), ptcle_user_vdist(i,k)   
                   ! check altitude sorting
                   IF ( ptcle_user_vdist_alt(i,k) .ge. ptcle_user_vdist_alt(i,k-1) ) THEN
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (in sub. read_data)  : ERROR '
                      WRITE(*,*) ' In file :', TRIM(ADJUSTL(filename_tmp))
                      WRITE(*,*) ' Altitudes must be unique and in decreasing order'
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ''
                      STOP 
                   ENDIF
                END DO
                CLOSE(UNIT = fvdist%unit)
                IF ( ptcle_user_vdist_alt(i,ptcle_user_vdist_nalt(i)).ne.0.0_dp ) THEN
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) ' In file :', TRIM(ADJUSTL(filename_tmp))
                   WRITE(*,*) ' Altitudes must start at 0.0 '
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                   STOP 
                ENDIF
                IF (ptcle_user_vdist_alt(i,1).lt.MAXVAL(alt_atm)) THEN
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) ' In file :', TRIM(ADJUSTL(filename_tmp))
                   WRITE(*,*) ' Maximum sampled altitude should be '
                   WRITE(*,*) ' greater or equal to maximum altitude in atmosphere '
                   WRITE(*,*) ' definition file '
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                   STOP 
                ENDIF
                ! we sort the alt in increasing order for further integration
                ptcle_user_vdist_alt(i,1:ptcle_user_vdist_nalt(i)) = &
                     REVERSE(ptcle_user_vdist_nalt(i), ptcle_user_vdist_alt(i,1:ptcle_user_vdist_nalt(i)))
                ptcle_user_vdist(i,1:ptcle_user_vdist_nalt(i)) = &
                     REVERSE(ptcle_user_vdist_nalt(i), ptcle_user_vdist(i,1:ptcle_user_vdist_nalt(i)))
             endif
          enddo
       endif

    endif ! do_rt

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !    SURFACE PROPERTIES
    if (do_rt) then

       SELECT CASE(surface_family)
       CASE('lambert')

          ALLOCATE(surface_albedo(nlambda))
          surface_albedo(:) = undef_dp

          if (surface_type.eq.'cste') then
             surface_albedo(:) = cste_alb
          else
             ! read and set the albedo for the lambertian surface
             filename_tmp = TRIMCAT(artdeco_path, dir_surf)
             filename_tmp = TRIMCAT(filename_tmp, fsurfalb%name)
             filename_tmp = TRIMCAT(filename_tmp, surface_type)
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             IF (verbose) THEN
                WRITE(*,*) ' '
                WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             INQUIRE(FILE=filename_tmp, EXIST=test_file)
             IF ( test_file .eqv. .FALSE. ) THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' You need to provide an albedo definition file '
                WRITE(*,*) ' for lambertian surface ', TRIM(ADJUSTL(surface_type))
                WRITE(*,*) ' Required file is :'
                WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             OPEN (UNIT = fsurfalb%unit, FILE = filename_tmp, STATUS = 'old')
             CALL SKIP_COM(fsurfalb%unit)
             READ (UNIT=fsurfalb%unit,FMT=*) nlamb
             allocate(lamb(nlamb), salb(nlamb))
             CALL SKIP_COM(fsurfalb%unit)
             ! read first value out of the loop
             READ (UNIT=fsurfalb%unit,FMT=*) lamb(1), salb(1)
             do i = 2 , nlamb
                READ (UNIT=fsurfalb%unit,FMT=*) lamb(i), salb(i)
                ! test whether wavelengths are increasing
                if (lamb(i).lt.lamb(i-1)) then
                   WRITE(*,*) ''
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                   WRITE(*,*) ' Wavelengths must be sorted in increasing order'
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                   STOP                    
                end if
             end do
             CLOSE(unit=fsurfalb%unit) 
             ! interpolate to the working wavelength
             IF (verbose) THEN
                WRITE(*,*) ' '
                WRITE(*,*) '                    Interpolate the surface albedo to working wavelengths'
             ENDIF
             do j = 1, nlambda
                if ((wlambda(j).lt.minval(lamb)) .or. (wlambda(j).gt.maxval(lamb))) then
                   WRITE(*,*) ''
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ' (in sub. read_data)  : ERROR '
                   WRITE(*,*) ' Working wavelength is out of the defined range in',TRIM(ADJUSTL(filename_tmp))
                   WRITE(*,FMT='(2x,A21,2x,F12.5,1x,A8)')  &
                        &'Working wavelenght = ', wlambda(j), 'microns'
                   WRITE(*,FMT='(2x,A21,2x,F12.5,1x,A8)')  &
                        &'Min defined        = ', minval(lamb), 'microns'
                   WRITE(*,FMT='(2x,A21,2x,F12.5,1x,A8)')  &
                        &'Max defined        = ', maxval(lamb), 'microns'
                   WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                   WRITE(*,*) ''
                   STOP                    
                endif
                surface_albedo(j) = LININTPOL(salb, lamb, nlamb, wlambda(j))
             end do
             deallocate(lamb, salb)
          endif

       CASE('brdf')

          IF (surface_type .eq. 'ocean') then

             surface_brdf_model = surface_type
             ! open ocean_spec.dat file
             filename_tmp = TRIMCAT (artdeco_path, dir_input)
             filename_tmp = TRIMCAT (filename_tmp, focean_spec%name)
             if (verbose) then
                WRITE(*,*) ''
                WRITE(*,*) '  (sub. read_data) read ', TRIM(filename_tmp)
             endif
             OPEN (UNIT = focean_spec%unit, FILE = filename_tmp, STATUS = 'old')
             CALL SKIP_COM(focean_spec%unit)
             READ (UNIT=focean_spec%unit,FMT=*) wind_spd
             CALL SKIP_COM(focean_spec%unit)
             READ (UNIT=focean_spec%unit,FMT=*) ocean_xsal
             CALL SKIP_COM(focean_spec%unit)
             READ (UNIT=focean_spec%unit,FMT=*) ocean_pcl
             CALL SKIP_COM(focean_spec%unit)
             READ (UNIT=focean_spec%unit,FMT=*) the_char
             the_char = TRIM(ADJUSTL(the_char))
             if (the_char.eq.'yes')  then
                ocean_shadow = .TRUE.
             elseif (the_char.eq.'no') then
                ocean_shadow = .FALSE.
             else
                write(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  :  ERROR          '
                write(*,*) '   In ocean_spec.dat :                  '
                write(*,*) '   Whether we turn on/off the shadowing '
                write(*,*) '   must be set to "yes" or "no"         '
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) ''
                STOP
             endif

             CLOSE(unit=focean_spec%unit)
             ocean_distr = 3 ! only isotropic surface implemented in RT models

             ALLOCATE(ocean_Rsw(nlambda))
             ocean_Rsw(:) = -1 

          else

             filename_tmp = TRIMCAT(artdeco_path, dir_surf)
             filename_tmp = TRIMCAT(filename_tmp, fsurfbrdf%name)
             filename_tmp = TRIMCAT(filename_tmp, surface_type)
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             IF (verbose) THEN
                WRITE(*,*) ' '
                WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
             ENDIF
             INQUIRE(FILE=filename_tmp, EXIST=test_file)
             IF ( test_file .eqv. .FALSE. ) THEN
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) '  (in sub. read_data)  : ERROR '
                WRITE(*,*) '  You need to provide an BRDF definition file '
                WRITE(*,*) '  for surface ', TRIM(ADJUSTL(surface_type))
                WRITE(*,*) '  Required file is :'
                WRITE(*,*)   TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP 
             ENDIF
             OPEN (UNIT = fsurfbrdf%unit, FILE = filename_tmp, STATUS = 'old')
             CALL SKIP_COM(fsurfbrdf%unit)
             READ (UNIT=fsurfbrdf%unit,FMT=*) surface_brdf_model
             surface_brdf_model = TRIM(ADJUSTL(surface_brdf_model))
             CALL SKIP_COM(fsurfbrdf%unit)
             READ (UNIT=fsurfbrdf%unit,FMT=*) nlamb
             CALL SKIP_COM(fsurfbrdf%unit)

             select case (surface_brdf_model)

             case('lisp_ross')

                ! local variables
                ALLOCATE(lamb(nlamb), &
                     klamb(nlamb), &
                     k2(nlamb), &
                     k1(nlamb), &
                     cv(nlamb), &
                     refind_r(nlamb), &
                     refind_i(nlamb))
                lamb(:) = undef_dp
                klamb(:) = undef_dp
                k2(:) = undef_dp
                k1(:) = undef_dp
                cv(:) = undef_dp
                refind_r(:) = undef_dp
                refind_i(:) = undef_dp
                READ (UNIT=fsurfbrdf%unit,FMT=*) lamb(1), &
                     klamb(1), &
                     k2(1), &
                     k1(1), &
                     cv(1), &
                     refind_r(1), &
                     refind_i(1)
                do j = 2, nlamb
                   READ (UNIT=fsurfbrdf%unit,FMT=*) lamb(j), &
                        klamb(j), &
                        k2(j), &
                        k1(j), &
                        cv(j), &
                        refind_r(j), &
                        refind_i(j)
                   if (lamb(j).lt.lamb(j-1)) then
                      WRITE(*,*) ''
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (in sub. read_data)  : ERROR '
                      WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                      WRITE(*,*) ' Wavelengths must be sorted in increasing order'
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ''
                      STOP                    
                   end if
                end do

                ! COMMON variables
                surface_n_par_brdf = 3
                surface_n_par_bpdf = 3
                ALLOCATE(surface_par_brdf(surface_n_par_brdf, nlambda),&
                     surface_par_bpdf(surface_n_par_bpdf, nlambda))
                surface_par_brdf(:, :) = undef_dp
                surface_par_bpdf(:, :) = undef_dp

                ! interpolate to the working wavelength
                IF (verbose) THEN
                   WRITE(*,*) ' '
                   WRITE(*,*) '                    Interpolate the surface brdf param. to working wavelengths'
                ENDIF
                do j = 1, nlambda
                   if ((wlambda(j).lt.minval(lamb)) .or. (wlambda(j).gt.maxval(lamb))) then
                      WRITE(*,*) ''
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ' (in sub. read_data)  : ERROR '
                      WRITE(*,*) ' Working wavelength is out of the defined range in',TRIM(ADJUSTL(filename_tmp))
                      WRITE(*,FMT='(2x,A21,2x,F12.5,1x,A8)')  &
                           &'Working wavelenght = ', wlambda(j), 'microns'
                      WRITE(*,FMT='(2x,A21,2x,F12.5,1x,A8)')  &
                           &'Min defined        = ', minval(lamb), 'microns'
                      WRITE(*,FMT='(2x,A21,2x,F12.5,1x,A8)')  &
                           &'Max defined        = ', maxval(lamb), 'microns'
                      WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                      WRITE(*,*) ''
                      STOP                    
                   endif
                   surface_par_brdf(1, j) = LININTPOL(klamb, lamb, nlamb, wlambda(j))
                   surface_par_brdf(2, j) = LININTPOL(k2, lamb, nlamb, wlambda(j))
                   surface_par_brdf(3, j) = LININTPOL(k1, lamb, nlamb, wlambda(j))
                   surface_par_bpdf(1, j) = LININTPOL(cv, lamb, nlamb, wlambda(j))
                   surface_par_bpdf(2, j) = LININTPOL(refind_r, lamb, nlamb, wlambda(j))
                   surface_par_bpdf(3, j) = LININTPOL(refind_i, lamb, nlamb, wlambda(j))
                end do

                DEALLOCATE(lamb, &
                     klamb, &
                     k2, &
                     k1, &
                     cv, &
                     refind_r, &
                     refind_i)

             case default

                write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                write(*,*) '  (sub. read_data) ERROR '
                write(*,*) '    The BRDF ', TRIM(ADJUSTL(surface_brdf_model))
                write(*,*) '    is not implemented'
                write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                STOP

             end select

             CLOSE(unit=fsurfbrdf%unit) 

          endif

       END SELECT

       ! implement surface temperature if needed
       IF (surface_temp .eq. -1) THEN
          IF (  alt_atm(nalt_atm) .ne. 0.0_dp ) THEN
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) "  The surface temperature can not be set to the atmosphere's bottom one"
             WRITE(*,*) '  T_atm(0km) is not sampled in the atmosphere definition file'
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP
          ELSE
             surface_temp = t_atm(nalt_atm)
          ENDIF
       ENDIF
    endif! do_rt

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! READ INCIDENT RADIATION FIELD
    if ((do_rt).and.(mode.ne.'mono'))then
       kdis_solrad = .FALSE.
       filename_tmp = TRIMCAT(artdeco_path, dir_solrad)
       filename_tmp = TRIMCAT(filename_tmp, fsolrad%name)
       filename_tmp = TRIMCAT(filename_tmp, 'kdis_')
       filename_tmp = TRIMCAT(filename_tmp, TRIM(ADJUSTL(kdis_model)))
       filename_tmp = TRIMCAT(filename_tmp, '_')
       filename_tmp = TRIMCAT(filename_tmp, solrad)
       filename_tmp = TRIMCAT(filename_tmp, '.dat')
       INQUIRE(FILE=filename_tmp, EXIST=test_file)
       IF ( test_file .eqv. .TRUE. ) THEN
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          kdis_solrad = .TRUE.
          OPEN (UNIT = fsolrad%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fsolrad%unit)
          READ (UNIT=fsolrad%unit,FMT=*) kdis_solrad_solcste
          CALL SKIP_COM(fsolrad%unit)
          READ (UNIT=fsolrad%unit,FMT=*) kdis_solrad_nwvl
          if (kdis_solrad_nwvl.ne.kdis_nwvl) then
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) ' Wavelength range should be the same than in'
             WRITE(*,*) ' file ',TRIM(ADJUSTL(filename_kdis))
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP                    
          end if
          allocate(kdis_solrad_wvl(kdis_solrad_nwvl), kdis_solrad_F0(kdis_solrad_nwvl))
          kdis_solrad_wvl(:) = undef_dp
          kdis_solrad_F0(:)  = undef_dp
          CALL SKIP_COM(fsolrad%unit)
          do i = 1 , kdis_solrad_nwvl
             READ (UNIT=fsolrad%unit,FMT=*) kdis_solrad_wvl(i), kdis_solrad_F0(i)
             if (kdis_solrad_wvl(i).ne.(kdis_wvlband(1,i))) then
                !if ( (kdis_solrad_wvl(i).gt.(kdis_wvlband(1,i)+EPS_LAMB)).or. &
                !     (kdis_solrad_wvl(i).lt.(kdis_wvlband(1,i)-EPS_LAMB)) ) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' Wavelength range should be the same than in'
                WRITE(*,*) ' file ',TRIM(ADJUSTL(filename_kdis))
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP                    
             end if
          end do
          CLOSE(unit=fsolrad%unit) 
       ELSE
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) No file ', TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) '                   Band integrated solar radiation will be computed'
          ENDIF
       ENDIF
       if  (test_file.eqv..false.) THEN
          filename_tmp = TRIMCAT(artdeco_path, dir_solrad)
          filename_tmp = TRIMCAT(filename_tmp, fsolrad%name)
          filename_tmp = TRIMCAT(filename_tmp, solrad)
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          IF (verbose) THEN
             WRITE(*,*) ' '
             WRITE(*,*) '  (sub. read_data) read ', TRIM(ADJUSTL(filename_tmp))
          ENDIF
          INQUIRE(FILE=filename_tmp, EXIST=test_file)
          IF ( test_file .eqv. .FALSE. ) THEN
             WRITE(*,*) ''
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ' (in sub. read_data)  : ERROR '
             WRITE(*,*) ' You need to provide a solar radiation definition file '
             WRITE(*,*) ' Required file is :'
             WRITE(*,*)  TRIM(ADJUSTL(filename_tmp))
             WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(*,*) ''
             STOP 
          ENDIF
          OPEN (UNIT = fsolrad%unit, FILE = filename_tmp, STATUS = 'old')
          CALL SKIP_COM(fsolrad%unit)
          READ (UNIT=fsolrad%unit,FMT=*) solrad_solcste
          CALL SKIP_COM(fsolrad%unit)
          READ (UNIT=fsolrad%unit,FMT=*) solrad_nwvl
          allocate(solrad_wvl(solrad_nwvl), solrad_F0(solrad_nwvl))
          solrad_wvl(:) = undef_dp
          solrad_F0(:)  = undef_dp
          CALL SKIP_COM(fsolrad%unit)
          ! read first value out of the loop
          READ (UNIT=fsolrad%unit,FMT=*) solrad_wvl(1), solrad_F0(1)
          do i = 2 , solrad_nwvl
             READ (UNIT=fsolrad%unit,FMT=*) solrad_wvl(i), solrad_F0(i)
             ! test whether wavelengths are increasing
             if (solrad_wvl(i).lt.solrad_wvl(i-1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' Wavelengths must be sorted in increasing order'
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP                    
             end if
          end do
          CLOSE(unit=fsolrad%unit) 
          if (mode.eq.'kdis') then
             ! check whether the "full" solrad wvl covers the kdis bands 
             if (kdis_wvlband(2,1).lt.solrad_wvl(1)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' The spectrum does not cover the k-dis wavelength range'
                WRITE(*,*) ' in file ',TRIM(ADJUSTL(filename_kdis)) 
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             endif
             if (kdis_wvlband(3,kdis_nwvl).gt.solrad_wvl(solrad_nwvl)) then
                WRITE(*,*) ''
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ' (in sub. read_data)  : ERROR '
                WRITE(*,*) ' In file ',TRIM(ADJUSTL(filename_tmp))
                WRITE(*,*) ' The spectrum does not cover the k-dis wavelength range'
                WRITE(*,*) ' in file ',TRIM(ADJUSTL(filename_kdis)) 
                WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
                WRITE(*,*) ''
                STOP
             endif
          endif
       ENDIF ! need to read solrad specified in artdeco_in.dat
    endif

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! allocate various arrays
    IF (do_rt) THEN

       ALLOCATE(rt_cputime(nlambda))
       rt_cputime(:) = undef_dp 

       ALLOCATE(rho_surf(nsza, nvza, nvaa, nlambda))
       ALLOCATE(SvR_levels(nsza, nvza, nvaa, nmat, nalt_atm, nlambda))
       ALLOCATE(SvR(nsza, nvza, nvaa, nmat, nlambda))
       ALLOCATE(SvR_sig(nsza, nvza, nvaa, nmat, nlambda))
       ALLOCATE(flux_out(nsza,3,nalt_atm,nlambda))
       rho_surf(:,:,:,:) = undef_dp
       SvR_levels(:,:,:,:,:,:) = undef_dp
       SvR(:,:,:,:,:)     = undef_dp
       SvR_sig(:,:,:,:,:) = undef_dp
       flux_out(:,:,:,:)  = undef_dp

       ! integrated values over kdis bands
       if (mode.eq.'kdis') then
          if (filter) then
             ALLOCATE(SvR_filter(nsza, nvza, nvaa, nmat, nfilter), & 
                  SvR_sig_filter(nsza, nvza, nvaa, nmat, nfilter), &
                  flux_filter(nsza, 3, nalt_atm, nfilter))
             SvR_filter(:,:,:,:,:)       = undef_dp
             SvR_sig_filter(:,:,:,:,:)   = undef_dp
             flux_filter(:,:,:,:)        = undef_dp
          else
             ALLOCATE(fluxint_out(nsza, 3, nalt_atm))
             fluxint_out(:,:,:) = undef_dp
          endif
       endif

       if (sinsca_corint .eqv. .true.) then
          ALLOCATE(SvR_levels_notms(nsza, nvza, nvaa, nmat, nalt_atm, nlambda))
          ALLOCATE(SvR_notms(nsza, nvza, nvaa, nmat, nlambda))
          ALLOCATE(SvR_sig_notms(nsza, nvza, nvaa, nmat, nlambda))
          SvR_levels_notms(:,:,:,:,:,:) = undef_dp
          SvR_notms(:,:,:,:,:)     = undef_dp
          SvR_sig_notms(:,:,:,:,:) = undef_dp
       endif
    ENDIF

  END SUBROUTINE READ_DATA

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

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

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

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

  SUBROUTINE WRITE_DATA

    USE MCOMMON
    USE MCONSTANTS
    USE MUTILITY, ONLY: TRIMCAT
    USE nr, only: indexx

    ! USE MSURFACE_BRDF
    USE MLAYERS, only : GET_MONO_LAYERS

    IMPLICIT NONE

    CHARACTER (LEN=max_len) :: filename_tmp

    INTEGER :: i, j, k, u, l
    INTEGER :: ifilter

    REAL(kind=dp) :: thetad, UMU0, UMU, polar

    character(len=8)  :: char8
    character(len=10) :: char10
    character(len=5)  :: char5

    TYPE (FICH) :: foutput        = FICH ('rad_'                , 11)
    TYPE (FICH) :: foutput_notms  = FICH ('rad_notms_'          , 12)
    TYPE (FICH) :: fopt           = FICH ('opt_'                , 13)
    TYPE (FICH) :: fopt_recomp    = FICH ('opt_recomp_'         , 14)
    TYPE (FICH) :: fbetal         = FICH ('Betal_'              , 15)
    TYPE (FICH) :: fsolrad        = FICH ('solrad_'             , 17)
    TYPE (FICH) :: fflux          = FICH ('flux_'               , 18)
    TYPE (FICH) :: ffluxkdis      = FICH ('flux_kdisband_'      , 19)
    ! TYPE (FICH) :: fopt_hg        = FICH ('opt_hg_'             , 20)
    TYPE (FICH) :: fheat          = FICH ('heat_'               , 21)
    TYPE (FICH) :: fkdisaitaui    = FICH ('kdis_aitaui_'         , 22)

    CHARACTER(LEN=*), PARAMETER  :: FMT_OUT    = '(200(2x,E18.9))'
    CHARACTER(LEN=*), PARAMETER  :: FMT_OUTII  = '(2x,I5,200(2x,E18.9))'
    CHARACTER(LEN=*), PARAMETER  :: FMT_OUTIII = '(2x,E18.9,2x,I5,200(2x,E18.9))'
    CHARACTER(LEN=*), PARAMETER  :: FMT_OUTIV  = '(2x,E18.9,2x,I5,2x,I5,200(2x,E18.9))'
    CHARACTER(LEN=*), PARAMETER  :: FMT_RADOUT = '(200(2x,F16.9))'
    CHARACTER(LEN=*), PARAMETER  :: FMT_FOUT   = FMT_RADOUT

    INTEGER :: nlamb
    INTEGER :: ialt
    REAL(kind=dp), ALLOCATABLE :: lamb(:)
    INTEGER, ALLOCATABLE :: ilamb(:)
    REAL(kind=dp), ALLOCATABLE :: cext(:)
    REAL(kind=dp), ALLOCATABLE :: ssa(:)
    REAL(kind=dp), ALLOCATABLE :: g(:)
    INTEGER, ALLOCATABLE :: nteta(:)
    REAL(kind=dp), ALLOCATABLE :: ptcle_u_phasemat_tmp(:,:)
    REAL(kind=dp), ALLOCATABLE :: ptcle_phasemat_tmp(:,:,:)

    REAL(kind=dp) :: rho_air
    REAL(kind=dp), ALLOCATABLE :: flux_net(:)
    REAL(kind=dp) :: net_rad_heat
    REAL(kind=dp) :: warm_rate

    integer :: iai
    integer :: ikdislamb

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

    IF (verbose) THEN
       WRITE (*,*) ''
       WRITE (*,*)'  (sub. write_data) write result files in directory ', TRIM(ADJUSTL(dir_root_out))
       WRITE (*,*) ''
    ENDIF

    IF (nptcle .gt. 0) THEN

       ptcle : DO i = 1, nptcle

          ! write optical properties if there are some new one
          if ((ptcle_flag_new_opt(i)).and.(ptcle_opt_model(i).ne.'baum')) then

             ! get teh total number of wavelengths to be in the new file
             if (ptcle_flag_store_opt(i)) then
                nlamb = ptcle_nlamb_phasemat_store(i) + ptcle_nlamb_phasemat_new(i)
             else
                nlamb = ptcle_nlamb_phasemat_new(i)
             endif
             ! create some arrays to contain new+stored optical properties
             allocate(lamb(nlamb), &
                  nteta(nlamb), &
                  Cext(nlamb), &
                  ssa(nlamb), &
                  g(nlamb), &
                  ptcle_u_phasemat_tmp(nang_max,nlamb), &
                  ptcle_phasemat_tmp(6, nang_max, nlamb), &
                  ilamb(nlamb))
             lamb(:)  = undef_dp
             nteta(:) = undef_i
             Cext(:)  = undef_dp
             ssa(:)   = undef_dp
             g(:)     = undef_dp
             ptcle_u_phasemat_tmp(:,:) = undef_dp
             ptcle_phasemat_tmp(:,:,:) = undef_dp
             ilamb(:) = undef_i
             ! implement new one first
             DO l = 1,  ptcle_nlamb_phasemat_new(i)
                lamb(l)  = ptcle_lamb_phasemat_new(i,l)
                nteta(l) = ptcle_nang_phasemat_new(i,l)
                Cext(l)  = ptcle_opt_new(i, 1, l)
                ssa(l)   = ptcle_opt_new(i, 2, l)
                g(l)     = ptcle_opt_new(i, 3, l)
                DO  k = 1, nteta(l)
                   ptcle_u_phasemat_tmp(k,l) = ptcle_u_phasemat_new(i, k, l)
                   ptcle_phasemat_tmp(:,k,l) = &
                        ptcle_phasemat_new(i,:, k, l)
                END DO
             END DO
             ! then old one
             IF (ptcle_flag_store_opt(i)) THEN
                DO l = 1,  ptcle_nlamb_phasemat_store(i)
                   lamb(l+ptcle_nlamb_phasemat_new(i))  = ptcle_lamb_phasemat_store(i,l)
                   nteta(l+ptcle_nlamb_phasemat_new(i)) = ptcle_nang_phasemat_store(i,l)
                   Cext(l+ptcle_nlamb_phasemat_new(i))  = ptcle_opt_store(i, 1, l)
                   ssa(l+ptcle_nlamb_phasemat_new(i))   = ptcle_opt_store(i, 2, l)
                   g(l+ptcle_nlamb_phasemat_new(i))     = ptcle_opt_store(i, 3, l)
                   DO  k = 1, nteta(l+ptcle_nlamb_phasemat_new(i))
                      ptcle_u_phasemat_tmp(k,l+ptcle_nlamb_phasemat_new(i)) = &
                           ptcle_u_phasemat_store(i, k, l)
                      ptcle_phasemat_tmp(:,k,l+ptcle_nlamb_phasemat_new(i)) = &
                           ptcle_phasemat_store(i, :, k, l)
                   END DO
                END DO
             ENDIF
             ! get the index for sorted wavelengths
             CALL INDEXX(lamb, ilamb)
             ! finally write the new optical properties file
             filename_tmp = TRIMCAT(artdeco_path, dir_opt)
             filename_tmp = TRIMCAT(filename_tmp, fopt%name)
             filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             OPEN (UNIT = fopt%unit, FILE = filename_tmp, STATUS = 'REPLACE')
             WRITE (UNIT=fopt%unit, FMT='(A01)')    '#'
             WRITE (UNIT=fopt%unit, FMT='(A42)')    '# Optical properties to be used by ARTDECO'
             WRITE (UNIT=fopt%unit, FMT='(A01)')    '#'
             WRITE (UNIT=fopt%unit, FMT='(A42)')    '# Used model to obtain that properties is:'
             WRITE (UNIT=fopt%unit, FMT=*)  TRIM(ADJUSTL(ptcle_opt_model(i)))
             WRITE (UNIT=fopt%unit, FMT='(A19)')    '# Used material is:'
             WRITE (UNIT=fopt%unit, FMT=*)  TRIM(ADJUSTL(ptcle_material(i)))
             WRITE (UNIT=fopt%unit, FMT='(A35)')    '# Number of phase matrix elements :'
             WRITE (UNIT=fopt%unit, FMT=*) ptcle_nelem(i)
             WRITE (UNIT=fopt%unit, FMT='(A23)')    '# number of wavelengths'
             WRITE (UNIT=fopt%unit, FMT=*) nlamb
             WRITE (UNIT=fopt%unit, FMT='(A)')    &
                  '#  wlambda(microns)   nteta   Cext (microns^2)       SSA              g '
             DO l = 1, nlamb
                WRITE (UNIT=fopt%unit, FMT=FMT_OUTIII) lamb(ilamb(l)), nteta(ilamb(l)), &
                     Cext(ilamb(l)), ssa(ilamb(l)), g(ilamb(l))
             enddo
             WRITE (UNIT=fopt%unit, FMT='(A14)') '# Phase matrix'
             if (ptcle_nelem(i).eq.6) then
                DO l = 1, nlamb
                   WRITE (UNIT=fopt%unit, FMT='(A,2x,F10.4)') '# lambda =', lamb(ilamb(l))
                   WRITE (UNIT=fopt%unit, FMT='(A121)')   '#         u                 F11   &
                        &             F22              F33               F44               F21               F34'
                   do k = 1, nteta(ilamb(l))
                      WRITE (UNIT=fopt%unit, FMT=FMT_OUT)  ptcle_u_phasemat_tmp(k,ilamb(l)), &
                           ptcle_phasemat_tmp(1,k,ilamb(l)), ptcle_phasemat_tmp(2,k,ilamb(l)), &
                           ptcle_phasemat_tmp(3,k,ilamb(l)), ptcle_phasemat_tmp(4,k,ilamb(l)), &
                           ptcle_phasemat_tmp(5,k,ilamb(l)), ptcle_phasemat_tmp(6,k,ilamb(l))
                   enddo
                ENDDO
             else
                DO l = 1, nlamb
                   WRITE (UNIT=fopt%unit, FMT='(A,2x,F10.4)') '# lambda =', lamb(ilamb(l))
                   WRITE (UNIT=fopt%unit, FMT='(A84)')   '#         u                 F11   &
                        &           F44               F21               F34'
                   do k = 1, nteta(ilamb(l))
                      WRITE (UNIT=fopt%unit, FMT=FMT_OUT)  ptcle_u_phasemat_tmp(k,ilamb(l)), &
                           ptcle_phasemat_tmp(1,k,ilamb(l)), ptcle_phasemat_tmp(4,k,ilamb(l)), &
                           ptcle_phasemat_tmp(5,k,ilamb(l)), ptcle_phasemat_tmp(6,k,ilamb(l))
                   enddo
                ENDDO
             endif
             CLOSE(UNIT = fopt%unit)

             deallocate(lamb, &
                  nteta, &
                  Cext, &
                  ssa, &
                  g, &
                  ptcle_u_phasemat_tmp, &
                  ptcle_phasemat_tmp, &
                  ilamb)

          ENDIF! new_opt

          ! WRITE COMPUTED BETAL & RECOMP. PHASE MAT
          IF (opt_only .eqv. .FALSE.) THEN

             if (nmaxbetal.gt.0) then

                if (print_betal) then
                   !========
                   ! Betal
                   filename_tmp = TRIMCAT(dir_root_out, fbetal%name)
                   if (ptcle_hg(i)) then
                      filename_tmp = TRIMCAT(filename_tmp, 'hg_' )
                   endif
                   filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
                   filename_tmp = TRIMCAT(filename_tmp, '_')
                   filename_tmp = TRIMCAT(filename_tmp, TRIM(ADJUSTL(trunc_method)) )
                   filename_tmp = TRIMCAT(filename_tmp, '.dat')
                   OPEN (UNIT = fbetal%unit, FILE = filename_tmp, STATUS = 'unknown')
                   WRITE (UNIT=fbetal%unit, FMT='(A01)')     '#'
                   WRITE (UNIT=fbetal%unit, FMT='(A18)')     '# ARTDECO : Beta_l'
                   WRITE (UNIT=fbetal%unit, FMT='(A01)')     '#'
                   WRITE (UNIT=fbetal%unit, FMT='(A32)')     '# Used model for truncation was:'
                   WRITE (UNIT=fbetal%unit, FMT=*)  TRIM(ADJUSTL(trunc_method)) 
                   WRITE (UNIT=fbetal%unit, FMT='(A23)')    '# number of wavelengths'
                   WRITE (UNIT=fbetal%unit, FMT=*) nlambda
                   WRITE (UNIT=fbetal%unit, FMT='(A)')    '# wlambda(microns)     nBetal     Trunc. Coeff.    &
                        & Cext (microns^2)       SSA              g '
                   DO l = 1, nlambda
                      WRITE (UNIT=fbetal%unit, FMT=FMT_OUTIII) wlambda(l), ptcle_nbetal(i,l), ptcle_trunccoeff(i,l), &
                           ptcle_opt(i, 1, l), ptcle_opt(i, 2, l), ptcle_opt(i, 3, l)
                   END DO
                   WRITE (UNIT=fbetal%unit, FMT='(A7)') '# Betal'
                   DO l = 1, nlambda
                      WRITE (UNIT=fbetal%unit, FMT='(A,2x,F10.4)') '# lambda =', wlambda(l)
                      WRITE (UNIT=fbetal%unit, FMT='(A107)')    '#     l     alpha1            alpha2            alpha3          &
                           &  alpha4             beta1            beta2'
                      do j = 0, ptcle_nbetal(i,l)
                         WRITE (UNIT=fbetal%unit, FMT=FMT_OUTII) j, ptcle_betal(i,1,j,l), ptcle_betal(i,2,j,l), &
                              ptcle_betal(i,3,j,l), ptcle_betal(i,4,j,l), ptcle_betal(i,5,j,l), ptcle_betal(i,6,j,l)
                      enddo
                   ENDDO
                   CLOSE(UNIT = fbetal%unit)
                endif

                if (print_recomp) then
                   ! write recomposed phase function (or matrix)
                   filename_tmp = TRIMCAT(dir_root_out, fopt_recomp%name)
                   if (ptcle_hg(i)) then
                      filename_tmp = TRIMCAT(filename_tmp, 'hg_' )
                   endif
                   filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
                   filename_tmp = TRIMCAT(filename_tmp, '_')
                   filename_tmp = TRIMCAT(filename_tmp, TRIM(ADJUSTL(trunc_method)) )
                   filename_tmp = TRIMCAT(filename_tmp, '.dat')
                   OPEN (UNIT = fopt_recomp%unit, FILE = filename_tmp, STATUS = 'unknown')
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A01)')    '#'
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A45)')    '# Recomposed optical properties out of ARTDECO'
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A01)')    '#'
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A32)')    '# Used model for truncation was:'
                   WRITE (UNIT=fopt_recomp%unit, FMT=*)  TRIM(ADJUSTL(trunc_method)) 
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A23)')    '# number of wavelengths'
                   WRITE (UNIT=fopt_recomp%unit, FMT=*) nlambda
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A)')     '#      lambda         nBl   &
                        &nangle   Trunc. Coeff.     sumF11    Cext (microns^2)       SSA              g  '
                   DO l = 1, nlambda
                      WRITE (UNIT=fopt_recomp%unit, FMT=FMT_OUTIV) wlambda(l), ptcle_nbetal(i,l), ptcle_nang_phasemat(i,l), &
                           ptcle_trunccoeff(i,l),ptcle_trunc_normphasemat(i,l), &
                           ptcle_opt(i, 1, l), ptcle_opt(i, 2, l), ptcle_opt(i, 3, l)
                   END DO
                   WRITE (UNIT=fopt_recomp%unit, FMT='(A14)') '# Phase matrix'
                   DO l = 1, nlambda
                      WRITE (UNIT=fopt_recomp%unit, FMT='(A,2x,F10.4)') '# lambda =', wlambda(l)
                      WRITE (UNIT=fopt_recomp%unit, FMT='(A121)')   '#         u                 F11   &
                           &             F22              F33               F44               F21               F34'
                      do j=1,ptcle_nang_phasemat(i,l)
                         WRITE (UNIT=fopt_recomp%unit, FMT=FMT_OUT) ptcle_u_phasemat(i, j, l), &
                              ptcle_trunc_phasemat(i, 1, j, l),ptcle_trunc_phasemat(i, 2, j, l),ptcle_trunc_phasemat(i, 3, j, l), &
                              ptcle_trunc_phasemat(i, 4, j, l),ptcle_trunc_phasemat(i, 5, j, l),ptcle_trunc_phasemat(i, 6, j, l)
                      enddo
                   ENDDO
                   CLOSE(UNIT = fopt_recomp%unit)
                endif

             endif

             ! if (ptcle_hg(i)) then
             !    ! write the H-G phase function (or matrix)
             !    filename_tmp = TRIMCAT(dir_root_out, fopt_hg%name)
             !    filename_tmp = TRIMCAT(filename_tmp, ptcle_type(i))
             !    filename_tmp = TRIMCAT(filename_tmp, '.dat')
             !    OPEN (UNIT = fopt_hg%unit, FILE = filename_tmp, STATUS = 'unknown')
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A)')      '#'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A)')      '# H-G approx. phase function obtained from the asymetry &
             !         &factor of the particle'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A01)')    '#'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A42)')    '# Used model to obtain that properties is:'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A)')      '  Henyey-Greenstein '
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A19)')    '# Used material is:'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A)')      '  unknown'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A35)')    '# Number of phase matrix elements :'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A)')      '   1'
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A23)')    '# number of wavelengths'
             !    WRITE (UNIT=fopt_hg%unit, FMT=*) nlambda
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A)')    &
             !         '#  wlambda(microns)   nteta   Cext (microns^2)       SSA              g '
             !    DO l = 1, nlambda
             !       WRITE (UNIT=fopt_hg%unit, FMT=FMT_OUTIII) wlambda(l), ptcle_nang_phasemat(i,l), &
             !            ptcle_opt(i,1,l),ptcle_opt(i,2,l),ptcle_opt(i,3,l)
             !    enddo
             !    WRITE (UNIT=fopt_hg%unit, FMT='(A14)') '# Phase matrix'
             !    DO l = 1, nlambda
             !       WRITE (UNIT=fopt_hg%unit, FMT='(A,2x,F10.4)') '# lambda =', wlambda(l)
             !       WRITE (UNIT=fopt_hg%unit, FMT='(A)')   '#         u                 F11   '
             !       do j=1,ptcle_nang_phasemat(i,l)
             !          WRITE (UNIT=fopt_hg%unit, FMT=FMT_OUT) ptcle_u_phasemat(i, j, l),ptcle_phasemat(i, 1, j, l)
             !       enddo
             !    ENDDO
             !    CLOSE(UNIT = fopt_hg%unit)
             ! endif

          END IF

       END DO ptcle

    ENDIF !(nptcle .gt 0)

    !=======================================================
    !                   SAVE RT PRODUCES
    ! radiances
    IF (do_rt) THEN

       if (mode.ne.'mono') then
          ! save solar irradiance used for computation
          filename_tmp = TRIMCAT(dir_root_out, fsolrad%name)
          filename_tmp = TRIMCAT(filename_tmp, out_root)
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, rt_model)
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          OPEN (UNIT = fsolrad%unit, FILE = filename_tmp, STATUS = 'NEW')
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# Solar irradiance on TOA (in W m-2 perpendicular to the light propagation)'
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# varsol is a correction factor applied to'
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# account for sun-earth distance evolution during a year'
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# The incoming flux a TOA is then F0 x varsol'
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# Number of wavelengths :'
          WRITE (UNIT=fsolrad%unit, FMT=*) nlambda
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# Number of solar configurations (sza / date / geographical coordinates)'
          WRITE (UNIT=fsolrad%unit, FMT=*) nsza 
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '################## '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '#  lambda (microns)    F0 (W m-2) '
          do i = 1, nlambda
             WRITE (UNIT=fsolrad%unit, FMT=FMT_RADOUT) wlambda(i), F0(i)
          end do
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '################## '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '#    sza         varsol       lon.   &
               &     lat.    day of year    hour (TU)  '
          if (get_possol) then
             DO j = 1, nsza
                WRITE (UNIT=fsolrad%unit, FMT='(4(2x,F10.6), 4x, I3, 8x, F10.6)')  sza(j), varsol_fact(j), &
                     lon, lat, day(j), h_tu(j)
             ENDDO
          else
             DO j = 1, nsza
                WRITE (UNIT=fsolrad%unit, FMT='(4(2x,F10.6), 4x, I3, 8x, F10.6)')   sza(j),  varsol_fact(j), &
                     -1.0, -1.0, -1, -1.0
             ENDDO
          end if
          CLOSE(fsolrad%unit)
       endif

       !        if (mode.eq.'mono') then
       !    WRITE (UNIT=foutput%unit, FMT='(A)')      '# F0 '
       ! else
       !    WRITE (UNIT=foutput%unit, FMT='(A)')      '# F0 (W m-2) '
       ! endif
       ! WRITE (UNIT=foutput%unit,FMT=FMT_RADOUT)        F0(l)

       if (mode.eq.'kdis') then

          ! if radiative quantities were integrated through ISRF :
          if (filter) then

             do ifilter = 1, nfilter

                ! filter integrated flux 
                filename_tmp = TRIMCAT (dir_root_out, fflux%name)
                filename_tmp = TRIMCAT(filename_tmp, out_root)
                filename_tmp = TRIMCAT(filename_tmp, '_')
                filename_tmp = TRIMCAT(filename_tmp, filter_name(ifilter))
                filename_tmp = TRIMCAT(filename_tmp, '_')
                filename_tmp = TRIMCAT(filename_tmp, rt_model)
                filename_tmp = TRIMCAT(filename_tmp, '.dat')
                OPEN (UNIT = fflux%unit, FILE = filename_tmp, STATUS = 'unknown')
                WRITE (UNIT=fflux%unit,FMT='(A)') '#  '
                WRITE (UNIT=fflux%unit,FMT='(A)') '#  ARTDECO   flux output '
                WRITE (UNIT=fflux%unit,FMT='(A)') '#  '
                WRITE (UNIT=fflux%unit,FMT='(A)') '# number of solar zenith angles '
                WRITE (UNIT=fflux%unit,FMT='(I4)')  nsza
                WRITE (UNIT=fflux%unit,FMT='(A)') '# number of sampled altitudes '
                WRITE (UNIT=fflux%unit,FMT='(I4)')  nalt_atm
                DO i = 1, nsza          
                   WRITE (UNIT=fflux%unit,FMT='(A)') '# Solar zenithal angle (deg) '
                   WRITE (UNIT=fflux%unit,FMT='(F6.3)')  sza(i)
                   WRITE (UNIT=fflux%unit,FMT='(A)') '#     ALTITUDE           PRESSION           FLUX DW   &
                        &       FLUX UP           FLUX DIR'  
                   WRITE (UNIT=fflux%unit,FMT='(A)') '#       (KM)               (hPa)             (W.M-2)  &
                        &        (W.M-2)          (W.M-2) '   
                   DO ialt =  1, nalt_atm
                      WRITE (UNIT=fflux%unit, FMT=FMT_FOUT) alt_atm(ialt), p_atm(ialt), (flux_filter(i, j, ialt, ifilter), j=1,3)
                   enddo
                enddo
                CLOSE(UNIT=fflux%unit)

                if (flux_only.eqv..FALSE.) then
                   ! filter integrated radiances
                   filename_tmp = TRIMCAT (dir_root_out, foutput%name)
                   filename_tmp = TRIMCAT(filename_tmp, out_root)
                   filename_tmp = TRIMCAT(filename_tmp, '_')
                   filename_tmp = TRIMCAT(filename_tmp, filter_name(ifilter))
                   filename_tmp = TRIMCAT(filename_tmp, '_')
                   filename_tmp = TRIMCAT(filename_tmp, rt_model)
                   filename_tmp = TRIMCAT(filename_tmp, '.dat')
                   OPEN (UNIT = foutput%unit, FILE = filename_tmp, STATUS = 'unknown')
                   WRITE (UNIT=foutput%unit, FMT='(A23)')      '# ARTDECO : output file'
                   WRITE (UNIT=foutput%unit, FMT='(A)')        '# Radiances at the top of the atmosphere (W/m2/sr)'
                   WRITE (UNIT=foutput%unit, FMT='(A04)')      '#   '
                   WRITE (UNIT=foutput%unit, FMT='(A32)')      '# Result given for each geometry :'
                   WRITE (UNIT=foutput%unit, FMT='(A38)')      '# sza : solar zenith angle(s) (degree)'
                   WRITE (UNIT=foutput%unit, FMT='(A38)')      '# vza : view  zenith angle(s) (degree)'
                   WRITE (UNIT=foutput%unit, FMT='(A40)')      '# vaa : view azimuthal angle(s) (degree)'
                   WRITE (UNIT=foutput%unit, FMT='(A1)')       '#'
                   WRITE (UNIT=foutput%unit, FMT='(A28)')      '# total number of geometries'
                   WRITE (UNIT=foutput%unit, FMT=*) nsza * nvza * nvaa
                   WRITE (UNIT=foutput%unit, FMT='(A18)')      '# nsza, nvza, nvaa'
                   WRITE (UNIT=foutput%unit, FMT=*) nsza, nvza, nvaa
                   WRITE (UNIT=foutput%unit, FMT='(A41)')      '# number of elements of the Stokes vector'
                   WRITE (UNIT=foutput%unit, FMT='(A58)')      '# 4 = full polarization, 3 = 3x3 approximation, 1 = scalar'
                   WRITE (UNIT=foutput%unit, FMT=*) nmat
                   WRITE (UNIT=foutput%unit, FMT='(A)')      '# CPU time for RT in sec(s) '
                   WRITE (UNIT=foutput%unit, FMT=FMT_OUT)      rt_cputime_filter(ifilter)
                   SELECT CASE (nmat)
                   CASE(1) 
                      WRITE (UNIT=foutput%unit, FMT='(A101)')      '#         sza               vza              vaa        &
                           &      theta              I                 dI'
                   CASE(3) 
                      WRITE (UNIT=foutput%unit, FMT='(A197)')    '#         sza               vza              vaa         &
                           &      theta              I                 Q                U                 dI               &
                           & dQ               &
                           & dU              Polar. (%)'
                   CASE(4)  
                      WRITE (UNIT=foutput%unit, FMT='(A233)')    '#         sza               vza              vaa         &
                           &      theta              I                 Q                U                  V               &
                           & dI               &
                           & dQ                dU                dV              Polar. (%)'
                   END SELECT
                   DO i = 1, nsza
                      UMU0 = COS(sza(i)*deg2rad)
                      DO j = nvza, 1, -1
                         UMU = COS(vza(j)*deg2rad)
                         DO k = 1, nvaa
                            thetad = ACOS(-UMU0 * UMU + SQRT( ( 1.0_dp-UMU0**2.0_dp ) * &
                                 ( 1.0_dp - UMU**2.0_dp ) ) * COS( vaa(k)*deg2rad )) / deg2rad
                            SELECT CASE (nmat)
                            CASE(1) 
                               WRITE (UNIT=foutput%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, &
                                    SvR_filter(i, j, k, 1, ifilter), &
                                    SvR_sig_filter(i, j, k, 1, ifilter)
                            CASE(3) 
                               polar = SQRT(SvR_filter(i, j, k, 2, ifilter)**2.0_dp + SvR_filter(i, j, k, 3, ifilter)**2.0_dp) / &
                                    SvR_filter(i, j, k, 1, ifilter) * 100.0_dp
                               if ((polar .gt. 100.0_dp) .or. (polar .lt. 0.0_dp)) polar = undef_dp
                               WRITE (UNIT=foutput%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, &
                                    (SvR_filter(i, j, k, u, ifilter), u=1,nmat), &
                                    (SvR_sig_filter(i, j, k, u, ifilter), u=1,nmat), polar
                            CASE(4) 
                               polar = SQRT(SvR_filter(i, j, k, 2, ifilter)**2.0_dp + SvR_filter(i, j, k, 3, ifilter)**2.0_dp + &
                                    SvR_filter(i, j, k, 4, ifilter)**2.0_dp) / &
                                    & SvR_filter(i, j, k, 1, ifilter) * 100.0_dp
                               if ((polar .gt. 100.0_dp) .or. (polar .lt. 0.0_dp)) polar = undef_dp
                               WRITE (UNIT=foutput%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, &
                                    (SvR_filter(i, j, k, u, ifilter), u=1,nmat), &
                                    (SvR_sig_filter(i, j, k, u, ifilter), u=1,nmat), polar 
                            END SELECT
                         ENDDO
                      ENDDO
                   ENDDO
                   CLOSE(UNIT = foutput%unit)
                endif! flux_only is FALSE
             enddo ! loop on filters

          else

             ! integrated flux
             filename_tmp = TRIMCAT (dir_root_out, fflux%name)
             filename_tmp = TRIMCAT(filename_tmp, out_root)
             filename_tmp = TRIMCAT(filename_tmp, '_')
             filename_tmp = TRIMCAT(filename_tmp, rt_model)
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             OPEN (UNIT = fflux%unit, FILE = filename_tmp, STATUS = 'unknown')
             WRITE (UNIT=fflux%unit,FMT='(A)') '#  '
             WRITE (UNIT=fflux%unit,FMT='(A)') '#  ARTDECO   flux output '
             WRITE (UNIT=fflux%unit,FMT='(A)') '#  '
             WRITE (UNIT=fflux%unit,FMT='(A)') '# number of solar zenith angles '
             WRITE (UNIT=fflux%unit,FMT='(I4)')  nsza
             WRITE (UNIT=fflux%unit,FMT='(A)') '# number of sampled altitudes '
             WRITE (UNIT=fflux%unit,FMT='(I4)')  nalt_atm
             DO i = 1, nsza          
                WRITE (UNIT=fflux%unit,FMT='(A)') '# Solar zenithal angle (deg) '
                WRITE (UNIT=fflux%unit,FMT='(F6.3)')  sza(i)
                WRITE (UNIT=fflux%unit,FMT='(A)') '#     ALTITUDE           PRESSION           FLUX DW   &
                     &       FLUX UP           FLUX DIR'  
                WRITE (UNIT=fflux%unit,FMT='(A)') '#       (KM)               (hPa)             (W.M-2)  &
                     &        (W.M-2)          (W.M-2) '   
                DO ialt = 1, nalt_atm
                   WRITE (UNIT=fflux%unit, FMT=FMT_FOUT) alt_atm(ialt), p_atm(ialt), (fluxint_out(i,j,ialt), j=1,3)
                enddo
             enddo
             CLOSE(UNIT=fflux%unit)

             ! Compute and write net radiative heating and  warnimg rate
             ALLOCATE(flux_net(nalt_atm))

             filename_tmp = TRIMCAT (dir_root_out, fheat%name)
             filename_tmp = TRIMCAT(filename_tmp, out_root)
             filename_tmp = TRIMCAT(filename_tmp, '_')
             filename_tmp = TRIMCAT(filename_tmp, rt_model)
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             OPEN (UNIT = fheat%unit, FILE = filename_tmp, STATUS = 'unknown')
             WRITE (UNIT=fheat%unit,FMT='(A)') '#  '
             WRITE (UNIT=fheat%unit,FMT='(A)') '#  ARTDECO heating rate '
             WRITE (UNIT=fheat%unit,FMT='(A)') '#  '
             WRITE (UNIT=fheat%unit,FMT='(A)') '# number of solar zenith angles '
             WRITE (UNIT=fheat%unit,FMT='(I4)')  nsza
             WRITE (UNIT=fheat%unit,FMT='(A)') '# number of sampled layers '
             WRITE (UNIT=fheat%unit,FMT='(I4)')  nalt_atm - 1
             DO i = 1, nsza 
                flux_net(:) = fluxint_out(i,2,:) - fluxint_out(i,1,:) ! FLUP - FLDW
                WRITE (UNIT=fheat%unit,FMT='(A)') '# Solar zenithal angle (deg) '
                WRITE (UNIT=fheat%unit,FMT='(F6.3)')  sza(i)
                WRITE (UNIT=fheat%unit,FMT='(A)') '#    LOWER BOUND.       UPPER BOUND.    NET RAD. HEAT. RATE   WARMING RATE'    
                WRITE (UNIT=fheat%unit,FMT='(A)') '#        (KM)               (KM)           (W M-2 M-1)       (K / DAY)    '
                DO ialt = 2, nalt_atm
                   rho_air      = (u_air_atm(ialt-1) + u_air_atm(ialt)) / 2.0_dp * mm_air * amu * 1D6 ! g m-3 
                   net_rad_heat = -(flux_net(ialt-1)-flux_net(ialt)) / (alt_atm(ialt-1)-alt_atm(ialt)) / 1.0D3 ! W m-2 m-1
                   warm_rate    = 86400_dp * net_rad_heat / (rho_air*cp_air)
                   WRITE (UNIT=fheat%unit, FMT=FMT_FOUT) alt_atm(ialt), alt_atm(ialt-1), net_rad_heat, warm_rate
                enddo
             enddo
             CLOSE(UNIT=fheat%unit)

             DEALLOCATE(flux_net)

          endif ! filter

          ! each kdis band flux
          filename_tmp = TRIMCAT (dir_root_out, ffluxkdis%name)
          filename_tmp = TRIMCAT(filename_tmp, out_root)
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, rt_model)
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          OPEN (UNIT = ffluxkdis%unit, FILE = filename_tmp, STATUS = 'unknown')
          WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '#  '
          WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '#  ARTDECO flux output for each k-dis band '
          WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '#  '
          WRITE (UNIT=ffluxkdis%unit, FMT='(A)') '# number of wavelengths'
          WRITE (UNIT=ffluxkdis%unit, FMT='(I4)')  nlambda
          WRITE (UNIT=ffluxkdis%unit, FMT='(A)') '# number of solar zenith angles'
          WRITE (UNIT=ffluxkdis%unit, FMT='(I4)')  nsza
          WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '# number of sampled altitudes '
          WRITE (UNIT=ffluxkdis%unit,FMT='(I4)')  nalt_atm
          WRITE (UNIT=ffluxkdis%unit, FMT='(A)')      '##################'
          DO l = 1, nlambda
             WRITE (UNIT=ffluxkdis%unit, FMT='(A23)')      '# Wavelength (microns) '
             WRITE (UNIT=ffluxkdis%unit,FMT=FMT_RADOUT)         wlambda(l)
             DO i = 1, nsza          
                WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '# Solar zenithal angle (deg) '
                WRITE (UNIT=ffluxkdis%unit,FMT='(F6.3)')  sza(i)
                WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '#     ALTITUDE           PRESSION           FLUX DW   &
                     &       FLUX UP           FLUX DIR'  
                WRITE (UNIT=ffluxkdis%unit,FMT='(A)') '#       (KM)               (hPa)             (W.M-2)  &
                     &        (W.M-2)          (W.M-2) '   
                DO ialt = 1, nalt_atm
                   WRITE (UNIT=ffluxkdis%unit, FMT=FMT_FOUT) alt_atm(ialt), p_atm(ialt), (flux_out(i,j,ialt,l), j=1,3)
                enddo
             enddo
          end DO
          CLOSE(UNIT=ffluxkdis%unit)

       endif !mode eq kdis

       if (flux_only.eqv..FALSE.) then
          ! radiance for each wavelength (or each kdis band)
          filename_tmp = TRIMCAT (dir_root_out, foutput%name)
          filename_tmp = TRIMCAT(filename_tmp, out_root)
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, rt_model)
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          OPEN (UNIT = foutput%unit, FILE = filename_tmp, STATUS = 'unknown')
          WRITE (UNIT=foutput%unit, FMT='(A23)')      '# ARTDECO : output file'
          if (mode.eq.'mono') then
             WRITE (UNIT=foutput%unit, FMT='(A)')      '# Radiances at the top of the atmosphere'
          else
             WRITE (UNIT=foutput%unit, FMT='(A)')      '# Radiances at the top of the atmosphere (W/m2/sr)'
          endif
          WRITE (UNIT=foutput%unit, FMT='(A04)')      '#   '
          WRITE (UNIT=foutput%unit, FMT='(A32)')      '# Result given for each geometry :'
          WRITE (UNIT=foutput%unit, FMT='(A38)')      '# sza : solar zenith angle(s) (degree)'
          WRITE (UNIT=foutput%unit, FMT='(A38)')      '# vza : view  zenith angle(s) (degree)'
          WRITE (UNIT=foutput%unit, FMT='(A40)')      '# vaa : view azimuthal angle(s) (degree)'
          WRITE (UNIT=foutput%unit, FMT='(A1)')       '#'
          WRITE (UNIT=foutput%unit, FMT='(A28)')      '# total number of geometries'
          WRITE (UNIT=foutput%unit, FMT=*) nsza * nvza * nvaa
          WRITE (UNIT=foutput%unit, FMT='(A18)')      '# nsza, nvza, nvaa'
          WRITE (UNIT=foutput%unit, FMT=*) nsza, nvza, nvaa
          WRITE (UNIT=foutput%unit, FMT='(A41)')      '# number of elements of the Stokes vector'
          WRITE (UNIT=foutput%unit, FMT='(A58)')      '# 4 = full polarization, 3 = 3x3 approximation, 1 = scalar'
          WRITE (UNIT=foutput%unit, FMT=*) nmat
          WRITE (UNIT=foutput%unit, FMT='(A23)')      '# number of wavelengths'
          WRITE (UNIT=foutput%unit, FMT=*) nlambda
          WRITE (UNIT=foutput%unit, FMT='(A18)')      '##################'
          DO l = 1, nlambda
             WRITE (UNIT=foutput%unit, FMT='(A23)')      '# Wavelength (microns) '
             WRITE (UNIT=foutput%unit,FMT=FMT_RADOUT)         wlambda(l)
             WRITE (UNIT=foutput%unit, FMT='(A)')      '# CPU time for RT in sec(s) '
             WRITE (UNIT=foutput%unit, FMT=FMT_OUT)      rt_cputime(l)
             SELECT CASE (nmat)
             CASE(1) 
                WRITE (UNIT=foutput%unit, FMT='(A101)')      '#         sza               vza              vaa        &
                     &      theta              I                 dI'
             CASE(3) 
                WRITE (UNIT=foutput%unit, FMT='(A197)')    '#         sza               vza              vaa         &
                     &      theta              I                 Q                U                 dI               &
                     & dQ               &
                     & dU              Polar. (%)'
             CASE(4)  
                WRITE (UNIT=foutput%unit, FMT='(A233)')    '#         sza               vza              vaa         &
                     &      theta              I                 Q                U                  V               &
                     & dI               &
                     & dQ                dU                dV              Polar. (%)'
             END SELECT
             DO i = 1, nsza
                UMU0 = COS(sza(i)*deg2rad)
                DO j = nvza, 1, -1
                   UMU = COS(vza(j)*deg2rad)
                   DO k = 1, nvaa
                      thetad = ACOS(-UMU0 * UMU + SQRT( ( 1.0_dp-UMU0**2.0_dp ) * &
                           ( 1.0_dp - UMU**2.0_dp ) ) * COS( vaa(k)*deg2rad )) / deg2rad
                      SELECT CASE (nmat)
                      CASE(1) 
                         
                         if (verbose.eqv..true.) then
                            if (UMU.gt.0.0) then
                               write(*,*) "sza, vza, vaa, lambda, I =", sza(i), vza(j), vaa(k),  wlambda(l), SvR(i, j, k, 1, l)
                            else
                               write(*,*) "sza, vza, vaa, lambda, I =", &
                                    sza(i), vza(j), vaa(k),  wlambda(l), SvR_levels(i, j, k, 1, nalt_atm, l) 
                            end if
                         end if

                         WRITE (UNIT=foutput%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, SvR(i, j, k, 1, l), &
                              SvR_sig(i, j, k, 1, l)
                      CASE(3)
                         
                         if (verbose.eqv..true.) then
                            write(*,*) "sza, vza, vaa, lambda, I , Q, U =", sza(i), vza(j), vaa(k),  &
                                 wlambda(l), SvR(i, j, k, 1, l), SvR(i, j, k, 2, l), SvR(i, j, k, 3, l)
                         end if
                         polar = SQRT(SvR(i, j, k, 2, l)**2.0_dp + SvR(i, j, k, 3, l)**2.0_dp) / SvR(i, j, k, 1, l) * 100.0_dp
                         if ((polar .gt. 100.0_dp) .or. (polar .lt. 0.0_dp)) polar = undef_dp
                         WRITE (UNIT=foutput%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, (SvR(i, j, k, u, l), u=1,nmat), &
                              (SvR_sig(i, j, k, u, l), u=1,nmat), polar
                      CASE(4) 
                         polar = SQRT(SvR(i, j, k, 2, l)**2.0_dp + SvR(i, j, k, 3, l)**2.0_dp + SvR(i, j, k, 4, l)**2.0_dp) / &
                              & SvR(i, j, k, 1, l) * 100.0_dp
                         if ((polar .gt. 100.0_dp) .or. (polar .lt. 0.0_dp)) polar = undef_dp
                         WRITE (UNIT=foutput%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, (SvR(i, j, k, u, l), u=1,nmat), &
                              (SvR_sig(i, j, k, u, l), u=1,nmat), polar 
                      END SELECT
                   ENDDO
                ENDDO
             ENDDO
          ENDDO ! lambda
          CLOSE(UNIT = foutput%unit)
          ! radiance without TMS correction
          if (sinsca_corint) then
             filename_tmp = TRIMCAT (dir_root_out, foutput_notms%name)
             filename_tmp = TRIMCAT(filename_tmp, out_root)
             filename_tmp = TRIMCAT(filename_tmp, '_')
             filename_tmp = TRIMCAT(filename_tmp, rt_model)
             filename_tmp = TRIMCAT(filename_tmp, '.dat')
             OPEN (UNIT = foutput_notms%unit, FILE = filename_tmp, STATUS = 'unknown')
             WRITE (UNIT=foutput_notms%unit, FMT='(A23)')      '# ARTDECO : output file'
             if (mode.eq.'mono') then
                WRITE (UNIT=foutput_notms%unit, FMT='(A)')      '# Radiances at the top of the atmosphere'
             else
                WRITE (UNIT=foutput_notms%unit, FMT='(A)')      '# Radiances at the top of the atmosphere (W/m2/sr)'
             endif
             WRITE (UNIT=foutput_notms%unit, FMT='(A1)')       '#'
             WRITE (UNIT=foutput_notms%unit, FMT='(A58)')      '# The single scattering intensity correction was not apply'
             WRITE (UNIT=foutput_notms%unit, FMT='(A1)')       '#'
             WRITE (UNIT=foutput_notms%unit, FMT='(A32)')      '# Result given for each geometry :'
             WRITE (UNIT=foutput_notms%unit, FMT='(A38)')      '# sza : solar zenith angle(s) (degree)'
             WRITE (UNIT=foutput_notms%unit, FMT='(A38)')      '# vza : view  zenith angle(s) (degree)'
             WRITE (UNIT=foutput_notms%unit, FMT='(A40)')      '# vaa : view azimuthal angle(s) (degree)'
             WRITE (UNIT=foutput_notms%unit, FMT='(A1)')       '#'
             WRITE (UNIT=foutput_notms%unit, FMT='(A28)')      '# total number of geometries'
             WRITE (UNIT=foutput_notms%unit, FMT=*) nsza * nvza * nvaa
             WRITE (UNIT=foutput_notms%unit, FMT='(A18)')      '# nsza, nvza, nvaa'
             WRITE (UNIT=foutput_notms%unit, FMT=*) nsza, nvza, nvaa
             WRITE (UNIT=foutput_notms%unit, FMT='(A41)')      '# number of elements of the Stokes vector'
             WRITE (UNIT=foutput_notms%unit, FMT='(A58)')      '# 4 = full polarization, 3 = 3x3 approximation, 1 = scalar'
             WRITE (UNIT=foutput_notms%unit, FMT=*) nmat
             WRITE (UNIT=foutput_notms%unit, FMT='(A23)')      '# number of wavelengths'
             WRITE (UNIT=foutput_notms%unit, FMT=*) nlambda
             WRITE (UNIT=foutput_notms%unit, FMT='(A18)')      '##################'
             DO l = 1, nlambda
                WRITE (UNIT=foutput_notms%unit, FMT='(A23)')      '# Wavelength (microns) '
                WRITE (UNIT=foutput_notms%unit,FMT=FMT_RADOUT)         wlambda(l)
                WRITE (UNIT=foutput_notms%unit, FMT='(A)')        '# CPU time for RT in sec(s) '
                WRITE (UNIT=foutput_notms%unit, FMT=FMT_OUT)      rt_cputime(l)
                SELECT CASE (nmat)
                CASE(1) 
                   WRITE (UNIT=foutput_notms%unit, FMT='(A101)')      '#         sza               vza              vaa        &
                        &      theta              I                 dI'
                CASE(3) 
                   WRITE (UNIT=foutput_notms%unit, FMT='(A197)')    '#         sza               vza              vaa         &
                        &      theta              I                 Q                U                 dI               &
                        & dQ               &
                        & dU              Polar. (%)'
                CASE(4)  
                   WRITE (UNIT=foutput_notms%unit, FMT='(A233)')    '#         sza               vza              vaa         &
                        &      theta              I                 Q                U                  V               &
                        & dI               &
                        & dQ                dU                dV              Polar. (%)'
                END SELECT
                DO i = 1, nsza
                   UMU0 = COS(sza(i)*deg2rad)
                   DO j = nvza, 1, -1
                      UMU = COS(vza(j)*deg2rad)
                      DO k = 1, nvaa
                         thetad = ACOS(-UMU0 * UMU + SQRT( ( 1.0_dp-UMU0**2.0_dp ) * &
                              ( 1.0_dp - UMU**2.0_dp ) ) * COS( vaa(k)*deg2rad )) / deg2rad
                         SELECT CASE (nmat)
                         CASE(1) 
                            WRITE (UNIT=foutput_notms%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, &
                                 Svr_notms(i, j, k, 1, l), &
                                 Svr_sig_notms(i, j, k, 1, l)
                         CASE(3) 
                            polar = SQRT(Svr_notms(i, j, k, 2, l)**2.0_dp + Svr_notms(i, j, k, 3, l)**2.0_dp) /&
                                 Svr_notms(i, j, k, 1, l) * 100.0_dp
                            if ((polar .gt. 100.0_dp) .or. (polar .lt. 0.0_dp)) polar = undef_dp
                            WRITE (UNIT=foutput_notms%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, &
                                 (Svr_notms(i, j, k, u, l), u=1,nmat), &
                                 (Svr_sig_notms(i, j, k, u, l), u=1,nmat), polar
                         CASE(4) 
                            polar = SQRT(Svr_notms(i, j, k, 2, l)**2.0_dp + Svr_notms(i, j, k, 3, l)**2.0_dp + &
                                 Svr_notms(i, j, k, 4, l)**2.0_dp) / &
                                 & Svr_notms(i, j, k, 1, l) * 100.0_dp
                            if ((polar .gt. 100.0_dp) .or. (polar .lt. 0.0_dp)) polar = undef_dp
                            WRITE (UNIT=foutput_notms%unit, FMT=FMT_RADOUT) sza(i), vza(j), vaa(k), thetad, &
                                 (Svr_notms(i, j, k, u, l), u=1,nmat), &
                                 (Svr_sig_notms(i, j, k, u, l), u=1,nmat), polar 

                         END SELECT
                      ENDDO
                   ENDDO
                ENDDO
             ENDDO ! lambda
             CLOSE(UNIT = foutput_notms%unit)
          ENDIF
       endif! flux_only FALSE

       ! SAVE SOLAR_KDIS_[...].DAT IF COMPUTED DURING THE RUN
       if ((mode.eq.'kdis') .and. (kdis_solrad.eqv..false.)) then
          filename_tmp = TRIMCAT(artdeco_path, dir_solrad)
          filename_tmp = TRIMCAT(filename_tmp, fsolrad%name)
          filename_tmp = TRIMCAT(filename_tmp, 'kdis_')
          filename_tmp = TRIMCAT(filename_tmp, TRIM(ADJUSTL(kdis_model)))
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, solrad)
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          OPEN (UNIT = fsolrad%unit, FILE = filename_tmp, STATUS = 'NEW')
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# Solar spectrum integrated over k-distribution bands'
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# Solar constant (W m-2) :'
          WRITE (UNIT=fsolrad%unit, FMT=FMT_RADOUT) kdis_solrad_solcste
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# Number of sampled values :'
          WRITE (UNIT=fsolrad%unit, FMT=*) kdis_solrad_nwvl
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# '
          WRITE (UNIT=fsolrad%unit, FMT='(A)') '# lambda (microns)         Spectr. (W m-2)'
          DO i = 1, kdis_solrad_nwvl
             WRITE (UNIT=fsolrad%unit, FMT='(F12.8, 3x, E18.9)') kdis_solrad_wvl(i), kdis_solrad_F0(i)
          ENDDO
          CLOSE(fsolrad%unit)
       endif


       if ((print_aitaui).and.(mode.eq.'kdis')) then
          IF (verbose) THEN
             WRITE (*, *)                   '  (sub. write_data) save ai,taui'
             WRITE (*,*) ''
          ENDIF
          filename_tmp = TRIMCAT (dir_root_out, fkdisaitaui%name)
          filename_tmp = TRIMCAT(filename_tmp, kdis_model)
          filename_tmp = TRIMCAT(filename_tmp, '_')
          filename_tmp = TRIMCAT(filename_tmp, atm_type)
          filename_tmp = TRIMCAT(filename_tmp, '.dat')
          OPEN (UNIT = fkdisaitaui%unit, FILE = filename_tmp, STATUS = 'unknown')
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '#' 
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# Effective (ai, taui) for each wavelength band and layer ' 
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# Taui is for gas absorption only (continuum included)'  
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '#' 
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# Number of layers ' 
          WRITE (UNIT=fkdisaitaui%unit, FMT=*) nlayers 
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# Number of wavelengths bands' 
          WRITE (UNIT=fkdisaitaui%unit, FMT=*) nlambda
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '#################################################' 
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# NOTE : layers are sorted in decreasing altitude'
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '#        (i.e. TOA first)'
          WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '#################################################' 
          DO j = 1, nlambda
             CALL GET_MONO_LAYERS(j)
             ikdislamb = lambda_ikdis(j)
             WRITE (UNIT=fkdisaitaui%unit, FMT='(A, 2x, I5)') '# Wavelengh band ', j
             WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '#  lamb_min,  lamb_center, lamb_max :'
             WRITE (UNIT=fkdisaitaui%unit, FMT='(3(2x,E18.9))') &
                  kdis_wvlband(2,ikdislamb), kdis_wvlband(1,ikdislamb), kdis_wvlband(3,ikdislamb)
             WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# nai :'
             WRITE (UNIT=fkdisaitaui%unit, FMT='(I5)')  n_aitaui_mono_layers
             WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# ai :'
             WRITE (UNIT=fkdisaitaui%unit, FMT='(200(2x,E18.9))') (ai_mono_layers(iai), iai=1, n_aitaui_mono_layers)
             WRITE (UNIT=fkdisaitaui%unit, FMT='(A)') '# taui :'
             do k = 1, nlayers
                WRITE (UNIT=fkdisaitaui%unit, FMT='(200(2x,E18.9))') &
                     (gas_dtau_mono_layers(k,iai)*(1._dp-gas_ssa_mono_layers(k,iai)), iai=1, n_aitaui_mono_layers)           
             end do
          ENDDO ! on lambda
          CLOSE (UNIT = fkdisaitaui%unit)
       endif


    ENDIF ! do_rt

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

    DEALLOCATE(wlambda)

    IF (nptcle.gt.0)then
       DEALLOCATE(ptcle_type,     &
            ptcle_hg,             &
            ptcle_opt_interp,     &
            ptcle_vdist_type,     &
            ptcle_h_vdist,        &
            ptcle_sd_gauss_vdist, &
            ptcle_tau_reflamb,    &
            ptcle_user_vdist_nalt,& 
            ptcle_user_vdist_alt, &
            ptcle_user_vdist )
       DEALLOCATE(ptcle_nelem,&
            ptcle_flag_store_opt,&
            ptcle_opt_exist_reflamb,&
            ptcle_opt_exist, &
            ptcle_ilambda,  &
            ptcle_ireflamb)
       IF (flag_store_opt) THEN
          DEALLOCATE(ptcle_opt_store,        &
               ptcle_nlamb_phasemat_store, & 
               ptcle_nang_phasemat_store,  &
               ptcle_lamb_phasemat_store,  &
               ptcle_u_phasemat_store,     &
               ptcle_phasemat_store)
       END IF
       DEALLOCATE(ptcle_material, &
            ptcle_opt_model,      &
            ptcle_param_mie3,     &
            ptcle_param_hm, &
            ptcle_param_baum)
       if (nmaterial .gt. 0) DEALLOCATE(material_name, material_refind, material_nlamb_refind)
       DEALLOCATE(ptcle_flag_new_opt)
       if (flag_new_opt) THEN
          DEALLOCATE(ptcle_opt_new, &
               ptcle_nlamb_phasemat_new, & 
               ptcle_nang_phasemat_new, &
               ptcle_lamb_phasemat_new, &
               ptcle_u_phasemat_new, &
               ptcle_phasemat_new)
       end if
       IF (OPT_ONLY.EQV..FALSE.) THEN
          DEALLOCATE(ptcle_Cext_reflamb,&
               ptcle_opt,           &
               ptcle_nang_phasemat, &
               ptcle_u_phasemat,    &
               ptcle_phasemat)
       END IF
       IF (OPT_ONLY.EQV..FALSE.) THEN
          DEALLOCATE(ptcle_nbetal,&
               ptcle_betal, ptcle_betal_flag)
          DEALLOCATE(ptcle_trunc_normphasemat,&
               ptcle_trunccoeff,&
               ptcle_trunc_phasemat)
       END IF
    endif ! nptcle >0

    if(do_rt) then

       if (mode.eq.'kdis') then
          if (filter) then
             DEALLOCATE(rt_cputime_filter, &
                  filter_name,&
                  nlamb_filter,&
                  lamb_filter,&
                  trans_filter, &
                  bound_lamb_filter)
          endif
          DEALLOCATE(lambda_ikdis)
          DEALLOCATE(kdis_solrad_F0, kdis_solrad_wvl)
          if (kdis_solrad.eqv..false.) then
             DEALLOCATE(solrad_F0, solrad_wvl)
          endif
       endif

       DEALLOCATE(varsol_fact)
       DEALLOCATE(F0)
       DEALLOCATE(sza, vza, vaa)
       DEALLOCATE(depol)

       SELECT CASE(surface_family)
       CASE('lambert')
          DEALLOCATE(surface_albedo)
       CASE('brdf')
          IF (surface_type .ne. 'ocean') then
             DEALLOCATE(surface_par_brdf,surface_par_bpdf)
          else
            DEALLOCATE(ocean_Rsw)
          end IF

       end SELECT
       DEALLOCATE(alt_atm, &
            p_atm, &
            t_atm)
       if (mode.eq.'kdis') then
          DEALLOCATE(kdis_wvlband,&
               kdis_p,&
               kdis_t)
          DEALLOCATE(u_air_atm, &
               gas_u_atm, &
               u_atm)
          if (kdis_nsp.gt.0) then
             DEALLOCATE(kdis_species,&
                  kdis_fcont,&
                  kdis_nai, &
                  kdis_ki,&                  
                  kdis_ai, &
                  kdis_xsect)
          endif
          if (kdis_nsp_c.gt.0) then
             DEALLOCATE(kdis_species_c,&
                  kdis_fcont_c, &
                  kdis_nai_c, &
                  kdis_ki_c,&                  
                  kdis_ai_c, &
                  kdis_xsect_c,&
                  kdis_c)
          endif
       endif
       DEALLOCATE(rt_cputime)
       DEALLOCATE(rho_surf)
       DEALLOCATE(SvR)
       DEALLOCATE(SvR_levels)
       DEALLOCATE(SvR_sig)
       DEALLOCATE(flux_out)

       if (mode.eq.'kdis') then
          if (filter) then
             DEALLOCATE(SvR_filter, & 
                  SvR_sig_filter, &
                  flux_filter)
          else
             DEALLOCATE(fluxint_out)
          endif
       endif

       if (sinsca_corint .eqv. .true.) then
          DEALLOCATE(SvR_levels_notms)
          DEALLOCATE(Svr_notms)
          DEALLOCATE(Svr_sig_notms)
       endif

       DEALLOCATE(nai_layers,&
            ai_layers)

       DEALLOCATE(gas_dtauabs_layers)

          DEALLOCATE(ptcle_dtau_layers, &
               gas_dtausca_layers, &
               ptcle_ssa_layers, & 
               gas_dtau_mono_layers,&
               gas_ssa_mono_layers)

       
       IF ( (rt_model .eq. 'mcrad1d') .or. &
            (rt_model .eq. 'sinsca') .or. &
            (sinsca_corint)) THEN
          DEALLOCATE(ptcle_nang_phasemat_layers, & 
               ptcle_phasemat_layers, &
               ptcle_u_phasemat_layers)
       END IF

       IF ( (rt_model.eq.'disort') .or. &
            (rt_model.eq.'addoub') .or. &
            (rt_model.eq.'doad') ) THEN
          DEALLOCATE(ptcle_dtauabs_layers,&
               tot_dtausca_layers,&
               nbetal_layers,&
               betal_layers,&
               tot_dtau_mono_layers,&
               ssa_mono_layers)
       ENDIF
       DEALLOCATE(ai_mono_layers)

       DEALLOCATE(ptcle_ssa_layers_unscaled, &
            ptcle_ssa_layers_tms, &
            ptcle_dtau_layers_unscaled)
       
       if (sinsca_corint) then

          if (allocated(ptcle_phasemat_actual_layers)) then
             DEALLOCATE(ptcle_phasemat_actual_layers)
          end if

          DEALLOCATE(ptcle_phasemat_tms_layers)
          if (is_secsca) then
             DEALLOCATE(ptcle_phasemat_delta_layers, &
                  ptcle_u_phasemat_delta_layers, &
                  ptcle_secsca_mu_eq, &
                  ptcle_secsca_neg_phas, &
                  ptcle_secsca_norm_phas)
          end if
       endif
    endif

    if (get_possol) DEALLOCATE(h_tu, day)

    if (flag_baum) then
       DEALLOCATE(lamb_baum, &
            de_baum, &
            Cext_baum, &
            ssa_baum, &
            g_baum, &
            mu_baum, &
            ptcle_phasemat_baum)
    endif

    CALL CPU_TIME(end_tcpu)
    CALL DATE_AND_TIME(char8, char10, char5, end_date)

    IF (verbose) THEN
       WRITE (*,*) ''
       WRITE (*,FMT='(1x,A43)')           '          Exit ARTDECO OK                  '
       WRITE (*,FMT='(1x,A17,E16.5, A9)') '      - CPU time ',  (end_tcpu-start_tcpu), ' sec(s) -'
       WRITE (*,*)                       '==========================================='
    ENDIF

  END SUBROUTINE WRITE_DATA

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

  FUNCTION COUNTCOL(file_input) RESULT(N)
    USE MCONSTANTS,only:undef_dp
    IMPLICIT NONE
    INTEGER, INTENT (IN) :: file_input
    ! to read column of unknown length
    INTEGER, PARAMETER :: MAXCOL = 500
    INTEGER :: N, IO
    REAL, DIMENSION(MAXCOL+1) :: FCOL
    !--------
    FCOL(:) = undef_dp
    READ(UNIT=file_input, FMT=*, IOSTAT = IO) ( FCOL(N), N = 1, MAXCOL+1 )
    BACKSPACE(UNIT=file_input)
    BACKSPACE(UNIT=file_input)
    IF (FCOL( MAXCOL+1).ne.undef_dp) THEN
       write(*,*) ''
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       WRITE(*,*) ' (in function.COUNTCOL) ERROR : '
       write(*,*) '                     NCOL > MAXCOL '
       write(*,*) '                     MAXCOL =', MAXCOL
       WRITE(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       write(*,*) ''
       STOP
    ENDIF
    N = 1
    DO WHILE((FCOL(N).ne.undef_dp))
       N=N+1
    ENDDO
    N=N-1
  END FUNCTION COUNTCOL

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

  FUNCTION COUNTLINE(file_input, filename) RESULT(N)
    USE MCONSTANTS, only : max_len
    IMPLICIT NONE
    INTEGER, INTENT (IN)    :: file_input
    CHARACTER (LEN=*) :: filename
    LOGICAL                 :: nocomment
    CHARACTER (LEN=max_len) :: line
    INTEGER :: IO
    INTEGER :: N, I
    !--------
    nocomment = .TRUE.
    N = 0
    io=0
    DO WHILE ((nocomment).and.(io.eq.0))
       READ (file_input, FMT='(A100)', iostat=io) line
       line = ADJUSTL(line)
       IF (line(1:1) /= '#') THEN
          nocomment = .TRUE.
          N = N + 1
       ELSE
          nocomment = .FALSE.
       END IF
    ENDDO
    if (io.lt.0) then
       write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       write(*,*) '  (in countline) ERROR   :    '
       write(*,*) '  End of file reached in file '
       write(*,*) TRIM(ADJUSTL(filename))
       write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
       STOP
    endif
    DO I=1,N+1
       BACKSPACE(file_input)
    ENDDO
  END FUNCTION COUNTLINE

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

  SUBROUTINE SKIP_COM(file_input)
    ! skips comment lines beginning with #
    USE MCONSTANTS, only : max_len
    IMPLICIT NONE
    INTEGER, INTENT (IN)    :: file_input
    LOGICAL                 :: comment_found
    CHARACTER (LEN=max_len) :: line
    ! Skip comments
    comment_found = .TRUE.
    DO WHILE (comment_found)
       READ (file_input, FMT='(A100)') line
       line = ADJUSTL(line)
       IF (line(1:1) /= '#') comment_found = .FALSE.
    ENDDO
    BACKSPACE(file_input)
  END SUBROUTINE SKIP_COM

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

  FUNCTION READ_COM(file_input) RESULT(line)
    ! reads comment lines beginning with #
    ! returns 1st line that is not a comment as
    ! CHARACTER
    USE MCONSTANTS, only : max_len
    IMPLICIT NONE
    INTEGER, INTENT (IN)    :: file_input
    LOGICAL                 :: comment_found
    CHARACTER (LEN=max_len) :: line
    ! Skip comments
    comment_found = .TRUE.
    DO WHILE (comment_found)
       READ (file_input, FMT='(A100)') line
       line = ADJUSTL(line)
       IF (line(1:1) /= '#') comment_found = .FALSE.
    ENDDO
  END FUNCTION READ_COM

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

  FUNCTION UPCASE(string_in, n) RESULT(upper)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: string_in
    INTEGER, INTENT(IN) :: n
    CHARACTER(LEN=n)   :: upper
    INTEGER                      :: j
    DO j = 1,len(string_in)
       IF(string_in(j:j) >= "a" .AND. string_in(j:j) <= "z") THEN
          upper(j:j) = achar(iachar(string_in(j:j)) - 32)
       ELSE
          upper(j:j) = string_in(j:j)
       ENDIF
    ENDDO
  END FUNCTION UPCASE

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

  FUNCTION LOWCASE(string_in, n) RESULT(lower)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: string_in
    INTEGER, INTENT(IN) :: n
    CHARACTER(LEN=n)   :: lower
    INTEGER                      :: j
    DO j = 1,len(string_in)
       IF(string_in(j:j) >= "A" .AND. string_in(j:j) <= "Z") THEN
          lower(j:j) = achar(iachar(string_in(j:j)) + 32)
       ELSE
          lower(j:j) = string_in(j:j)
       ENDIF
    ENDDO
  END FUNCTION LOWCASE

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

  SUBROUTINE CHECK_NSTR(NSTR, UMU0, MAXCMU, CHNSTR)

    USE MCONSTANTS, only : dp
    USE nr, ONLY : gauleg

    INTEGER, INTENT(INOUT) :: NSTR
    REAL(KIND=DP), INTENT(IN) :: UMU0
    INTEGER, INTENT(IN) :: MAXCMU
    LOGICAL, INTENT(OUT) :: CHNSTR

    INTEGER :: NN
    INTEGER :: IPASS
    INTEGER :: IQ

    REAL(KIND=DP) :: CMU(MAXCMU)
    REAL(KIND=DP) :: CWT(MAXCMU)

    LOGICAL TEST

    IF (NSTR > MAXCMU) THEN
       PRINT*, ""
       PRINT*, " ( CHECK_NSTR ) ERROR"
       PRINT*, " ( CHECK_NSTR ) NSTR > MAXCMU"
       PRINT*, ""
    END IF

    CHNSTR = .FALSE.

    IPASS=0
2   NN = NSTR / 2
    TEST=.TRUE.
    CALL gauleg( DBLE(0.0), DBLE(1.0), CMU(1:NN),CWT(1:NN) )

    ! print*, UMU0
    ! print*, NN
    ! print*, CMU(1:NN)
    ! print*, CWT(1:NN)

    DO IQ = 1, NN
       IF ( ABS(UMU0-CMU(IQ))/UMU0 .LT. 1.E-4)THEN
          NSTR=NSTR+2
          TEST=.FALSE.
          IPASS=IPASS+1
          CHNSTR = .TRUE.
          GOTO 1
       ENDIF
    ENDDO
1   IF(.NOT.TEST.AND.IPASS.LT.4)GOTO 2

    DO IQ = 1, NN
       IF ( ABS(UMU0-CMU(IQ))/UMU0 .LT. 1.E-4 ) THEN
          PRINT*, '(CHECK_NSTR) BEAM ANGLE=COMPUTATIONAL ANGLE; CHANGE NSTR'
          STOP
       ENDIF
    ENDDO
    RETURN

  END SUBROUTINE CHECK_NSTR


END MODULE MINOUT
