module gray_params

  use const_and_precisions, only : wp_
  implicit none
  integer, parameter :: lenfnm = 256
  integer, parameter :: headw = 132, headl = 21

  ! Antenna/wave launcher parameters
  type antenna_parameters
    ! From gray_params.data:
    real(wp_) :: alpha, beta        ! Launching angles
    real(wp_) :: power              ! Initial power
    real(wp_) :: psi, chi           ! Initial polarisation ellipse parameters
    integer :: iox                  ! Initial wave mode
    integer :: ibeam                ! Beam kind
    character(len=lenfnm) :: filenm ! beamdata.txt filename

    ! From beamdata.txt:
    real(wp_) :: fghz   ! EC frequency
    real(wp_) :: pos(3) ! Launcher position (tokamak frame)
    real(wp_) :: w(2)   ! Beam waists w(z) for the amplitude           (local frame)
    real(wp_) :: ri(2)  ! Beam inverse radii 1/R(z) for the phase      (local frame)
    real(wp_) :: phi(2) ! Axes orientation angles for amplitude, phase (local frame)
  end type

  ! MHD equilibrium parameters
  type equilibrium_parameters
    real(wp_) :: ssplps, ssplf      ! Spline tensions for ψ(R,Z), I(ψ)
    real(wp_) :: factb              ! Rescaling factor for B
    integer :: sgnb, sgni           ! Sign of B, I
    integer :: ixp                  ! Position of X point
    integer :: iequil               ! Equilibrium kind
    integer :: icocos               ! COCOS index
    integer :: ipsinorm             ! Normalisation of ψ
    integer :: idesc, ifreefmt      ! G-EQDSK quirks
    character(len=lenfnm) :: filenm ! Equilibrium data filepath
  end type

  ! Kinetic plasma profiles
  type profiles_parameters
    real(wp_) :: psnbnd             ! Plasma boundary (ψ: ne(ψ)=0)
    real(wp_) :: sspld              ! Density spline tension
    real(wp_) :: factne, factte     ! Rescale factors for n_e, T_e
    integer :: iscal                ! Rescaling model
    integer :: irho                 ! Radial coordinate
    integer :: iprof                ! Profile kind
    character(len=lenfnm) :: filenm ! Profiles data filepath
  end type

  ! Raytracing parameters
  type raytracing_parameters
    real(wp_) :: dst                ! Integration step size
    real(wp_) :: rwmax              ! Normalized maximum radius of beam power
    integer :: nray, nrayr, nrayth  ! Total/radial/angular number of rays
    integer :: igrad                ! Complex eikonal switch
    integer :: nstep                ! Max number of integration steps
    integer :: idst                 ! Choice of the integration variable
    integer :: ipass                ! Number of plasma passes
    integer :: ipol                 ! Whether to compute wave polarisation
  end type

  ! EC resonant heating & Current Drive parameters
  type ecrh_cd_parameters
    integer :: iwarm  ! Choice of power absorption model
    integer :: ieccd  ! Choice of current drive model
    integer :: ilarm  ! Larmor-radius expansion order
    integer :: imx    ! Max iterations for solving the dispersion relation
  end type

  ! Output data parameters
  type output_parameters
    integer :: ipec  ! Grid spacing for profiles (profili EC)
    integer :: nrho  ! Number of grid steps for EC profiles + 1
    integer :: istpr ! Subsampling factor for beam cross-section (units 8, 12)
    integer :: istpl ! Subsampling factor for outer rays (unit 33)
  end type

  ! Other parameters
  type misc_parameters
    real(wp_) :: rwall ! R of the limiter (fallback)
  end type

  ! MHD equilibrium data
  type equilibrium_data
    real(wp_), allocatable :: rv(:)     ! R of the uniform grid
    real(wp_), allocatable :: zv(:)     ! Z of the uniform grid
    real(wp_), allocatable :: rlim(:)   ! R of the limiter contour (wall)
    real(wp_), allocatable :: zlim(:)   ! Z of the limiter contour
    real(wp_), allocatable :: rbnd(:)   ! R of the boundary contour (plasma)
    real(wp_), allocatable :: zbnd(:)   ! Z of the boundary contour
    real(wp_), allocatable :: fpol(:)   ! Poloidal current function
    real(wp_), allocatable :: qpsi(:)   ! Safety factor on the flux grid
    real(wp_), allocatable :: psin(:,:) ! Poloidal flux on a uniform grid
    real(wp_), allocatable :: psinr(:)  ! Poloidal flux
    real(wp_) :: psia                   ! Poloidal flux at edge - flux at magnetic axis
    real(wp_) :: rvac                   ! Reference R₀ (B = B₀R₀/R without the plasma)
    real(wp_) :: rax                    ! R of the magnetic axis
    real(wp_) :: zax                    ! Z of the magnetic axis
  end type

  ! Kinetic plasma profiles data
  type profiles_data
    real(wp_), allocatable :: psrad(:) ! Radial coordinate
    real(wp_), allocatable :: terad(:) ! Electron temperature profile
    real(wp_), allocatable :: derad(:) ! Electron density profile
    real(wp_), allocatable :: zfc(:)   ! Effective charge profile
  end type

  ! All GRAY parameters
  type gray_parameters
    type(antenna_parameters)     :: antenna
    type(equilibrium_parameters) :: equilibrium
    type(profiles_parameters)    :: profiles
    type(raytracing_parameters)  :: raytracing
    type(ecrh_cd_parameters)     :: ecrh_cd
    type(output_parameters)      :: output
    type(misc_parameters)        :: misc
  end type

  ! All GRAY input data
  type gray_data
    type(equilibrium_data) :: equilibrium
    type(profiles_data)    :: profiles
  end type

  ! GRAY final results (only some)
  type gray_results
    real(wp_) :: pabs                 ! Total absorbed power
    real(wp_) :: icd                  ! Total driven current
    real(wp_), allocatable :: dpdv(:) ! Absorbed power density
    real(wp_), allocatable :: jcd(:)  ! Driven current density
  end type

  ! Global variables exposed for gray_core
  integer, save :: iequil, iprof, ipol
  integer, save :: iwarm, ilarm, imx, ieccd
  integer, save :: igrad, idst, ipass
  integer, save :: istpr0, istpl0
  integer, save :: ipec, nnd

contains

  subroutine print_parameters(params, strout)
    implicit none

    ! subroutine arguments
    type(gray_parameters), intent(in) :: params
    character(len=*), dimension(:), intent(out) :: strout  ! min len=110, dimension(21)

    ! local variables
    character(len=8)  :: date
    character(len=10) :: time
#ifndef REVISION
    character(len=*), parameter :: REVISION="unknown"
#endif

    ! Date and time
    call date_and_time(date, time)
    write(strout(1), '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
                      date(1:4), date(5:6), date(7:8), &
                      time(1:2), time(3:4), time(5:10)

    ! Git revision
    write(strout(2), '("# GRAY Git revision: ",a)') REVISION

    ! MHD equilibrium parameters
    if (params%equilibrium%iequil > 0) then
      write(strout(3), '("# EQL input: ",a)') trim(params%equilibrium%filenm(1:95))
      ! TODO add missing values
      write(strout(7), '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
        0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_
    else
      write(strout(3), '("# EQL input: N/A (vacuum)")')
      write(strout(7), '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
    end if

    write(strout(4), '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
      params%equilibrium%iequil, params%equilibrium%sgnb, &
      params%equilibrium%sgni, params%equilibrium%factb
    if (params%equilibrium%iequil > 1) then
      write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
        params%equilibrium%icocos, params%equilibrium%ipsinorm, &
        params%equilibrium%idesc, params%equilibrium%ifreefmt
      write(strout(6), '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
        params%equilibrium%ssplps, params%equilibrium%ssplf, &
        params%equilibrium%ixp
    else
      write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
      write(strout(6), '("# EQL ssplps ssplf ixp: N/A (analytical)")')
    end if

    ! Profiles parameters
    if (params%equilibrium%iequil > 0) then
      write(strout(8), '("# PRF input: ",a)') trim(params%profiles%filenm(1:95))
      write(strout(9), '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
        params%profiles%iprof, params%profiles%iscal, &
        params%profiles%factne, params%profiles%factte
      if (params%profiles%iprof > 0) then
        write(strout(10), '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
          params%profiles%irho, params%profiles%psnbnd, params%profiles%sspld
      else
        write(strout(10), '("# PRF irho psnbnd sspld: N/A (analytical)")')
      end if
      ! TODO: add missing values
      write(strout(11), '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') &
        0._wp_, 0._wp_, 0._wp_
    else
      write(strout(8), '("# PRF input: N/A (vacuum)")')
      write(strout(9), '("# PRF iprof iscal factne factte: N/A (vacuum)")')
      write(strout(10), '("# PRF irho psnbnd sspld: N/A (vacuum)")')
      write(strout(11), '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
    end if

    ! Antenna parameters
    write(strout(12), '("# ANT input: ",a)') trim(params%antenna%filenm(1:95))
    write(strout(13), '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
      params%antenna%ibeam, params%antenna%iox, &
      params%antenna%psi, params%antenna%chi
    write(strout(14), '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
      params%antenna%alpha, params%antenna%beta, &
      params%antenna%power, params%antenna%fghz
    write(strout(15), '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos
    write(strout(16), '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') &
      params%antenna%w, params%antenna%ri, params%antenna%phi

    ! Other parameters
    write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall

    ! code parameters
    write(strout(18), '("# COD igrad idst ipass ipol:",4(1x,i4))') &
      params%raytracing%igrad, params%raytracing%idst, &
      params%raytracing%ipass, params%raytracing%ipol
    write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') &
      params%raytracing%nrayr, params%raytracing%nrayth, &
      params%raytracing%nstep, params%raytracing%rwmax, &
      params%raytracing%dst
    write(strout(20), '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') &
      params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, &
      params%ecrh_cd%imx, params%ecrh_cd%ieccd
    write(strout(21), '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
      params%output%ipec, params%output%nrho, &
      params%output%istpr, params%output%istpl
  end subroutine print_parameters


  function update_parameter(params, name, value) result(err)
    ! Updates the value of a parameter, addressed by a string
    ! The return error is:
    !   ERR_SUCCESS on success;
    !   ERR_VALUE   on invalid parameter value;
    !   ERR_UNKNOWN on unknown parameter name.
    !
    ! Ex. update_parameter(params, 'raytracing.nrayr', '10')
    use ini_parser, only : ini_error, ERR_SUCCESS, ERR_VALUE, ERR_UNKNOWN

    implicit none

    ! function arguments
    type(gray_parameters), intent(inout) :: params
    character(*),          intent(in)    :: name, value
    integer(kind(ini_error)) :: err

    ! hope for the best...
    err = ERR_SUCCESS

    select case (name)
#include "gray_params.inc"
      case default
        ! unknown parameter
        err = ERR_UNKNOWN
    end select

  end function update_parameter


  subroutine read_gray_config(filename, params, err)
    ! Reads the GRAY parameters from the gray.ini configuration file
    use ini_parser, only : parse_ini, property_handler, ini_error

    implicit none

    ! subroutine arguments
    character(len=*),         intent(in)  :: filename
    type(gray_parameters),    intent(out) :: params
    integer(kind(ini_error)), intent(out) :: err

    call parse_ini(filename, ini_handler, err)

    contains

      function ini_handler(section, name, value) result(err)
        ! This function handles a single INI property and updates
        ! the `params` structure

        implicit none

        ! function arguments
        character(*), intent(in) :: section, name, value
        integer(kind(ini_error)) :: err

        err = update_parameter(params, section // "." // name, value)
      end function ini_handler

  end subroutine read_gray_config


  subroutine read_gray_params(filename, params, err)
    ! Reads the GRAY parameters from the legacy gray_params.data file
    use utils,  only : get_free_unit
    use logger, only : log_error

    implicit none

    ! subrouting arguments
    character(len=*),      intent(in)  :: filename
    type(gray_parameters), intent(out) :: params
    integer,               intent(out) :: err

    ! local variables
    integer :: u

    u = get_free_unit()

    open(u, file=filename, status='old', action='read', iostat=err)
    if (err /= 0) then
      call log_error('opening gray_params file ('//filename//') failed!', &
                     mod='gray_params', proc='read_gray_params')
      return
    end if

    read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, &
               params%raytracing%rwmax
    read(u, *) params%raytracing%igrad, params%raytracing%ipass, &
               params%raytracing%ipol
    read(u, *) params%raytracing%dst, params%raytracing%nstep, &
               params%raytracing%idst

    read(u, *) params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, params%ecrh_cd%imx
    read(u, *) params%ecrh_cd%ieccd

    read(u, *) params%antenna%alpha, params%antenna%beta
    read(u, *) params%antenna%power
    read(u, *) params%antenna%iox, params%antenna%psi, params%antenna%chi
    read(u, *) params%antenna%ibeam
    read(u, *) params%antenna%filenm

    read(u, *) params%equilibrium%iequil
    read(u, *) params%equilibrium%filenm
    read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, &
               params%equilibrium%idesc, params%equilibrium%ifreefmt
    read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps, &
               params%equilibrium%ssplf
    read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, &
               params%equilibrium%factb

    read(u, *) params%misc%rwall

    read(u, *) params%profiles%iprof, params%profiles%irho
    read(u, *) params%profiles%filenm
    read(u, *) params%profiles%psnbnd, params%profiles%sspld
    read(u, *) params%profiles%factte, params%profiles%factne, &
               params%profiles%iscal

    read(u, *) params%output%ipec, params%output%nrho
    read(u, *) params%output%istpr, params%output%istpl

    close(u)
  end subroutine read_gray_params


  subroutine set_globals(params)
    ! Set global variables exposed by this module.

    use logger, only : log_warning

    implicit none

    ! subroutine arguments
    type(gray_parameters), intent(in) :: params

    iequil = params%equilibrium%iequil
    iprof  = params%profiles%iprof

    ipec   = params%output%ipec
    nnd    = params%output%nrho
    istpr0 = params%output%istpr
    istpl0 = params%output%istpl

    ipol   = params%raytracing%ipol
    igrad  = params%raytracing%igrad
    idst   = params%raytracing%idst
    ipass  = params%raytracing%ipass

    if (params%raytracing%nrayr < 5) then
      igrad = 0
      call log_warning('nrayr < 5 ⇒ optical case only', &
                       mod="gray_params", proc="set_globals")
    end if

    iwarm = params%ecrh_cd%iwarm
    ilarm = params%ecrh_cd%ilarm
    imx   = params%ecrh_cd%imx
    ieccd = params%ecrh_cd%ieccd

  end subroutine set_globals

end module gray_params