!     $Id$
!***************************************************************
!
!   Set constants and initial parameters
!
!***************************************************************

SUBROUTINE TXINIT
  use tx_commons

  implicit none

  !   ***** Configuration parameters *****

  !   Plasma minor radius (m)
  RA = 0.35D0

  !   Wall radius (m)
  RB = 0.38D0

  !   Mesh accumulation radius (m)
  !     RACCUM coincides with RA when RACCUM is minus. (-1: default)
  !     RACCUM is valid when RACCUM is plus.
  RACCUM = - 1.d0

  !   Plasma major radius (m)
  RR = 1.3D0

  !   Toroidal magnetic field (T)
  BB = 1.3D0

  !   Plasma current start (MA)
  rIPs= 0.15D0

  !   Plasma current end (MA)
  rIPe= 0.15D0

  !   ***** Plasma components *****

  !   Atomic number of ion
  PA = 1.D0

  !   Charge number of ion
  PZ = 1.D0

  !   Effective charge
  Zeff = 2.D0

  !   ***** Initial plasma parameters *****

  !   Initial electron density at r = 0 (10^20 m^-3)
  PN0 = 0.4D0

  !   Initial electron density at r = a (10^20 m^-3)
  PNa = 0.05D0

  !   Electron density in diverter region (Minimum density in SOL)
  PNeDIV = 10.D-3

  !   Initial electron temperature at r = 0 (keV)
  PTe0 = 700.D-3

  !   Initial electron temperature at r = a (keV)
  PTea =  50.D-3

  !   Electron temperature in diverter region (Minimum Te in SOL)
  PTeDIV = 10.D-3

  !   Initial ion temperature  at r = 0 (keV)
  PTi0 = 700.D-3

  !   Initial ion temperature  at r = a (keV)
  PTia =  50.D-3

  !   Ion temperature in diverter region (Minimum Ti in SOL)
  PTiDIV = 10.D-3

  !   Initial current profile parameter
  PROFJ = 2.D0

  !   Initial density profile parameters
  PROFN1 = 3.D0
  PROFN2 = 1.D0

  !   Initial temperature profile parameters
  PROFT1 = 2.D0
  PROFT2 = 2.D0

  !   ***** Particle diffusivity and viscosity parameters *****

  !   Turbulent pinch velocity parameter
  VWpch0 = 0.D0

  !   Electron-driven diffusion parameter
  De0 = 0.05D0

  !   Ion-driven diffusion parameter
  Di0 = 0.D0

  !   Electron viscosity parameter
  rMue0 = 0.3D0

  !   Ion viscosity parameter
  rMui0 = 0.3D0

  !   Drift frequency parameter (omega/omega*e)
  WPM0 = 0.D0

  !  ***** Thermal diffusivity parameters *****

  !   Electron thermal diffusivity parameter (Chie/D)
  !     0 for fixed temperature profile
  Chie0 = 0.3D0

  !   Ion thermal diffusivity parameter (Chie/D)
  !     0 for fixed temperature profile
  Chii0 = 0.3D0

  !   ***** Turbulent transport control parameters *****

  !   Fixed transport coefficient parameter
  !     1 : particle transport
  !     2 : momentum transport
  !     3 : thermal transport
  FSDFIX(1:3) = 1.D0

  !   Anomalous turbulent transport models
  !     1 : CDBM model (Current diffusive ballooning mode)
  !     2 : CDIM model (Current diffusive interchange mode)
  !     3 : MMM95 (Multi-mode model)
  !   If MDANOM is negative, smoothing of the pressure gradient in the direction of time is off.
  MDANOM = 1

  !   Switch and amplification of turbulence coefficients (typically 0 or 1)
  !     1 : particle transport
  !     2 : momentum transport
  !     3 : thermal transport
  FSANOM(1:3) = 0.D0

  !   Effect of ExB shear stabilization
  FSCBSH = 0.D0

  !   Position of ETB shoulder (valid only when MDLETB /= 0)
  !     1 : particle transport
  !     2 : momentum transport
  !     3 : thermal transport
  RhoETB(1:3) = 0.9D0

  !   ==== CDBM transport coefficient parameters ============================
  !      The finer time step size, typically less than or equal to DT=5.D-4,
  !         is anticipated when using CDBM.

  !   Effect of magnetic curvature
  !     (It could destabilize a numerical robustness when FS2 > FS1 (see TRCOFS).)
  FSCBKP = 0.D0

  !   Effect of elongation (CDBM05 model)
  !      [M.Honda and A.Fukuyama NF 46 (2006) 580]
  FSCBEL = 1.D0

  !   Factor of E x B rotation shear
  rG1 = 24.D0

  !   =======================================================================

  !   Bohm transport coefficient parameter in SOL
  FSBOHM = 0.D0

  !   Pseud-classical particle transport coefficient parameter in SOL
  FSPCLD = 0.D0

  !   Pseud-classical mom. transport coefficient parameter in SOL
  FSPCLM = 0.D0

  !   Pseud-classical heat transport coefficient parameter in SOL
  FSPCLC = 0.D0

  !   Controller for the degree of the effect of Te'
  !     on the turbulent particle transport
  !   Only valid if MDVAHL = 1 or 2.
  FSVAHL = 0.D0

  !   Particle diffusion coefficient profile parameter (D(r=a)/D(r=0))
  PROFD =  3.D0

  !   Exponent of particle diffusion coefficient profile
  PROFD1 = 3.D0

  !   Gaussian modification of particle diffusion coefficient profile
  PROFD2 = 0.D0

  !   Mom. diffusion coefficient profile parameter (D(r=a)/D(r=0))
  PROFM = 10.D0

  !   Exponent of mom. diffusion coefficient profile
  PROFM1 = 2.D0

  !   Heat diffusion coefficient profile parameter (D(r=a)/D(r=0))
  PROFC = 10.D0

  !   Exponent of heat diffusion coefficient profile
  PROFC1 = 2.D0

  !   ***** Other transport parameters *****

  !   Charge exchange parameter
  FSCX = 1.D0

  !   Orbit loss parameter
  !     FSLC = 0 : No orbit loss included.
  !            1 : Loss term is expressed as damping term.
  !            2 : Loss term is expressed as source and sink term.
  FSLC = 0.D0

  !   Orbit loss model
  !     MDLC = 1 : Shaing model
  !            2 : Itoh model
  MDLC = 1

  !   Ripple loss parameter
  !     FSRP = 0.D0 : No ripple effect
  !            1.D0 : Full ripple effect
  !            2.D0 : Ripple effect but no ripple trapped ions
  FSRP = 0.D0

  !   Alpha heating paramter
  FSNF = 0.D0

  !   Neoclassical thermal diffusivity parameter
  ChiNC = 0.D0

  !   Poloidal neoclassical viscosity parameter
  FSNC = 1.D0

  !   Helical neoclassical viscosity parameter
  FSHL = 0.D0

  !   Particle loss to divertor parameter (default = 0.3)
  FSLP = 0.3D0

  !   Heat loss to divertor parameter (default = 1.0)
  FSLTE = 1.D0
  FSLTI = 1.D0

  !   Heat limiting model
  FSLIM = 0.D0

  !   Ionization parameter
  FSION = 1.D0

  !   Slow neutral diffusion factor
  FSD01 = 1.D0

  !   Thermal neutral diffusion factor
  FSD02 = 1.D0

  !   Halo neutral diffusion factor
  FSD03 = 1.D0

  !   Poloidal torque paramter due to neoclassical heat flux
  FSNCPL = 0.D0

  !   ***** initial parameters *****

  !   Initial Density scale length in SOL (m), valid if MDITSN /= 0
  rLn = 0.03D0

  !   Initail Temperature scale length in SOL (m), valid if MDITST /= 0
  rLT = 0.03D0

  !   ***** Heating parameters *****

  !   NBI beam energy (keV)
  Eb = 32.D0

  !   Heating radius of perp NBI heating (m)
  RNBP  = 0.175D0

  !   Heating center of perp NBI heating (m)
  RNBP0 = 0.D0

  !   Heating radius of first tangential NBI heating (m)
  RNBT1  = 0.175D0

  !   Heating radius of second tangential NBI heating (m)
  RNBT2  = 0.175D0

  !   Heating center of first tangential NBI heating (m)
  RNBT10 = 0.D0

  !   Heating center of second tangential NBI heating (m)
  RNBT20 = 0.D0

  !   Perpendicular NBI input power (MW)
  PNBHP = 0.D0

  !   First tangential NBI input power (MW)
  PNBHT1 = 0.D0

  !   Second tangential NBI input power (MW)
  PNBHT2 = 0.D0

  !   NBI input power from input file (MW)
  PNBHex = 0.D0

  !   NBI current drive parameter
  PNBCD = 1.D0

  !   Rate of the collisional slowing down part of the perpendicular NBI
  PNBMPD = 0.D0

  !   Different NBI deposition profiles for electrons and ions due to banana orbit effect
  !     MDLNBD = 0 : No charge separation
  !              1 : Orbit effect for banana particles only
  !              2 : Orbit effect for all beam ions including passing particles
  MDLNBD = 0

  !   Momentum input model for beam ions
  !     MDLMOM = 0 : Parallel velocity = initial beam velocity, Toroidal torque only
  !              1 : Toroidal and poloidal torques are given distinctly from "Vbpara".
  MDLMOM = 0 

  !   Refractive index of RF waves
  rNRFe = 0.D0
  rNRFi = 0.D0

  !   Heating width of RF heating (m)
  RRFew = 0.175D0
  RRFiw = 0.175D0

  !   Heating center radius of RF heating (m)
  RRFe0 = 0.D0
  RRFi0 = 0.D0

  !   RF input power (MW)
  PRFHe = 0.D0
  PRFHi = 0.D0

  !   Additional torque input (N m)
  Tqt0  = 0.D0  ! Toroidal
  Tqp0  = 0.D0  ! Poloidal

  !   ***** Neutral parameters *****

  !   Initial Neutral density (10^20 m^-3)
  PN0s = 1.D-8

  !   Neutral thermal velocity (m/s)
  !     V0 = SQRT(2.D0*X[eV]*AEE/(PA*AMP))
  V0 = 1.5D3

  !   Recycling rate in SOL
  rGamm0 = 0.8D0

  !   Gas-puff particle flux (10^20 m^-2 1/s)
  !      If you input Gamma_0 [particles/sec], you translate it into
  !      rGASPF [1/(m^2 s)]:
  !        rGASPF = Gamma_0 / (2.D0*PI*RR*2.D0*PI*RB)
  rGASPF = 0.1D0

  !   ***** Ripple loss parameters *****

  ! Number of toroidal field coils in JT-60U
  NTCOIL = 18

  ! Minimum ripple amplitude in the system (R = Rmag0 = 2.4m on JT-60U)
  DltRPn = 0.0002D0 * 1.D-2

  ! Effective ellipticity for ripple amplitude estimation
  kappa = 1.2D0

  !   ***** Parameters for toroidal neoclassical viscosity *****

  ! Maximum poloidal mode number of error field (should be power of two)
  m_pol  = 32

  ! Toroidal mode number of error field
  n_tor  = NTCOIL

  ! Magnetic braiding parameters ***AF (2008-06-08)
  DMAG0   = 0.D0  ! Magnetic field line diffusivity : (Delta r)**2/Delta z [m]
  RMAGMN  = RA    ! minimum radius of the braiding region [m]
  RMAGMX  = RA    ! maximum radius of the braiding region [m]

  !   Helical ripple amplitude at r=a, linear in (r/a)
!  EpsH = 0.1D0     ! amplitude 
  EpsHM(1:NHFMmx,0:3) = 0.d0
!  EpsHM(1,0:3) = 0.D0, 0.D0, 0.1D0, 0.D0     ! amplitude 
  EpsHM(1,2) = 0.1D0
!  NCph = 5        ! toroidal pitch number
!  NCth = 2        ! poloidal pitch number
  HPN(1:NHFMmx, 1:2) = 0
  HPN(1,1) = 2    ! poloidal pitch number
  HPN(1,2) = 5    ! toroidal pitch number
  Q0 = 3.D0       ! q(0) by external coils
  QA = 2.D0       ! q(a) by external coils

  !   ***** Numerical parameters *****

  !   Time step size(s)
!!!!  DT = 1.D-4
  DT = 1.D-3

  !   Convergence parameter
  !      Convergence calculation is going on until IC = ICMAX when EPS is negative.
  EPS = 1.D-2

  !   Iteration
  ICMAX = 10

  !   Time-advancing method
  !     ADV = 0     : Explicit scheme       (Not usable)
  !           0.5   : Crank-Nicolson scheme (Not usable)
  !           2/3   : Galerkin scheme       (Not usable)
  !           0.878 : Liniger scheme
  !           1     : Implicit scheme       (Recommended)
  ADV = 1.D0

  !   Mode of Backward Differential (Gear's) Formula (Second-order accuracy)
  !   (http://www.scholarpedia.org/article/Backward_differentiation_formulas)
  !   (Bibun Houteishiki no Suuchi Kaihou I, Taketomo Mitsui, Iwanami, 1993, p29,30)
  !   (P.M.Gresho and R.L.Sani, "Incompressible Flow and the Finite Element Method,
  !    Vol.1 Advection-Diffusion", John Wiley and Sons, 2000, p.263)
  !     IGBDF = 0   : not used
  !             1   : Use BDF2
  IGBDF = 0

  !   Lower bound of dependent variables
  tiny_cap = 1.d-14

  !   Permittivity switch for numerical convergence
  rMUb1 = rMU0
  rMUb2 = 1.d0

  !   Convergence model
  !     (Definition: ERROR means |X^j - X^{j-1}|, where j denotes the number of iteration.)
  !     0: Relative error 
  !          ERROR divided by root-mean-square,
  !          evaluated for each equation and each grid point.
  !     1: Relative error
  !          L2 norm of ERROR divided by L2 norm of X^{j-1},
  !          evaluated for each equation.
  MODECV = 0

  !   ***** Mesh number parameters *****

  !   Magnitude of mesh peakness
  CMESH0 =  2.D0
  CMESH  = 10.D0

  !   Width of mesh peakness
  WMESH0 = 0.2D0
  WMESH  = 5.D-2

  !   Number of nodes
  NRMAX  = 60

  !   Number of elements
  NEMAX = NRMAX

  !   ***** Time parameter *****

  !   Number of time step
  NTMAX = 10

  !   ***** Graphic parameters *****

  !   Time step interval between print output
  NTSTEP = 10

  !   Time step interval between lines in f(r) graph
  NGRSTP = 1

  !   Time step interval between points in f(t) graph
  NGTSTP = 1

  !   Time step interval between points in f(t) graph
  NGVSTP = 1

  !   Mode of Graph
  !   0 : for Display (with grid, w/o power)
  !   1 : for Display (with grid and power)
  !   2 : for Print Out (w/o grid, with power)
  MODEG = 1

  !   MODE of Graph Line
  !   0 : Change Line Color (Last Color Fixed)
  !   1 : Change Line Color and Style
  !   2 : Change Line Color and Style (With Legend)
  !   3 : Change Line Color, Style and Mark
  !   4 : Change Line Color, Style and Mark (With Legend)
  MODEGL = 1

  !   Mode of AV (Diagnostic message in terms of convergence)
  !   0 : OFF (recommended)
  !   n : Interval of displaying diagnostic message
  MODEAV = 0

  !   Diagnostic parameter
  !   0 : OFF
  !   1 : debug message output (ntstep == 1)
  !   2 : debug message output (few)
  !   3 : debug message output (many)
  !   4 : debug message output (many,each prof)
  ! +10 : convergence output for each equation
  !  -1 : message for steady state check at NGTSTP intervals
  IDIAG = 0

  !   Minimum value of radial axis of TXGRFR/S (normalized radius)
  GRMIN = 0.0
  !   Maximum value of radial axis of TXGRFR/S (IF 0.0, RB/RA is used)
  GRMAX = 0.0

  !   ***** Model parameters *****

  !   Mode of LAPACK
  !   0    : Use BANDRD
  !   else : Use LAPACK_DGBSV or LA_GBSV
  MDLPCK = 0

  !   Mode of orbit squeezing effect in NCLASS
  !   0    : No orbit squeezing effect
  !   1    : Orbit squeezing effect
  !   2    : Orbit squeezing effect, fixed during iteration
  MDOSQZ = 2

  !   Mode of neoclassical resistivity model
  !   0    : original
  !   1    : NCLASS
  !   2    : Sauter
  !   3    : Hirshman, Hawryluk and Birge
  MDLETA = 0

  !   Mode of fixed temperature profile
  !   0    : not fixed
  !   1    : fixed
  MDFIXT = 0

  !   Mode of nonlinear iteration method
  !   0    : fixed point method (or Picard method), linear convergence
  !   1    : Secant method, golden ratio (1.618...) convergence (NOT working now)
  !   2    : Steffensen's method, quadratic convergence (Validity unconfirmed)
  MDSOLV = 0

  !   Mode of initial density profiles
  !   -2   : read from file and smooth
  !   -1   : read from file
  !   0    : original
  MDINTN = 0

  !   Mode of initial temperature profiles
  !   -2   : read from file and smooth
  !   -1   : read from file
  !   0    : original
  !   1    : pedestal model
  !   2    : empirical steady state temperature profile
  MDINTT = 0

  !   Mode of initial current density profiles
  !   -2   : read from file and smooth
  !   -1   : read from file
  !   0    : original
  MDINTC = 0

  !   Mode of initial density profile in the SOL
  !   0    : polynominal model
  !   1    : exponential decay model
  MDITSN = 1

  !   Mode of initial temperature profile in the SOL
  !   0    : polynominal model
  !   1    : exponential decay model
  !   2    : exponential decay model 2
  !          This should be chosen if MDINTT=2.
  MDITST = 1

  !   Mode of initial plasma profils
  !   0    : minimal profiles are calculated
  !   1    : many profiles are analytically calculated 
  MDINIT = 1

  !   Mode of inherent convection annihilator
  !   0    : Nothing to do
  !   1    : Use parameter FSVAHL controlling the contribution of d(lnT)/dln(n)
  !          (Te gradient term is only changed according to FSVAHL.)
  !   2    : Annihilate the contribution of Er to the flux
  !   3    : Annihilate the inherent convection (= obtain pure diffusion)
  !          when FSNC = 0
  MDVAHL = 0

  !   Mode of Edge Transport barrier
  !   0    : Nothing to do
  !   1    : Turn off anomalous effect outside RhoETB
  MDLETB = 0

  !   Multiplication factor for graphic in the radial direction
  !   default : 1.0
  !
  gDIV(1:NGYRM) = 1.0
  gDIV(1)  = 1.E20
  gDIV(2)  = 1.E14
  gDIV(4)  = 1.E3
  gDIV(5)  = 1.E3
  gDIV(7)  = 1.E3
  gDIV(8)  = 1.E3
  gDIV(9)  = 1.E3
  gDIV(12) = 1.E18
  gDIV(13) = 1.E3
  gDIV(16) = 1.E16
  gDIV(17) = 1.E-3
  gDIV(19) = 1.E3
  gDIV(21) = 1.E6
  gDIV(22) = 1.E6
  gDIV(23) = 1.E3
  gDIV(24) = 1.E3
  gDIV(25) = 1.E3
  gDIV(26) = 1.E3
  gDIV(27) = 1.E3
  gDIV(28) = 1.E3
  gDIV(35) = 1.E15
  gDIV(36) = 1.E13
  gDIV(37) = 1.E20
  gDIV(38) = 1.E3
  gDIV(39) = 1.E-3
  gDIV(41) = 1.E3
  gDIV(42) = 1.E3
  gDIV(45) = 1.E3
  gDIV(46) = 1.E3
  gDIV(50) = 1.E3
  gDIV(51) = 1.E3
  gDIV(53) = 1.E-2
  gDIV(55) = 1.E3
  gDIV(56) = 1.E-4
  gDIV(57) = 1.E-6
  gDIV(60) = 1.E6
  gDIV(64) = 1.E-20
  gDIV(65) = 1.E-20
  gDIV(72) = 1.E3
  gDIV(73) = 1.E3
  gDIV(74) = 1.E6
  gDIV(75) = 1.E6
  gDIV(85) = 1.E6
  gDIV(86) = 1.E6
  gDIV(89) = 1.E3
  gDIV(90) = 1.E20
  gDIV(92) = 1.E3
  gDIV(96) = 1.E6
  gDIV(97) = 1.E6
  gDIV(98) = 1.E6
  gDIV(100) = 1.E6
  gDIV(101) = 1.E6
  gDIV(102) = 1.E6
  gDIV(103) = 1.E6
  gDIV(104) = 1.E6
  gDIV(109) = 1.E15
  gDIV(110) = 1.E-4
!  gDIV(115) = 1.E15
  gDIV(121) = 1.E6
  gDIV(123) = 1.E6
  gDIV(124) = 1.E6
  gDIV(133) = 1.E20
  gDIV(136) = 1.E3
  gDIV(137) = 1.E-6
  gDIV(138) = 1.E13
  gDIV(139) = 1.E3
  gDIV(142) = 1.E27
  gDIV(143) = 1.E27

  !   *** Density perturbation technique ***

  !   Normalized radius where density increase
  DelRho = 0.2D0

  !   Amount of increase in density (10^20 m^-3)
  DelN = 1.D-1

  !   ***

  !   Index for graphic save interval

  NGR=-1

  RETURN
END SUBROUTINE TXINIT

!***************************************************************
!
!        Calculate mesh, etc.
!
!***************************************************************

SUBROUTINE TXCALM

  use tx_commons
  use tx_interface, only : LORENTZ, LORENTZ_PART, BISECTION

  implicit none
  INTEGER(4) :: NR, NRL, NR_RC_NEAR
  real(8)    :: MAXAMP, C1L, C2L, W1L, W2L, RL, RC, RCL, CLNEW

  !   Ion mass number
  AMI   = PA * AMP
  !   Beam ion mass number
  AMB   = AMI
  !   Number of equations
  NQMAX = NQM

  !   Square root permittivity for LQm1
  !     for the sake of acceleration of convergence
  sqeps0 = sqrt(EPS0)

  !  Mesh

  if(RACCUM < 0.d0) then
     RC = RA
  else
     RC = RACCUM
  end if

  !  As a trial, generate mesh using given CL and WL and seek the position in the
  !  original coordinate, which becomes the nearest mesh of separatrix after mapping.
  C1L = CMESH0
  C2L = CMESH
  W1L = WMESH0
  W2L = WMESH
  MAXAMP = LORENTZ(RB,C1L,C2L,W1L,W2L,0.D0,RC) / RB
  NR_RC_NEAR = 0
  R(0) = 0.D0
  DO NR = 1, NRMAX - 1
     RL = DBLE(NR) / DBLE(NRMAX) * RB
     CALL BISECTION(LORENTZ,C1L,C2L,W1L,W2L,0.D0,RC,MAXAMP,RL,RB,R(NR))
     IF(ABS(R(NR)-RC) <= ABS(R(NR)-R(NR_RC_NEAR))) NR_RC_NEAR = NR
  END DO
  R(NRMAX) = RB

  !  Construct new CL value that separatrix is just on mesh.
  !  New CL is chosen in order not to be settle so far from given CL.
  !  The mesh finally obtained is well-defined.
  RCL = DBLE(NR_RC_NEAR) / DBLE(NRMAX) * RB
  CLNEW = ( (RC - RCL) * RB - RCL * C1L * LORENTZ_PART(RB,W1L,W2L,0.D0,RC,0) &
       &   + RB * C1L * LORENTZ_PART(RC,W1L,W2L,0.D0,RC,0)) &
       &  / (  RCL * LORENTZ_PART(RB,W1L,W2L,0.D0,RC,1) &
       &     - RB  * LORENTZ_PART(RC,W1L,W2L,0.D0,RC,1))
  MAXAMP = LORENTZ(RB,C1L,CLNEW,W1L,W2L,0.D0,RC) / RB
  R(0) = 0.D0
  DO NR = 1, NRMAX - 1
     RL = DBLE(NR) / DBLE(NRMAX) * RB
     CALL BISECTION(LORENTZ,C1L,CLNEW,W1L,W2L,0.D0,RC,MAXAMP,RL,RB,R(NR))
  END DO
  R(NRMAX) = RB

!!$  !  Maximum NR till RA
!!$
!!$  DO NR = 0, NRMAX-1
!!$     IF(R(NR) <= RA .AND. R(NR+1) >= RA) THEN
!!$        NRL = NR
!!$        EXIT
!!$     END IF
!!$  END DO

  !  Maximum NR till RC

  DO NR = 0, NRMAX-1
     IF(R(NR) <= RC .AND. R(NR+1) >= RC) THEN
        NRL = NR
        EXIT
     END IF
  END DO

  !  Adjust RC on mesh (for the time being RC is assumed to be equivalent to RA)

  IF(ABS(R(NRL)-RC) < ABS(R(NRL+1)-RC)) THEN
     NRA = NRL
     RL = 0.5d0 * ABS(R(NRA) - RC)
     R(NRA  ) = RC
     R(NRA-1) = R(NRA-1) + RL
  ELSE
     NRA = NRL + 1
     RL = 0.5d0 * ABS(R(NRA) - RC)
     R(NRA  ) = RC
     R(NRA+1) = R(NRA+1) - RL
  END IF

  !  RACCUM doesn't coincide with the plasma surface RA

  IF(RACCUM >= 0.D0 .AND. RACCUM /= RA) THEN

     !  Maximum NR till RA

     DO NR = 0, NRMAX-1
        IF(R(NR) <= RA .AND. R(NR+1) >= RA) THEN
           NRL = NR
           EXIT
        END IF
     END DO

     !  Adjust RA on mesh

     IF(ABS(R(NRL)-RA) < ABS(R(NRL+1)-RA)) THEN
        NRA = NRL
        RL = 0.5d0 * ABS(R(NRA) - RA)
        R(NRA  ) = RA
        R(NRA-1) = R(NRA-1) + RL
     ELSE
        NRA = NRL + 1
        RL = 0.5d0 * ABS(R(NRA) - RA)
        R(NRA  ) = RA
        R(NRA+1) = R(NRA+1) - RL
     END IF
     
  END IF

  !  Mesh number at the center of the core plasma

  DO NR = 0, NRMAX-1
     IF(R(NR) <= 0.5d0*RA.AND.R(NR+1) >= 0.5d0*RA) THEN
        IF(ABS(R(NR)-0.5d0*RA) < ABS(R(NR+1)-0.5d0*RA)) THEN
           NRC = NR
        ELSE
           NRC = NR + 1
        END IF
        EXIT
     END IF
  END DO

  !  Mesh coordinate

  RHO(0:NRMAX) = R(0:NRMAX) / RA
  PSI(0:NRMAX) = R(0:NRMAX)**2

  !  Mesh interval

  H   (1:NEMAX) = R  (1:NRMAX) - R  (0:NRMAX-1)
  HPSI(1:NEMAX) = PSI(1:NRMAX) - PSI(0:NRMAX-1)

  RETURN
END SUBROUTINE TXCALM

!***************************************************************
!
!   Initialize profiles
!
!***************************************************************

SUBROUTINE TXPROF

  use tx_commons
  use tx_variables
  use tx_interface, only : INTG_P, INTDERIV3, detect_datatype, INTG_F, dfdx, &
       &                   initprof_input, moving_average, CORR, coulog

  implicit none
  INTEGER(4) :: NR, IER, ifile, NHFM, NR_smt, NR_smt_start = 10
  REAL(8) :: RL, PROF, PROFN, PROFT, PTePROF, PTiPROF!, QL, dRIP
  REAL(8) :: AJFCT, SUM_INT
  REAL(8) :: ALP, dPe, dPi, DR1, DR2
  REAL(8) :: EpsL, Vte, Wte, rNuAsE_inv, FTL, EFT, CR
  real(8) :: FACT, PBA, dPN, CfN1, CfN2, pea, pia, pediv, pidiv, dpea, dpia, &
       &     Cfpe1, Cfpe2, Cfpi1, Cfpi2, sigma, fexp, PN0L, PNaL, PNeDIVL, &
       &     PTe0L, PTi0L, PTeaL, PTiaL, PTeDIVL, PTiDIVL
  REAL(8) :: FCTR!, DERIV4
  real(8), dimension(:), allocatable :: AJPHL, TMP, RHSV, dPedr, dPidr, Prof1, Prof2
  real(8), dimension(:,:), allocatable :: CMTX

  !  Read spline table for neoclassical toroidal viscosity
  IF(FSRP /= 0.D0) CALL Wnm_spline(fmnq, wnm, umnq, nmnqm)

  NEMAX = NRMAX

  !  Define basic quantities like mass of particles, mesh, etc.

  CALL TXCALM

  !  Contribution of perturbed magnetic field
!  IF(DltRPn /= 0.D0) CALL perturb_mag

  !  Initialize variable vector

  X(1:NQMAX,0:NRMAX) = 0.D0

  !  Variables

  if(MDINTN < 0 .or. MDINTT < 0 .or. ABS(MDINTC) /= 0) call initprof_input
  if(MDINTN < 0) then ! density at the boundaries
     call initprof_input(  0,1,PN0L)
     call initprof_input(NRA,1,PNaL)
     PNeDIVL = 0.25d0 * PNaL
  else
     PN0L = PN0
     PNaL = PNa
     PNeDIVL = PNeDIV
  end if
  if(MDINTT < 0) then ! temperature at the boundaries
     call initprof_input(  0,2,PTe0L)
     call initprof_input(  0,3,PTi0L)
     call initprof_input(NRA,2,PTeaL)
     call initprof_input(NRA,3,PTiaL)
     PTeDIVL = 0.25d0 * PTeaL
     PTiDIVL = 0.25d0 * PTiaL
  else
     PTe0L   = PTe0
     PTi0L   = PTi0
     PTeaL   = PTea
     PTiaL   = PTia
     PTeDIVL = PTeDIV
     PTiDIVL = PTiDIV
  end if
  PBA   = RB - RA
  dPN   = - 3.D0 * (PN0L - PNaL) / RA
  CfN1  = - (3.D0 * PBA * dPN + 4.D0 * (PNaL - PNeDIVL)) / PBA**3
  CfN2  =   (2.D0 * PBA * dPN + 3.D0 * (PNaL - PNeDIVL)) / PBA**4
  IF(MDFIXT == 0) THEN
     pea   = PNaL    * PTeaL   ; pia   = PNaL    / PZ * PTiaL
     dpea  = dPN     * PTeaL   ; dpia  = dPN     / PZ * PTiaL
     pediv = PNeDIVL * PTeDIVL ; pidiv = PNeDIVL / PZ * PTiDIVL
  ELSE
     pea   = PTeaL   ; pia   = PTiaL
     dpea  = 0.d0    ; dpia  = 0.d0
     pediv = PTeDIVL ; pidiv = PTiDIVL
  END IF
  Cfpe1 = - (3.D0 * PBA * dpea + 4.D0 * (pea - pediv)) / PBA**3
  Cfpe2 =   (2.D0 * PBA * dpea + 3.D0 * (pea - pediv)) / PBA**4
  Cfpi1 = - (3.D0 * PBA * dpia + 4.D0 * (pia - pidiv)) / PBA**3
  Cfpi2 =   (2.D0 * PBA * dpia + 3.D0 * (pia - pidiv)) / PBA**4

  DO NR = 0, NRMAX
     RL = R(NR)
     IF (RL <= RA) THEN ! +++ Core +++
        IF(MDINTN < 0) THEN
           call initprof_input(NR,1,X(LQe1,NR)) ! Ne
           X(LQi1,NR) = X(LQe1,NR) / PZ         ! Ni
        ELSE
           PROFN = (1.D0 - RHO(NR)**PROFN1)**PROFN2
           X(LQe1,NR) = (PN0L - PNaL) * PROFN + PNaL ! Ne
           X(LQi1,NR) = X(LQe1,NR) / PZ           ! Ni
        END IF
        IF(MDINTT < 0) THEN
           call initprof_input(NR,2,PTePROF) ! Te
           call initprof_input(NR,3,PTiPROF) ! Ti
        ELSE IF(MDINTT == 0) THEN
           PROFT = (1.D0 - RHO(NR)**PROFT1)**PROFT2
           PTePROF = (PTe0L - PTeaL) * PROFT + PTeaL
           PTiPROF = (PTi0L - PTiaL) * PROFT + PTiaL
        ELSE IF(MDINTT == 1) THEN
           PTePROF = (PTe0L - PTeaL) * (0.8263d0 * (1.d0 - RHO(NR)**2 )**1.5d0 &
                &                     + 0.167d0  * (1.d0 - RHO(NR)**30)**1.25d0) + PTeaL
           PTiPROF = (PTi0L - PTiaL) * (0.8263d0 * (1.d0 - RHO(NR)**2 )**1.5d0 &
                &                     + 0.167d0  * (1.d0 - RHO(NR)**30)**1.25d0) + PTiaL
        ELSE
           PTePROF = PTe0L * exp(- log(PTe0L / PTeaL) * RHO(NR)**2)
           PTiPROF = PTi0L * exp(- log(PTi0L / PTiaL) * RHO(NR)**2)
        END IF
        IF(MDFIXT == 0) THEN
           X(LQe5,NR) = PTePROF * X(LQe1,NR) ! Ne*Te
           X(LQi5,NR) = PTiPROF * X(LQi1,NR) ! Ni*Ti
        ELSE 
           X(LQe5,NR) = PTePROF ! Te
           X(LQi5,NR) = PTiPROF ! Ti
        END IF
     ELSE ! +++ SOL +++
        ! density
        IF(MDITSN == 0) THEN
           X(LQe1,NR) = PNaL + dPN * (RL - RA) + CfN1 * (RL - RA)**3 &
                &                              + CfN2 * (RL - RA)**4
           X(LQi1,NR) = X(LQe1,NR) / PZ
        ELSE
           X(LQe1,NR) = PNaL * EXP(- (RL - RA) / rLn)
           X(LQi1,NR) = X(LQe1,NR) / PZ
        END IF
        IF(MDITST == 0) THEN
           X(LQe5,NR) = pea + dpea * (RL - RA) + Cfpe1 * (RL - RA)**3 &
                &                              + Cfpe2 * (RL - RA)**4
           X(LQi5,NR) = pia + dpia * (RL - RA) + Cfpi1 * (RL - RA)**3 &
                &                              + Cfpi2 * (RL - RA)**4
        ELSE
           IF(MDITST == 1) THEN
              PTePROF = PTeaL * EXP(- (RL - RA) / rLT)
              PTiPROF = PTiaL * EXP(- (RL - RA) / rLT)
           ELSE
              sigma = 0.3d0 * (RHO(NRMAX) - 1.d0)
              fexp  = exp(- (RHO(NR) - 1.d0)**2 / (2.d0 * sigma**2))
              PTePROF = PTe0L   * exp(- log(PTe0L / PTeaL) * RHO(NR)**2) *         fexp &
                   &  + PTeDIVL                                          * (1.d0 - fexp)
              PTiPROF = PTi0L   * exp(- log(PTi0L / PTiaL) * RHO(NR)**2) *         fexp &
                   &  + PTiDIVL                                          * (1.d0 - fexp)
           END IF
           ! pressure (MDFIXT=0) or temperature (MDFIXT=1)
           IF(MDFIXT == 0) THEN
              X(LQe5,NR) = PTePROF * X(LQe1,NR)
              X(LQi5,NR) = PTiPROF * X(LQi1,NR)
           ELSE
              X(LQe5,NR) = PTePROF
              X(LQi5,NR) = PTiPROF
           END IF
        END IF
     END IF
     ! N0_1 (slow neutrals)
     X(LQn1,NR) = PN0s
     ! N0_2 (thermal neutrals)
!     X(LQn2,NR) = 0.D0
     X(LQn2,NR) = 1.D-20 ! when ThntSW = 0.D0
     ! N0_3 (halo neutrals)
     X(LQn3,NR) = 0.D0
     ! Bphi
     X(LQm5,NR) = 0.5D0 * PSI(NR) * BB / rMU0
     BphV(NR)   = BB
     ! Fixed densities to keep them constant during iterations
     PNeV_FIX(NR) = X(LQe1,NR)
     PNiV_FIX(NR) = X(LQi1,NR)
     IF(MDFIXT == 0) THEN
        PTeV_FIX(NR) = X(LQe5,NR) / X(LQe1,NR)
        PTiV_FIX(NR) = X(LQi5,NR) / X(LQi1,NR)
     ELSE
        PTeV_FIX(NR) = X(LQe5,NR)
        PTiV_FIX(NR) = X(LQi5,NR)
     END IF
  END DO

  IF(MDINTT == -2) THEN ! Smoothing temperatures
     allocate(Prof1(0:NRMAX),Prof2(0:NRMAX))
     Prof1(0:NRMAX) = X(LQe5,0:NRMAX) / X(LQe1,0:NRMAX)
     Prof2(0:NRMAX) = X(LQi5,0:NRMAX) / X(LQi1,0:NRMAX)
     NR_smt = NRA - NR_smt_start ! smoothing data only in the edge region
     DO NR = NR_smt, NRMAX
        X(LQe5,NR) = moving_average(NR,Prof1,NRMAX) * X(LQe1,NR)
        X(LQi5,NR) = moving_average(NR,Prof2,NRMAX) * X(LQi1,NR)
     END DO
     deallocate(Prof1,Prof2)
  END IF
  IF(MDINTN == -2) THEN ! Smoothing densities
     allocate(Prof1(0:NRMAX))
     X(LQe5,0:NRMAX) = X(LQe5,0:NRMAX) / X(LQe1,0:NRMAX)
     X(LQi5,0:NRMAX) = X(LQi5,0:NRMAX) / X(LQi1,0:NRMAX)
     Prof1(0:NRMAX) = X(LQe1,0:NRMAX)
     NR_smt = NRA - NR_smt_start ! smoothing data only in the edge region
     DO NR = NR_smt, NRMAX
        X(LQe1,NR) = moving_average(NR,Prof1,NRMAX)
        X(LQi1,NR) = X(LQe1,NR) / PZ         ! Ni
     END DO
     X(LQe5,0:NRMAX) = X(LQe5,0:NRMAX) * X(LQe1,0:NRMAX)
     X(LQi5,0:NRMAX) = X(LQi5,0:NRMAX) * X(LQi1,0:NRMAX)
     deallocate(Prof1)
  END IF

  ! Poloidal magnetic field

  IF(FSHL == 0.D0) THEN
     BthV(0) = 0.D0
     DO NR = 1, NRMAX
        RL = R(NR)
        IF(RL < RA) THEN
           PROF = 1.D0 - (RL / RA)**2
           ! Btheta
           BthV(NR) = rMUb1 * rIPs * 1.D6 / (2.D0 * PI * RL) * (1.D0 - PROF**(PROFJ+1))
        ELSE
           BthV(NR) = rMUb1 * rIPs * 1.D6 / (2.D0 * PI * RL)
        END IF
     END DO
  ELSE
     BthV(0) = 0.D0
     Q(0) = Q0
     DO NR = 1, NRMAX
        RL = R(NR)
!        QL = (Q0 - QA) * (1.D0 - (RL / RA)**2) + QA
        Q(NR) = (Q0 - QA) * (1.D0 - (RL / RA)**2) + QA
        BthV(NR) = BB * RL / (Q(NR) * RR)
     END DO
  END IF
  Bthb = BthV(NRMAX)

  ! Toroidal electron current

  allocate(AJPHL(0:NRMAX))
  ifile = detect_datatype('LQe4')
  if(ifile == 0) then
     IF(MDINTC <= -1) THEN
        DO NR = 0, NRA
           call initprof_input(NR,4,AJPHL(NR)) ! NeUeph
        END DO
        AJPHL(NRA+1:NRMAX) = 0.D0
        AJFCT = rIPs * 1.D6 / (2.D0 * PI * INTG_F(AJPHL))
        ! Artificially extrapolate a current density in the SOL for numerical stability
        DO NR = NRA+1, NRMAX
           AJPHL(NR) = AJPHL(NRA) * EXP(- (R(NR) - RA) / (0.5d0 * rLn))
        END DO
        
        IF(MDINTC == -2) THEN ! Smoothing current density
           allocate(Prof1(0:NRMAX))
           Prof1(0:NRMAX) = AJPHL(0:NRMAX)
           NR_smt = NRA - NR_smt_start ! smoothing data only in the edge region
           DO NR = NR_smt, NRMAX
              AJPHL(NR) = moving_average(NR,Prof1,NRMAX)
           END DO
           deallocate(Prof1)
        END IF

        DO NR = 0, NRMAX
           AJPHL(NR) = AJFCT * AJPHL(NR)
           X(LQe4,NR) = - AJPHL(NR) / (AEE * 1.D20)
           AJOH(NR) = AJPHL(NR)
        END DO
        BthV(0) = 0.d0
        SUM_INT = 0.d0
        DO NR = 1, NRMAX
           SUM_INT = SUM_INT + INTG_P(AJPHL,NR,0)
           BthV(NR) = rMUb1 * SUM_INT / R(NR)
        END DO
     ELSE ! (MDINTC == 0)
        DO NR = 0, NRMAX
!!$           AJPHL(NR) = 2.D0 / rMUv1 * DERIV4(NR,PSI,R(0:NRMAX)*BthV(0:NRMAX),NRMAX,0)
!!$           X(LQe4,NR) = - AJPHL(NR) / (AEE * 1.D20)

           IF((1.D0-(R(NR)/RA)**2) <= 0.D0) THEN
              PROF= 0.D0    
           ELSE             
              PROF= 1.D0-(R(NR)/RA)**2
           END IF
           IF(FSHL == 0.D0) THEN
              ! Ne*UePhi
              AJPHL(NR)  = rIPs * 1.D6 / (PI * RA**2) * (PROFJ + 1.D0) * PROF**PROFJ
              X(LQe4,NR) = - AJPHL(NR) / (AEE * 1.D20)
              AJOH(NR)   = AJPHL(NR)
!              AJOH(NR)   = PROF
           ELSE
              AJPHL(NR)  = 0.D0
              X(LQe4,NR) = 0.D0
              AJOH(NR)   = 0.D0
           END IF
        END DO
     END IF
  else
     call inexpolate(infiles(ifile)%nol,infiles(ifile)%r,infiles(ifile)%data,NRMAX,RHO,5,AJPHL)
     AJFCT = rIPs * 1.D6 / (2.D0 * PI * INTG_F(AJPHL))
     DO NR = 0, NRMAX
        AJPHL(NR) = AJFCT * AJPHL(NR)
        X(LQe4,NR) = - AJPHL(NR) / (AEE * 1.D20)
        AJOH(NR) = AJPHL(NR)
     END DO
     BthV(0) = 0.d0
     SUM_INT = 0.d0
     DO NR = 1, NRMAX
        SUM_INT = SUM_INT + INTG_P(AJPHL,NR,0)
        BthV(NR) = rMUb1 * SUM_INT / R(NR)
     END DO
  end if

  if(MDINTN < 0 .or. MDINTT < 0 .or. ABS(MDINTC) /= 0) call initprof_input(idx = 0)

  ! Inverse matrix of derivative formula for integration

  allocate(CMTX(1:NRMAX,1:NRMAX),RHSV(1:NRMAX))
  CMTX(:,:) = 0.D0
  DO NR = 1, NRMAX
     IF(NR == 1) THEN
        DR1 = PSI(NR-1) - PSI(NR)
        DR2 = PSI(NR+1) - PSI(NR)
        CMTX(NR,NR  ) = - (DR1 + DR2) / (DR1 * DR2)
        CMTX(NR,NR+1) = - DR1**2 / (DR1 * DR2 * (DR2 - DR1))
     ELSEIF(NR == NRMAX) THEN
        DR1 = PSI(NR-1) - PSI(NR)
        DR2 = PSI(NR-2) - PSI(NR)
        CMTX(NR,NR-2) = - DR1**2 / (DR1 * DR2 * (DR2 - DR1))
        CMTX(NR,NR-1) =   DR2**2 / (DR1 * DR2 * (DR2 - DR1))
        CMTX(NR,NR  ) = - (DR1 + DR2) / (DR1 * DR2)
     ELSE
        DR1 = PSI(NR-1) - PSI(NR)
        DR2 = PSI(NR+1) - PSI(NR)
        CMTX(NR,NR-1) =   DR2**2 / (DR1 * DR2 * (DR2 - DR1))
        CMTX(NR,NR  ) = - (DR1 + DR2) / (DR1 * DR2)
        CMTX(NR,NR+1) = - DR1**2 / (DR1 * DR2 * (DR2 - DR1))
     END IF
  END DO
  CALL INVMRD(CMTX,NRMAX,NRMAX,IER)

  ! LQm4

  IF((FSHL == 0.0) .AND. (ifile == 0) .AND. (MDINTC == 0) .AND. &
     (PROFJ == 1 .OR. PROFJ == 2 .OR. PROFJ == 3 .OR. &
      PROFJ == 4 .OR. PROFJ == 5)) THEN

  ! Analytic solution for AphV (Valid for PROFJ=1,2,3,4,5, otherwise use below.)

     DO NR = 0, NRMAX
        IF(R(NR) < RA) THEN
           RL = R(NR) / RA
           IF(PROFJ == 1) THEN
              FACT = RL**2*(2.d0-0.5d0*RL**2)
           ELSE IF(PROFJ == 2) THEN
              FACT = RL**2*(3.d0-1.5D0*RL**2+RL**4/3.D0)
           ELSE IF(PROFJ == 3) THEN
              FACT = RL**2*(4.d0-3.d0*RL**2+(4.d0/3.d0)*RL**4-0.25d0*RL**6)
           ELSE IF(PROFJ == 4) THEN
              FACT = RL**2*(5.d0-5.d0*RL**2+(10.d0/3.d0)*RL**4-1.25d0*RL**6+0.2d0*RL**8)
           ELSE IF(PROFJ == 5) THEN
              FACT = RL**2*(6.d0-7.5d0*RL**2+(20.d0/3.d0)*RL**4-3.75d0*RL**6+1.2d0*RL**8 &
                   &-RL**10/6.d0)
           END IF
           X(LQm4,NR) = - rMUb1 * rIPs * 1.D6 / (4.D0 * PI) * FACT
        ELSE
           RL = 1.D0
           IF(PROFJ == 1) THEN
              FACT = RL**2*(2.d0-0.5d0*RL**2)
           ELSE IF(PROFJ == 2) THEN
              FACT = RL**2*(3.d0-1.5D0*RL**2+RL**4/3.D0)
           ELSE IF(PROFJ == 3) THEN
              FACT = RL**2*(4.d0-3.d0*RL**2+(4.d0/3.d0)*RL**4-0.25d0*RL**6)
           ELSE IF(PROFJ == 4) THEN
              FACT = RL**2*(5.d0-5.d0*RL**2+(10.d0/3.d0)*RL**4-1.25d0*RL**6+0.2d0*RL**8)
           ELSE IF(PROFJ == 5) THEN
              FACT = RL**2*(6.d0-7.5d0*RL**2+(20.d0/3.d0)*RL**4-3.75d0*RL**6+1.2d0*RL**8 &
                   &-RL**10/6.d0)
           END IF
           X(LQm4,NR) = - rMUb1 * rIPs * 1.D6 / (4.D0 * PI) * FACT &
                &       - rMUb1 * rIPs * 1.D6 / (4.D0 * PI) * LOG(PSI(NR)/RA**2)
        END IF
     END DO
  ELSE

  ! Numerical solution for AphV

     RHSV(1:NRMAX) = - 0.5D0 * BthV(1:NRMAX) / R(1:NRMAX)
     X(LQm4,0) = 0.D0
     X(LQm4,1:NRMAX) = matmul(CMTX,RHSV)
  END IF

  ! Poloidal current density (Virtual current for helical system)

  IF(FSHL == 0.D0) THEN
     AJV(0:NRMAX)=0.D0

     Q(1:NRMAX) = ABS(R(1:NRMAX) * BB / (RR * BthV(1:NRMAX)))
     Q(0) = FCTR(R(1),R(2),Q(1),Q(2))
  ELSE
     ! Integrate 1 / (r * rMU0) * d/dr (r * BthV) to obtain AJV
!!$     allocate(TMP(0:NRMAX))
!!$     TMP(0:NRMAX) = R(0:NRMAX) * BthV(0:NRMAX)
!!$     DO NR = 1, NRMAX
!!$        dRIP = DERIV4(NR,R,TMP,NRMAX,0) * 2.D0 * PI / rMUb1
!!$        AJV(NR)=dRIP / (2.D0 * PI * R(NR))
!!$     END DO
!!$     AJV(0)=FCTR(R(1),R(2),AJV(1),AJV(2))
!!$     deallocate(TMP)
     AJV(0:NRMAX) = BB / (RR * rMU0) * 2.d0 * Q0 / Q(0:NRMAX)**2
  END IF

  ! Toroidal electric field for initial NCLASS calculation

  DO NR = 0, NRMAX
     ! +++ Hirshman, Hawryluk and Birge model +++
     PNeV(NR) = X(LQe1,NR)
     PNiV(NR) = X(LQi1,NR)
     IF(MDFIXT == 0) THEN
        PTeV(NR) = X(LQe5,NR) / X(LQe1,NR)
        PTiV(NR) = X(LQi5,NR) / X(LQi1,NR)
     ELSE
        PTeV(NR) = X(LQe5,NR)
        PTiV(NR) = X(LQi5,NR)
     END IF
     ! Inverse aspect ratio
     EpsL = R(NR) / RR
     Vte = SQRT(2.D0 * ABS(PTeV(NR)) * rKeV / AME) ! Thermal velocity for ions
     Wte = Vte / (Q(NR) * RR) ! Omega_te; transit frequency for electrons
     rNuei(NR) = PNiV(NR) * 1.D20 * PZ**2 * AEE**4 &
          &     * coulog(1.d0,PNeV(NR),PTeV(NR),PTiV(NR),AEP,1.d0,PA,PZ) &
          &     / (6.D0 * PI * SQRT(2.D0 * PI) * EPS0**2 * SQRT(AME) &
          &     * (ABS(PTeV(NR)) * rKeV)**1.5D0)
     rNuAsE_inv = EpsL**1.5D0 * Wte / (SQRT(2.D0) * rNuei(NR))
     ! Trapped particle fraction
     FTL  = 1.46D0 * SQRT(EpsL) - 0.46D0 * EpsL**1.5D0
     EFT  = FTL * rNuAsE_inv / (rNuAsE_inv + (0.58D0 + 0.20D0 * Zeff))
     ! Spitzer resistivity for hydrogen plasma
     ETAS(NR) = CORR(1.D0) * AME * rNuei(NR) / (PNeV(NR) * 1.D20 * AEE**2)
     CR   = 0.56D0 * (3.D0 - Zeff) / ((3.D0 + Zeff) * Zeff)
     IF(ABS(FSNC) > 0.D0) THEN
        ETA(NR) = ETAS(NR) * Zeff * (1.D0 + 0.27D0 * (Zeff - 1.D0)) &
             &  /((1.D0 - EFT) * (1.D0 - CR * EFT) * (1.D0 + 0.47D0 * (Zeff - 1.D0)))
     ELSE
        ETA(NR) = ETAS(NR)
     END IF
     IF(FSHL == 0.D0) THEN
        ! Ephi
        X(LQm3,NR) = - ETA(NR) *  AJPHL(NR)
     ELSE
        X(LQm3,NR) = 0.D0
     END IF
     IF(X(LQm3,NR) == 0.D0) X(LQm3,NR) = -1.D-4
  END DO
  deallocate(AJPHL)

  !   NBI total input power (MW)
  PNBH = PNBHP + PNBHT1 + PNBHT2

  T_TX=0.D0
  NGT=-1
  NGR=-1
  NGVV=-1
  rIP=rIPs

  !  Check whether (m=0, n>0) Fourier component exists or not.   miki_m 10-09-07
  UHphSwitch = 0
  do NHFM = 1, NHFMmx 
     if (HPN(NHFM,1) == 0 .and. HPN(NHFM,2) > 0) UHphSwitch = 1
  enddo

  !  Define physical variables from X

  CALL TXCALV(X,0) ! Set variables as well as pres0 and ErV0
  !  --- Fixed Er to keep it constant during iterations ---
  ErV_FIX(0:NRMAX) = ErV(0:NRMAX)

  !  Calculate various physical quantities

  CALL TXCALC(0)

  !  Initial condition Part II

  IF(MDINIT == 1) THEN

     allocate(dPedr(0:NRMAX),dPidr(0:NRMAX))
     dPedr(0:NRMAX) = 2.D0 * R(0:NRMAX) * dfdx(PSI,PeV,NRMAX,0) * rKeV
     dPidr(0:NRMAX) = 2.D0 * R(0:NRMAX) * dfdx(PSI,PiV,NRMAX,0) * rKeV

     DO NR = 1, NRMAX
        dPe = dPedr(NR)
        dPi = dPidr(NR)
        IF(NR == NRMAX) THEN
           dPe = 0.D0
           dPi = 0.D0
        END IF
        IF(rNueNC(NR) == 0.D0) THEN
           X(LQe3,NR) = 0.D0
           X(LQi3,NR) = 0.D0
        ELSE
           ALP = (AMI / AME) * (rNuiNC(NR) / rNueNC(NR))
           X(LQi3,NR) =( (- BthV(NR) * X(LQe4,NR) + (dPe + dPi) / AEE) / BphV(NR) &
                &       +(FQeth(NR) + AMI / AME * FQith(NR)) /(rNueNC(NR) * 1.D20)) &
                &     / (PZ + ALP) * R(NR)
           X(LQe3,NR) =- ALP * X(LQi3,NR) + (FQeth(NR) + AMI / AME * FQith(NR)) &
                &                         /(rNueNC(NR) * 1.D20) * R(NR)
        END IF
        ErV(NR)    =- BphV(NR) / PNiV(NR) * X(LQi3,NR) / R(NR) + dPi / (PZ * AEE * PNiV(NR))
        X(LQe2,NR) = (rNueNC(NR) * X(LQe3,NR) - FQeth(NR) * 1.D-20 * R(NR)) &
             &     * AME / (AEE * BphV(NR))
        X(LQi2,NR) = X(LQe2,NR) / PZ
        X(LQm3,NR) = BthV(NR) / PNeV(NR) * X(LQe2,NR) / R(NR) &
             &     + AME * rNuei3(NR) / (AEE * PNeV(NR)) * X(LQe4,NR)
     END DO
     X(LQe3,0) = 0.D0 ; X(LQi3,0) = 0.D0
     ErV(0)    = 0.D0
     X(LQe2,0) = 0.D0 ; X(LQi2,0) = 0.D0
     X(LQm3,0) = AME * rNuei3(0) / (AEE * PNeV(0)) * X(LQe4,0)

     deallocate(dPedr,dPidr)

     ! Scalar potential

     allocate(TMP(0:NRMAX))
     ! TMP(0) is an arbitrary value (INTDERIV3 does not require the value at axis node
     !                               in case of (LAST ARGUMENT)=1.)
     TMP(1:NRMAX) = - 0.5D0 * ErV(1:NRMAX) / R(1:NRMAX)
     TMP(0) = FCTR(PSI(1),PSI(2),TMP(1),TMP(2))
     CALL INTDERIV3(TMP,PSI,X(LQm1,0:NRMAX),0.D0,NRMAX,1)
     
     ! AthV

     TMP(1:NRMAX) = 0.5D0 * rMU0 * AEE * (X(LQe3,1:NRMAX) - PZ * X(LQi3,1:NRMAX)) * 1.D20 &
          &       / PSI(1:NRMAX)
     TMP(0) = FCTR(PSI(1),PSI(2),TMP(1),TMP(2))
     CALL INTDERIV3(TMP,PSI,BphV,BB,NRMAX,1)

     RHSV(1:NRMAX) = 0.5D0 * BphV(1:NRMAX)
     X(LQm5,0) = 0.D0
     X(LQm5,1:NRMAX) = matmul(CMTX,RHSV) / rMU0
     deallocate(CMTX,RHSV,TMP)

     CALL TXCALV(X,0) ! Set variables as well as ErV0
     ErV_FIX(0:NRMAX) = ErV(0:NRMAX)

     CALL TXCALC(0)

  END IF

  !  Calculate global quantities for storing and showing initial status

  CALL TXGLOB

  RETURN
END SUBROUTINE TXPROF

!*****************************************************************************************

module tx_parameter_control
  use tx_commons
  implicit none
  public
  NAMELIST /TX/ &
       & RA,RB,RACCUM,RR,BB, &
       & PA,PZ,Zeff, &
       & PN0,PNa,PTe0,PTea,PTi0,PTia,PROFJ,PROFN1,PROFN2,PROFT1,PROFT2, &
       & De0,Di0,VWpch0,rMue0,rMui0,WPM0, &
       & Chie0,Chii0,ChiNC, &
       & FSDFIX,FSANOM,FSCBKP,FSCBEL,FSCBSH,FSBOHM,FSPCLD,FSPCLM,FSPCLC,FSVAHL,MDANOM,RhoETB, &
       & PROFD,PROFD1,PROFD2,PROFM,PROFM1,PROFC,PROFC1, &
       & FSCX,FSLC,FSRP,FSNF,FSNC,FSLP,FSLTE,FSLTI,FSION, &
       & FSD01,FSD02,FSD03,FSNCPL,MDLC,FSLIM, &
       & rLn,rLT, &
       & Eb,RNBP,RNBP0,RNBT1,RNBT2,RNBT10,RNBT20,PNBHP,PNBHT1,PNBHT2,PNBCD,PNBMPD, &
       & rNRFe,RRFew,RRFe0,PRFHe,Tqt0,Tqp0,rNRFi,RRFiw,RRFi0,PRFHi, &
       & PN0s,V0,rGamm0,rGASPF,PNeDIV,PTeDIV,PTiDIV, &
       & NTCOIL,DltRPn,kappa,m_pol,n_tor, &
       & DT,EPS,ICMAX,ADV,tiny_cap,CMESH0,CMESH,WMESH0,WMESH, &
       & NRMAX,NTMAX,NTSTEP,NGRSTP,NGTSTP,NGVSTP, &
       & DelRho,DelN, &
!       & DMAG0,RMAGMN,RMAGMX,EpsH,NCph,NCth, &
       & DMAG0,RMAGMN,RMAGMX, &
       & rG1,FSHL,Q0,QA, &
       & rIPs,rIPe, &
       & MODEG,gDIV,MODEAV,MODEGL,MDLPCK,MODECV, &
       & MDOSQZ,MDLETA,MDFIXT,MDITSN,MDITST,MDINTN,MDINTT,MDINTC,MDINIT,MDVAHL,MDLETB, &
       & IDIAG,IGBDF,MDSOLV,MDLNBD,MDLMOM, & ! 09/06/17~ miki_m
       & EpsHM,HPN, & ! 10/08/06 miki_m
       & GRMAX,GRMIN  ! 13/08/22 AF
  private :: TXPLST

contains
!***************************************************************
!
!   Change input parameters
!
!***************************************************************

  SUBROUTINE TXPARM(KID)

    INTEGER(4) :: IST
    character(len=*)  :: KID

    DO 
       WRITE(6,*) '# INPUT &TX :'
       READ(5,TX,IOSTAT=IST)
       IF(IST > 0) THEN
          CALL TXPLST
          CYCLE
       ELSE IF(IST < 0) THEN
          KID='Q'
          EXIT
       ELSE
          KID=' '
          EXIT
       END IF
    END DO
    IERR=0

  END SUBROUTINE TXPARM

  SUBROUTINE TXPARL(KLINE)

    integer(4) :: IST
    character(len=*) :: KLINE
    character(len=90) :: KNAME

    KNAME=' &TX '//KLINE//' &END'
    WRITE(7,'(A90)') KNAME
    REWIND(7)
    READ(7,TX,IOSTAT=IST)
    IF(IST /= 0) THEN
       CALL TXPLST
    ELSE
       WRITE(6,'(A)') ' ## PARM INPUT ACCEPTED.'
    END IF
    REWIND(7)
    IERR=0

  END SUBROUTINE TXPARL

  SUBROUTINE TXPARF(KPNAME)

    integer(4) :: IST, KL
    logical :: LEX
    character(len=*) :: KPNAME

    INQUIRE(FILE=KPNAME,EXIST=LEX)
    IF(.NOT.LEX) RETURN

    OPEN(25,FILE=KPNAME,IOSTAT=IST,STATUS='OLD')
    IF(IST > 0) THEN
       WRITE(6,*) 'XX PARM FILE OPEN ERROR : IOSTAT = ',IST
       RETURN
    END IF
    READ(25,TX,IOSTAT=IST)
    IF(IST > 0) THEN
       WRITE(6,*) 'XX PARM FILE READ ERROR'
       RETURN
    ELSE IF(IST < 0) THEN
       WRITE(6,*) 'XX PARM FILE EOF ERROR'
       RETURN
    END IF
    CALL KTRIM(KPNAME,KL)
    WRITE(6,*) &
         &     '## FILE (',KPNAME(1:KL),') IS ASSIGNED FOR PARM INPUT'
    IERR=0

  END SUBROUTINE TXPARF

  SUBROUTINE TXPARM_CHECK

    integer(4) :: idx

    idx = 1
    DO 
       ! /// idx = 1 - 10 ///
       ! System integers
       IF(NRMAX < 0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(NTMAX < 0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(NTSTEP < 0 .OR. NGRSTP < 0 .OR. NGTSTP < 0 .OR. NGVSTP < 0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       ! Physical variables
       IF(RA >= RB) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(RACCUM > RB) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(RR <= RB) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rIPs < 0.D0 .OR. rIPe < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(PA < 0.D0 .OR. PZ < 0.D0 .OR. Zeff < 1.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(PN0 < 0.D0 .OR. PNa < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(PTe0 < 0.D0 .OR. PTea < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       ! /// idx = 11 - 20 ///
       IF(PTi0 < 0.D0 .OR. PTia < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(De0 < 0.D0 .OR. Di0 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rMue0 < 0.D0 .OR. rMui0 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(Chie0 < 0.D0 .OR. Chii0 < 0.D0 .OR. ChiNC < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(minval(FSDFIX) < 0.D0 .OR. minval(FSANOM) < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(minval(RhoETB) < 0.D0 .OR. maxval(RhoETB) > 1.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(FSCBKP < 0.D0 .OR. FSCBEL < 0.D0 .OR. FSCBSH < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(FSBOHM < 0.D0 .OR. FSPCLD < 0.D0 .OR. FSPCLM < 0.D0 .OR. FSPCLC < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(FSLC < 0.D0 .OR. FSRP < 0.D0 .OR. FSNC < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(FSHL < 0.D0 .OR. FSNF < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       ! /// idx = 21 - 30 ///
       IF(FSLP < 0.D0 .OR. FSLTE < 0.D0 .OR. FSLTI < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(FSION < 0.D0)  THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(FSD01 < 0.D0 .OR. FSD02 < 0.D0 .OR. FSD03 < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(FSNCPL < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(MDLC /= 1 .AND. MDLC /= 2) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rG1 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(Eb < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(PNBHP < 0.D0 .OR. PNBHT1 < 0.D0 .OR. PNBHT2 < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(ABS(PNBCD) > 1.D0 .OR. ABS(PNBMPD) > 1.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       IF(RNBP0 > RB .OR. RNBP0 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       ! /// idx = 31 - 40 ///
       IF(RNBT10 > RB .OR. RNBT10 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(RNBT20 > RB .OR. RNBT20 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rNRFe < 0.D0 .OR. PRFHe < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rNRFi < 0.D0 .OR. PRFHi < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(RRFe0 > RB .OR. RRFe0 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(RRFi0 > RB .OR. RRFi0 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(PN0s < 0.D0 .OR. V0 < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rGamm0 < 0.D0 .OR. rGamm0 > 1.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(rGASPF < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(PNeDIV < 0.D0 .OR. PTeDIV < 0.D0 .OR. PTiDIV < 0.D0) THEN ; EXIT ; ELSE
          idx = idx + 1 ; ENDIF
       ! /// idx = 41 - 46 ///
       IF(NTCOIL <= 0 .OR. DltRPn < 0.D0 .OR. DltRPn > 1.D0 .OR. kappa < 0.D0) THEN
          EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(DT < 0.D0 .OR. EPS == 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(ICMAX < 0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(ADV < 0.D0 .OR. ADV > 1.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(tiny_cap < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(CMESH0 < 0.D0 .OR. CMESH < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF
       IF(WMESH0 < 0.D0 .OR. WMESH < 0.D0) THEN ; EXIT ; ELSE ; idx = idx + 1 ; ENDIF

       RETURN
    END DO

    WRITE(6,'(A,I3)') 'XX INPUT ERROR: Please check consistency of INPUT PARAMETERS. idx =',idx
    STOP
    
  END SUBROUTINE TXPARM_CHECK

  !***** INPUT PARAMETER LIST *****

  SUBROUTINE TXPLST

    WRITE(6,601)
    RETURN

601 FORMAT(' ','# &TX : RA,RB,RACCUM,RR,BB,PA,PZ,Zeff,'/ &
         &       ' ',8X,'PN0,PNa,PTe0,PTea,PTi0,PTia,PROFJ,,PROFN1,PROFN2,PROFT1,PROFT2,'/ &
         &       ' ',8X,'De0,Di0,VWpch0,rMue0,rMui0,WPM0,'/ &
         &       ' ',8X,'Chie0,Chii0,ChiNC,'/ &
         &       ' ',8X,'FSDFIX,FSANOM,FSCBKP,FSCBEL,FSCBSH,FSBOHM,FSPCLD,FSPCLM,FSPCLC,FSVAHL,'/ &
         &       ' ',8X,'MDANOM,RhoETB,'/ &
         &       ' ',8X,'PROFD,PROFD1,PROFD2,PROFM,PROFM1,PROFC,PROFC1,'/ &
         &       ' ',8X,'FSCX,FSLC,FSRP,FSNF,FSNC,FSLP,FSLTE,FSLTI,FSION,'/ &
         &       ' ',8X,'FSD01,FSD02,FSD03,FSNCPL,FSLIM,'/&
         &       ' ',8X,'MDLC,rLn,rLT,'/ &
         &       ' ',8X,'Eb,RNBP,RNBP0,RNBT1,RNBT2,RNBT10,RNBT20,PNBHP,PNBHT1,PNBHT2,'/ &
         &       ' ',8X,'PNBCD,PNBMPD,rNRFe,RRFew,RRFe0,PRFHe,Tqt0,Tqp0,rNRFe,RRFew,RRFe0,PRFHe,'/&
         &       ' ',8X,'PN0s,V0,rGamm0,rGASPF,PNeDIV,PTeDIV,PTiDIV,'/ &
         &       ' ',8X,'NTCOIL,DltRPn,kappa,m_pol,n_tor,'/ &
         &       ' ',8X,'DT,EPS,ICMAX,ADV,tiny_cap,CMESH0,CMESH,WMESH0,WMESH,'/ &
         &       ' ',8X,'NRMAX,NTMAX,NTSTEP,NGRSTP,NGTSTP,NGVSTP,'/ &
         &       ' ',8X,'DelRho,DelN,'/ &
!         &       ' ',8X,'Dmag0,RMAGMN,RMAGMX,EpsH,NCph,NCth,'/ & 
         &       ' ',8X,'Dmag0,RMAGMN,RMAGMX,EpsHM,HPN,'/ &   ! 10/08/06 miki_m
         &       ' ',8X,'rG1,FSHL,Q0,QA,'/ &
         &       ' ',8X,'rIPs,rIPe,'/ &
         &       ' ',8X,'MODEG,gDIV,MODEAV,MODEGL,MDLPCK,MODECV,'/ &
         &       ' ',8X,'MDOSQZ,MDLETA,MDFIXT,MDITSN,MDITST,MDINTN,MDINTT,MDINIT,MDVAHL,MDLETB,' / & 
         &       ' ',8X,'IDIAG,IGBDF,MDSOLV,MDLNBD,MDLMOM,GRMIN,GRMAX')
  END SUBROUTINE TXPLST

!***************************************************************
!
!   View input parameters
!
!***************************************************************

  SUBROUTINE TXVIEW

    WRITE(6,'((1X,A6," =",1PD9.2,3(2X,A6," =",1PD9.2)))') &
         &   'RA    ', RA    ,  'RB    ', RB    ,  &
         &   'RR    ', RR    ,  'BB    ', BB    ,  &
         &   'PA    ', PA    ,  'PZ    ', PZ    ,  &
         &   'PN0   ', PN0   ,  'PNa   ', PNa   ,  &
         &   'PTe0  ', PTe0  ,  'PTea  ', PTea  ,  &
         &   'PTi0  ', PTi0  ,  'PTia  ', PTia  ,  &
         &   'rIP   ', rIP   ,  'Zeff  ', Zeff  ,  &
         &   'PROFJ ', PROFJ ,  'PROFN1', PROFN1,  &
         &   'PROFN2', PROFN2,  'PROFT1', PROFT1,  &
         &   'PROFT2', PROFT2,  'CMESH0', CMESH0,  &
         &   'WMESH0', WMESH0,  'CMESH ', CMESH ,  &
         &   'WMESH ', WMESH ,  'ADV   ', ADV   ,  &
         &   'De0   ', De0   ,  'Di0   ', Di0   ,  &
         &   'rMue0 ', rMue0 ,  'rMui0 ', rMui0 ,  &
         &   'VWpch0', VWpch0,  'WPM0  ', WPM0  ,  &
         &   'PROFD ', PROFD ,  'PROFD1', PROFD1,  &
         &   'PROFD2', PROFD2,  'PROFM ', PROFM ,  &
         &   'PROFM1', PROFM1,  'PROFC ', PROFC ,  &
         &   'PROFC1', PROFC1,  'ChiNC ', ChiNC ,  &
         &   'Chie0 ', Chie0 ,  'Chii0 ', Chii0 ,  &
         &   'FSDFX1', FSDFIX(1),  'FSDFX2', FSDFIX(2),  &
         &   'FSDFX3', FSDFIX(3),  'FANOM1', FSANOM(1),  &
         &   'FANOM2', FSANOM(2),  'FANOM3', FSANOM(3),  &
         &   'RoETB1', RhoETB(1),  'RoETB2', RhoETB(2),  &
         &   'RoETB3', RhoETB(3),  'FSCBKP', FSCBKP,  &
         &   'FSCBEL', FSCBEL,  'FSCBSH', FSCBSH,  &
         &   'FSBOHM', FSBOHM,  'FSPCLD', FSPCLD,  &
         &   'FSPCLM', FSPCLM,  'FSPCLC', FSPCLC,  &
         &   'FSVAHL', FSVAHL,  'FSCX  ', FSCX  ,  &
         &   'FSLC  ', FSLC  ,  'FSRP  ', FSRP  ,  &
         &   'FSNF  ', FSNF  ,  'FSNC  ', FSNC  ,  &
         &   'FSLP  ', FSLP  ,  'FSLTE ', FSLTE ,  &
         &   'FSLTI ', FSLTI ,  'FSION ', FSION ,  &
         &   'FSD01 ', FSD01 ,  'FSD02 ', FSD02 ,  &
         &   'FSD03 ', FSD03 ,  'FSNCPL', FSNCPL,  &
         &   'FSLIM ', FSLIM ,  &
         &   'rLn   ', rLn   ,  'rLT   ', rLT   ,  &
         &   'Eb    ', Eb    ,  &
         &   'RNBP  ', RNBP  ,  'RNBP0 ', RNBP0 ,  &
         &   'RNBT1 ', RNBT1 ,  'RNBT10', RNBT10,  &
         &   'RNBT2 ', RNBT2 ,  'RNBT20', RNBT20,  &
         &   'PNBHP ', PNBHP ,  'PNBHT1', PNBHT1,  &
         &   'PNBHT2', PNBHT2,  'PNBHex', PNBHex,  &
         &   'PNBMPD', PNBMPD,  'rNRFe ', rNRFe ,  &
         &   'RRFew ', RRFew ,  'RRFe0 ', RRFe0 ,  &
         &   'PRFHe ', PRFHe ,  'rNRFi ', rNRFe ,  &
         &   'RRFiw ', RRFiw ,  'RRFi0 ', RRFi0 ,  &
         &   'PRFHi ', PRFHi ,  &
         &   'Tqt0  ', Tqt0  ,  'Tqp0  ', Tqp0  ,  &
         &   'rGamm0', rGamm0,  'V0    ', V0    ,  &
         &   'rGASPF', rGASPF,  'PNeDIV', PNeDIV,  &
         &   'PTeDIV', PTeDIV,  'PTiDIV', PTiDIV,  &
         &   'DltRPn', DltRPn,  'kappa ', kappa ,  &
         &   'PN0s  ', PN0s  ,  'EPS   ', EPS   ,  &
         &   'tiny  ', tiny_cap,'DT    ', DT    ,  &
         &   'rG1   ', rG1   ,  'Zeff  ', Zeff  ,  &
         &   'rIPs  ', rIPs  ,  'rIPe  ', rIPe  ,  &
         &   'DMAG0 ', DMAG0 ,  'RMAGMN', RMAGMN,  &
         &   'RMAGMX', RMAGMX,  &
         &   'FSHL  ', FSHL  ,    &  ! Too many elements of EpsHM to show miki_m 10-08-11
         &   'Q0    ', Q0    ,  'QA    ', QA    ,  &
!!$         &   'EpsHM(1,:) ', EpsHM(1,0:3) ,         &
!!$         &   'EpsHM(2,:) ', EpsHM(2,0:3) ,         &
!!$         &   'EpsHM(3,:) ', EpsHM(3,0:3)
         &   'EpsH10 ', EpsHM(1,0), 'EpsH11 ', EpsHM(1,1), &
         &   'EpsH12 ', EpsHM(1,2), 'EpsH13 ', EpsHM(1,3), &
         &   'EpsH20 ', EpsHM(2,0), 'EpsH21 ', EpsHM(2,1), &
         &   'EpsH22 ', EpsHM(2,2), 'EpsH23 ', EpsHM(2,3), &
         &   'EpsH30 ', EpsHM(3,0), 'EpsH31 ', EpsHM(3,1), &
         &   'EpsH32 ', EpsHM(3,2), 'EpsH33 ', EpsHM(3,3), &
         &   'EpsH40 ', EpsHM(4,0), 'EpsH41 ', EpsHM(4,1), &
         &   'EpsH42 ', EpsHM(4,2), 'EpsH43 ', EpsHM(4,3), &
         &   'GRMIN  ', GRMIN,      'GRMAX  ', GRMAX
    WRITE(6,'((" ",A6," =",I5,3(6X,A6," =",I5)))') &
         &   'NRMAX ', NRMAX ,  &
         &   'NTMAX ', NTMAX ,  'NTSTEP', NTSTEP,  &
         &   'NGRSTP', NGRSTP,  'NGTSTP', NGTSTP,  &
         &   'NGVSTP', NGVSTP,  'ICMAX ', ICMAX ,  &
         &   'MODEG ', MODEG ,  'MODEAV', MODEAV,  &
         &   'MODEGL', MODEGL,  'MDLPCK', MDLPCK,  &
         &   'MODECV', MODECV,  &
         &   'MDOSQZ', MDOSQZ,  'MDLETA', MDLETA,  &
         &   'MDANOM', MDANOM,  'MDFIXT', MDFIXT,  &
         &   'MDITSN', MDITSN,  'MDITST', MDITST,  &
         &   'MDINTN', MDINTN,  'MDINTT', MDINTT,  &
         &   'MDINTC', MDINTC,  'MDINIT', MDINIT,  &
         &   'MDVAHL', MDVAHL,  'MDLETB', MDLETB,  &
         &   'IDIAG ', IDIAG ,  'IGBDF ', IGBDF,   &
         &   'MDSOLV', MDSOLV,  &
         &   'NTCOIL', NTCOIL,  'MDLC  ', MDLC,    &
         &   'm_pol ', m_pol ,  'n_tor ', n_tor,   &
         &   'MDLNBD', MDLNBD,  'MDLMOM', MDLMOM ,  &
!         &   'NCph  ', NCph  ,  'NCth  ', NCth,    &
         &   'HPNth1   ', HPN(1,1), 'HPNph1   ', HPN(1,2), &
         &   'HPNth2   ', HPN(2,1), 'HPNph2   ', HPN(2,2), &
         &   'HPNth3   ', HPN(3,1), 'HPNph3   ', HPN(3,2), &
         &   'HPNth4   ', HPN(4,1), 'HPNph4   ', HPN(4,2)
    RETURN
  END SUBROUTINE TXVIEW
end module tx_parameter_control
