!> Create a mapping between GTC internal varibles and their physically
!> meaningful names, full description, and normalization constants to CGS unit
!>
!> The Convention is:
!>     Every variable in this module is named as
!>     {subroutine_name}_{internal_variable_name}_[name, dscr or norm]
!>     for example, a GTC variable "pdf" in the subroutine "snapshot" will have
!>     the following quantities associated:
!>     CHARACTER(20) snapshot_pdf_name
!>     CHARACTER(1000) snapshot_pdf_dscr
!>     real(lk) snapshot_pdf_norm

!>     The types of variables in this module will be CHARACTER(len=20) for
!>     names, CHARACTER(len=1000) for dscriptions, and real(lk) for 
!>     normalization constants.

!>     In some cases, an internal array is used to contain multiple physical
!>     quantities in GTC, corresponding naming and normalization arrays will
!>     also be dimensional and the indexes determines the correspondance.

!> \author Lei Shi
!> \date 05/24/2018
module physical_meaning

  use precision, only: lk
  use utility, only: check_allocation, mstat

  implicit none
  public

  ! number of fields written out in snapshot
  integer, parameter :: snapshot_nfield=3

  ! Variables in Snapshot
  CHARACTER(len=20), dimension(:,:), allocatable :: snapshot_profile_name
  CHARACTER(len=1000), dimension(:,:), allocatable :: snapshot_profile_dscr
  real(lk), dimension(:,:), allocatable :: snapshot_profile_norm

  CHARACTER(len=20), dimension(:,:), allocatable :: snapshot_pdf_name
  CHARACTER(len=1000), dimension(:,:), allocatable :: snapshot_pdf_dscr
  real(lk), dimension(:,:), allocatable :: snapshot_pdf_norm

  CHARACTER(len=20) :: snapshot_poloidata_name(snapshot_nfield+2)
  CHARACTER(len=1000) :: snapshot_poloidata_dscr(snapshot_nfield+2)
  real(lk) :: snapshot_poloidata_norm(snapshot_nfield+2)

  CHARACTER(len=20) :: snapshot_fluxdata_name(snapshot_nfield)
  CHARACTER(len=1000) :: snapshot_fluxdata_dscr(snapshot_nfield)
  real(lk) :: snapshot_fluxdata_norm(snapshot_nfield)

  save
  contains
    subroutine create_name
      call dynamic_name
    end subroutine create_name

    subroutine static_name
    end subroutine static_name

    subroutine dynamic_name
      call snapshot_dynamic_name
    end subroutine dynamic_name

    subroutine snapshot_dynamic_name
      use global_parameters, only: nspecies, r0, nhybrid, magnetic, cs_method
      use particle_array, only: iload, fload, feload
      integer :: ns

      allocate(snapshot_profile_name(6,nspecies), &
               snapshot_profile_dscr(6,nspecies), &
               snapshot_profile_norm(6,nspecies), &
               snapshot_pdf_name(4,nspecies), &
               snapshot_pdf_dscr(4,nspecies), &
               snapshot_pdf_norm(4,nspecies), STAT=mstat)
      call check_allocation(mstat, "snapshot name creation")

      ! assign names according to GTC particle species
      ns = 1 ! First entry is always ion
      if(iload == 1 .or. iload == 2 .or. iload == 3) then
        snapshot_profile_name(1, ns) = "ni_bar"
        snapshot_profile_dscr(1, ns) = "total ion marker weight normalized to"&
          //"initial local ion marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(1, ns) = 1.0_lk

        snapshot_profile_name(2, ns) = "dni_bar"
        snapshot_profile_dscr(2, ns) = "perturbed ion marker weight normalized to"&
          //"total ion marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(2, ns) = 1.0_lk

        snapshot_profile_name(3, ns) = "upara_i_bar"
        snapshot_profile_dscr(3, ns) = "flux-averaged ion marker parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal ion velocity"
        snapshot_profile_norm(3, ns) = 1.0_lk

        snapshot_profile_name(4, ns) = "dupara_i_bar"
        snapshot_profile_dscr(4, ns) = "flux-averaged ion perturbed parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal ion velocity"
        snapshot_profile_norm(4, ns) = 1.0_lk

        snapshot_profile_name(5, ns) = "Ei_bar"
        snapshot_profile_dscr(5, ns) = "flux-averaged ion full energy;"&
          //"1D function of psi; unitless, normalized to on-axis electron temperature Te_a"
        snapshot_profile_norm(5, ns) = 1.0_lk

        snapshot_profile_name(6, ns) = "dEi_bar"
        snapshot_profile_dscr(6, ns) = "flux-averaged ion perturbed energy;"&
          //"1D function of psi; unitless, normalized to on-axis electron temperature Te_a"
        snapshot_profile_norm(6, ns) = 1.0_lk

        snapshot_pdf_name(1, ns) = "f_E_ion"
        snapshot_pdf_dscr(1, ns) = "volume integrated ion energy distribution; "&
          //"1D function of energy; unit: super particle number"
        snapshot_pdf_norm(1, ns) = 1.0_lk

        snapshot_pdf_name(2, ns) = "df_E_ion_bar"
        snapshot_pdf_dscr(2, ns) = "volume integrated ion perturbed energy distribution; "&
          //"1D function of energy; unitless, normalized to super particle energy distribution"
        snapshot_pdf_norm(2, ns) = 1.0_lk

        snapshot_pdf_name(3, ns) = "f_lam_ion_bar"
        snapshot_pdf_dscr(3, ns) = "volume integrated ion pitch angle distribution; "&
          //"1D function of pitch angle; unit: super particle number"
        snapshot_pdf_norm(3, ns) = 1.0_lk

        snapshot_pdf_name(4, ns) = "df_lam_ion_bar"
        snapshot_pdf_dscr(4, ns) = "volume integrated ion perturbed pitch angle distribution; "&
          //"1D function of pitch angle; unitless, normalized to super particle pitch angle distribution"
        snapshot_pdf_norm(4, ns) = 1.0_lk
      else ! inconventional run without normal kinetic ions, NOT IMPLEMENTED
        snapshot_profile_name(1, ns) = "unknown_profile_1_1"
        snapshot_profile_name(2, ns) = "unknown_profile_1_2"
        snapshot_profile_name(3, ns) = "unknown_profile_1_3"
        snapshot_profile_name(4, ns) = "unknown_profile_1_4"
        snapshot_profile_name(5, ns) = "unknown_profile_1_5"
        snapshot_profile_name(6, ns) = "unknown_profile_1_6"
        snapshot_profile_dscr(:, ns) = "UNKNOWN"
        snapshot_profile_norm(:, ns) = 1.0_lk

        snapshot_pdf_name(1, ns) = "unknown_pdf_1_1"
        snapshot_pdf_name(2, ns) = "unknown_pdf_1_2"
        snapshot_pdf_name(3, ns) = "unknown_pdf_1_3"
        snapshot_pdf_name(4, ns) = "unknown_pdf_1_4"
        snapshot_pdf_dscr(:, ns) = "UNKNOWN"
        snapshot_pdf_norm(:, ns) = 1.0_lk
      endif

      ! Now take care of the other species

      ! First, if kinetic electrons exist
      if(nhybrid > 0) then
        ns = ns+1
        
        snapshot_profile_name(1, ns) = "ne_bar"
        snapshot_profile_dscr(1, ns) = "total electron marker weight normalized to"&
          //"initial local electron marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(1, ns) = 1.0_lk

        snapshot_profile_name(2, ns) = "dne_bar"
        snapshot_profile_dscr(2, ns) = "perturbed electron marker weight normalized to"&
          //"total electron marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(2, ns) = 1.0_lk

        snapshot_profile_name(3, ns) = "upara_e_bar"
        snapshot_profile_dscr(3, ns) = "flux-averaged electron marker parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal electron velocity"
        snapshot_profile_norm(3, ns) = 1.0_lk

        snapshot_profile_name(4, ns) = "dupara_e_bar"
        snapshot_profile_dscr(4, ns) = "flux-averaged electron perturbed parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal electron velocity"
        snapshot_profile_norm(4, ns) = 1.0_lk

        snapshot_profile_name(5, ns) = "Ee_bar"
        snapshot_profile_dscr(5, ns) = "flux-averaged electron full energy;"&
          //"1D function of psi; unitless, normalized to on-axis electron temperature Te_a"
        snapshot_profile_norm(5, ns) = 1.0_lk

        snapshot_profile_name(6, ns) = "dEe_bar"
        snapshot_profile_dscr(6, ns) = "flux-averaged electron perturbed energy;"&
          //"1D function of psi; unitless, normalized to on-axis electron temperature Te_a"
        snapshot_profile_norm(6, ns) = 1.0_lk

        snapshot_pdf_name(1, ns) = "f_E_electron"
        snapshot_pdf_dscr(1, ns) = "volume integrated electron energy distribution; "&
          //"1D function of energy; unit: super particle number"
        snapshot_pdf_norm(1, ns) = 1.0_lk

        snapshot_pdf_name(2, ns) = "df_E_electron_bar"
        snapshot_pdf_dscr(2, ns) = "volume integrated electron perturbed energy distribution; "&
          //"1D function of energy; unitless, normalized to super particle energy distribution"
        snapshot_pdf_norm(2, ns) = 1.0_lk

        snapshot_pdf_name(3, ns) = "f_lam_electron_bar"
        snapshot_pdf_dscr(3, ns) = "volume integrated electron pitch angle distribution; "&
          //"1D function of pitch angle; unit: super particle number"
        snapshot_pdf_norm(3, ns) = 1.0_lk

        snapshot_pdf_name(4, ns) = "df_lam_electron_bar"
        snapshot_pdf_dscr(4, ns) = "volume integrated electron perturbed pitch angle distribution; "&
          //"1D function of pitch angle; unitless, normalized to super particle pitch angle distribution"
        snapshot_pdf_norm(4, ns) = 1.0_lk
      endif

      ! second, if kinetic fast ions exist
      if(fload > 0) then
        ns = ns+1
        
        snapshot_profile_name(1, ns) = "nf_bar"
        snapshot_profile_dscr(1, ns) = "total fast ion marker weight normalized to"&
          //"initial local fast ion marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(1, ns) = 1.0_lk

        snapshot_profile_name(2, ns) = "dnf_bar"
        snapshot_profile_dscr(2, ns) = "perturbed fast ion marker weight normalized to"&
          //"total fast ion marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(2, ns) = 1.0_lk

        snapshot_profile_name(3, ns) = "upara_f_bar"
        snapshot_profile_dscr(3, ns) = "flux-averaged fast ion marker parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal fast ion velocity"
        snapshot_profile_norm(3, ns) = 1.0_lk

        snapshot_profile_name(4, ns) = "dupara_f_bar"
        snapshot_profile_dscr(4, ns) = "flux-averaged fast ion perturbed parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal fast ion velocity"
        snapshot_profile_norm(4, ns) = 1.0_lk

        snapshot_profile_name(5, ns) = "Ef_bar"
        snapshot_profile_dscr(5, ns) = "flux-averaged fast ion full energy;"&
          //"1D function of psi; unitless, normalized to on-axis fast ion temperature Te_a"
        snapshot_profile_norm(5, ns) = 1.0_lk

        snapshot_profile_name(6, ns) = "dEf_bar"
        snapshot_profile_dscr(6, ns) = "flux-averaged fast ion perturbed energy;"&
          //"1D function of psi; unitless, normalized to on-axis fast ion temperature Te_a"
        snapshot_profile_norm(6, ns) = 1.0_lk

        snapshot_pdf_name(1, ns) = "f_E_fast"
        snapshot_pdf_dscr(1, ns) = "volume integrated fast ion energy distribution; "&
          //"1D function of energy; unit: super particle number"
        snapshot_pdf_norm(1, ns) = 1.0_lk

        snapshot_pdf_name(2, ns) = "df_E_fast_bar"
        snapshot_pdf_dscr(2, ns) = "volume integrated fast ion perturbed energy distribution; "&
          //"1D function of energy; unitless, normalized to super particle energy distribution"
        snapshot_pdf_norm(2, ns) = 1.0_lk

        snapshot_pdf_name(3, ns) = "f_lam_fast_bar"
        snapshot_pdf_dscr(3, ns) = "volume integrated fast ion pitch angle distribution; "&
          //"1D function of pitch angle; unit: super particle number"
        snapshot_pdf_norm(3, ns) = 1.0_lk

        snapshot_pdf_name(4, ns) = "df_lam_fast_bar"
        snapshot_pdf_dscr(4, ns) = "volume integrated fast ion perturbed pitch angle distribution; "&
          //"1D function of pitch angle; unitless, normalized to super particle pitch angle distribution"
        snapshot_pdf_norm(4, ns) = 1.0_lk
      endif

      ! finally, if kinetic fast electrons exist
      if(feload > 0) then
        ns = ns+1
        
        snapshot_profile_name(1, ns) = "nfe_bar"
        snapshot_profile_dscr(1, ns) = "total fast electron marker weight normalized to"&
          //"initial local fast electron marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(1, ns) = 1.0_lk

        snapshot_profile_name(2, ns) = "dnfe_bar"
        snapshot_profile_dscr(2, ns) = "perturbed fast electron marker weight normalized to"&
          //"total fast electron marker weight; 1D function of psi; unitless"
        snapshot_profile_norm(2, ns) = 1.0_lk

        snapshot_profile_name(3, ns) = "upara_fe_bar"
        snapshot_profile_dscr(3, ns) = "flux-averaged fast electron marker parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal fast electron velocity"
        snapshot_profile_norm(3, ns) = 1.0_lk

        snapshot_profile_name(4, ns) = "dupara_fe_bar"
        snapshot_profile_dscr(4, ns) = "flux-averaged fast electron perturbed parallel velocity;"&
          //"1D function of psi; unitless, normalized to on-axis thermal fast electron velocity"
        snapshot_profile_norm(4, ns) = 1.0_lk

        snapshot_profile_name(5, ns) = "Efe_bar"
        snapshot_profile_dscr(5, ns) = "flux-averaged fast electron full energy;"&
          //"1D function of psi; unitless, normalized to on-axis fast electron temperature Te_a"
        snapshot_profile_norm(5, ns) = 1.0_lk

        snapshot_profile_name(6, ns) = "dEfe_bar"
        snapshot_profile_dscr(6, ns) = "flux-averaged fast electron perturbed energy;"&
          //"1D function of psi; unitless, normalized to on-axis fast electron temperature Te_a"
        snapshot_profile_norm(6, ns) = 1.0_lk

        snapshot_pdf_name(1, ns) = "f_E_fast_electron"
        snapshot_pdf_dscr(1, ns) = "volume integrated fast electron energy distribution; "&
          //"1D function of energy; unit: super particle number"
        snapshot_pdf_norm(1, ns) = 1.0_lk

        snapshot_pdf_name(2, ns) = "df_E_fast_electron_bar"
        snapshot_pdf_dscr(2, ns) = "volume integrated fast electron perturbed energy distribution; "&
          //"1D function of energy; unitless, normalized to super particle energy distribution"
        snapshot_pdf_norm(2, ns) = 1.0_lk

        snapshot_pdf_name(3, ns) = "f_lam_fast_electron_bar"
        snapshot_pdf_dscr(3, ns) = "volume integrated fast electron pitch angle distribution; "&
          //"1D function of pitch angle; unit: super particle number"
        snapshot_pdf_norm(3, ns) = 1.0_lk

        snapshot_pdf_name(4, ns) = "df_lam_fast_electron_bar"
        snapshot_pdf_dscr(4, ns) = "volume integrated fast electron perturbed pitch angle distribution; "&
          //"1D function of pitch angle; unitless, normalized to super particle pitch angle distribution"
        snapshot_pdf_norm(4, ns) = 1.0_lk
      endif

      ! Now the field related names
      ! Field names are currently relatively static, only dynamic part is the
      ! third field, in normal hybrid model it's total perturbed electron density,
      ! while in CS model, it's the non-adiabatic A_para

      ! First the static parts
      ! 1st entry is always electric potential
      snapshot_fluxdata_name(1) = "phi_flux"
      snapshot_fluxdata_dscr(1) = "Electrostatic potential on diag_flux, (ntheta, mtoroidal) in Fortran, with duplicated grid in theta. unit: GTC unit"
      ! TODO: Normalization factor to CGS is currently not implemented
      snapshot_fluxdata_norm(1) = -1.0_lk

      snapshot_poloidata_name(1) = "phi_poloidal"
      snapshot_poloidata_dscr(1) = "Electrostatic potential on zeta=0 plane, (ntheta, mpsi+1) in Fortran, with duplicated grid in theta. unitless, normalized to Te on axis"
      snapshot_poloidata_norm(1) = 1.0_lk

      ! second entry depends on magnetic feature on/off
      if(magnetic == 0) then
        ! electrostatic run, no magnetic quantities
        snapshot_fluxdata_name(2) = "UNDEFINED_FLUX_2"
        snapshot_fluxdata_dscr(2) = "UNDEFINED FLUX"
        snapshot_fluxdata_norm(2) = -1.0_lk

        snapshot_fluxdata_name(3) = "UNDEFINED_FLUX_3"
        snapshot_fluxdata_dscr(3) = "UNDEFINED FLUX"
        snapshot_fluxdata_norm(3) = -1.0_lk

        snapshot_poloidata_name(2) = "UNDEFINED_POLOID_2"
        snapshot_poloidata_dscr(2) = "UNDEFINED POLOID"
        snapshot_poloidata_norm(2) = -1.0_lk

        snapshot_poloidata_name(3) = "UNDEFINED_POLOID_3"
        snapshot_poloidata_dscr(3) = "UNDEFINED POLOID"
        snapshot_poloidata_norm(3) = -1.0_lk
      else
        ! electromagnetic run has these data
        snapshot_fluxdata_name(2) = "apara_flux"
        snapshot_fluxdata_dscr(2) = "Parallel vector potential on diag_flux, (ntheta, mtoroidal) in Fortran, with duplicated grid in theta. unit: GTC unit"
        !TODO: Apara normalization factor not implemented
        snapshot_fluxdata_norm(2) = -1.0_lk

        snapshot_poloidata_name(2) = "apara_poloidal"
        snapshot_poloidata_dscr(2) = "Parallel vector potential on zeta=0 plane, (ntheta, mpsi+1) in Fortran, with duplicated grid in theta. unitless, normalized to proper magnitude comparable to phi_bar in ideal shear Alfven wave"
        snapshot_poloidata_norm(2) = 1.0_lk

        if (cs_method == 0) then
          ! normal hybrid case, nfield=3 corresponds to electron density
          snapshot_fluxdata_name(3) = "dne_flux"
          snapshot_fluxdata_dscr(3) = "Perturbed electron density on diag_flux, (ntheta, mtoroidal) in Fortran, with duplicated grid in theta. Unitless, normalized to local electron equilibrium density "
          !TODO: dne_flux normalization factor not implemented
          snapshot_fluxdata_norm(3) = -1.0_lk

          snapshot_poloidata_name(3) = "dne_poloidal"
          snapshot_poloidata_dscr(3) = "Perturbed electron density on zeta=0 plane, (ntheta, mpsi+1) in Fortran, with duplicated grid in theta. Unitless, normalized to local electron equilibrium density"
          !TODO: dne_poloidal normalization factor not implemented
          snapshot_poloidata_norm(3) = -1.0_lk
        else
          ! CS model, nfield=3 corresponds to Apara_NA
          snapshot_fluxdata_name(3) = "apara_na_flux"
          snapshot_fluxdata_dscr(3) = "Non-adiabatic Apara on diag_flux, (ntheta, mtoroidal) in Fortran, with duplicated grid in theta. Unit: GTC unit "
          !TODO: apara_na_flux normalization factor not implemented
          snapshot_fluxdata_norm(3) = -1.0_lk

          snapshot_poloidata_name(3) = "apara_na_poloidal"
          snapshot_poloidata_dscr(3) = "Non-adiabatic Apara on zeta=0 plane, (ntheta, mpsi+1) in Fortran, with duplicated grid in theta. Unitless, normalized to proper magnitude comparable to phi_bar in ideal shear Alfven wave"
          !TODO: apara_na_poloidal normalization factor not implemented
          snapshot_poloidata_norm(3) = -1.0_lk
        endif
      endif

      !Finally, the coordinates for poloidata
      snapshot_poloidata_name(snapshot_nfield + 1) = "R_poloidal"
      snapshot_poloidata_dscr(snapshot_nfield + 1) = "Major Radius coordinate corresponds to poloidal fields mesh points. Unit: GTC unit"
      snapshot_poloidata_norm(snapshot_nfield + 1) = r0

      snapshot_poloidata_name(snapshot_nfield + 2) = "Z_poloidal"
      snapshot_poloidata_dscr(snapshot_nfield + 2) = "Vertical coordinate corresponds to poloidal fields mesh points. Unit: GTC unit"
      snapshot_poloidata_norm(snapshot_nfield + 2) = r0
    endsubroutine snapshot_dynamic_name

end module physical_meaning
