gray/src/main.f90

program main
  use const_and_precisions, only : wp_, one, zero
  use logger,               only : INFO, ERROR, set_log_level, log_message
  use units,                only : set_active_units, close_units
  use utils,                only : dirname
  use gray_cli,             only : cli_options, parse_cli_options, &
                                   deinit_cli_options, parse_param_overrides
  use gray_core,            only : gray_main
  use gray_params,          only : gray_parameters, gray_data, gray_results, &
                                   read_gray_params, read_gray_config, &
                                   params_set_globals => set_globals

  implicit none

  ! CLI options
  type(cli_options) :: opts

  ! gray_main subroutine arguments
  type(gray_parameters) :: params  ! Inputs
  type(gray_data)       :: data    !
  type(gray_results)    :: results ! Outputs
  integer               :: err     ! Exit code

  ! Store the original working directory
  character(len=256) :: cwd
  call getcwd(cwd)

  ! Parse the command-line options
  call parse_cli_options(opts)

  ! Initialise logging
  if (opts%quiet) opts%verbose = ERROR
  call set_log_level(opts%verbose)

  ! Activate the given output units
  call set_active_units(opts%units)

  ! Load the parameters from file and move to its directory
  ! (all other filepaths are assumed relative to it)
  if (allocated(opts%config_file)) then
    ! Use the newer gray.ini configuration file
    call log_message(level=INFO, mod='main', msg='using GRAY INI config')
    call read_gray_config(opts%config_file, params)
    err = chdir(dirname(opts%config_file))
  else
    ! Use the legacy gray_params.data file
    call log_message(level=INFO, mod='main', msg='using GRAY params file')
    call read_gray_params(opts%params_file, params)
    err = chdir(dirname(opts%params_file))
  end if

  if (err /= 0) then
    call log_message(level=ERROR, mod='main', &
      msg='chdir to configuration file directory failed!')
    call exit(1)
  end if

  ! Apply CLI overrides and copy the parameters into
  ! global variables exported by the gray_params module
  call parse_param_overrides(params)
  call params_set_globals(params)

  ! Read the input data and set the global variables
  ! of the respective module. Note: order matters.
  call init_equilibrium(params, data)
  call init_antenna(params%antenna)
  call init_profiles(params%profiles, params%equilibrium%factb, &
                     params%antenna%pos, data%profiles)
  call init_misc(params, data)

  ! Change the current directory to output files there
  if (allocated(opts%output_dir)) then
    if (chdir(trim(cwd) // '/' // opts%output_dir) /= 0) then
      call log_message(level=ERROR, mod='main', &
        msg='chdir to output_dir ('//opts%output_dir//') failed!')
      call exit(1)
    end if
  end if

  if (allocated(opts%sum_filelist)) then
    call log_message(level=INFO, mod='main', msg='summing profiles')

    sum: block
      real(wp_), dimension(:), allocatable :: jphi
      real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
      real(wp_), dimension(:), allocatable :: extracol
      integer :: i, j, k, l, n, ngam, nextracol, irt
      character(len=255) :: filename, headerline, fmtstr
      real(wp_), dimension(5) :: f48v
      real(wp_) :: jphip,dpdvp, &
                   rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
      allocate(jphi(params%output%nrho), currins(params%output%nrho), &
               pins(params%output%nrho), rtin(params%output%nrho),    &
               rpin(params%output%nrho))

      open(100, file=opts%sum_filelist, action='read', status='old', iostat=err)
      if (err /= 0) then
        call log_message(level=ERROR, mod='main', &
          msg='opening file list ('//opts%sum_filelist//') failed!')
        call exit(1)
      end if

      open(unit=100 -1, file='f48sum.txt', action='write', status='unknown')
      open(unit=100 -2, file='f7sum.txt', action='write', status='unknown')

      read(100, *) n, ngam, nextracol
      allocate(extracol(nextracol))
      do i=1,n
        read(100, *) filename
        open(100 + i, file=filename, action='read', status='old', iostat=err)
        if (err /= 0) then
          call log_message(level=ERROR, mod='main', &
            msg='opening summand file ('//trim(filename)//') failed!')
          call exit(1)
        end if

        do j=1,22
          read(100 + i, '(a255)') headerline
          if (i == 1) then
            write(100 -1, '(a)') trim(headerline)
            write(100 -2, '(a)') trim(headerline)
          end if
        end do
      end do

      allocate(results%jcd(params%output%nrho), results%dpdv(params%output%nrho))
      do k=1,ngam
        jphi         = zero
        results%jcd  = zero
        results%dpdv = zero
        currins      = zero
        pins         = zero
        do j=1,params%output%nrho
          do i=1,n
            read(100+i, *) (extracol(l), l=1,nextracol),          &
                rpin(j), rtin(j), f48v(1:5), irt
            jphi(j)         = f48v(1) + jphi(j)
            results%jcd(j)  = f48v(2) + results%jcd(j)
            results%dpdv(j) = f48v(3) + results%dpdv(j)
            currins(j)      = f48v(4) + currins(j)
            pins(j)         = f48v(5) + pins(j)
          end do
          write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7
          write(100 -1,trim(fmtstr))                           &
              (extracol(l), l=1,nextracol), rpin(j), rtin(j),     &
              jphi(j), results%jcd(j), results%dpdv(j),           &
              currins(j), pins(j), irt
        end do
        results%pabs = pins(params%output%nrho)
        results%icd = currins(params%output%nrho)
        write(100 -1, *)
        call sum_profiles(params, jphi, results%jcd, results%dpdv,       &
                          currins, pins, results%pabs, results%icd,      &
                          jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
                          drhotjava, drhotpav, ratjamx, ratjbmx)
        write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12
        write(100 -2, trim(fmtstr))                                   &
             (extracol(l), l=1,nextracol), results%icd, results%pabs,    &
             jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav,              &
             drhotjava, drhotpav, ratjamx, ratjbmx
      end do
      do i=-2,n
        close(100 + i)
      end do
     deallocate(jphi, currins, pins, rtin, rpin)
     deallocate(extracol)
     deallocate(opts%params_file)
   end block sum
  else
    call gray_main(params, data, results, err)
  end if

  print_res: block
    character(256) :: msg
    write(msg, '(a,g0.3," MW")') 'total absoption: P=', results%pabs
    call log_message(msg, level=INFO, mod='main')
    write(msg, '(a,g0.3," kA")') 'total current drive: I=', results%icd * 1.0e3_wp_
    call log_message(msg, level=INFO, mod='main')
  end block print_res

  ! Free memory
  call deinit_equilibrium(data%equilibrium)
  call deinit_profiles(data%profiles)
  call deinit_misc
  call deinit_cli_options(opts)
  deallocate(results%dpdv, results%jcd)
  call close_units

contains

  subroutine init_equilibrium(params, data)
    ! Reads the MHD equilibrium file (either in the G-EQDSK format
    ! or an analytical description) and initialises the respective
    ! GRAY parameters and data.
    use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, &
                            set_equil_an, set_equil_spline, eq_scal

    implicit none

    ! subroutine arguments
    type(gray_parameters), intent(inout) :: params
    type(gray_data),       intent(out)   :: data

    if (params%equilibrium%iequil < 2) then
      ! Analytical equilibrium
      call read_equil_an(params%equilibrium%filenm, &
                         params%raytracing%ipass,   &
                         data%equilibrium)

      ! Set psia sign to give the correct sign to Iphi
      ! (COCOS=3: psia<0 for Iphi>0)
      data%equilibrium%psia = sign(one, data%equilibrium%qpsi(2) &
                                      * data%equilibrium%fpol(1))
    else
      ! Numerical equilibrium
      call read_eqdsk(params%equilibrium, data%equilibrium)
      call change_cocos(data%equilibrium, params%equilibrium%icocos, 3)
    end if

    ! Rescale B, I and/or force their signs
    call eq_scal(params%equilibrium, data%equilibrium)

    ! Set global variables (for splines)
    if (params%equilibrium%iequil < 2) then
      call set_equil_an(data%equilibrium)
    else
      call set_equil_spline(params%equilibrium, data%equilibrium)
    end if
  end subroutine init_equilibrium


  subroutine deinit_equilibrium(data)
    ! Free all memory allocated by the init_equilibrium subroutine.
    use gray_params, only : equilibrium_data
    use equilibrium, only : unset_equil_spline, unset_rho_spline, unset_q

    implicit none

    ! subroutine arguments
    type(equilibrium_data), intent(inout) :: data

    ! Free the MHD equilibrium arrays
    if (allocated(data%rv))   deallocate(data%rv, data%zv, data%fpol, data%qpsi)
    if (allocated(data%psin)) deallocate(data%psin, data%psinr)
    if (allocated(data%rbnd)) deallocate(data%rbnd, data%zbnd)
    if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)

    ! Unset global variables of the `equilibrium` module
    call unset_equil_spline
    call unset_rho_spline
    call unset_q
  end subroutine deinit_equilibrium


  subroutine init_profiles(params, factb, launch_pos, data)
    ! Reads the plasma kinetic profiles file (containing the elecron
    ! temperature, density and plasma effective charge) and initialises
    ! the respective GRAY data structure.
    use gray_params,  only : profiles_parameters, profiles_data
    use equilibrium,  only : frhopolv
    use coreprofiles, only : read_profiles_an, read_profiles, &
                             scale_profiles, set_profiles_an, &
                             set_profiles_spline
    implicit none

    ! subroutine arguments
    type(profiles_parameters), intent(inout) :: params
    real(wp_),                 intent(in)    :: factb
    real(wp_),                 intent(in)    :: launch_pos(3)
    type(profiles_data),       intent(out)   :: data

    if (params%iprof == 0) then
      ! Analytical profiles
      call read_profiles_an(params%filenm, data)
    else
      ! Numerical profiles
      call read_profiles(params%filenm, data)

      ! Convert psrad to ψ
      select case (params%irho)
        case (0) ! psrad is ρ_t = √Φ (toroidal flux)
          data%psrad = frhopolv(data%psrad)**2
        case (1) ! psrad is ρ_p = √ψ (poloidal flux)
          data%psrad = data%psrad**2
        case default ! psrad is already ψ
      end select
    end if

    ! Rescale input data
    call scale_profiles(params, factb, data)

    ! Set global variables
    if (params%iprof == 0) then
      ! Analytical profiles
      call set_profiles_an(params, data)
    else
      ! Numerical profiles
      call set_profiles_spline(params, data, launch_pos)
    end if
  end subroutine init_profiles


  subroutine deinit_profiles(data)
    ! Free all memory allocated by the init_profiles subroutine.
    use gray_params,  only : profiles_data
    use coreprofiles, only : unset_profiles_spline

    implicit none

    ! subroutine arguments
    type(profiles_data), intent(inout) :: data

    ! Free the plasma kinetic profiles arrays
    if (allocated(data%psrad)) deallocate(data%psrad)
    if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)

    ! Unset global variables of the `coreprofiles` module
    call unset_profiles_spline
  end subroutine deinit_profiles


  subroutine init_antenna(params)
    ! Reads the wave launcher file (containing the wave frequency, launcher
    ! position, direction and beam description) and initialises the respective
    ! GRAY parameters.
    use beams,       only : read_beam0, read_beam1, read_beam2
    use gray_params, only : antenna_parameters

    implicit none

    ! subroutine arguments
    type(antenna_parameters), intent(inout) :: params

    ! Note: α, β are loaded from gray_params.data
    select case (params%ibeam)
      case (2)
        ! 2 degrees of freedom
        ! w(z, α, β), 1/R(z, α, β)
        ! FIXME: 1st beam is always selected, iox read from table
        call read_beam2(params, beamid=1)
      case (1)
        ! 1 degree of freedom
        ! w(z, α), 1/R(z, α)
        call read_beam1(params)
      case default
        ! fixed w(z), 1/R(z)
        call read_beam0(params)
    end select
  end subroutine init_antenna


  subroutine init_misc(params, data)
    ! Performs miscellanous initial tasks, before the gray_main subroutine.
    use reflections, only : range2rect
    use limiter,     only : limiter_set_globals=>set_globals

    implicit none

    ! subroutine arguments
    type(gray_parameters), intent(inout) :: params
    type(gray_data),       intent(inout) :: data

    ! Build a basic limiter when one is not provided by the EQDSK
    if (.not. allocated(data%equilibrium%rlim) &
        .or. params%raytracing%ipass < 0) then
      block
        real(wp_) :: rmxm, r0m, z0m, dzmx
        r0m  = norm2(params%antenna%pos(1:2)) * 0.01_wp_
        dzmx = abs(params%raytracing%ipass) * &
               params%raytracing%dst * params%raytracing%nstep * 0.01_wp_
        z0m  = params%antenna%pos(3) * 0.01_wp_

        allocate(data%equilibrium%rlim(5))
        allocate(data%equilibrium%zlim(5))
        params%raytracing%ipass = abs(params%raytracing%ipass)
        if (params%equilibrium%iequil < 2) then
          rmxm = (data%equilibrium%rv(1) + data%equilibrium%rv(2)) * 0.01_wp_
        else
          rmxm = data%equilibrium%rv(size(data%equilibrium%rv))
        end if
        call range2rect(params%misc%rwall, max(r0m, rmxm), &
                        z0m - dzmx, z0m + dzmx,            &
                        data%equilibrium%rlim, data%equilibrium%zlim)
      end block
    end if

    ! Set the global variables of the `limiter` module
    call limiter_set_globals(data%equilibrium)
  end subroutine init_misc


  subroutine deinit_misc
    ! Free all memory allocated by the init_misc subroutine.
    use limiter, only : limiter_unset_globals=>unset_globals

    implicit none

    ! Unset the global variables of the `limiter` module
    call limiter_unset_globals
  end subroutine deinit_misc


  subroutine sum_profiles(params, jphi, jcd, dpdv, currins, pins, pabs, icd, &
                          jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav,     &
                          drhotjava, drhotpav, ratjamx, ratjbmx)
    use const_and_precisions, only : zero, degree
    use coreprofiles,         only : set_profiles_an, set_profiles_spline, &
                                     temp, fzeff
    use dispersion,           only : expinit
    use gray_params,          only : gray_parameters, print_parameters
    use beams,                only : launchangles2n, xgygcoeff
    use magsurf_data,         only : flux_average, dealloc_surfvec
    use beamdata,             only : init_btr, dealloc_beam
    use pec,                  only : pec_init, postproc_profiles, dealloc_pec, &
                                     rhop_tab, rhot_tab
    use gray_core,            only : print_headers, print_finals, print_pec, &
                                     print_bres, print_prof, print_maps, &
                                     print_surfq
    implicit none

    ! subroutine arguments
    type(gray_parameters),   intent(inout) :: params
    real(wp_),               intent(in)    :: pabs, icd
    real(wp_), dimension(:), intent(in)    :: jphi, jcd, dpdv, currins, pins
    real(wp_), intent(out) :: jphip, dpdvp, rhotj, rhotjava,       &
                              rhotp, rhotpav, drhotjava, drhotpav, &
                              ratjamx, ratjbmx

    ! local variables
    real(wp_) :: ak0, bres, xgcn
    real(wp_) :: chipol, psipol, st
    real(wp_) :: drhotp, drhotj, dpdvmx, jphimx

    real(wp_), dimension(3) :: anv0
    real(wp_), dimension(:, :),    pointer :: yw=>null(), ypw=>null(), gri=>null()
    real(wp_), dimension(:, :, :), pointer :: xc=>null(), du1=>null(), ggri=>null()

    real(wp_), dimension(:, :), pointer :: psjki=>null(), ppabs=>null(), ccci=>null()
    real(wp_), dimension(:),    pointer :: tau0=>null(), alphaabs0=>null(), &
                                           dids0=>null(), ccci0=>null()
    real(wp_),    dimension(:), pointer :: p0jk=>null()
    complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
    integer,      dimension(:), pointer :: iiv=>null()

    ! ======== set environment BEGIN ========
    ! Compute X=(ω_pe/ω)² and Y=ω_ce/ω (with B=1)
    call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)

    ! Compute the initial cartesian wavevector (anv0)
    call launchangles2n(params%antenna, anv0)

    ! Initialise the ray variables (beamtracing)
    call init_btr(params%raytracing, yw, ypw, xc, du1, &
                  gri, ggri, psjki, ppabs, ccci,       &
                  tau0, alphaabs0, dids0, ccci0,       &
                  p0jk, ext, eyt, iiv)

    ! Initialise the dispersion module
    if (params%ecrh_cd%iwarm > 1) call expinit

    ! Initialise the magsurf_data module
    call flux_average ! requires frhotor for dadrhot,dvdrhot

    ! Initialise the output profiles
    call pec_init(params%output%ipec)
    ! ======= set environment  END  ======

    ! ======== pre-proc prints BEGIN ========
    call print_headers(params)

    ! Print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
    call print_surfq([1.5_wp_, 2.0_wp_])

    ! Print ne, Te, q, Jphi versus psi, rhop, rhot
    call print_bres(bres)
    call print_prof
    call print_maps(bres, xgcn,                                &
                    norm2(params%antenna%pos(1:2)) *0.01_wp_ , &
                    sin(params%antenna%beta*degree))
    ! ========= pre-proc prints END =========

    ! Print power and current density profiles
    call print_pec(rhop_tab, rhot_tab, jphi, jcd, &
                   dpdv, currins, pins, index_rt=1)
    ! Compute profiles width
    call postproc_profiles(pabs, icd, rhot_tab, dpdv, jphi,            &
                           rhotpav, drhotpav, rhotjava, drhotjava,     &
                           dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, &
                           dpdvmx, jphimx, ratjamx, ratjbmx)
    ! Print 0D results
    call print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, &
                      drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp,      &
                      drhotj, ratjamx, ratjbmx, st, psipol, chipol,         &
                      1, params%antenna%power, cpl1=zero, cpl2=zero)

    ! Free memory
    call dealloc_surfvec ! for fluxval
    call dealloc_beam(yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci,  &
                      tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
    call dealloc_pec
  end subroutine sum_profiles

end program main