subroutine get_haze_and_cloud_opacity(FTYPE,CTYPE,m0,m3,ich,tauext,wbar,gbar)

  use datafile_mod, only : datadir
  use radinc_h, only : L_NSPECTI, L_NSPECTV

  implicit none

  !============================================================================
  !
  !     Purpose
  !     -------
  !     Calculates the opacity of a GCM cell from the m0, m3
  !     and particle type information.
  !
  !
  !     INPUT: 
  !     ------
  !     FTYPE : which fractal particle type will be used in the GCM.
  !             This is used at the first call, and defined once for all.
  !             FTYPE = 1 : Df = 2.00
  !             FTYPE = 2 : Df = 2.15
  !             FTYPE = 3 : Df = 2.30
  !             FTYPE = 4 : Df = 2.45
  !
  !     CTYPE : which particle type is currently called.
  !             CTYPE = 0 : CLOUD DROPLETS
  !             CTYPE = 1 : FRACTAL AEROSOLS (Df = 2.00)
  !             CTYPE = 2 : FRACTAL AEROSOLS (Df = 2.15)
  !             CTYPE = 3 : FRACTAL AEROSOLS (Df = 2.30)
  !             CTYPE = 4 : FRACTAL AEROSOLS (Df = 2.45)
  !             CTYPE = 5 : SPHERICAL AEROSOLS (MODE OF SMALL AEROSOLS)
  !
  !     m0  = 0th order moment of particule.
  !     m3  = 3rd order moment of particule.
  !     ich = spectral channel [1,...,46].
  !
  !
  !     OUPUT:
  !     ------
  !     tauext = opacity.
  !     wbar   = average single scattering albedo.
  !     gbar   = average asymmetry parameter.
  !
  !
  !     Author
  !     ------
  !     B. de Batz de Trenquelleon (08/2022).
  !
  !============================================================================


  !-------------------------
  ! 0. Declarations
  !-------------------------

  ! INPUTS :
  !---------
  integer,intent(in) :: FTYPE
  integer,intent(in) :: CTYPE

  real*8,intent(in) :: m0
  real*8,intent(in) :: m3

  integer,intent(in) :: ich

  ! OUTPUTS :
  !----------
  real*8,intent(out) :: tauext
  real*8,intent(out) :: wbar
  real*8,intent(out) :: gbar
  
  ! PARAMETERS :
  !-------------
  integer :: nwl
  integer :: nrad
  integer :: ii,iii,kk,kkk

  parameter(nwl = L_NSPECTI + L_NSPECTV, nrad = 33)

  real*8 :: xwlnv(nwl)
  
  ! Spherical aerosols :
  real*8, save :: rc_as(nrad)
  real*8, save :: ssca_as(nwl,nrad),sabs_as(nwl,nrad),sext_as(nwl,nrad)
  real*8, save :: wbar_as(nwl,nrad),gbar_as(nwl,nrad)
  ! Fractal aerosols :
  real*8, save :: rc_f(nrad)
  real*8, save :: ssca_f(nwl,nrad),sabs_f(nwl,nrad),sext_f(nwl,nrad)
  real*8, save :: wbar_f(nwl,nrad),gbar_f(nwl,nrad)
  ! Cloud droplets :
  real*8, save :: rc_s(nrad)
  real*8, save :: ssca_s(nwl,nrad),sabs_s(nwl,nrad),sext_s(nwl,nrad)
  real*8, save :: wbar_s(nwl,nrad),gbar_s(nwl,nrad)
  
  ! LOCAL VARIABLES FOR INTERPOLATION :
  !------------------------------------
  integer :: idx
  real*8  :: rinit,rc,step,xfact
  character(len=100) :: aer_sph_file,aer_fra_file,cloud_file
  logical            :: file_ok
  logical,save       :: firstcall = .true.

  ! INITIALISATION :
  !-----------------
  tauext = 0.0
  wbar   = 0.0
  gbar   = 0.0

  !---------------------------------------------------
  ! 1. First call : open and read the first time
  !---------------------------------------------------
  
  ! Check spectral channel :
  if(ich.lt.1 .or. ich.gt.nwl) then
    print*, "ERROR : ich.lt.1 .or. ich.gt.nwl in get_haze_and_cloud_opacity"
    call abort
  endif

  if(firstcall) then

    ! Spherical aerosols :
    aer_sph_file = trim(datadir)//'/optical_tables/table_sigma_int2_sph3.00_rm50nm.dat'
    inquire(file=trim(aer_sph_file),exist=file_ok)
    if(.not.file_ok) then
      write(*,*) 'The file ',TRIM(aer_sph_file),' was not found by get_haze_and_cloud_opacity.F90 !'
      write(*,*) 'Check in ', TRIM(datadir),' that you have the one corresponding !!'
      call abort
    endif
    open(unit = 27, file = aer_sph_file)
    
    ! Fractal aerosols :
    if(FTYPE.eq.1) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.00_rm50nm.dat'
    if(FTYPE.eq.2) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.15_rm50nm.dat'
    if(FTYPE.eq.3) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.30_rm50nm.dat'
    if(FTYPE.eq.4) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.45_rm50nm.dat'
    inquire(file=trim(aer_fra_file),exist=file_ok)
    if(.not.file_ok) then
      write(*,*) 'The file ',TRIM(aer_fra_file),' was not found by get_haze_and_cloud_opacity.F90 !'
      write(*,*) 'Check in ', TRIM(datadir),' that you have the one corresponding !!'
      call abort
    endif
    open(unit = 26, file = aer_fra_file)
    
    ! And clouds :
    cloud_file = trim(datadir)//'/optical_tables/table_sigma_int2_sphere.dat'
    inquire(file=trim(cloud_file),exist=file_ok)
    if(.not.file_ok) then
      write(*,*) 'The file ',TRIM(cloud_file),' was not found by get_haze_and_cloud_opacity.F90 !'
      write(*,*) 'Check in ', TRIM(datadir),' that you have the one corresponding !!'
      call abort
    endif
    open(unit = 25, file = cloud_file)

    do ii = 1, nwl    ! 46 loops on wlnv
      do kk = 1, nrad ! 33 loops on rc 

        read(25,*) iii, kkk, xwlnv(ii), rc_s(kk), &
                   ssca_s(ii,kk), sabs_s(ii,kk), sext_s(ii,kk), &
                   gbar_s(ii,kk), wbar_s(ii,kk)
        if(ii.ne.iii .or. kk.ne.kkk) then
          write(*,*) 'lecture unit = 25'
          write(*,*) 'ii, iii : ', ii, iii
          write(*,*) 'kk, kkk : ', kk, kkk
          stop
        endif

        read(26,*) iii, kkk, xwlnv(ii), rc_f(kk), &
                   ssca_f(ii,kk), sabs_f(ii,kk), sext_f(ii,kk), &
                   gbar_f(ii,kk), wbar_f(ii,kk)
        if(ii.ne.iii .or. kk.ne.kkk) then
          write(*,*) 'lecture unit = 26'
          write(*,*) 'ii, iii : ', ii, iii
          write(*,*) 'kk, kkk : ', kk, kkk
          stop
        endif

        read(27,*) iii, kkk, xwlnv(ii), rc_as(kk), &
                   ssca_as(ii,kk), sabs_as(ii,kk), sext_as(ii,kk), &
                   gbar_as(ii,kk), wbar_as(ii,kk)
        if(ii.ne.iii .or. kk.ne.kkk) then
          write(*,*) 'lecture unit = 27'
          write(*,*) 'ii, iii : ', ii, iii
          write(*,*) 'kk, kkk : ', kk, kkk
          stop
       endif

      enddo ! End loop on kk
    enddo   ! End loop on ii

    if(firstcall) firstcall = .false.

    close(25)
    close(26)
    close(27)
  endif ! End firstcall


  !-------------------------
  ! 2. Interpolate
  !-------------------------

  rinit = 1.e-9                ! fin = rinit*step**(33-1) = 1.e-5
  if(CTYPE.eq.0) rinit = 1.e-7 ! fin = rinit*step**(33-1) = 1.e-3

  ! No alpha(k) !
  ! Because it is taken into account in the optical tables ...
  rc = (m3/m0)**(1./3.)
  
  step = 10.**(1./8.)
  idx = int(log(rc/rinit) / log(step)) + 1

  ! Spherical (cloud) :
  !--------------------
  IF(CTYPE.EQ.0) THEN

    if(abs(rc_s(1)-rinit)/rc_s(1).gt.1.e-5 .or. &
       abs(rc_s(nrad)-rinit*step**32)/rc_s(nrad).gt.1.e-5) then
        write(*,*) 'rinit          = ',rinit
        write(*,*) 'rc_s(1)        = ',rc_s(1)
        write(*,*) 'rinit*step**32 = ',rinit*step**32
        write(*,*) 'rc_s(nrad)     = ',rc_s(nrad)
        stop
    endif

    if(rc.le.rc_s(1)) then
      wbar   = wbar_s(ich,1)
      gbar   = gbar_s(ich,1)
      tauext = sext_s(ich,1)*(rc/rc_s(1))**3
    endif

    if(rc.ge.rc_s(nrad)) then
      wbar   = wbar_s(ich,nrad)
      gbar   = gbar_s(ich,nrad)
      tauext = sext_s(ich,nrad)*(rc/rc_s(nrad))**3.
    endif

    if(rc.gt.rc_s(1) .and. rc.lt.rc_s(nrad)) then
      xfact  = (log(rc)-log(rc_s(idx))) / (log(rc_s(idx+1))-log(rc_s(idx)))
      tauext = log(sext_s(ich,idx))*(1.-xfact) + log(sext_s(ich,idx+1))*xfact 
      tauext = exp(tauext)
      wbar   = wbar_s(ich,idx)*(1.-xfact) + wbar_s(ich,idx+1)*xfact 
      gbar   = gbar_s(ich,idx)*(1.-xfact) + gbar_s(ich,idx+1)*xfact 
    endif

  ENDIF ! (CTYPE.EQ.0)

  ! Spherical (haze) :
  !-------------------
  IF(CTYPE.EQ.5) THEN

    if(abs(rc_as(1)-rinit)/rc_as(1).gt.1.e-5 .or. &
       abs(rc_as(nrad)-rinit*step**32)/rc_as(nrad).gt.1.e-5) then
        write(*,*) 'rinit          = ',rinit
        write(*,*) 'rc_as(1)       = ',rc_as(1)
        write(*,*) 'rinit*step**32 = ',rinit*step**32
        write(*,*) 'rc_as(nrad)    = ',rc_as(nrad)
        stop
    endif

    if(rc.le.rc_as(1)) then
      wbar   = wbar_as(ich,1)
      gbar   = gbar_as(ich,1)
      tauext = sext_as(ich,1)*(rc/rc_as(1))**3
    endif

    if(rc.ge.rc_as(nrad)) then
      wbar   = wbar_as(ich,nrad)
      gbar   = gbar_as(ich,nrad)
      tauext = sext_as(ich,nrad)*(rc/rc_as(nrad))**3.
    endif

    if(rc.gt.rc_as(1) .and. rc.lt.rc_as(nrad)) then
      xfact  = (log(rc)-log(rc_as(idx))) / (log(rc_as(idx+1))-log(rc_as(idx)))
      tauext = log(sext_as(ich,idx))*(1.-xfact) + log(sext_as(ich,idx+1))*xfact
      tauext = exp(tauext)
      wbar   = wbar_as(ich,idx)*(1.-xfact) + wbar_as(ich,idx+1)*xfact
      gbar   = gbar_as(ich,idx)*(1.-xfact) + gbar_as(ich,idx+1)*xfact
    endif
  
  ENDIF ! (CTYPE.EQ.5)

  ! Fractal (haze) :
  !-----------------
  IF(CTYPE.NE.0 .AND. CTYPE.NE.5) THEN

    if(abs(rc_f(1)-rinit)/rc_f(1).gt.1.e-5 .or. &
       abs(rc_f(nrad)-rinit*step**32)/rc_f(nrad).gt.1.e-5) then
        write(*,*) 'rinit          = ',rinit
        write(*,*) 'rc_f(1)        = ',rc_f(1)
        write(*,*) 'rinit*step**32 = ',rinit*step**32
        write(*,*) 'rc_f(nrad)     = ',rc_f(nrad)
        stop
    endif

    if(rc.le.rc_f(1)) then
      wbar   =  wbar_f(ich,1)
      gbar   =  gbar_f(ich,1)
      tauext = sext_f(ich,1)*(rc/rc_f(1))**3
    endif

    if(rc.ge.rc_f(nrad)) then
      wbar   = wbar_f(ich,nrad)
      gbar   = gbar_f(ich,nrad)
      tauext = sext_f(ich,nrad)*(rc/rc_f(nrad))**3.
    endif
    
    if(rc.gt.rc_f(1) .and. rc.lt.rc_f(nrad)) then
      xfact  = (log(rc)-log(rc_f(idx))) / (log(rc_f(idx+1))-log(rc_f(idx)))
      tauext = log(sext_f(ich,idx))*(1.-xfact) + log(sext_f(ich,idx+1))*xfact 
      tauext = exp(tauext)
      wbar   = wbar_f(ich,idx)*(1.-xfact) + wbar_f(ich,idx+1)*xfact 
      gbar   = gbar_f(ich,idx)*(1.-xfact) + gbar_f(ich,idx+1)*xfact 
    endif

  ENDIF ! (CTYPE.NE.0 .AND. CTYPE.NE.5)

  ! Sanity check :
  !---------------
  IF (CTYPE.LT.0 .AND. CTYPE.GT.5) THEN
    write(*,*) 'CTYPE < 0 ET > 5  = ',CTYPE
    call abort
  ENDIF

  tauext = tauext * m0

  return

end subroutine get_haze_and_cloud_opacity