2008/01/28

TASK コードの利用説明書
1

TASK コードとは

TASK (Transport Analyzing System for takamaK) コードは，主にトカマクプラズマにおける平衡，輸送，波動伝播，速度分布を解析するコード群である．
1.1

TASK コードの特色• トカマクの時間発展シミュレーション– モジュール構造の統合シミュレーション– 様々な加熱・電流駆動機構の解析– 高い移植性– ヘリカル系への拡張– MPI ライブラリを用いた並列分散処理– 実験データベースの利用• 核燃焼プラズマ統合コード構想のコアコード– 最小限の統合コード：各モジュールは交換可能– インターフェースの標準化：実装の検証– 利用者の拡大：マニュアル等の整備

1.2

TASK コードのモジュール構成TASK/EQ

２次元平衡解析

固定境界，トロイダル回転効果

TR

１次元輸送解析

拡散型輸送方程式，輸送モデル

WR

幾何光学的波動解析

EC, LH: 光線追跡法，ビーム追跡法

WM

波動光学的波動解析

IC, AW: アンテナ励起，固有モード

FP

速度分布解析

相対論的，軌道平均，３次元

DP
PL

波動分散解析
データ交換

局所誘電率テンソル，任意速度分布
座標変換，標準データ，分布データベース

LIB
MTX
MPI
TOT

共通ライブラリ
行列解法
並列化
一体化

行列解法，特殊関数
直接法／反復法，並列化
並列化ライブラリインターフェース
一体化コード

1

1.3

仕様

使用言語は FORTRAN95 であり，allocatable 変数を derived type の component とすることを可能にする拡張を含むコンパイラで，コンパイルすることができる．
使用するグラフィックライブラリは GSAF である．このライブラリは，
http://p-grp.nucleng.kyoto-u.ac.jp/gsaf/

から入手することができる．このライブラリを用いて，図形を X-window に表示し，図形デー
タをファイルに保存し，このデータを postscript ファイルに変換することができる．数値計算ライブラリとして LAPACK を使用する場合があるが，その機能が不必要であれば省略することができる．並列化ライブラリとして MPI を使用しているが，MPI ライブラリがなくてもコンパイ
ルできる．

2

TASK コードのインストール

2.1

TASK コードの入手
TASK コードは次のいずれかの方法で入手することができる．

1. TASK と GSAF の最新ソースファイルは，CVS サーバーから以下の手順で download することができる．この場合に必要な password は fukuyama@nucleng.kyoto-u.ac.jp まで問い合わせること．なお，この方法では修正したファイルを upload することはできない．
export CVS RSH=ssh
export CVSROOT=:pserver:anonymous@p-grp.nucleng.kyoto-u.ac.jp:/home/fukuyama/cvs
cvs login
password:

XXXXXXXX

cvs co bpsd
cvs co task
cvs co gsaf

2. CVS サーバーに account を有する場合には，修正したファイルを upload することができる．この場合の download 手順はexport CVS RSH=ssh
export CVSROOT=p-grp.nucleng.kyoto-u.ac.jp:/home/fukuyama/cvs
cvs co bpsd
cvs co task
cvs co gsaf

3. 圧縮されたファイルツリー taskXXX.tar.gz ファイルは
http://bpsi.nucleng.kyoto-u.ac.jp/task/から入手できる．最新のファイルツリーである保証はないが，比較的安定なバージョンである．得られた task.tar.gz は gnu tar の場合
tar xvzf taskXXX.tar.gz

で，あるいは通常の tar の場合
gzcat taskXXX.tar.gz | tar xvf

等によって解凍することができる．

2

2.2

TASK コードの更新

CVS を利用している場合には，次のいずれかの方法で TASK コードを最新バージョンに更新することができる．ただし，最新バージョンにはバグが含まれている可能性もあるので，注意すること．
1. pserver を利用して download した場合：task ディレクトリにおいて
export CVS RSH=ssh
cvs login
password:

XXXXXXXX

cvs update -Pdを実行する．

2. CVS サーバーに account を有する場合：task ディレクトリにおいて
export CVS RSH=ssh
cvs update -Pdを実行する．

2.3
2.3.1

コードのコンパイル

GSAF のコンパイル

1.

cd gsaf/src：ソースディレクトリに移動する．

2.

cp ../arch/XXXX-XXXX/Makefile.arch .：必要な設定ファイルをコピーする．

• linux-g77: intel86 系 Linux において g77, gcc を利用してコンパイル．
• linux-pgf: intel86 系 Linux において PGI 製 pgf77, gcc を利用してコンパイル．
• linux-pgf: intel86 系 Linux において Intel 製 ifort, icc を利用してコンパイル．
• macosx-g77: Mac OSX において g77, gcc を利用してコンパイル．
• macosx-xlf: Mac OSX において IBM/Absoft 製 xlf, xlc を利用してコンパイル．
• sx-nifs: NIFS 汎用計算機 において sxf90, sxcc を利用してコンパイル．
• hp-nifs: NIFS Application server において f90, cc を利用してコンパイル．
3. Makefile.arch の中の directory BINPATH, LIBPATH を適切に設定する．
• /usr/local が利用できる場合には，/usr/local/bin, /usr/local/lib を利用• /usr/local が利用できない場合には $HOME/bin,$HOME/lib を利用
4.

make：ライブラリを作成する．

5. BINPATH, LIBPATH に書込権限をもつ user として，
make install：コマンドを作成し，ライブラリとともにインストールする．
6. LIBPATH がライブラリ検索 path に含まれているようにする．

3

• LIBPATH を /etc/ld.so.conf に登録し，ldconfig を実行しておく．これらの操作には root 権限が必要．
• 環境変数 LD LIBRARY PATH に LIBPATH を含める．
例：export LD LIBRARY PATH=/usr/local/lib
7.

cd test：テストプログラムを作成する場合には

8. 環境変数 LD LIBRARY PATH に ../lib を含める．
export LD LIBRARY PATH=../lib

9.

make

10.

./bsctest：基本テストプログラムを実行する．

11.

cd ../../..

2.3.2

TASK のコンパイル

まず，task ディレクトリの make.header.org ファイルを make.header という名前でコピーし，その内容を環境に合わせて変更する．基本的には，利用するコンパイラに関連した script 行のコ
メント文字 # を削除する．
次に利用するコードのディレクトリに移動し，
make libs makeと入力すれば，コンパイルリンクが実行される．
原則的に，task ディレクトリで update した場合や他のサブディレクトリのファイルを修正した場合には，make libs が必要であり，現在のサブディレクトリ内のファイルだけを修正した場合には make だけでよい．

2.4

コードの各プログラムの起動

各プログラムはグラフィックライブラリに GSAF を使用しており，プログラム起動時（正確には CALL GSOPEN が実行されたとき）に，原則としてグラフィック出力設定を問い合わせる．

1. 最初の問い合わせは解像度の指定であり，１文字で指定する．’0’ が指定された場合には，画面には出力されず，必要に応じてファイルに図形データが出力さえる．
2. ２番目の問い合わせは，図形ファイルへの出力の指定であり，やはり１文字で指定する．ファイルに出力せずに続行する場合は ’C’ を指定する．表示内容をファイルに保存する場合
は，’Y’，’F’ を指定し，出力ファイル名の問い合わせがある．’Y’ が指定された場合には常に保存，’F’ が指定された場合にはファイルは指定されたが保存は可能になっていない．最初の問い合わせへの回答が無出力 ’0’ の場合には，ファイルが指定されると，特に指定しな
い限り，各ページを保存する．ファイルが指定されなければ，保存はしない．最初の問い合
わせへの回答が ’0’ 以外の場合には，各ページ出力後に問い合わせがあり，単に改行の場合
は標準値，’ Y’ は連続保存，’S’ は単発保存，’ N’ は保存せず，となる．
3. グラフィック出力設定の入力は，環境変数 GSGDP を指定することで，省略することができる．例えばexport GSGDP=3c

4

2.5

コードの各プログラムの入力

入力行の解釈１．行中に "=" を含む場合は namelist 変数の入力
２．先頭の１文字が英字の場合はコマンドの入力
３．それ以外は要求されている変数の入力

namelist 変数の入力
・すべての namelist 変数はファイル xxinit.f のサブルーチン XXINIT で初期化されている．
・実行時のディレクトリにファイル xxparm があれば，namelist 入力として読み込まれる．
・入力行が "=" を含む場合は，先頭に " &xx"，末尾に " &end" を付加し，
namelist 入力として読み込まれる．
・コマンド "P" を入力して，namelist 入力を行うことも可能．
コマンド入力・先頭の１文字が英字の入力行はコマンドの入力
・コマンドの種類を表す英字は大文字と小文字を区別しない．
Ｐ：namelist 変数の入力
Ｖ：namelist 変数等の表示
Ｒ：プログラムの実行開始
Ｃ：プログラムの実行継続
Ｇ：計算結果の図形表示
Ｓ：計算結果の保存
Ｌ：計算結果の読込
Ｑ：実行の終了

3

コンパイル・パラメータ

XXcomn.inc の中で配列の大きさを指定するパラメータが設定されており，それらを変更する
ことにより，計算パラメータ領域を拡大あるいは縮小できる．
plcom1.inc
NSM

[I]

5 Maximum number of particle species

plcom2.inc
NHM
NXM
NPM

[I] 100 Maximum number of cyclotron harmonics[I] 1001 Maximum number of 1D graph mesh points[I] 100 Maximum number of momentum amplitude mesh
5

NTHM
NRM
NGXM
NGYM

[I]
[I]
[I]
[I]

100
25
101
101

Maximum number of pitch angle mesh
Maximum number of radial mesh
Maximum number of graphic x mesh
Maximum number of graphic y mesh

eqcom1.inc
NRGM
NZGM
NPSM

[I]
[I]
[I]

131 Maximum number of R mesh (display)
131 Maximum number of Z mesh (display)
131 Maximum number of psi mesh (display)

eqcom2.inc
NSGM
NTGM
NUGM

[I]
[I]
[I]

64 Maximum number of radial mesh (solver)
64 Maximum number of poloidal mesh (solver)
128 Maximum number of poloidal mesh of radial boundary (solver)

eqcom3.inc
NRM
NTHM
NSUM
NNM

[I] 201 Maximum number of radial mesh (interpolation)[I] 129 Maximum number of poloidal mesh (interpolation)[I] 1025 Maximum number of poloidal mesh of radial boundary (interpolation)[I] 1000 Maximum number of poloidal mesh to follow field line

eqcom4.inc
NTRM

[I]

200 Maximum number of radial mesh (interpolation)

trcom0.inc
NRM
NTM
NGM
NSM
NSZM
N0M
NFM
NCGM
NCTM
NTUM

[I]
80
[I] 10001[I] 1001
[I]
4
[I]
2
[I]
2
[I]
2
[I]
22
[I] 100[I] 1001

Maximum number of radial mesh
Maximum number of time mesh for global quantities
Maximum number of time mesh for profiles
Maximum number of particle species (bulk ions and electrons)
Maximum number of impurity species
Maximum number of neutral species
Maximum number of fast ion species
Maximum number of profiles
Maximum number of global quantities
Maximum number of time mesh in UFILE

6

wrcom1.inc
NEQ
NBEQ
NBVAR
NRAYM
NITM
NRADM

[I]
8
[I]
19
[I]
53
[I]
9
[I] 10000[I] 1000

Number of equations in ray tracing
Number of equations in beam tracing
Number of variables in beam tracing
Maximum number of rays and beams
Maximum number of iterations in tracing
Maximum number of radial division in power deposition

wmcomm.inc
NRM
NDPM
MDPM
NAM
NSUM
NGZM
NTHGM

[I] 201 Maximum number of radial mesh[I]
2 Maximum power of 2 for toroidal mode number[I]
4 Maximum power of 2 for poloidal mode number[I]
4 Maximum number of antennas[I] 2049 Maximum number of poloidal mesh on radial boundary[I] 401 Maximum number of graphic mesh[I]
64 Maximum number of poloidal mesh for graphics

vmcomm.inc
Under reorganization
fpcom1.inc
NRM
NPM
NTHM
NTG1M
NTG2M
NCRM
NLM

[I]
[I]
[I]
[I]
[I]
[I]
[I]

25
50
50
21
501
5
13

Maximum number of radial mesh
Maximum number of momentum amplitude mesh
Maximum number of pitch angle mesh
Maximum number of detail time mesh
Maximum number of time mesh
Maximum number of cyclotron harmonics
Maximum number of Legendre harmonics

dpcom1.inc
NHM
NXM
NPM
NTHM
NRM
NGXM
NGYM

[I] 100 Maximum number of cyclotron harmonics[I] 1001 Maximum number of one-dimensional graphic mesh[I] 100 Maximum number of momentum amplitude mesh[I] 100 Maximum number of pitch angle mesh[I]
25 Maximum number of radial mesh[I] 101 Maximum number of x-mesh for two-dimensional graphics[I] 1011 Maximum number of y-mesh for two-dimensional graphics

7

4
C
C
C
C
C
C
C
C
C
C
C
C

pl コード======( PHYSICAL CONSTANTS )======
PI
AEE
AME
AMP
VC
RMU0
EPS0
CI
PI
AEE
AME
AMP
VC
RMU0
EPS0
CI

C
C
C
C
C
C
C
C
C
C
C
C
C
C

: Pi
: Elementaty charge: Electron mass
: Proton mass
: Speed of light in vacuum: Permeability of free space: Permittivity of free space: Imaginary unit

= 2.D0*ACOS(0.D0)= 1.6021892 D-19
= 9.109534
D-31= 1.6726485 D-27
= 2.99792458 D 8= 4.D0*PI*1.D-7= 1.D0/(VC*VC*RMU0)= (0.D0,1.D0)======( DEVICE PARAMETERS )======
RR: Plasma major radius(m)
RA: Plasma minor radius(m)RB: Wall minor radius(m)
RKAP : Plasma shape elongation
RDEL : Plasma shape triangularity *BB: Magnetic field at center(T)
Q0: Safety factor at center
QA: Safety factor on plasma surface
RIP: Plasma current(MA)
PROFJ : Curren density profile parameter (power of (1 - rho^2))
RR
RA
RB
RKAP
RDLT

= 3.D0= 1.D0= 1.2D0= 1.D0= 0.D0BB
Q0QA= 3.D0= 1.D0= 3.D0C

8

RIP= 3.D0PROFJ = 2.D0
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

======( PLASMA PARAMETERS )======NSMAX : Number of particle species
PA: Mass number
PZ: Charge number
PN: Density at center(1.0E20/m**3)
PNS: Density on plasma surface(1.0E20/m**3)
PZCL : Ratio of collision frequency to wave frequency
PTPR : Parallel temperature at center(keV)
PTPP : Perpendicular temperature at center
(keV)
PTS
: Temperature on surface(keV)
PU: Toroidal rotation velocity at center
(m/s)
PUS
: Toroidal rotation velocity on surface(m/s)
PNITB : Density increment at ITB(1.0E20/Mm*3)
PTITB : Temperature increment at ITB(keV)
PUITB : Toroidal rotation velocity increment at ITB(m/s)
NSMAX = MIN(2,NSM)

C
PA(1)= AME/AMP
PZ(1)=-1.0D0
PN(1)= 1.0D0
PNS(1) = 0.0D0
PZCL(1) = 0.00D0
PTPR(1) = 5.0D0
PTPP(1) = 5.0D0
PTS(1) = 0.05D0PU(1)= 0.D0
PUS(1) = 0.D0
PNITB(1)= 0.D0
PTITB(1)= 0.D0
PUITB(1)= 0.D0C
IF(NSM.GE.2) THEN
PA(2)= 1.0D0
PZ(2)= 1.0D0
PN(2)= 1.0D0
PNS(2) = 0.0D0
PZCL(2) = 0.00D0
PTPR(2) = 5.0D0
PTPP(2) = 5.0D0
9

PTS(2) = 0.05D0PU(2)= 0.D0
PUS(2) = 0.D0
PNITB(2)= 0.D0
PTITB(2)= 0.D0
PUITB(2)= 0.D0
ENDIF
C
DO NS=3,NSM
PA(NS)= 1.0D0PZ(NS)= 1.0D0
PN(NS)= 0.0D0PNS(NS) = 0.0D0
PZCL(NS) = 0.0D0
PTPR(NS) = 5.0D0
PTPP(NS) = 5.0D0
PTS(NS) = 0.0D0
PU(NS)= 0.D0PUS(NS) = 0.D0
PNITB(NS)= 0.D0
PTITB(NS)= 0.D0
PUITB(NS)= 0.D0
ENDDO
C
C
C
C
C
C
C
C
C
C
C

======( PROFILE PARAMETERS )======PROFN1: Density profile parameter (power of rho)
PROFN2: Density profile parameter (power of (1 - rho^PROFN1))
PROFT1: Temperature profile parameter (power of rho)
PROFT2: Temperature profile parameter (power of (1 - rho^PROFN1))PROFU1: Rotation profile parameter (power of rho)
PROFU2: Rotation profile parameter (power of (1 - rho^PROFN1))PROFN1= 2.D0
PROFN2= 0.5D0
PROFT1= 2.D0
PROFT2= 1.D0
PROFU1= 2.D0
PROFU2= 1.D0

C
C
C
C

======( MODEL PARAMETERS )======
MODELG: Control plasma geometry model
10

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

0: Slab geometry
1: Cylindrical geometry
2: Toroidal geometry
3: TASK/EQ output geometry
4: VMEC output geometry
5: EQDSK output geometry
6: Boozer output geometry
MODELN: Control plasma profile
0: Calculated from PN,PNS,PTPR,PTPP,PTS,PU,PUS
1: PT calculated from TASK/EQ pressure profile
2: PN*PT proportional to TASK/EQ pressure profile
8: Read from file by means of WMDPRF routine (DIII-D)
9: Read from file by means of WMXPRF routine (JT-60)
MODELQ: Control safety factor profile (for MODELG=0,1,2)
0: Parabolic q profile (Q0,QA,RHOMIN,RHOITB)
1: Given current profile (RIP,PROFJ)MODELG= 2
MODELN= 0
MODELQ= 0

C
C
C
C
C
C

RHOMIN: rho at minimum q (0 for positive shear)
QMIN : q minimum for reversed shear
RHOITB: rho at ITB (0 for no ITB)RHOEDG: rho at EDGE for smoothing (1 for no smooth)
RHOMIN = 0.D0
QMIN
= 1.5D0
RHOITB = 0.D0
RHOEDG = 1.D0

C
C
C
C
C
C

======( GRAPHIC PARAMETERS )======RHOGMN: minimum rho in radial profile
RHOGMX: maximum rho in radial profile
RHOGMN = 0.D0
RHOGMX = 1.D0

C
C
C
C
C
C

======( MODEL PARAMETERS )======KNAMEQ: Filename of equilibrium data
KNAMWR: Filename of ray tracing data
KNAMWM: Filename of full wave data
11

C
C
C
C

KNAMFP: Filename of Fokker-Planck data
KNAMFO: Filename of File output
KNAMPF: Filename of profile data
KNAMEQ = ’eqdata’
KNAMWR = ’wrdata’
KNAMWM = ’wmdata’
KNAMFP = ’fpdata’
KNAMFO = ’fodata’
KNAMPF = ’pfdata’

5
C
C
C
C
C
C
C
C

EQ コード*** CONSTANTS ****
PI
RMU0
AMP
AEE
PI
RMU0
AMP
AEE

C
C
C
C
C
C
C
C
C
C
C

: Pi
: Permeability of free space: Proton mass
: Electron charge= 2.D0*ASIN(1.D0)= 4.D0*PI*1.D-7= 1.6726231D-27= 1.60217733D-19*** CONFIGURATION PARAMETERS ***
RR
RA
RB
RKAP
RDLT
BB
RIP
RR
RA
RB
RKAP
RDLT
BB
RIP

: Plasma major radius: Plasma minor radius: Wall minor radius: Plasma shape elongation: Plasma shape triangularity: Magnetic field at center: Plasma current= 3.D0= 1.D0= RA*1.1D0= 1.6D0= 0.25D0= 3.D0= 3.D0

C
12

(m)
(m)
(m)

(T)
(MA)

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

*** PROFILE PARAMETERS ***PP0: Plasma pressure (main component)
PP1: Plasma pressure (sub component)
PP2: Plasma pressure (increment within ITB)PROFP0: Pressure profile parameter
PROFP1: Pressure profile parameter
PROFP2: Pressure profile parameter

(MPa)
(MPa)
(MPa)

PPSI=PP0*(1.D0-PSIN**PROFR0)**PROFP0+PP1*(1.D0-PSIN**PROFR1)**PROFP1+PP2*(1.D0-(PSIN/PSIITB)**PROFR2)**PROFP2&&

The third term exists for RHO < RHOITB
PP0= 0.001D0
PP1= 0.0D0
PP2= 0.0D0
PROFP0 = 1.5D0
PROFP1 = 1.5D0
PROFP2 = 2.0D0

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

PJ0
: Current density at R=RR (main component) : Fixed to 1
PJ1
: Current density at R=RR (sub component)
(arb)
PJ2
: Current density at R=RR (sub component)(arb)
PROFJ0: Current density profile parameter
PROFJ1: Current density profile parameter
PROFJ2: Current density profile parameter
HJPSI=-PJ0*(1.D0-PSIN**PROFR0)**PROFJ0*PSIN**(PROFR0-1.D0)-PJ1*(1.D0-PSIN**PROFR1)**PROFJ1*PSIN**(PROFR1-1.D0)-PJ2*(1.D0-PSIN**PROFR2)**PROFJ2*PSIN**(PROFR2-1.D0)

&
&
&
&
&

The third term exists for RHO < RHOITB
PJ0= 1.00D0
PJ1= 0.0D0
PJ2= 0.0D0
PROFJ0 = 1.5D0
PROFJ1 = 1.5D0
PROFJ2 = 1.5D0
13

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

FF0
: Current density at R=RR (main component) : Fixed to 1
FF1
: Current density at R=RR (sub component)
(arb)
FF2
: Current density at R=RR (sub component)(arb)
PROFF0: Current density profile parameter
PROFF1: Current density profile parameter
PROFF2: Current density profile parameter
FPSI=BB*RR+FF0*(1.D0-PSIN**PROFR0)**PROFF0+FF1*(1.D0-PSIN**PROFR1)**PROFF1+FF2*(1.D0-PSIN**PROFR2)**PROFF2&&&

The third term exists for RHO < RHOITB
FF0= 1.0D0
FF1= 0.0D0
FF2= 0.0D0
PROFF0 = 1.5D0
PROFF1 = 1.5D0
PROFF2 = 1.5D0

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

QQ0
: Safety factor on axis for QQ1=QQ2=0
QQS: Safety factor on surface
QQ0: Safety factor
QQ1: Safety factor (sub component)
QQ2: Safety factor (increment within ITB)PROFQ0: Safety factor profile parameter
PROFQ1: Safety factor profile parameter
PROFP2: Pressure profile parameter&&&

QPSI=QQS+(QQ0-QQS)*(1.D0-PSIN**PROFR0)**PROFQ0+QQ1*(1.D0-PSIN**PROFR1)**PROFQ1+QQ2*(1.D0-(PSIN/PSIITB)**PROFR2)**PROFQ2
The third term exists for RHO < RHOITB
QQ0= 1.D0QQS= 3.D0
QQ1= 0.0D0
QQ2= 0.0D0
PROFQ0 = 1.0D0
PROFQ1 = 1.0D0
14

PROFQ2 = 1.0D0
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

PT0
: Plasma temperature (main component)
PT1: Plasma temperature (sub component)
PT2: Plasma temperature (increment within ITB)PTS: Plasma temperature (at surface)
PROFT0: Temperature profile parameter
PROFT1: Temperature profile parameter
PROFT2: Temperature profile parameter

&
&
&

(keV)
(keV)
(keV)
(keV)

TPSI=PTS+(PT0-PTS)*(1.D0-PSIN**PROFR0)**PROFT0+PT1*(1.D0-PSIN**PROFR1)**PROFT1+PT2*(1.D0-PSIN/PSIITB)**PROFR2)**PROFT2+PTS
The third term exits for RHO < RHOITB
PT0= 1.0D0
PT1= 0.0D0
PT2= 0.0D0
PTS= 0.05D0
PROFT0 = 1.5D0
PROFT1 = 1.5D0
PROFT2 = 2.0D0

C
C
C
C
C
C
C
C
C
C
C
C
C
C

PV0
: Toroidal rotation (main component)
PV1: Toroidal rotation (sub component)
PV2: Toroidal rotation (increment within ITB)PROFV0: Velocity profile parameter
PROFV1: Velocity profile parameter
PROFV2: Velocity profile parameter&&

PVSI=PV0*(1.D0-PSIN**PROFR0)**PROFV0+PV1*(1.D0-PSIN**PROFR1)**PROFV1+PV2*(1.D0-(PSIN/PSIITB)**PROFR2)**PROFV2
The third term exits for RHO < RHOITB
PV0= 0.0D0
PV1= 0.0D0
PV2= 0.0D0
PROFV0 = 1.5D0
PROFV1 = 1.5D0
PROFV2 = 2.0D0
15

(m/s)
(m/s)
(m/s)

C
C
C

PN0 : Plasma number density(constant)
PN0

C
C
C
C
C
C

= 1.D20

PROFR0: Profile parameter
PROFR1: Profile parameter
PROFR2: Profile parameter
RHOITB: Normalized radius SQRT(PSI/PSIA) at ITB
PROFR0 = 1.D0
PROFR1 = 2.D0
PROFR2 = 2.D0
RHOITB = 0.5D0

C
C
C
C

OTC
HM

: Constant OMEGA**2/TPSI: Constant(Am)

OTC = 0.15D0
HM = 1.D6
C
C
C
C
C
C
C
C
C
C
C
C
C

*** MESH PARAMETERS ***NSGMAX: Number of radial mesh points for Grad-Shafranov eq.
NTGMAX: Number of poloidal mesh points for Grad-Shafranov eq.
NUGMAX: Number of radial mesh points for flux-average quantities
NRGMAX: Number of horizontal mesh points in R-Z plane
NZGMAX: Number of vertical mesh points in R-Z plane
NPSMAX: Number of flux surfaces
NRMAX : Number of radial mesh points for flux coordinates
NTHMAX: Number of poloidal mesh points for flux coordinates
NSUMAX: Number of boundary points
NSGMAX = 32
NTGMAX = 32
NUGMAX = 32C
NRGMAX = 32
NZGMAX = 32
NPSMAX = 21C
NRMAX = 50
NTHMAX = 64
NSUMAX = 65
16

C
C
C
C
C

*** CONTROL PARAMETERS ***EPSEQ : Convergence criterion for equilibrium
EPSEQ

C
C
C
C
C
C
C
C
C
C
C

= 1.D-6

MDLEQF : Profile parameter
0: given analytic profile P,Jtoroidal,T,Vph
1: given analytic profile P,F,T,Vph
2: given analytic profile P,Jparallel,T,Vph
3: given analytic profile P,q,T,Vph
5: given spline profile P,Jtoroial,T,Vph
6: given spline profile P,F,T,Vph
7: given spline profile P,Jparapllel,T,Vph
8: given spline profile P,q,T,Vph
MDLEQF = 0

C
C
C
C
C

MDLEQC : Poloidal coordinate parameter
0: Poloidal length coordinate
1: Boozer coordinate
MDLEQC = 0

CCC

NPRINT: Level print out
NPRINT= 0

C
C
C
C
C

*** FILE NAME ***KNAMEQ: Filename of equilibrium data
KNAMEQ = ’eqdata’

6
C
C
C
C
C
C
C

WR コード*** CONSTANTS ****
PI
AEE
AME
AMM

: Pi
: Elementaty charge: Electron mass
: Proton mass
17

C
C
C
C
C

VC
RMU0
EPS0
VOID
PI
AEE
AME
AMM
VC
RMU0
EPS0
RKEV
VOID

C
C
C
C
C
C
C
C
C
C
C
C
C

= ASIN(1.D0)*2.D0= 1.60217733D-19= 9.1093897D-31= 1.6726231D-27= 2.99792458D8= 4.D0*PI*1.D-7= 1.D0/(VC*VC*RMU0)= AEE*1.D3= 0.D0==== DEVICE PARAMETERS ====
RR
RA
RKAP
RDLT
BB
RIPS
RIPE
RIPSS
RHOA
RR
RA
RKAP
RDLT
BB
RIPS
RIPE
RIPSS
RHOA

C
C
C
C
C
C
C
C

: Speed of light in vacuum: Permeability of free space: Permittivity of free space: 0.D0: PLASMA MAJOR RADIUS (M): PLASMA MINOR RADIUS (M): ELIPTICITY OF POLOIDAL CROSS SECTION: TRIANGULARITY OF POLOIDAL CROSS SECTION: TOROIDAL MAGNETIC FIELD ON PLASMA AXIS (T): INITIAL VALUE OF PLASMA CURRENT (MA): FINAL VALUE OF PLASMA CURRENT (MA): VALUE OF PLASMA CURRENT FOR INITIAL CONVERGENCE: EDGE OF CALCULATE REGION (NORMALIZED SMALL RADIUS)= 3.0D0= 1.2D0= 1.5D0= 0.0D0= 3.D0= 3.D0= 3.D0= 3.D0= 1.D0==== PLASMA PARAMETERS ====NSMAX : NUMBER OF MAIN PARTICLE SPECIES (NS=1:ELECTRON)NSZMAX : NUMBER OF IMPURITIES SPECIESNSNMAX : NUMBER OF NEUTRAL SPECIESPA(NS) : ATOMIC NUMBER
18

C
C
C
C
C
C

PZ(NS) : CHARGE NUMBERPN(NS) : INITIAL NUMBER DENSITY ON AXIS (1.E20 M**-3)
PNS(NS): INITIAL NUMBER DENSITY ON SURFACE (1.E20 M**-3)
PT(NS) : INITIAL TEMPERATURE ON AXIS (KEV)PTS(IS): INITIAL TEMPERATURE ON SURFACE (KEV)NSMAX=2
NSZMAX=0
NSNMAX=2! the number of impurities! the number of neutrals, 0 or 2 fixed

C
PA(1)
PZ(1)
PN(1)
PT(1)
PTS(1)PNS(1)= AME/AMM=-1.D0= 0.5D0= 1.5D0= 0.05D0= 0.05D0

PA(2)
PZ(2)
PN(2)
PT(2)
PTS(2)PNS(2)= 2.D0= 1.D0= 0.5D0-2.D-7= 1.5D0= 0.05D0= 0.05D0-2.D-8

PA(3)
PZ(3)
PN(3)
PT(3)
PTS(3)PNS(3)= 3.D0= 1.D0= 1.D-7= 1.5D0= 0.05D0= 1.D-8

PA(4)
PZ(4)
PN(4)
PT(4)
PTS(4)PNS(4)= 4.D0= 2.D0= 1.D-7= 1.5D0= 0.05D0= 1.D-8

PA(5)
PZ(5)
PN(5)
PT(5)
PTS(5)PNS(5)= 12.D0= 2.D0= VOID= 0.D0= 0.D0= VOID

C

C

C

C

19

C
PA(6)
PZ(6)
PN(6)
PT(6)
PTS(6)PNS(6)= 12.D0= 4.D0= VOID= 0.D0= 0.D0= VOID

PA(7)
PZ(7)
PN(7)
PT(7)
PTS(7)PNS(7)= 2.D0= 0.D0= 1.D-15= 0.D0= 0.D0= 2.D-4

PA(8)
PZ(8)
PN(8)
PT(8)
PTS(8)PNS(8)= 2.D0= 0.D0= 1.D-15= 0.D0= 0.D0= 1.D-15

C

C

C
C
C
C
C
C
C
C
C

==== IMPURITY PARAMETERS ====
PNC
PNFE
PNNU
PNNUS
PNC
PNFE
PNNU
PNNUS

C
C
C
C
C
C
C
C
C
C

: CARBON DENSITY FACTOR: IRON DENSITY FACTORCOMPARED WITH ITER PHYSICS DESIGN GUIDELINE: NEUTRAL NUMBER DENSITY ON AXIS (1.E20 M**-3):
ON SURFACE (1.E20 M**-3)= 0.D0= 0.D0= 0.D0= 0.D0==== PROFILE PARAMETERS ====PROFN*: PROFILE PARAMETER OF INITIAL DENSITYPROFT*: PROFILE PARAMETER OF INITIAL TEMPERATURE
PROFU*: PROFILE PARAMETER OF INITIAL NEUTRAL DENSITYPROFJ*: PROFILE PARAMETER OF INITIAL CURRENT DENSITY(X0-XS)(1-RHO**PROFX1)**PROFX2+XS
ALP: ADDITIONAL PARAMETERS
20

C
C
C
C

ALP(1): RADIUS REDUCTION FACTOR
ALP(2): MASS WEIGHTING FACTOR FOR NC
ALP(3): CHARGE WEIGHTING FACTOR FOR NC
PROFN1 = 2.D0
PROFN2 = 0.5D0
PROFT1 = 2.D0
PROFT2 = 1.D0
PROFU1 =12.D0
PROFU2 = 1.D0
PROFJ1 =-2.D0
PROFJ2 = 1.D0

C
ALP(1) = 1.0D0
ALP(2) = 0.D0
ALP(3) = 0.D0
C
C
C
C
C
C
C
C

==== TRANSPORT PARAMETERS ====AV0
: INWARD PARTICLE PINCH FACTOR
AD0: PARTICLE DIFFUSION FACTOR
CNC: COEFFICIENT FOR NEOCLASICAL DIFFUSION
CDW(8) : COEFFICIENTS FOR DW MODEL
AV0
AD0= 0.5D0= 0.5D0C
CNC= 1.D0
CDW(1) = 0.04D0
CDW(2) = 0.04D0
CDW(3) = 0.04D0
CDW(4) = 0.04D0
CDW(5) = 0.04D0
CDW(6) = 0.04D0
CDW(7) = 0.04D0
CDW(8) = 0.04D0
C
C
C
C
C
C
C
C

==== TRANSPORT MODEL ====MDLKAI: TURBULENT TRANSPORT MODEL***** 0.GE.MDLKAI.LT.10 : CONSTANT COEFFICIENT MODEL ********** 10.GE.MDLKAI.LT.20 : DRIFT WAVE (+ITG +ETG) MODEL ********** 20.GE.MDLKAI.LT.30 : REBU-LALLA MODEL *****
21

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

***** 30.GE.MDLKAI.LT.40 : CURRENT-DIFFUSIVITY DRIVEN MODEL ********** 40.GE.MDLKAI.LT.60 : DRIFT WAVE BALLOONING MODEL **********MDLKAI.GE.60 : ITG(/TEM, ETG) MODEL ETC *****
*****
*****
*****
*****

MDLKAI.EQ. 0
MDLKAI.EQ. 1
MDLKAI.EQ. 2
MDLKAI.EQ. 3: CONSTANT *****: CONSTANT/(1-A*rho^2) *****: CONSTANT*(dTi/drho)^B/(1-A*rho^2) *****: CONSTANT*(dTi/drho)^B*Ti^C *****

*****
*****
*****
*****
*****
*****

MDLKAI.EQ. 10
MDLKAI.EQ. 11
MDLKAI.EQ. 12
MDLKAI.EQ. 13
MDLKAI.EQ. 14
MDLKAI.EQ. 15: etac=1 *****: etac=1 1/(1+exp) *****: etac=1 1/(1+exp) *q *****: etac=1 1/(1+exp) *(1+q^2) *****: etac=1+2.5*(Ln/RR-0.2) 1/(1+exp) *****: etac=1 1/(1+exp) func(q,eps,Ln) **********MDLKAI.EQ. 20: Rebu-Lalla model *****

*****
*****
*****
*****
*****
*****
*****
*****
*****
*****
*****

MDLKAI.EQ. 30
MDLKAI.EQ. 31
MDLKAI.EQ. 32
MDLKAI.EQ. 33
MDLKAI.EQ. 34
MDLKAI.EQ. 35
MDLKAI.EQ. 36
MDLKAI.EQ. 37
MDLKAI.EQ. 38
MDLKAI.EQ. 39
MDLKAI.EQ. 40: CDBM 1/(1+s) *****: CDBM F(s,alpha,kappaq) *****: CDBM F(s,alpha,kappaq)/(1+WE1^2) *****: CDBM F(s,0,kappaq) *****: CDBM F(s,0,kappaq)/(1+WE1^2) *****: CDBM (s-alpha)^2/(1+s^2.5) *****: CDBM (s-alpha)^2/(1+s^2.5)/(1+WE1^2) *****: CDBM s^2/(1+s^2.5) *****: CDBM s^2/(1+s^2.5)/(1+WE1^2) *****: CDBM F2(s,alpha,kappaq,a/R) *****: CDBM F3(s,alpha,kappaq,a/R)/(1+WS1^2) *****

*****
*****
*****
*****
*****

MDLKAI.EQ. 60
MDLKAI.EQ. 61
MDLKAI.EQ. 62
MDLKAI.EQ. 63
MDLKAI.EQ. 64: GLF23 model *****: GLF23 (stability enhanced version) *****: IFS/PPPL model *****: Weiland model *****: Bohm/Gyro-Bohm model *****MDLETA: RESISTIVITY MODEL
0: CLASSICAL
1: NEOCLASSICALMDLAD : PARTICLE DIFFUSION MODEL
0: NO PARTICL TRANSPORT
1: CONSTANT DMDLAVK: HEAT PINCH MODEL
0: NO HEAT PINCH
22

C
C
C

MDLJBS: BOOTSTRAP CURRENT MODELMDLKNS: NEOCLASSICAL TRANSPORT MODEL
MDLKAI = 31
MDLETA = 3
MDLAD = 3
MDLAVK = 3
MDLJBS = 5
MDLKNC = 1

C
C
C
C
C

MDLWLD : Weiland model mode selector0: using effective transport coefficientselse : using transport coefficients’ vectors
MDLWLD=0

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

MDDW : mode selector for anomalous particle transport coefficient.
you must NOT modify this parameter.
0: if MDDW=0 from start to finish when you choosea certain transport model (MDLKAI),
you could control a ratio of anomalous particle
transport to total particle transport to manipulate
the factor of AD0.
else : this is because you chose MDLKAI=60, 61, or 63
which assign the transport models that can calculate
an anomalous particle transport coefficienton their own.
MDDW=0==== Semi-Empirical Parameter for Anomalous Transport ====
CHP
= 0.D0CK0= 12.D0 ! for electron
CK1= 12.D0 ! for ionsCWEB
= 0.D0 ! for omega ExB
CALF= 1.D0 ! for s-alpha
CKALFA = 0.D0
CKBETA = 0.D0
CKGUMA = 0.D0

C
C
C
C
C

==== CONTROL PARAMETERS ====
DT
NRMAX

: SIZE OF TIME STEP
: NUMBER OF RADIAL MESH POINTS
23

C
C
C
C
C
C
C

NTMAX : NUMBER OF TIME STEPNTSTEP : INTERVAL OF SNAP DATA PRINTNGRSTP : INTERVAL OF RADIAL PROFILE SAVENGTSTP : INTERVAL OF TIME EVOLUTION SAVENGPST : ???
TSST
: ???
DT
= 0.01D0
NRMAX = 50
NTMAX = 100
NTSTEP = 10
NGRSTP = 100
NGTSTP = 2
NGPST = 4
TSST
= 1.D9

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

==== Convergence Parameter ====EPSLTR : CONVERGENCE CRITERION OF ITERATIONLMAXTR : MAXIMUM COUNT OF ITERATIONEPSLTR = 0.001D0EPSLTR = 1.D99
LMAXTR = 10==== SAWTOOTH PARAMETERS ====
TPRST
MDLST

IZERO
TPRST
MDLST
IZERO

C
C
C
C
C
C
C

: SAWTOOTH PERIOD (S): SAWTOOTH MODEL TYPE
0:OFF
1:ON
: SAWTOOTH CRASH TYPE

= 0.1D0= 0= 3==== FUSION REACTION PARAMETERS ====
MDLNF

MDLNF

: FUSION REACTION MODEL TYPE
0:OFF
1:ON

= 0
24

C
C
C
C
C
C
C
C
C
C
C
C
C

==== NBI HEATING PARAMETERS ====PNBTOT : NBI TOTAL INPUT POWER (MW)PNBR0 : RADIAL POSITION OF NBI POWER DEPOSITION (M)PNBRW : RADIAL WIDTH OF NBI POWER DEPOSITION (M)
PNBENG : NBI BEAM ENERGY (keV)PNBRTG : TANGENTIAL RADIUS OF NBI BEAM (M)PNBCD : CURRENT DRIVE FACTORMDLNB : NBI MODEL TYPE
0:OFF
1:GAUSSIAN
2:PENCIL BEAMPNBTOT = 0.D0
PNBR0 = 0.D0
PNBRW = 0.5D0
PNBENG = 80.D0
PNBRTG = 3.D0
PNBCD = 1.D0
MDLNB = 1

C
C
C
C
C
C
C
C
C
C
C

==== ECRF PARAMETERS ====PECTOT : ECRF INPUT POWER (MW)PECR0 : RADIAL POSITION OF POWER DEPOSITION (M)PECRW : RADIAL WIDTH OF POWER DEPOSITION (M)PECTOE : POWER PARTITION TO ELECTRONPECNPR : PARALLEL REFRACTIVE INDEXPECCD : CURRENT DRIVE FACTOR
MDLEC : ECRF MODELPECTOT = 0.D0
PECR0 = 0.D0
PECRW = 0.2D0
PECTOE = 1.D0
PECNPR = 0.D0
PECCD = 0.D0
MDLEC = 0

C
C
C
C
C
C

==== LHRF PARAMETERS ====PLHTOT : LHRF INPUT POWER (MW)
PLHR0 : RADIAL POSITION OF POWER DEPOSITION (M)
PLHRW : RADIAL WIDTH OF POWER DEPOSITION (M)
25

C
C
C
C
C

PLHTOE : POWER PARTITION TO ELECTRONPLHNPR : PARALLEL REFRACTIVE INDEXPLHCD : CURRENT DRIVE FACTORMDLLH : LHRF MODEL
PLHTOT = 0.D0
PLHR0 = 0.D0
PLHRW = 0.2D0
PLHTOE = 1.D0
PLHNPR = 2.D0
MDLLH = 0

C
C
C
C
C
C
C
C
C
C
C

==== ICRF PARAMETERS ====PICTOT : ICRF INPUT POWER (MW)PICR0 : RADIAL POSITION OF POWER DEPOSITION (M)
PICRW : RADIAL WIDTH OF POWER DEPOSITION (M)
PICTOE : POWER PARTITION TO ELECTRONPICNPR : PARALLEL REFRACTIVE INDEXPICCD : CURRENT DRIVE FACTORMDLIC : ICRF MODEL
PICTOT = 0.D0
PICR0 = 0.D0
PICRW = 0.5D0
PICTOE = 0.5D0
PICNPR = 2.D0
PICCD = 0.D0
MDLIC = 0

C
C
C
C
C
C
C
C

==== CURRENT DRIVE PARAMETERS ====PBSCD : BOOTSTRAP CURRENT DRIVE FACTORMDLCD : CURRENT DRIVE OPERATION MODEL
0: TOTAL PLASMA CURRENT FIXED
1: TOTAL PLASMA CURRENT VARIABLE
PBSCD
MDLCD

C
C
C
C
C

= 1.D0= 0==== PELLET INJECTION PARAMETERS ====MDLPEL : PELLET INJECTION MODEL TYPE
0:OFF 1:GAUSSIAN 2:NAKAMURA
26

3:HO

C
C
C
C
C
C
C
C

PELTOT : TOTAL NUMBER OF PARTICLES IN PELLETPELR0 : RADIAL POSITION OF PELLET DEPOSITION (M)PELRW : RADIAL WIDTH OF PELLET DEPOSITION (M)PELRAD : RADIUS OF PELLET (M)PELVEL : PELLET INJECTION VELOCITY (M/S)
PELTIM : TIME FOR PELLET TO BE INJECTEDPELPAT : PARTICLE RATIO IN PELLET’
MDLPEL = 1
PELTOT = 0.D0
PELR0 = 0.D0
PELRW = 0.5D0PELRAD = 0.D0
PELVEL = 0.D0
PELTIM = -10.D0

C
DO NS=1,NSMAXPELPAT(NS) = 1.0D0
ENDDO
C
C
C
C
C

==== DEVICE NAME AND SHOT NUMBER IN UFILE DATA ====KUFDEV : DEVICE NAMEKUFDCG : DISCHARGE NUMBERKUFDEV=’jet’
KUFDCG=’19649’

C
C
C
C

==== LOG FILE NAME ====KFNLOG : LOG FILE NAMEKFNLOG=’trf.log’

C
C
C
C
C
C
C
C
C

==== INTERACTION WITH EQ ====
MODELG: 0 : TR ONLY
3 : USING GEOMETRIC FACTORS FROM EQ FOR INITIAL PROFILEMODELQ: 0 : TR ONLY
3 : TR/EQ COUPLEDNTEQIT: STEP INTERVAL OF EQ CALCULATIONMODELG=0
MODELQ=0
NTEQIT=10

C
27

C
C
C
C
C
C
C
C

==== INPUT FROM UFILE ====MDLUF :
0 : not used
1 : time evolution
2 : steady state
3 : compared with TOPICS
MDLUF=0

C
C
C
C
C
C
C
C
C

==== IMPURITY TREATMENT ====MDNI :
0 : NSMAX=2, ne=ni
1 : calculate nimp and zeff profiles from NE, ZIMP and NM1
2 : calculate nimp and ni profiles from NE, ZIMP and ZEFFR
3 : calculate zeff and ni profiles from NE, ZIMP and NIMPMDNI=0CCCCC==== INITIAL PROFILE SWITCH ====MODEP : initial profile selector for steady-state simulation
MODEP=3

C
C
C
C
C
C
C
C

==== INITIAL CURRENT PROFILE SWITCH ====MDLJQ :
0 : create AJ(NR) profile from experimental Q profile
1 : create QP(NR) profile from experimental CURTOT profile
MDLJQ=0C
C
C*** Eqs. Selection Parameter ***MDLEQB=1
MDLEQN=0
MDLEQT=1
MDLEQU=0
MDLEQZ=0
MDLEQ0=0
MDLEQE=0! 0/1 for B_theta! 0/1 for density! 0/1 for heat! 0/1 for rotation! 0/1 for impurity! 0/1 for neutral! 0/1 for electron density
28

C
MDLEOI=0
C
C
C
C
C
C
C
C
C

7

! 0/1/2 for electron only or bulk ion only if NSMAX=1*** NCLASS SWITCH ***0: offelse : on
MDNCLS=0*** MODERATE TIME EVOLUTION FOR ANOMALOUS TRANSPORT COEFFICIENTS ***0: offelse : multiplyer for TAUK (which is the required time ofaveraging magnetic surface)
MDTC=0WR コードWR コードは光線追跡法あるいはビーム追跡法による伝播解析，分散関係表示などができる．

7.1

WR コードの実行
P : パラメータを変更する．
V : パラメータを表示する．
R : 光線追跡法を実行する．
B : ビーム追跡法を実行する．
G : 伝播解析結果のグラフを表示する．
S : 伝播解析結果のデータを保存する．
1 : 分散関係を表示する．
2 : 分散関係を表示する．
3 : 分散関係を表示する．
F : 分散式の解を求める．
Q : 終了する．

一連の作業の流れとしては P でパラメータを変更した後，R で光線追跡を実行した後，G でグラフを見たり，S で伝播解析の結果を保存したりする．

7.2

P,V : 入力パラメータ

入力パラメータは，namelist を用いて，任意のパラメータを変更することができる．入力行は，まず空白１文字の後，’&wr ’ に引き続いて，「パラメータ名＝新しい値」の形式で設定を繰29

り返し，最後に，’&end ’ を入力して終了する．入力例は&wr RR=8.14,RA=2.8,BB=5.68 &end

&wr PN=0.2,0.1,0.1 &end

ここで PN= は PN(1),PN(2),PN(3)=に対応する．
以下に入力パラメータの説明と標準値を示す．

RR

プラズマ主半径 : 3.0 [m]RA

プラズマ小半径 : 1.0 [m]RB

壁小半径 : 1.2 [m]RKAP

楕円率 : 1.0

RDLT

三角形度 : 0.0BB

中心での磁場 : 3.0 [T]Q0r=0 での安全係数 : 1.0

QA

r=RA での安全係数 : 3.0

RIP

全電流 : 3.0 [MA]

PROFJ

電流分布パラメータ : 2.0

PROFN1

密度分布形状パラメータ : 2.0D0

PROFN2密度分布形状パラメータ : 0.5D0

PROFT1

温度分布形状パラメータ : 2.0D0

PROFT2温度分布形状パラメータ : 1.0D0

PROFU1平行速度分布形状パラメータ : 2.0D0

PROFU2平行速度分布形状パラメータ : 1.0D0

NSMAX

粒子種の数 : 2PA

原子の質量 [陽子質量]: PA(1)=5.4462e-4, PA(2)=1.0

PZ

電荷の数 [素電荷]: PZ(1)=-1.0, PZ(2)=1.0

PN

中心密度 [1020 /m3 ] : PN(1)=1.0, PN(2)=1.0

PNS

周辺密度 [1020 /m3 ] : PNS(1)=0.0, PNS(2)=0.0

PZCL

衝突周波数 [ν/ω] : PZCL(1)=0.0, PZCL(2)=0.0

PTPR

中心平行方向温度 [keV] : PTPR(1)=5.0, PTPR(2)=5.0

PTPP

中心垂直方向温度 [keV] : PTPP(1)=5.0, PTPP(2)=5.0

PTS

周辺温度 [keV] : PTS(1)=0.05, PTS(2)=0.05

PU

中心平行方向速度 [m/s] : PU(1)=0.0, PU(2)=0.0

PUS

周辺平行方向速度 [m/s] : PUS(1)=0.0, PUS(2)=0.0

PNITB

ITB での密度増分 [1020 /m3 ] : PNITB(1)=0.0

PTITB

ITB での温度増分 [keV] : PTITB(1)=0.0

PUITB

ITB での速度増分 [m/s] : PUITB(1)=0.0（1 : 電子）30

MODELG

配位モデル

0 : 平板モデル
1 : 円柱モデル
2 : トカマクモデル
3 : T ASK/EQ 平衡配位
4 : V M EC 平衡配位

MODELN

径方向分布モデル

0 : P N, P N S, P T, P T S 等で指定
1 : P N は指定，P T は平衡圧力から計算
2 : P N ∗ P T が平衡圧力に比例
9 : 分布データ読み込み

MODELQ

安全係数分布モデル (MODELG=0,1,2)

0 : QO, QA を指定
1 : RIP, P ROF J を指定

RHOMIN

安全係数が極小となる規格化半径: 0.D0

QMIN

極小安全係数RHOITBITB を与える規格化半径: 0.D0

RHOEDG

プラズマ表面での分布の平滑化を与える規格化半径: 1.D0

RHOGMN

径方向分布グラフの規格化半径の下限: 0.00

RHOGMX

径方向分布グラフの規格化半径の上限: 1.D0

KNAMEQ

平衡データファイル名: eqdata

KNAMWR

波動伝播データファイル名: wrdata

KNAMFP

速度分布データファイル名: fpdata

KNAMFO

数値データファイル名: fodata

MODELP

誘電率テンソル : MODELP(1)=5,MODELP(2)=0

0 : 無衝突の冷たいプラズマ
1 : 衝突のある冷たいプラズマ
2 : 理想電磁流体プラズマ
3 : 抵抗性電磁流体プラズマ31

4 : 有限ラーモア半径効果を無視した運動論的プラズマ
5 : 有限ラーモア半径効果を取り入れた運動論的プラズマ
6 : 相対論効果を取り入れた運動論的プラズマ
7 : 速度分布を与えた運動論的プラズマ
8 : ジャイロ運動論的プラズマ
9 : 速度分布を与えたジャイロ運動論的プラズマ
0 − 9 : 伝播＝与えられたモデル
偏波＝与えられたモデル
吸収＝与えられたモデル

10 − 19 : 伝播＝冷たいプラズマモデル
偏波＝与えられたモデル
吸収＝与えられたモデル

20 − 29 : 伝播＝冷たいプラズマモデル
偏波＝冷たいプラズマモデル
吸収＝与えられたモデル

NDISP1

最小サイクロトロン高調波番号：NDISP1(1)=-2, NDISP1(2)=-2

NDISP2

最大サイクロトロン高調波番号：NDISP2(1)=2, NDISP2(2)=2MODELV

速度分布モデル0 : 非相対論的マクスウェル速度分布
1 : 非相対論的任意速度分布（ファイルから読み込み）
2 : 相対論的マクスウェル速度分布
3 : 相対論的任意速度分布（ファイルから読み込み）RF

周波数 [MHz]

RPI

初期大半径位置 [m]ZPI

初期垂直位置 [m]PHII

初期トロイダル角 [Rad]RNZI

初期垂直方向屈折率

RNPHII

初期トロイダル方向屈折率

RKR0

径方向波数の初期推定値（Newton 法の初期値）

UUI

規格化パワー初期値

SMAX

光線長の最大値 [m] : 5.00

DELS

光線の刻み幅 [m] : 1.00e-2

32

UUMIN

光線を追跡する最小パワー : 1.00e-4

NRAYMX

光線本数

EPSRAY

常微分方程式の収束判定条件

DELRAY

常微分方程式のステップ幅の下限

DELDER

数値微分のステップ幅

DELKR

ニュートン法における数値微分のステップ幅

EPSNW

ニュートン法における収束判定条件

LMAXNW

ニュートン法における反復回数の上限

INTYPE

計算出発パラメータの入力形式

0 : RF, RP, ZP, P HI, RKR0, RN Z, RN P HI, U U
1 : RF, RP, ZP, P HI, RKR0, AN GZ, AN GP H, U U
2 : RF, RP, ZP, P HI, M ODE, AN GZ, AN GP H, U U
IGTYPE

グラフの表示形式

0:
1:
IQTYPE

光線追跡法の常微分方程式解法

0:
1:
2:

7.3

全トーラス
部分トーラス

Runge-Kutta, 固定幅
Runge-Kutta, 自動幅Runge-Kutta-Fahlberg, 自動幅

NRZMAX

吸収パワーの径方向分布を求めるための分割数

NRADMX

ビーム追跡法における吸収パワー分布分割数

RCURVA

ビーム波面の初期曲率半径（k と B に垂直）RCURVB

ビーム波面の初期曲率半径（k と k × B に垂直）

RBRADAビームの初期半径（k と B に垂直）

RBRADBビームの初期半径（k と k × B に垂直）R : 光線追跡法を実行する
R を入力した場合は光線追跡法を実行する．

170e3,10.8,,,-2000,,0.8/など．現在値のままの場合は省略可能であり，/ は以降を省略する．光線の数だけ繰り返す．
パラメータの物理的意味を以下に示す．

• INTYPE=033

入力データ
RF
RPI
ZPI
PHII
RKR0
RNZI
RNPHIIUUI

周波数 [MHz]
初期主半径 R [m]初期垂直方向位置 Z [m]初期トロイダル角 [radian]
径方向波数の初期推測値（ニュートン法の初期値）[1/m]
初期垂直方向屈折率
初期トロイダル方向屈折率
規格化パワー初期値• INTYPE=1
入力データ
RF
RPI
ZPI
PHII
RKR0
ANGZ
ANGPH
UUI

周波数 [MHz]
初期主半径 R [m]初期垂直方向位置 Z [m]初期トロイダル角 [radian]
径方向波数の初期推測値（ニュートン法の初期値）[1/m]
初期ポロイダル方向入射角
初期トロイダル方向入射角
規格化パワー初期値• INTYPE=2（未サポート）
入力データ
RF
RPI
ZPI
PHII
MODEW
ANGZ
ANGPH
UUI

7.4

周波数 [MHz]
初期主半径 R [m]初期垂直方向位置 Z [m]初期トロイダル角 [radian]
モード選択（0:slow wave, 1:fast wave）
初期ポロイダル方向入射角
初期トロイダル方向入射角
規格化パワー初期値

G : グラフ表示’1’ ：ポロイダル軌跡とパワー分布’2’ ：径方向依存性１’3’ ：径方向依存性２’4’ ：ビーム軌跡とパワー分布
34

’5’ ：偏光面角度と s’6’ ：波数方向’X’ ：終了
8

FP コードFP コードは波動による電流駆動を解析するために，相対論効果や捕捉粒子の寄与を含めて速度分布の時間発展を記述することができる．
8.1

FP コードを使う
R : FP 方程式の時間発展計算を開始する．
C : FP 方程式の時間発展計算を続行する．
P : パラメータを変更する．
V : パラメータを表示する．
G : 結果のグラフを表示する．
F : 結果を ascii 形式でファイルに出力する．
I : 過去の履歴データをクリアする．
W : 結果を再表示する．
Y : FP 方程式の係数を計算する．
S : 速度分布関数をファイルに保存する．
L : 速度分布関数をファイルから読み込む．
．
Q : 終了する．

作業の流れの例としては，P でパラメータを変更した後，R でフォッカープランク方程式を解き，，G でグラフを見たり，S で速度分布関数を保存したりする．
コンパイル・パラメータ：
XXcomn.inc の中で配列の大きさを指定するパラメータが設定されており，それらを変更する
ことにより，計算パラメータ領域を拡大あるいは縮小できる．

fpcom1.inc
NRM
NPM
NTHM
NTG1M
NTG2M

[I]
[I]
[I]
[I]
[I]

25
50
50
21
501

Maximum number of radial mesh
Maximum number of momentum amplitude mesh
Maximum number of pitch angle mesh
Maximum number of detail time mesh
Maximum number of time mesh
35

NCRM
NLM

8.2

[I]
[I]

5 Maximum number of cyclotron harmonics
13 Maximum number of Legendre harmonics

P,V : 入力パラメータ

入力パラメータは，namelist を用いて，任意のパラメータを変更することができる．入力行は，
まず空白１文字の後，’&fp ’ に引き続いて，「パラメータ名＝新しい値」の形式で設定を繰り返し，
最後に，’&end ’ を入力して終了する．入力例は&fp DELT=0.1,NTMAX=5, &end
以下に入力パラメータの説明と標準値を示す．

R1

NRMAX=1 の場合の半径方向の位置 [m]DELR1

NRMAX=1 の場合の半径方向の仮想的間隔 [m]RMIN

NRMAX6=1 の場合の最小半径

RMAX

NRMAX6=1 の場合の最大半径E0

トロイダル電場 [V/m] :DRR0

半径方向の拡散係数 [m2 /s]

DEC

規格化された電子サイクロトロン波による拡散係数

PEC1

電子サイクロトロン波の N パラレルスペクトルの中心

PEC2電子サイクロトロン波の N パラレルスペクトルの幅

RFEC

電子サイクロトロン周波数 [MHz]DELYEC

電子サイクロトロン波ビームの垂直幅 [m]DLH

規格化された低域混成波の拡散係数

PLH1

低域混成波スペクトル（最小速度，またはスペクトルの中心）

PLH2

低域混成波スペクトル（最大速度，またはスペクトルの幅）

RLH

低域混成波の最小近接小半径 [m]DFW

規格化された速波の拡散係数

RFW

速波の最小近接短軸半径 [m]PFW1

速波スペクトル（最小速度，またはスペクトルの中心）

PFW2

速波スペクトル（最大速度，またはスペクトルの幅）

RFDW

周波数 [MHz] :DELNPR トロイダル方向屈折率のスペクトル幅
NCMIN

サイクロトロン高調波番号の下限

NCMAX

サイクロトロン高調波番号の上限

CEWR

波動電界の径方向成分36

CEWTH

波動電界のポロイダル方向成分

CEWPH

波動電界のトロイダル方向成分

RKWR

波数の径方向成分

RKWTH

波数のポロイダル方向成分

RKWPH

波数のトロイダル方向成分

REWY

光線の垂直方向位置

DREWY

光線の垂直方向幅

PMAX

中心熱運動量で正規化された最大運動量

DELT

時間ステップ幅 [s]RIMPL

計算の implicit パラメーター

EPSM

行列方程式を解く場合の収束限界

EPSE

電場計算での収束限界

LMAXE

電場計算での最大繰り返し数

EPSDE

2 重指数積分法での収束限界

H0DE

2 重指数積分法での初期ステップ幅

NGLINE

等高線の最大本数

LLMAX

Legendre 展開の最大次数NPMAX

運動量の大きさ方向の分割数

NTHMAX 運動量の角度方向の分割数
NRMAX

半径方向の分割数

NAVMAX 波の拡散係数を計算するための軌道平均の分割数
NTMAX

最大時間ステップ幅

NTSTP1

半径方向分布形状データを保存するための時間ステップ幅

NTSTP2

全体的なデータを保存するための時間ステップ幅

NTSTPC

係数を再計算するための時間ステップ幅

MODELE 電場の計算をする場合には 1 とする
MODELR 相対効果を含む場合には 1 とする
MODELA バウンス平均をする場合には 1 とする
MODELC 非線形衝突演算子を用いる場合には 1 とする
MODELW 拡散係数の計算モデル
0 : 近似モデル式
1 : 近似電磁界を用いて計算
2 : 波動伝播解析 W R の結果を読み込んで計算

37

PWAVE

入力パワー

LMAXNWR 光線と磁気面の交点を求める Newton 法の反復回数上限
EPSNWR 光線と磁気面の交点を求める Newton 法の収束判定条件
8.3

G, F : グラフ表示およびファイル出力
各種グラフ出力およびファイル出力の説明を以下に示す．

F1
F2
FX2
FS11
FS12
FS21
FS22
DPP
DPT
DTP
DTT
DRR
DCPP
DCPT
DCTP
DCTT
DCRR
DWPP
DWPT
DWTP
DWTT
DWRR
FP
FT
FR
FCP
FCT
FCR
FEP
FET
FER
RN
RI
RW
RPC
RPW

１次元速度分布２次元速度分布２次元速度分布のピッチ角差分
１次元速度分布の内側境界値
２次元速度分布の内側境界値
１次元速度分布の外側境界値
２次元速度分布の外側境界値
拡散係数 Dpp の運動量依存性
拡散係数 Dpθ の運動量依存性
拡散係数 Dθp の運動量依存性
拡散係数 Dθθ の運動量依存性
拡散係数 Drr の運動量依存性
衝突拡散係数 DCpp の運動量依存性
衝突拡散係数 DCpθ の運動量依存性
衝突拡散係数 DCθp の運動量依存性
衝突拡散係数 DCθθ の運動量依存性
衝突拡散係数 DCrr の運動量依存性
波動による拡散係数 DWpp の運動量依存性
波動による拡散係数 DWpθ の運動量依存性
波動による拡散係数 DWθp の運動量依存性
波動による拡散係数 DWθθ の運動量依存性
波動による拡散係数 DWrr の運動量依存性
摩擦係数 Fp の運動量依存性
摩擦係数 Fθ の運動量依存性
摩擦係数 Fr の運動量依存性
衝突摩擦係数 FCp の運動量依存性
衝突摩擦係数 FCθ の運動量依存性
衝突摩擦係数 FCr の運動量依存性
電界加速係数 FEp の運動量依存性
電界加速係数 FEθ の運動量依存性
電界加速係数 FEr の運動量依存性
電子密度の径方向依存性
電流密度の径方向依存性
エネルギー密度の径方向依存性
衝突吸収パワー密度の径方向依存性
波動吸収パワー密度の径方向依存性

38

RPE
RT
RQ
RE
TN
TI
TW
TPC
TPW
TPE
TT
TQ
TE

9

WM コード9.1

利用説明

電界吸収パワー密度の径方向依存性
温度の径方向依存性
安全係数の径方向依存性
トロイダル電界の径方向依存性
電子密度の時間依存性
電流の時間依存性
エネルギーの時間依存性
衝突吸収パワーの時間依存性
波動吸収パワーの時間依存性
電界吸収パワーの時間依存性
温度の時間依存性
安全係数の時間依存性
トロイダル電界の時間依存性

コマンド入力Ｐ：namelist 変数の入力
Ｖ：namelist 変数等の表示
Ｗ：アンテナ励起による波動伝播の計算
Ａ：与えられた複素周波数に対する振幅の計算
Ｆ：振幅パラメータの実周波数依存性の計算
Ｃ：振幅パラメータの複素周波数平面における等高線表示
Ｅ：与えられた複素周波数から出発して固有周波数の計算
Ｓ：固有周波数のパラメータ依存性の計算
Ｇ：固有関数の図形表示
Ｑ：実行の終了
コードにおける振幅パラメータの定義：
プラズマ中の電子密度に比例する分布をもつ励起電流に対して，
電界振幅の絶対値の自乗を空間積分した値の逆数
コードにおける固有周波数の定義：
複素周波数平面において振幅が極小をもち，極大値の逆数が EPSNW よりも小さい周波数
図形指定コマンド入力

ABC/DEF/GHI は A,B,C の中から１文字，D,E,F の中から１文字，
G,H,I の中から１文字の計３文字を順に並べた入力を表す．

39

大文字と小文字は区別しない．
注）現在，電流密度分布は正しく計算されていない．
電磁界の径方向分布
R/EB/ATMN
R: 径方向分布
E: 波動電界と吸収パワー密度
B: 波動磁界と駆動電流密度
A: ポロイダル角・トロイダル角指定（パワーは角度積分値）
T: ポロイダル角・トロイダル角指定（パワーは角度指定）
M: ポロイダルモード・トロイダルモード指定
N: 複数ポロイダルモード・単一トロイダルモード指定例）REN は波動電界と吸収パワー密度の径方向分布を，
複数のポロイダルモード成分を重ねて表示する．
電磁界の２次元分布
CPM/EB/RTZ/RIA
C: 円形断面等高線表示
P: ポロイダル断面等高線表示
M: ポロイダルモード数・径方向面等高線表示
E: 波動電界
B: 波動磁界
R: 径方向成分
T: ポロイダル成分
Z: トロイダル成分
R: 電磁界の実数成分（励起電流と同位相の成分）
I: 電磁界の虚数成分（励起電流と位相が 90 度遅れた成分）
A: 電磁界の絶対値
吸収パワー密度の２次元分布
RCPM/P/123
R: 径方向分布
C: 円形断面等高線表示
P: ポロイダル断面等高線表示
M: ポロイダルモード数・径方向面等高線表示
P: 吸収パワー密度1: 粒子種１＝電子
2: 粒子種２＝通常は多数イオン
3: 粒子種３＝通常は少数イオン
駆動電流密度の２次元分布
CP/J

40

C: 円形断面等高線表示
P: ポロイダル断面等高線表示
プラズマの空間分布
P/F/TPQB
T: 温度分布P: 圧力分布
Q: 安全係数
B: 磁界強度
測度テンソル成分の径方向分布
P/G
ヘリカル系における磁気面S
２次元表示の種類設定
G/1234
1: 等高線表示
2: 塗り分け表示
3: 鳥瞰図格子表示
4: 鳥瞰図等高線表示
図形表示の終了
X

9.2
C
C
C
C
C
C
C

入力パラメータ**** ALPHA PARTICLE PARAMETERS ****
PNA
PNAL
PTA

: Alpha density at center: Density scale length: Effective temperature

PNA = 0.02D0
PNAL = 0.5D0
PTA = 3.5D3
C
C
C
C
C

**** ZEFF PARAMETERS ****ZEFF

: Effective Z (sum n Z^2 / sum n Z)

ZEFF

= 2.D041(1.0E20/Mm*3)
(m)
(keV)

C
C
C
C
C
C
C
C
C
C
C

*** WAVE PARAMETER ***CRF: Wave frequency(MHz)RD: Antenna minor radius(m)
BETAJ : Antenna current profile parameter
NTH0 : Central value of poloidal mode number
NPH0 : Central value of toroidal mode number
NHC: Number of helical coils
PRFIN : Input Power (0 for given antenna current) (W)
CRF
RD
BETAJ
NTH0
NPH0
NHC
PRFIN

C
C
C
C
C
C
C
C
C
C
C

= (50.0D0,0.D0)= 1.1D0= 0.D0= 0= 8= 10= 0.D0*** ANTENNA PARAMETERS ***NAMAX : Number of antennae
AJ: Antenna current density
APH: Antenna phase
THJ1 : Start poloidal angle of antenna
THJ2 : End poloidal angle of antenna
PHJ1 : Start toroidal angle of antenna
PHJ2 : End toroidal angle of antenna
NAMAX=1
DO NA=1,NAM
AJ(NA)= 1.D0APH(NA) = 0.D0
THJ1(NA) =-45.D0
THJ2(NA) = 45.D0
PHJ1(NA) = 0.D0
PHJ2(NA) = 0.D0
ENDDO

C
C
C
C
C
C

*** MESH PARAMETERS ***NRMAX : Number of radial mesh points
NTHMAX: Number of poloidal mesh points
NPHMAX : Number of toroidal mesh points
42

(A/m)
(degree)(degree)(degree)(degree)(degree)C
NRMAX
NTHMAXNPHMAX
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

= 50
= 1
= 1

*** CONTROL PARAMETERS ***NPRINT: Control print output
0: No print out
1: Minimum print out (without input data)
2: Minimum print out (with input data)
3: Standard print out
4: More print out
NGRAPH: Control graphic output
0: No graphic out
1: Standard graphic out (2D: Coutour)
2: Standard graphic out (2D: Paint)
3: Standard graphic out (2D: Bird’s eye)
MODELJ: Control antenna current model
0: Real antenna
1: Real antenna
2: Poloidal current
3: Toroidal current
2X: Vacuum eigen mode, poloidal current
3X: Vacuum eigen mode, toroidal current
0: Vacuum-1: MHD plasma-2: Cold plasma-3: Hot plasma (No FLR)
-4: Hot plasma (Cold FLR)
-5: Hot plasma (FLR)
8: TASK/DP (No FLR)
9: TAKS/DP (Cold FLR)
MODELA: Control alpha particle contribution
0: No alpha effect
1: Precession of alpha particles
2: Precession of electrons
3: Precession of both alpha particles and electrons
4: Calculate alpha particle density using slowing down
MODELK: Control mode number cutoff
0: No cutoff
1: With cutoff (this should not be used)
MODELM: Control matrix solver
0: BANDCD
43

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

1: BCGCDB
2: CGSCDB
3: BCGSTAB
4: BANDCDM
5: BCGCDBM
6: CGSCDBM
7: BSTABCDBM
8: BANDCDBM
9: BCGCDBMA
10: CGSCDBMA
11: BSTABCDBMA
12: BANDCDBMODELW: Control writting a data of absorped power
0: Not writting
1: Writting
RHOMIN: rho at minimum q for reversed shear
QMIN : q minimum for reversed shear
RHOITB: rho at ITBNPRINT = 2
NGRAPH = 1
MODELJ = 0
MODELA = 0
MODELK = 0
MODELM = 2
MODELW = 0

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

*** EIGEN VALUE PARAMETERS ***FRMIN : Minimum real part of frequency in amplitude scan
FRMAX : Maximum real part of frequency in amplitude scan
FIMIN : Minimum imag part of frequency in amplitude scan
FIMAX : Maximum imag part of frequency in amplitude scan
FI0
: Imag part of frequency in 1D amplitude scan
NGFMAX: Number of real freq mesh in 1D amplitude scan
NGXMAX: Number of real freq mesh in 2D amplitude scan
NGYMAX: Number of imag freq mesh in 2D amplitude scan
SCMIN : Minimum value in parameter scan
SCMAX : Maximum value in parameter scan
NSCMAX: Number of mesh in parameter scan

44

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

LISTEG: Listing in parameter scan
FRINI : Initial real part of frequency in Newton method
FIINI : Initial imag part of frequency in Newton method
DLTNW : Step size in evaluating derivatives in Newton method
EPSNW : Convergence criterion in Newton method
LMAXNW: Maximum iteration count in Newton method
LISTNW: Listing in Newton method
MODENW: Type of Newton method
NCONT : Number of contour lines
ILN1 : Line type of lower contours
IBL1 : Line boldness of lower contours
ICL1 : Line color of lower contours
ILN2 : Line type of higher contours
IBL2 : Line boldness of higher contours
ICL2 : Line color of higher contours
FRMIN = 0.1D0FRMAX = 1.D0
FIMIN =-0.1D0
FIMAX = 0.1D0
FI0= 0.D0C
FRINI = DBLE(CRF)
FIINI = DIMAG(CRF)
C
NGFMAX= 11
NGXMAX= 11
NGYMAX= 11C
SCMIN = 0.1D0
SCMAX = 1.D0
NSCMAX= 11C
LISTEG= 1C
DLTNW = 1.D-6
EPSNW = 1.D-6
LMAXNW= 10
LISTNW= 1
MODENW= 0
C
45

C
C
C
C
C

*** ALFVEN FREQUENCY PARAMETERS ***WAEMIN : Minimum frequency in Alfven frequency scan
WAEMAX : Maximum frequency in Alfven frequency scan
WAEMIN = 0.001D0
WAEMAX = 0.200D046