! Fortran module for computation of aerosol opacities
! author : Antoine Bierjon, November 2022

MODULE aeropt_mod

 IMPLICIT NONE
      
 integer,save :: nwvl                            ! Number of wavelengths in the domain (VIS or IR)
 integer,save :: nsize                           ! Number of particle sizes stored in the file
 real,dimension(:),save,allocatable :: wvl       ! Wavelength axis [m]
 real,dimension(:),save,allocatable :: radiusdyn ! Particle size axis [m]
 real,dimension(:,:),save,allocatable :: Qext    ! Extinction efficiency coefficient [/]
 real,dimension(:,:),save,allocatable :: omeg    ! Single scattering albedo [/]
 
 integer,save :: aeroptlogfileID = 100           ! Log file's unit ID

 CONTAINS 
 
!*******************************************************************************
 SUBROUTINE read_optpropfile(datadir,optpropfile)
  !==============================================================================
  ! Purpose:
  ! Read optical properties from ASCII file given as input
  ! and store them in Fortran arrays
  !==============================================================================

  use netcdf

  implicit none

  !==============================================================================
  ! Arguments:
  !==============================================================================
  character(len=256), intent(in) :: datadir     ! directory containing the ASCII file
  character(len=128), intent(in) :: optpropfile ! name of the ASCII optical properties file
  
  ! No intent(out) arguments, since only aeropt_mod module variable are updated
  ! by this subroutine : wvl,radiusdyn,Qext,omeg
  
  !==============================================================================
  ! Local variables:
  !==============================================================================
  logical :: file_ok                         ! to know if the file can be opened
  integer :: file_unit = 60                  ! ASCII file ID
  integer :: jfile                           ! ASCII file scan index
  logical :: endwhile                        ! Reading loop boolean
  character(len=132) :: scanline             ! ASCII file scanning line
  integer :: read_ok                         ! to know if the line is well read
  integer :: iwvl,isize                      ! Wavelength and Particle size loop indices
  
  !==============================================================================
  ! aeropt_mod module variables used here:
  !==============================================================================
  ! nwvl,nsize,wvl,radiusdyn,Qext,omeg

  !==========================================================================
  ! 1. Open the ASCII file
  !==========================================================================
  INQUIRE(FILE=trim(datadir)//'/'//trim(optpropfile),EXIST=file_ok)
  if(.not.file_ok) then
    write(*,*)'Problem opening ',trim(optpropfile)
    write(*,*)'It should be in: ',trim(datadir)
    stop 
  endif
  OPEN(UNIT=file_unit,FILE=trim(datadir)//'/'//trim(optpropfile),FORM='formatted')

  !==========================================================================
  ! 2. Retrieve the optical properties' dimensions
  !==========================================================================
  jfile = 1
  endwhile = .false.
  DO WHILE (.NOT.endwhile)
    READ(file_unit,*,iostat=read_ok) scanline
    if (read_ok.ne.0) then
      write(*,*)'Error reading file ',&
                trim(datadir)//'/'//trim(optpropfile)
      stop
    endif
    IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
      BACKSPACE(file_unit)
reading1_seq: SELECT CASE (jfile) ! FIRST READING SEQUENCE
      CASE(1) reading1_seq ! nwvl ----------------------------
          read(file_unit,*,iostat=read_ok) nwvl
          if (read_ok.ne.0) then
            write(*,*)'reading1_seq: Error while reading line: ',&
                      trim(scanline)
            write(*,*)'   of file ',&
                      trim(datadir)//'/'//trim(optpropfile)
            stop
          endif
          jfile = jfile+1
      CASE(2) reading1_seq ! nsize ---------------------------
          read(file_unit,*,iostat=read_ok) nsize
          if (read_ok.ne.0) then
            write(*,*)'reading1_seq: Error while reading line: ',&
                      trim(scanline)
            write(*,*)'   of file ',&
                      trim(datadir)//'/'//trim(optpropfile)
            stop
          endif
          endwhile = .true.
      CASE DEFAULT reading1_seq ! default --------------------
          write(*,*) 'reading1_seq: ',& 
                     'Error while loading optical properties.' 
          stop 
      END SELECT reading1_seq ! ==============================
    ENDIF
  ENDDO ! DO WHILE(.NOT.endwhile)

  write(*,*) "nwvl=",nwvl," ; nsize=",nsize

  ALLOCATE(wvl(nwvl))             ! Wavelength axis
  ALLOCATE(radiusdyn(nsize))      ! Aerosol radius axis
  ALLOCATE(Qext(nwvl,nsize))      ! Extinction efficiency coeff
  ALLOCATE(omeg(nwvl,nsize))      ! Single scattering albedo

  !==========================================================================
  ! 3. Retrieve the optical properties
  !==========================================================================
  jfile = 1
  endwhile = .false.
  DO WHILE (.NOT.endwhile)
    READ(file_unit,*) scanline
    IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
      BACKSPACE(file_unit)
reading2_seq: SELECT CASE (jfile) ! SECOND READING SEQUENCE
      CASE(1) reading2_seq ! wvl -----------------------------
          read(file_unit,*,iostat=read_ok) wvl
          if (read_ok.ne.0) then
            write(*,*)'reading2_seq: Error while reading line: ',&
                      trim(scanline)
            write(*,*)'   of file ',&
                      trim(datadir)//'/'//trim(optpropfile)
            stop
          endif
          jfile = jfile+1
      CASE(2) reading2_seq ! radiusdyn -----------------------
          read(file_unit,*,iostat=read_ok) radiusdyn
          if (read_ok.ne.0) then
            write(*,*)'reading2_seq: Error while reading line: ',&
                      trim(scanline)
            write(*,*)'   of file ',&
                      trim(datadir)//'/'//trim(optpropfile)
            stop
          endif
          jfile = jfile+1
      CASE(3) reading2_seq ! Qext ----------------------------
          isize = 1
          DO WHILE (isize .le. nsize)
            READ(file_unit,*,iostat=read_ok) scanline
            if (read_ok.ne.0) then
              write(*,*)'reading2_seq: Error while reading line: ',&
                        trim(scanline)
              write(*,*)'   of file ',&
                        trim(datadir)//'/'//trim(optpropfile)
              stop
            endif
            IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
              BACKSPACE(file_unit)
              read(file_unit,*) Qext(:,isize)
              isize = isize + 1
            ENDIF
          ENDDO
          jfile = jfile+1
      CASE(4) reading2_seq ! omeg ----------------------------
          isize = 1
          DO WHILE (isize .le. nsize)
            READ(file_unit,*,iostat=read_ok) scanline
            if (read_ok.ne.0) then
              write(*,*)'reading2_seq: Error while reading line: ',&
                        trim(scanline)
              write(*,*)'   of file ',&
                        trim(datadir)//'/'//trim(optpropfile)
              stop
            endif
            IF ((scanline(1:1).ne.'#').and.(scanline(1:1).ne.' ')) THEN
              BACKSPACE(file_unit)
              read(file_unit,*) omeg(:,isize)
              isize = isize + 1
            ENDIF
          ENDDO
          endwhile = .true.
      CASE DEFAULT reading2_seq ! default --------------------
          write(*,*) 'reading2_seq: ',&
                     'Error while loading optical properties.'
          stop
      END SELECT reading2_seq ! ==============================
    ENDIF
  ENDDO

  ! Close the file
  CLOSE(file_unit)

  write(*,*)""
  write(*,*) "Wavelength array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
  write(*,*) "ranging from ",wvl(1)," to ",wvl(nwvl)," meters"
  write(*,*) "Effective radius array loaded from file ",trim(datadir)//'/'//trim(optpropfile)
  write(*,*) "ranging from ",radiusdyn(1)," to ",radiusdyn(nsize)," meters"
 
 
 END SUBROUTINE read_optpropfile
 
  
!*******************************************************************************
 
 SUBROUTINE interp_wvl_reff(wvl_val,reff_val,missval,&
                            Qext_val,omeg_val)
  !==============================================================================
  ! Purpose:
  ! Interpolate the optical properties to a given wavelength and effective radius
  ! from properties tables
  !==============================================================================

  use netcdf

  implicit none

  !==============================================================================
  ! Arguments:
  !==============================================================================
  real, intent(in) :: wvl_val                    ! Reference wavelength [m]
  real, intent(in) :: reff_val                   ! Input reff [m]
  real, intent(in) :: missval                    ! Value to put in outfile when the reff is out of bounds
  
  real, intent(out) :: Qext_val                  ! Qext after both interpolations
  real, intent(out) :: omeg_val                  ! omega after both interpolations
  
  
  !==============================================================================
  ! Local variables:
  !==============================================================================
  integer :: iwvl,isize                 ! Wavelength and Particle size loop indices
  integer :: iwvl1,iwvl2,isize1,isize2  ! Wavelength and Particle size indices for the interpolation
  real :: coeff                         ! Interpolation coefficient
  real,dimension(2) :: reffQext         ! Qext after reff interpolation
  real,dimension(2) :: reffomeg         ! omeg after reff interpolation
  
  !==============================================================================
  ! aeropt_mod module variables used here:
  !==============================================================================
  ! nwvl,nsize,wvl,radiusdyn,Qext,omeg,aeroptlogfileID

  !==========================================================================
  ! 1. Get the optpropfile wavelength values encompassing the ref wavelength
  !==========================================================================
  iwvl=1
  DO WHILE ((iwvl.lt.nwvl).and.(wvl(iwvl).lt.wvl_val))
    iwvl=iwvl+1
  ENDDO
  if ((wvl(nwvl).lt.wvl_val).or.(wvl(1).gt.wvl_val)) then
    write(*,*) "ERROR: the reference wavelength (",wvl_val,"m) is out of the bounds"
    write(*,*) "of the optical properties' tables"
    write(*,*) "(wvl_inf=",wvl(1),"m ; wvl_sup=",wvl(nwvl),"m)"
    write(*,*) "Ensure you wrote it well (in meters),"
    write(*,*) "or supply new optical properties' tables."
    stop
  endif
  if (wvl(iwvl).eq.wvl_val) then
    ! if the ref wavelength is in the optpropfile
    iwvl1 = iwvl
    iwvl2 = iwvl
  else ! wvl(iwvl)>wvl_val and wvl(iwvl-1)<wvl_val
    iwvl1 = iwvl-1
    iwvl2 = iwvl
  endif
    
  !==========================================================================
  ! 2. Get the optpropfile reff values encompassing the input reff
  !==========================================================================
      
  isize=1
  DO WHILE ((isize.lt.nsize).and.(radiusdyn(isize).lt.reff_val))
    isize=isize+1
  ENDDO

  if (reff_val.eq.missval) then
    ! if the effective radius is a missing value,
    ! put a missing value for the optical properties
    Qext_val = missval
    omeg_val = missval
!    write(*,*) "reff missing value detected"
    
  else if ((radiusdyn(nsize).lt.reff_val).or.(radiusdyn(1).gt.reff_val)) then
    write(aeroptlogfileID,*) "WARNING: the effective radius (",reff_val,"m) is out of the bounds"
    write(aeroptlogfileID,*) "of the optical properties' tables"
    write(aeroptlogfileID,*) "(reff_inf=",radiusdyn(1),"m ; reff_sup=",radiusdyn(nsize),"m)"
    write(aeroptlogfileID,*) "A missing value will be written for optical properties."

    Qext_val = missval
    omeg_val = missval

  else
    
  !==========================================================================
  ! 3. Interpolate the properties
  !==========================================================================
    if (radiusdyn(isize).eq.reff_val) then
      ! if the user reff is exactly in the optpropfile
      isize1 = isize
      isize2 = isize
    else ! radius(isize)>reff and radiusdyn(isize-1)<reff
      isize1 = isize-1
      isize2 = isize
    endif

    ! Qext COMPUTATION
    ! Linear interpolation in effective radius
    if (isize2.ne.isize1) then
      coeff = (reff_val-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
    else
      coeff = 0.
    endif
    reffQext(1) = Qext(iwvl1,isize1)+coeff*(Qext(iwvl1,isize2)-Qext(iwvl1,isize1))
    reffQext(2) = Qext(iwvl2,isize1)+coeff*(Qext(iwvl2,isize2)-Qext(iwvl2,isize1))
    ! Linear interpolation in wavelength
    if (iwvl2.ne.iwvl1) then
      coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
    else
      coeff = 0.
    endif
    Qext_val = reffQext(1)+coeff*(reffQext(2)-reffQext(1))             


    ! omeg COMPUTATION
    ! Linear interpolation in effective radius
    if (isize2.ne.isize1) then
      coeff = (reff_val-radiusdyn(isize1))/(radiusdyn(isize2)-radiusdyn(isize1))
    else
      coeff = 0.
    endif
    reffomeg(1) = omeg(iwvl1,isize1)+coeff*(omeg(iwvl1,isize2)-omeg(iwvl1,isize1))
    reffomeg(2) = omeg(iwvl2,isize1)+coeff*(omeg(iwvl2,isize2)-omeg(iwvl2,isize1))
    ! Linear interpolation in wavelength
    if (iwvl2.ne.iwvl1) then
      coeff = (wvl_val-wvl(iwvl1))/(wvl(iwvl2)-wvl(iwvl1))
    else
      coeff = 0.
    endif
    omeg_val = reffomeg(1)+coeff*(reffomeg(2)-reffomeg(1))

  endif
 
 END SUBROUTINE interp_wvl_reff
   
!*******************************************************************************
 SUBROUTINE end_aeropt_mod
  !==============================================================================
  ! Purpose:
  ! Deallocate module variables
  !==============================================================================

  implicit none
  
  !==============================================================================
  ! aeropt_mod module variables used here:
  !==============================================================================
  ! wvl,radiusdyn,Qext,omeg

  if (allocated(wvl)) deallocate(wvl)
  if (allocated(radiusdyn)) deallocate(radiusdyn)
  if (allocated(Qext)) deallocate(Qext)
  if (allocated(omeg)) deallocate(omeg)

 END SUBROUTINE end_aeropt_mod 
 
!*******************************************************************************
 SUBROUTINE create_logfile(filename)
  !==============================================================================
  ! Purpose:
  ! Deallocate module variables
  !==============================================================================

  implicit none
  
  !==============================================================================
  ! Arguments:
  !==============================================================================
  character(len=128),intent(in) :: filename ! name of the aeropt log file
  !==============================================================================
  ! aeropt_mod module variables used here:
  !==============================================================================
  ! aeroptlogfileID
  
  ! Create output file
  write(*,*) "Creating aeropt log file: ",trim(filename)
  write(*,*) "It will overwrite any existing file with the same name."
  OPEN(UNIT=aeroptlogfileID,FILE=trim(filename),FORM='formatted',STATUS='replace')
  ! NB: STATUS='replace' will overwrite any existing file with the same name
  
 END SUBROUTINE create_logfile 
 
!*******************************************************************************

END MODULE aeropt_mod