tee_mpi_module.f90
Go to the documentation of this file.
1 !--
2 !----------------------------------------------------------------------
3 ! Copyright(c) 2009-2022 SPMDODEL Development Group. All rights reserved.
4 !----------------------------------------------------------------------
5 !表題 tee_mpi_module
6 !
7 ! spml/tee_module モジュールは平行平板間での 3 次元流体運動を
8 ! スペクトル法によって数値計算するための Fortran90 関数を提供する
9 ! ものである.
10 !
11 ! 水平方向にフーリエ数変換および上下の境界壁を扱うための
12 ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
13 ! 関数を提供する.
14 !
15 ! 内部で ee_module, at_module を用いている. 最下部ではフーリエ
16 ! およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
17 ! サブルーチンを用いている.
18 !
19 !
20 !履歴 2009/12/19 竹広真一 新規作成
21 ! 2010/03/10 佐々木洋平 threadprivate 削除(コンパイラ依存)
22 ! 2010/03/27 竹広真一 goto 文を削除
23 ! 2011/12/05 竹広真一 tee_Div_zyx_zyx_zyx を public 変数に
24 ! 2012/07/08 竹広真一 tee_LaplaPol2PolTau_tee を追加
25 ! 2012/08/24 竹広真一 tee_LaplaPol2Pol_tee public に追加
26 ! 2013/08/20 竹広真一 gnu fortran 対応
27 ! 2016/09/22 竹広真一 tee_module を mpi 並列化
28 ! 2020/08/16 竹広真一 変数 ip を save し忘れ
29 ! 2020/12/19 竹広真一 zef_PoloidalEnergySpectrum_tef bug fix
30 ! 2021/02/18 竹広真一 zef_ToroidalEnstropySpectrum_tef,
31 ! zef_PoloidalEnstropySpectrum_tef を追加
32 ! 2022/03/29 竹広真一 tef_LaplaInvDD_tef, tef_LaplaInvNN_tef を追加
33 !
34 !凡例
35 ! データ種類と index
36 ! x : 水平(X) y : 水平(Y) v : 水平(Y,MPI 分割) z : 鉛直
37 ! e : フーリエ変換スペクトル
38 ! f : フーリエ変換スペクトル(MPI 分割
39 ! l : フーリエ関数スペクトル(X 方向波数)
40 ! m : フーリエ関数スペクトル(Y 方向波数)
41 ! t : チェビシェフ関数スペクトル
42 ! a : 任意の次元
43 !
44 ! zvx : 3 次元格子点データ(MPI 分割)
45 ! vx : 水平 2 次元格子点データ(MPI 分割)
46 ! zv : 鉛直 2 次元格子点データ(MPI 分割)
47 ! zx : 鉛直 2 次元格子点データ
48 !
49 ! zef : 水平スペクトル鉛直格子点データ(MPI 分割)
50 ! tef : スペクトルデータ鉛直チェビシェフデータ(MPI 分割)
51 !
52 !++
53 module tee_mpi_module
54  !
55  != tee_mpi_module
56  !
57  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
58  ! Version:: $Id$
59  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
60  !
61  !== 概要
62  !
63  ! spml/tee_mpi_module モジュールは平行平板間での 3 次元流体運動を
64  ! スペクトル法によって数値計算するための Fortran90 関数を提供する
65  ! ものである.
66  !
67  ! 水平方向にフーリエ数変換および上下の境界壁を扱うための
68  ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
69  ! 関数を提供する.
70  !
71  ! 内部で ee_mpi_module, at_module を用いている. 最下部ではフーリエ
72  ! およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
73  ! サブルーチンを用いている.
74  !
75  !== 関数・変数の名前と型について
76  !
77  !=== 命名法
78  !
79  ! * 関数名の先頭 (tef_, zvx_, zef_, ef_, vx_, x_, v_, z_, a_) は,
80  ! 返す値の形を示している.
81  ! tef_ :: スペクトルデータ(2 重フーリエ・チェビシェフ変換)
82  ! zvx_ :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
83  ! zef_ :: 水平スペクトル, 鉛直格子点データ
84  ! e2a_ :: 1 次元化した水平スペクトル, 任意座標データ
85  ! aef_ :: 任意座標データ, 水平スペクトルデータ
86  !
87  ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla,..)
88  ! は, その関数の作用を表している.
89  !
90  ! * 関数名の最後 (tef_, zvx_, zef_, ef_, vx_, x_, v_, z_, a_) は, 入力変数の
91  ! 形がスペクトルデータおよび格子点データであることを示している.
92  ! _tef :: スペクトルデータ(2 重フーリエ・チェビシェフ変換)
93  ! _zvx :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
94  ! _zef :: 水平スペクトル, 鉛直格子点データ
95  ! _e2a :: 1 次元化した水平スペクトル, 任意座標データ
96  ! _aef :: 任意座標データ, 水平スペクトルデータ
97  !
98  !=== 各データの種類の説明
99  !
100  ! * zvx : 3 次元格子点データ(鉛直, 水平 2 次元, MPI 分割)
101  ! * 変数の種類と次元は real(8), dimension(0:km,jc(ip),0:im-1).
102  ! * im, km はそれぞれ水平 X, 鉛直 Z 座標の格子点数である.
103  ! サブルーチン tee_mpi_Initial にてあらかじめ設定しておく.
104  ! * jc(ip) は tee_mpi_initial にて全 Y 格子点数 jm を設定後に使用できる.
105  !
106  ! * tef : スペクトルデータ
107  ! * 変数の種類と次元は real(8), dimension(0:nm,-mm:mm,2*lc(ip)).
108  ! * mm は Y 方向最大波数, nm はチェビシェフ多項式の最大次数
109  ! であり, サブルーチン tee_Initial にてあらかじめ設定しておく.
110  ! * lc(ip) は tee_mpi_initial にて全 Y 方向最大波数 lm を設定後に使用できる.
111  !
112  ! * zef : 水平スペクトル, 鉛直格子点データ.
113  ! * 変数の種類と次元は real(8), dimension(0:km,-mm:mm,2*lc(ip)).
114  !
115  ! * e2a : 1 次元化した水平スペクトル, 任意座標データ
116  ! * 変数の種類と次元は real(8), dimension((2*mm+1)*2*lc(ip),:)
117  !
118  ! * aef : 任意座標データ, 水平スペクトルデータ
119  ! * 変数の種類と次元は real(8), dimension(:,-mm:mm,2*lc(ip)).
120  !
121  ! * tef_ で始まる関数が返す値はスペクトルデータに同じ.
122  !
123  ! * zvx_ で始まる関数が返す値は 3 次元格子点データに同じ.
124  !
125  ! * zef_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ.
126  !
127  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
128  ! 微分などを作用させたデータをスペクトル変換したものことである.
129  !
130  !
131  !== 変数・手続き群の要約
132  !
133  !==== 初期化
134  !
135  ! tee_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
136  !
137  !==== 座標変数
138  !
139  ! x_X, v_Y, z_Z :: 格子点座標(水平 X,Y, 鉛直 Z 座標)を
140  ! 格納した1 次元配列
141  ! x_X_Weight, v_X_Weight, z_Z_Weight :: 重み座標を格納した 1 次元配列
142  ! zvx_X, zvx_Y, zyx_Z :: 格子点データの水平鉛直座標(X,Y,Z)
143  ! (格子点データ型 3 次元配列)
144  ! vx_X, vx_Y :: 格子点データの水平座標(X,Y)
145  ! zv_Z, zv_Y :: 格子点データの鉛直水平座標(Z,Y)
146  ! zx_Z, zx_X :: 格子点データの鉛直水平座標(Z,X)
147  ! (格子点データ型 2 次元配列)
148  !
149  !==== 基本変換
150  !
151  ! zvx_tee, tee_zvx :: スペクトルデータと 3 次元格子データの間の変換
152  ! (2 重フーリエ, チェビシェフ変換)
153  ! zvx_zee, zee_zvx :: 3 次元格子データと水平スペクトル・鉛直格子データとの間
154  ! の変換 (2 重フーリエ変換)
155  ! zef_tee, tee_zef :: スペクトルデータと水平スペクトル・鉛直格子データとの間
156  ! の変換 (チェビシェフ変換)
157  ! ee_vx, vx_ee :: スペクトルデータと 2 次元水平格子データの間の変換
158  ! (2 重フーリエ変換)
159  ! az_at, at_az :: 同時に複数個行う (チェビシェフ変換)格子データと
160  ! チェビシェフデータの間の変換を
161  ! e2a_aef, aef_e2a :: 水平スペクトル軸を 1 次元化し転置, 転置 ２次元化する.
162  !
163  !==== 微分
164  !
165  ! tef_Dx_tef :: スペクトルデータに動径微分∂/∂x を作用させる
166  ! tef_Dy_tef :: スペクトルデータに動径微分∂/∂y を作用させる
167  ! tef_Dz_tef :: スペクトルデータに動径微分∂/∂z を作用させる
168  !
169  ! tef_Lapla_tef :: スペクトルデータにラプラシアンを作用させる
170  ! tef_LaplaH_tef :: スペクトルデータに水平ラプラシアンを作用させる
171  ! tef_LaplaHInv_tef :: スペクトルデータに逆水平ラプラシアンを作用させる
172  ! tef_LaplaInvDD_tef :: スペクトルデータに逆ラプラシアンを作用させる(ディリクレ)
173  ! tef_LaplaInvNN_tef :: スペクトルデータに逆ラプラシアンを作用させる(ノイマン)
174  !
175  ! tef_Div_zvx_zvx_zvx :: ベクトル成分である 3 つの格子データに
176  ! 発散を作用させる
177  !
178  !==== トロイダルポロイダル計算用微分
179  !
180  ! tef_ZRot_zvx_zvx :: ベクトル v の渦度と動径ベクトル r の内積
181  ! z・(▽×v) を計算する
182  ! tef_ZRotRot_zvx_zvx_zvx :: ベクトルの v の z・(▽×▽×v) を計算する
183  ! tef_Potential2Vector :: トロイダルポロイダルポテンシャルから
184  ! ベクトル場を計算する
185  ! tef_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される
186  ! 非発散ベクトル場の回転の各成分を計算する
187  !
188  !==== ポロイダル/トロイダルモデル用スペクトル解析
189  !
190  ! zef_ToroidalEnergySpectrum_tef, zk_ToroidalEnergySpectrum_tef ::
191  ! トロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
192  ! zef_PoloidalEnergySpectrum_tef, zk_PoloidalEnergySpectrum_tef ::
193  ! ポロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
194  ! zef_ToroidalEnstrophySpectrum_tef, zk_ToroidalEnstrophySpectrum_tef ::
195  ! トロイダルポテンシャルからエンストロフィーのフーリエ各成分を計算する
196  ! zef_PoloidalEnstrophySpectrum_tef, zk_PoloidalEnstrophySpectrum_tef ::
197  ! ポロイダルポテンシャルからエンストロフィーのフーリエ各成分を計算する
198  !
199  !==== 境界値問題
200  !
201  ! tef_BoundariesTau, tef_BoundariesGrid, tef_Boundaries ::
202  ! ディリクレ, ノイマン境界条件を適用する(タウ法, 選点法)
203  !
204  ! tef_TorBoundariesTau, tef_TorBoundariesGrid, tef_TorBoundaries ::
205  ! 速度トロイダルポテンシャルの境界条件を適用する(タウ法,選点法)
206  !
207  ! tef_LaplaPol2PolTau_tef, zee_LaplaPol2Pol_zee, tef_LaplaPol2PolGrid_tef ::
208  ! 速度ポロイダルポテンシャルΦを▽^2Φから求める
209  ! (入出力がそれぞれチェビシェフ格子点,チェビシェフ係数)
210  !
211  ! tef_TorMagBoundariesTau, tef_TorMagBoundariesGrid, tef_TorMagBoundaries ::
212  ! 磁場トロイダルポテンシャルの境界条件を適用する(タウ法, 選点法)
213  !
214  ! tef_PolMagBoundariesTau, tef_PolMagBoundariesGrid, tef_PolMagBoundaries ::
215  ! 磁場トロイダルポテンシャル境界の境界条件を適用する(タウ法, 選点法)
216  !
217  !==== 積分・平均(3 次元データ)
218  !
219  ! IntZYX_zvx, AvrZYX_zvx :: 3 次元格子点データの全領域積分および平均
220  ! z_IntYX_zvx, z_AvrYX_zvx :: 3 次元格子点データの水平積分および平均
221  ! v_IntZX_zvx, v_AvrZX_zvx :: 3 次元格子点データのZX積分および平均
222  ! z_IntZY_zvx, z_AvrZY_zvx :: 3 次元格子点データのZY積分および平均
223  ! zv_IntX_zvx, zv_AvrX_zvx :: 3 次元格子点データの水平X方向積分および平均
224  ! zx_IntY_zvx, zx_AvrY_zvx :: 3 次元格子点データの水平Y方向積分および平均
225  ! zx_IntZ_zvx, zx_AvrZ_zvx :: 3 次元格子点データの鉛直方向積分および平均
226  !
227  !==== 積分・平均(2 次元データ)
228  !
229  ! IntYX_vx, AvrYX_vx :: 2 次元格子点データの水平積分および平均
230  ! IntZX_zx, AvrZX_zx :: 2 次元(ZX)格子点データのZX積分および平均
231  ! IntZY_zv, AvrZY_zv :: 2 次元(ZY)格子点データのZY積分および平均
232  ! y_IntX_vx, y_AvrX_vx :: 水平 2 次元格子点データのX方向積分および平均
233  ! x_IntY_vx, x_AvrY_vx :: 水平2 次元格子点データのY方向積分および平均
234  ! z_IntX_zx, z_AvrX_zx :: 2 次元(ZX)格子点データのX方向積分および平均
235  ! x_IntZ_zx, x_AvrZ_zx :: 2 次元(ZX)格子点データのZ方向積分および平均
236  ! z_IntY_zv, z_AvrY_zv :: 2 次元(ZY)格子点データのY方向積分および平均
237  ! y_IntZ_zv, y_AvrZ_zv :: 2 次元(ZY)格子点データのZ方向積分および平均
238  !
239  !==== 積分・平均(1 次元データ)
240  !
241  ! IntX_x, AvrX_x :: 1 次元(X)格子点データのX方向積分および平均
242  ! IntY_v, AvrY_v :: 1 次元(Y)格子点データのY方向積分および平均
243  ! IntZ_z, AvrZ_z :: 1 次元(Z)格子点データのZ方向積分および平均
244  !
245  !==== 補間計算
246  !
247  ! Interpolate_tee :: スペクトルデータから任意の点の値を補間する.
248  !
249  use dc_message
250  use mpi
251  use lumatrix
252  use ee_mpi_module, ls=>ks, le=>ke, lc=>kc, lf_l=>kf_k
253  use at_module, only : z_z => g_x, z_z_weight => g_x_weight, &
254  at_initial, at_az => at_ag, az_at => ag_at, &
255  t_z => t_g, z_t => g_t, &
256  t_dz_t => t_dx_t, at_dz_at => at_dx_at, &
257  a_interpolate_at, &
258  at_boundaries, t_boundaries, &
259  at_boundariestau_dd, at_boundariestau_dn, &
260  at_boundariestau_nd, at_boundariestau_nn, &
261  at_boundariesgrid_dd, at_boundariesgrid_dn, &
262  at_boundariesgrid_nd, at_boundariesgrid_nn
263 
264  implicit none
265  private
266 
267  public tee_mpi_initial
268 
269  public ls, le, lc
270  public js, je, jc
271  public lf_l
272 
273  public x_x, x_x_weight
274  public v_y, v_y_weight
275  public z_z, z_z_weight
276 
277  public vx_x, vx_y, zv_z, zv_y, zx_x, zx_z
278  public zvx_x, zvx_y, zvx_z
279  public zef_z
280 
281  public tee_vmiss
282 
283  public ef_vx, vx_ef
284  public at_dz_at, t_dz_t, az_at, at_az, z_t, t_z
285  public zvx_tef, tef_zvx, zvx_zef, zef_zvx, zef_tef, tef_zef
286  public e2z_e2t, e2t_e2z
287  public aef_e2a, e2a_aef, ef_e2, e2_ef
288  public tef_dx_tef, tef_dy_tef, tef_dz_tef
289  public tef_lapla_tef, tef_laplah_tef, tef_laplahinv_tef
290  public tef_laplainvdd_tef, tef_laplainvnn_tef
291  public tef_laplainvtaudd_tef, tef_laplainvtaunn_tef
292  public tef_div_zvx_zvx_zvx
293 
294  public zv_intx_zvx, zx_inty_zvx, vx_intz_zvx
295  public x_intzy_zvx, v_intzx_zvx, z_intyx_zvx
296  public intzyx_zvx
297 
298  public x_inty_vx, v_intx_vx, intyx_vx
299  public z_inty_zv, v_intz_zv, intzy_zv
300  public z_intx_zx, x_intz_zx, intzx_zx
301  public intx_x, inty_v, intz_z
302 
303  public zv_avrx_zvx, zx_avry_zvx, vx_avrz_zvx
304  public x_avrzy_zvx, v_avrzx_zvx, z_avryx_zvx
305  public avrzyx_zvx
306 
307  public x_avry_vx, v_avrx_vx, avryx_vx
308  public z_avry_zv, v_avrz_zv, avrzy_zv
309  public z_avrx_zx, x_avrz_zx, avrzx_zx
310  public avrx_x, avry_v, avrz_z
311 
312  public tef_zrot_zvx_zvx, tef_zrotrot_zvx_zvx_zvx
313  public tef_potential2vector, tef_potential2rotation
314 
315  public interpolate_tef
316 
317  public zef_toroidalenergyspectrum_tef ! nz_ToroidalEnergySpectrum_wt
318  public zef_poloidalenergyspectrum_tef ! nz_PoloidalEnergySpectrum_wt
319  public zef_toroidalenstrophyspectrum_tef ! nz_ToroidalEnstrophySpectrum_wt
320  public zef_poloidalenstrophyspectrum_tef ! nz_PoloidalEnstrophySpectrum_wt
321 
322  public tef_boundaries, zef_laplapol2pol_zef
323  public tef_torboundaries, tef_laplapol2pol_tef
324  public tef_tormagboundaries, tef_polmagboundaries
325 
326  public tef_boundariestau, tef_laplapol2poltau_tef, tef_torboundariestau
327  public tef_tormagboundariestau, tef_polmagboundariestau
328 
329  public tef_boundariesgrid, tef_laplapol2polgrid_tef, tef_torboundariesgrid
330  public tef_tormagboundariesgrid, tef_polmagboundariesgrid
331 
332  public at_boundaries, t_boundaries
333 
334  interface tef_laplainvdd_tef
335  module procedure tef_laplainvtaudd_tef
336  end interface
337 
338  interface tef_laplainvnn_tef
339  module procedure tef_laplainvtaunn_tef
340  end interface
341 
342  interface tef_boundaries
343  module procedure tef_boundariestau
344  end interface
345 
346  interface tef_torboundaries
347  module procedure tef_torboundariestau
348  end interface
349 
350  interface tef_laplapol2pol_tef
351  module procedure tef_laplapol2poltau_tef
352  end interface
353 
354  interface tef_tormagboundaries
355  module procedure tef_tormagboundariestau
356  end interface
357 
358  interface tef_polmagboundaries
359  module procedure tef_polmagboundariestau
360  end interface
361 
362  integer :: im=32, jm=32, km=16 ! 格子点の設定(X,Y,Z)
363  integer :: lm=10, mm=10, nm=16 ! 切断波数の設定(水平X,Y, 鉛直Z)
364  real(8) :: xl, yl, zl ! 領域の大きさ(水平X,Y, 鉛直Z)
365 
366  real(8), parameter :: pi=3.1415926535897932385d0
367 
368  real(8), dimension(:,:,:), allocatable :: zvx_x, zvx_y, zvx_z ! 座標
369  real(8), dimension(:,:), allocatable :: zv_z, zv_y ! 座標
370  real(8), dimension(:,:), allocatable :: zx_z, zx_x ! 座標
371  real(8), dimension(:,:,:), allocatable :: zef_z ! 座標
372 
373  real(8) :: tee_vmiss = -999.0 ! 欠損値
374 
375 ! integer :: np ! MPI プロセス数
376  integer :: ip ! MPI プロセス ID
377 
378  save im, jm, km, lm, mm, nm, xl, yl, zl, ip
379 
380  contains
381  !--------------- 初期化 -----------------
382  subroutine tee_mpi_initial(i,j,k,l,m,n,xmin,xmax,ymin,ymax,zmin,zmax)
383  !
384  ! スペクトル変換の格子点数, 波数, 各座標の範囲を設定する.
385  !
386  ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
387  ! しなければならない.
388  !
389  integer,intent(in) :: i ! 格子点数(水平X)
390  integer,intent(in) :: j ! 格子点数(水平Y)
391  integer,intent(in) :: k ! 格子点数(鉛直Z)
392  integer,intent(in) :: l ! 切断波数(水平X波数)
393  integer,intent(in) :: m ! 切断波数(水平Y波数)
394  integer,intent(in) :: n ! 切断波数(鉛直Z波数)
395 
396  real(8),intent(in) :: xmin, xmax ! X 方向領域
397  real(8),intent(in) :: ymin, ymax ! Y 方向領域
398  real(8),intent(in) :: zmin, zmax ! Z 方向領域
399 
400  integer :: ierr
401 
402  im = i ; jm = j ; km = k
403  lm = l ; mm = m ; nm = n
404  xl = xmax - xmin
405  yl = ymax - ymin
406  zl = zmax - zmin
407 
408  CALL mpi_comm_rank(mpi_comm_world,ip,ierr)
409 
410  call ee_mpi_initial(im,jm,lm,mm,xmin,xmax,ymin,ymax)
411 
412  call at_initial(km,nm,zmin,zmax)
413 
414  allocate(zvx_x(0:km,jc(ip),0:im-1))
415  allocate(zvx_y(0:km,jc(ip),0:im-1))
416  allocate(zvx_z(0:km,jc(ip),0:im-1))
417 
418  allocate(zv_z(0:km,jc(ip)))
419  allocate(zv_y(0:km,jc(ip)))
420  allocate(zx_z(0:km,0:im-1))
421  allocate(zx_x(0:km,0:im-1))
422 
423  allocate(zef_z(0:km,-mm:mm,2*lc(ip)))
424 
425  zvx_x = spread(vx_x,1,km+1)
426  zvx_y = spread(vx_y,1,km+1)
427  zvx_z = spread(spread(z_z,2,jc(ip)),3,im)
428 
429  zv_z = spread(z_z,2,jc(ip))
430  zv_y = spread(v_y,1,km+1)
431  zx_z = spread(z_z,2,im)
432  zx_x = spread(x_x,1,km+1)
433 
434  zef_z = spread(spread(z_z,2,2*mm+1),3,2*lc(ip))
435 
436  call messagenotify('M','tee_mpi_initial', &
437  'tee_mpi_module (2022/03/29) is initialized')
438 
439  end subroutine tee_mpi_initial
440 
441  !--------------- 基本変換 -----------------
442 
443  function zvx_tef(tef)
444  !
445  ! スペクトルデータから 3 次元格子点データへ(逆)変換する.
446  !
447  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
448  !(in) 2 重フーリエチェビシェフスペクトルデータ
449  real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_tef
450  !(out) 3 次元格子点データ
451 
452  zvx_tef = zvx_zef(zef_tef(tef))
453 
454  end function zvx_tef
455 
456  function tef_zvx(zvx)
457  !
458  ! 3 次元格子点データからスペクトルデータへ(正)変換する.
459  !
460  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
461  !(in) 3 次元格子点データ
462  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_zvx
463  !(out) 2 重フーリエチェビシェフスペクトルデータ
464 
465  tef_zvx = tef_zef(zef_zvx(zvx))
466 
467  end function tef_zvx
468 
469  function zvx_zef(zef)
470  !
471  ! 水平スペクトル・鉛直格子点データから 3 次元格子点データへ(逆)変換する.
472  !
473  real(8), dimension(0:km,-mm:mm,2*lc(ip)), intent(in) :: zef
474  !(in) 2 次元水平スペクトル・鉛直格子点データ
475  real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_zef
476  !(out) 3 次元格子点データ
477  real(8), dimension(-mm:mm,2*lc(ip)) :: ef
478  real(8), dimension(jc(ip),0:im-1) :: vx
479 
480  integer :: k
481 
482  do k = 0, km
483  !zyx_zef(k,:,:) = yx_ee(zef(k,:,:))
484  !
485  ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
486  !
487  ef = zef(k,-mm:mm,1:2*lc(ip))
488  vx = vx_ef(ef)
489  zvx_zef(k,:,:) = vx
490  enddo
491 
492  end function zvx_zef
493 
494  function zef_zvx(zvx)
495  !
496  ! 3 次元格子データから水平スペクトル・鉛直格子点データへ(正)変換する.
497  !
498  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
499  !(in) 3 次元格子点データ
500  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_zvx
501  !(out) 2 次元スペクトル・鉛直格子点データ
502  real(8), dimension(-mm:mm,2*lc(ip)) :: ef
503  real(8), dimension(jc(ip),0:im-1) :: vx
504 
505  integer :: k
506 
507  do k = 0, km
508  !zef_zyx(k,:,:) = ee_yx(zyx(k,:,:))
509  !
510  ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
511  !
512  vx = zvx(k,1:jc(ip),0:im-1)
513  ef = ef_vx(vx)
514  zef_zvx(k,:,:) = ef
515  enddo
516 
517  end function zef_zvx
518 
519  function zef_tef(tef)
520  !
521  ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する.
522  !
523  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
524  !(in) 2 次元水平スペクトル鉛直チェビシェフスペクトルデータ
525  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_tef
526  !(out) 2 次元水平スペクトル・鉛直格子点データ
527 
528  zef_tef = aef_e2a(e2z_e2t(e2a_aef(tef)))
529 
530  end function zef_tef
531 
532  function tef_zef(zef)
533  !
534  ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する.
535  !
536  real(8), dimension(0:km,-mm:mm,2*lc(ip)), intent(in) :: zef
537  !(in) 2 次元水平スペクトル・鉛直格子点データ
538  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_zef
539  !(out) 2 次元水平鉛直チェビシェフスペクトルデータ
540 
541  tef_zef = aef_e2a(e2t_e2z(e2a_aef(zef)))
542 
543  end function tef_zef
544 
545  function e2z_e2t(e2t)
546  !
547  ! 鉛直スペクトルを格子点に変換する
548  !
549  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t
550  !(in) 1 次元化された水平スペクトル・鉛直スペクトル
551  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_e2t
552  !(out) 1 次元化された水平スペクトル・鉛直格子座標
553 
554  e2z_e2t = az_at(e2t)
555 
556  end function e2z_e2t
557 
558  function e2t_e2z(e2z)
559  !
560  ! 鉛直格子点をスペクトルに変換する
561  !
562  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z
563  !(in) 1 次元化された水平スペクトル・鉛直格子座標
564  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_e2z
565  !(out) 1 次元化された水平スペクトル・鉛直スペクトル
566 
567  e2t_e2z = at_az(e2z)
568 
569  end function e2t_e2z
570 
571  function e2a_aef(aef)
572  !
573  ! 水平スペクトルを 1 次元化転置する.
574  !
575  real(8), dimension(:,-mm:,:), intent(in) :: aef
576  !(in) 任意座標・水平 2 次元スペクトルデータ dimension(:,-mm:mm,2*lc(ip))
577  real(8), dimension(2*lc(ip)*(2*mm+1),size(aef,1)) :: e2a_aef
578  !(out) 1 次元化された水平スペクトル・任意座標データ
579 
580  integer :: j,k,l,m
581 
582  if ( size(aef,2) /= 2*mm+1 ) &
583  call messagenotify('E','e2a_aef',&
584  '2nd dimension of input data invalid')
585  if ( size(aef,3) /= 2*lc(ip) ) &
586  call messagenotify('E','e2a_aef',&
587  '3rd dimension of input data invalid')
588 
589  !e2a_aee = transpose(reshape(aee,(/size(aee,1),2*lc(ip)*(2*mm+1)/)))
590  !
591  ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
592  !
593 
594  do k=1,size(aef,1)
595  do l=1,2*lc(ip)
596  do m=-mm,mm
597  j = (m+mm+1) + (2*mm+1)*(l-1)
598  e2a_aef(j,k) = aef(k,m,l)
599  enddo
600  enddo
601  enddo
602 
603  end function e2a_aef
604 
605  function aef_e2a(e2a)
606  !
607  ! 水平スペクトルを転置展開する.
608  !
609  real(8), dimension(:,:),intent(in) :: e2a
610  !(in) 1 次元化された水平スペクトル・任意座標データ
611  ! dimmension((2*mm*1)*(2*lm*1),:)
612  real(8), dimension(size(e2a,2),-mm:mm,2*lc(ip)) :: aef_e2a
613  !(out) 任意座標・水平 2 次元スペクトルデータ
614 
615  integer :: j,k,l,m
616 
617  if ( size(e2a,1) /= (2*mm+1)*2*lc(ip) ) &
618  call messagenotify('E','aef_e2a',&
619  '1st dimension of input data invalid')
620 
621  !aee_e2a = reshape(transpose(e2a),(/size(e2a,2),2*mm+1,2*lc(ip)/))
622  !
623  ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
624  !
625 
626  do k=1,size(e2a,2)
627  do l=1,2*lc(ip)
628  do m=-mm,mm
629  j = (m+mm+1) + (2*mm+1)*(l-1)
630  aef_e2a(k,m,l) = e2a(j,k)
631  enddo
632  enddo
633  enddo
634 
635  end function aef_e2a
636 
637  function e2_ef(ef)
638  !
639  ! 水平スペクトルを 1 次元化転置する.
640  !
641  real(8), dimension(-mm:,:), intent(in) :: ef
642  !(in) 任意座標・水平 2 次元スペクトルデータ dimension(-mm:mm,2*lc(ip))
643  real(8), dimension(2*lc(ip)*(2*mm+1)) :: e2_ef
644  !(out) 1 次元化された水平スペクトル・任意座標データ
645 
646  integer :: j,l,m
647 
648  !e2_ee = reshape(e2a_aee(reshape(ee,(/1,2*mm+1,2*lc(ip)/))), &
649  ! (/lm*(2*mm+1)/))
650  !
651  ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
652  !
653 
654  do l=1,2*lc(ip)
655  do m=-mm,mm
656  j = (m+mm+1)+ (2*mm+1)*(l-1)
657  e2_ef(j) = ef(m,l)
658  enddo
659  enddo
660 
661  end function e2_ef
662 
663  function ef_e2(e2)
664  !
665  ! 水平スペクトルを転置展開する.
666  !
667  real(8), dimension(:),intent(in) :: e2
668  !(in) 1 次元化された水平スペクトル・任意座標データ
669  ! dimmension((2*mm*1)*2*lc(ip),:)
670  real(8), dimension(-mm:mm,2*lc(ip)) :: ef_e2
671  !(out) 任意座標・水平 2 次元スペクトルデータ
672 
673  integer :: j,l,m
674 
675  !ee_e2 = reshape(aee_e2a(reshape(e2,(/(2*mm+1)*2*lc(ip),1/))), &
676  ! (/2*mm+1,2*lc(ip)/))
677  !
678  ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
679  !
680 
681  do l=1,2*lc(ip)
682  do m=-mm,mm
683  j = (m+mm+1) + (2*mm+1)*(l-1)
684  ef_e2(m,l) = e2(j)
685  enddo
686  enddo
687 
688  end function ef_e2
689 
690  !--------------- 微分計算 -----------------
691  function tef_dx_tef(tef)
692  !
693  ! 入力スペクトルデータに水平微分 ∂/∂x を作用する.
694  !
695  ! スペクトルデータのX微分とは, 対応する格子点データにX微分を
696  ! 作用させたデータのスペクトル変換のことである.
697  !
698  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
699  !(in) 入力スペクトルデータ
700 
701  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Dx_tef
702  !(in) X微分されたスペクトルデータ
703 
704  integer :: n
705 
706  do n=0,nm
707  tef_dx_tef(n,:,:) = ef_dx_ef(tef(n,:,:))
708  enddo
709 
710  end function tef_dx_tef
711 
712  function tef_dy_tef(tef)
713  !
714  ! 入力スペクトルデータに水平微分 ∂/∂y を作用する.
715  !
716  ! スペクトルデータのY微分とは, 対応する格子点データにY微分を
717  ! 作用させたデータのスペクトル変換のことである.
718  !
719  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
720  !(in) 入力スペクトルデータ
721 
722  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Dy_tef
723  !(in) X微分されたスペクトルデータ
724 
725  integer :: n
726 
727  do n=0,nm
728  tef_dy_tef(n,:,:) = ef_dy_ef(tef(n,:,:))
729  enddo
730 
731  end function tef_dy_tef
732 
733  function tef_dz_tef(tef)
734  !
735  ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する.
736  !
737  ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を
738  ! 作用させたデータのスペクトル変換のことである.
739  !
740  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
741  !(in) 入力スペクトルデータ
742 
743  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Dz_tef
744  !(in) 鉛直微分されたスペクトルデータ
745 
746  tef_dz_tef = aef_e2a(at_dz_at(e2a_aef(tef)))
747 
748  end function tef_dz_tef
749 
750  function tef_lapla_tef(tef)
751  !
752  ! 入力スペクトルデータにラプラシアン
753  !
754  ! ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2
755  !
756  ! を作用する.
757  !
758  ! スペクトルデータのラプラシアンとは, 対応する格子点データに
759  ! ラプラシアンを作用させたデータのスペクトル変換のことである.
760  !
761  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
762  !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
763 
764  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Lapla_tef
765  !(out) ラプラシアンを作用された 2 次元スペクトルデータ
766 
767  tef_lapla_tef = tef_laplah_tef(tef) + tef_dz_tef(tef_dz_tef(tef))
768 
769  end function tef_lapla_tef
770 
771  function tef_laplah_tef(tef)
772  !
773  ! 入力スペクトルデータに水平ラプラシアン
774  !
775  ! ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2
776  !
777  ! を作用する.
778  !
779  ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに
780  ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである.
781  !
782  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
783  !(in) 入力スペクトルデータ
784 
785  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaH_tef
786  !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ
787 
788  integer :: n
789 
790  do n=0,nm
791  tef_laplah_tef(n,:,:) = ef_lapla_ef(tef(n,:,:))
792  enddo
793 
794  end function tef_laplah_tef
795 
796  function tef_laplahinv_tef(tef)
797  !
798  ! 入力スペクトルデータに逆水平ラプラシアン
799  !
800  ! ▽^-2_H = (∂^2/∂X^2 + ∂^2/∂Y^2)^-1
801  !
802  ! を作用する.
803  !
804  ! スペクトルデータの逆水平ラプラシアンとは, 対応する格子点データに
805  ! 逆水平ラプラシアンを作用させたデータのスペクトル変換のことである.
806  !
807  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
808  !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
809 
810  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaHInv_tef
811  !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ
812 
813  integer :: n
814 
815  do n=0,nm
816  tef_laplahinv_tef(n,:,:) = ef_laplainv_ef(tef(n,:,:))
817  enddo
818 
819  end function tef_laplahinv_tef
820 
821  function tef_div_zvx_zvx_zvx(zvx_VX,zvx_VY,zvx_VZ)
822  !
823  ! べクトル成分である 3 つの格子データに発散を作用させた
824  ! スペクトルデータを返す.
825  !
826  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX
827  !(in) ベクトル場のX成分
828  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY
829  !(in) ベクトル場のY成分
830 
831  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VZ
832  !(in) ベクトル場のZ成分
833 
834  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Div_zvx_zvx_zvx
835  !(out) ベクトル場の発散
836 
837  tef_div_zvx_zvx_zvx = tef_dx_tef(tef_zvx(zvx_vx)) &
838  + tef_dy_tef(tef_zvx(zvx_vy)) &
839  + tef_dz_tef(tef_zvx(zvx_vz))
840 
841  end function tef_div_zvx_zvx_zvx
842 
843  function tef_laplainvtaudd_tef(tef,new)
844  !
845  ! ディリクレ境界条件でΦを▽^2Φから計算する.
846  ! チェビシェフ空間でタウ法による境界条件を適用している.
847  !
848  ! Φを f = ▽^2Φから定める式は
849  !
850  ! ▽^2Φ = f, Φ = 0. at boundaries.
851  !
852  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
853  !
854  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(in) :: tef
855  !(in) 入力▽^2φ分布
856 
857  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaInvTauDD_tef
858  !(out) 出力ポロイダルポテンシャル分布
859 
860  logical, intent(IN), optional :: new
861  !(in) true だと境界条件計算用行列を強制的に新たに作る.
862  ! default は false.
863 
864  real(8), dimension(:,:,:), allocatable :: alu
865  integer, dimension(:,:), allocatable :: kp
866 
867  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
868  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
869  logical :: rigid1 = .true.
870  logical :: rigid2 = .true.
871 
872  logical :: first = .true.
873  logical :: new_matrix = .false.
874  integer :: n
875  save :: alu, kp, first
876 
877  if (.not. present(new)) then
878  new_matrix=.false.
879  else
880  new_matrix=new
881  endif
882 
883  if ( first .OR. new_matrix ) then
884  first = .false.
885 
886  if ( allocated(alu) ) deallocate(alu)
887  if ( allocated(kp) ) deallocate(kp)
888  allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
889 
890  alu = 0.0d0
891  do n=0,nm
892  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
893 
894  ! 各水平波数に関して独立の式
895  alu(:,:,n) = e2a_aef(tef_lapla_tef(aef_e2a(e2t_work)))
896  enddo
897 
898  do n=0,nm
899  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
900  e2z_work = az_at(e2t_work)
901 
902  ! ディリクレ条件. 境界で値が 0
903  alu(:,nm,n) = e2z_work(:,0)
904  alu(:,nm-1,n) = e2z_work(:,km)
905  enddo
906 
907  call ludecomp(alu,kp)
908 
909  call messagenotify('M','tef_LaplaInvTauDD_tef',&
910  'Matrix to apply b.c. newly produced.')
911  endif
912 
913  e2t_work = e2a_aef(tef)
914  e2t_work(:,nm) = 0.0d0 ! ディリクレ条件
915  e2t_work(:,nm-1) = 0.0d0 ! ディリクレ条件
916 
917  tef_laplainvtaudd_tef = aef_e2a(lusolve(alu,kp,e2t_work))
918 
919  end function tef_laplainvtaudd_tef
920 
921  function tef_laplainvtaunn_tef(tef,new)
922  !
923  ! ディリクレ境界条件でΦを▽^2Φから計算する.
924  ! チェビシェフ空間でタウ法による境界条件を適用している.
925  !
926  ! Φを f = ▽^2Φから定める式は
927  !
928  ! ▽^2Φ = f, dΦ/dz = 0. at boundaries.
929  !
930  ! 直流成分は z 平均が 0 である制約を適用する.
931  !
932  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
933  !
934  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(in) :: tef
935  !(in) 入力▽^2φ分布
936 
937  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaInvTauNN_tef
938  !(out) 出力ポロイダルポテンシャル分布
939 
940  logical, intent(IN), optional :: new
941  !(in) true だと境界条件計算用行列を強制的に新たに作る.
942  ! default は false.
943 
944  real(8), dimension(:,:,:), allocatable :: alu
945  integer, dimension(:,:), allocatable :: kp
946 
947  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
948  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
949  logical :: rigid1 = .true.
950  logical :: rigid2 = .true.
951 
952  logical :: first = .true.
953  logical :: new_matrix = .false.
954  integer :: n,k
955  save :: alu, kp, first
956 
957  if (.not. present(new)) then
958  new_matrix=.false.
959  else
960  new_matrix=new
961  endif
962 
963  if ( first .OR. new_matrix ) then
964  first = .false.
965 
966  if ( allocated(alu) ) deallocate(alu)
967  if ( allocated(kp) ) deallocate(kp)
968  allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
969 
970  alu = 0.0d0
971  do n=0,nm
972  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
973 
974  ! 各水平波数に関して独立の式
975  alu(:,:,n) = e2a_aef(tef_lapla_tef(aef_e2a(e2t_work)))
976  enddo
977 
978  do n=0,nm
979  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
980  e2z_work = az_at(at_dz_at(e2t_work))
981 
982  ! ノイマン条件. 境界で微分値が 0
983  alu(:,nm,n) = e2z_work(:,0)
984  alu(:,nm-1,n) = e2z_work(:,km)
985 
986  ! 直流成分が決まらないので 1 点値を与える
987  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
988  e2z_work = az_at(e2t_work)
989  alu(:,0,n) = 0.0d0
990  do k=0,km
991  alu(:,0,n) = alu(:,0,n) + e2z_work(:,k)*z_z_weight(k)
992  end do
993  enddo
994 
995  call ludecomp(alu,kp)
996 
997  call messagenotify('M','tef_LaplaInvTauDD_tef',&
998  'Matrix to apply b.c. newly produced.')
999  endif
1000 
1001  e2t_work = e2a_aef(tef)
1002  e2t_work(:,nm) = 0.0d0 ! ノイマン条件
1003  e2t_work(:,nm-1) = 0.0d0 ! ノイマン条件
1004  e2t_work(:,0) = 0.0d0 ! 直流成分
1005 
1006  tef_laplainvtaunn_tef = aef_e2a(lusolve(alu,kp,e2t_work))
1007 
1008  end function tef_laplainvtaunn_tef
1009 
1010 
1011  !--------------- 積分計算 -----------------
1012  !----(入力データ zvx)---
1013  function zv_intx_zvx(zvx) ! X積分
1014  !
1015  ! 3 次元格子点データのX方向積分.
1016  !
1017  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1018  !(in) 3 次元格子点データ
1019 
1020  real(8), dimension(0:km,jc(ip)) :: zv_IntX_zvx
1021  !(out) X方向積分された 2 次元ZV格子点データ
1022 
1023  integer :: i
1024 
1025  zv_intx_zvx = 0.0d0
1026  do i=0,im-1
1027  zv_intx_zvx(:,:) = zv_intx_zvx(:,:) &
1028  + zvx(:,:,i) * x_x_weight(i)
1029  enddo
1030  end function zv_intx_zvx
1031 
1032  function zx_inty_zvx(zvx)
1033  !
1034  ! 3 次元格子点データのY方向域積分.
1035  !
1036  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1037  !(in) 3 次元格子点データ
1038 
1039  real(8), dimension(0:km,0:im-1) :: zx_IntY_zvx
1040  !(out) Y積分された 2 次元ZV格子点データ.
1041 
1042  real(8), dimension(0:km,0:im-1) :: zx_IntYTmp
1043  integer :: j, ierr
1044 
1045  zx_intytmp = 0.0d0
1046  do j=1,jc(ip)
1047  zx_intytmp(:,:) = zx_intytmp(:,:) &
1048  + zvx(:,j,:) * v_y_weight(j)
1049  enddo
1050 
1051  CALL mpi_allreduce(zx_intytmp,zx_inty_zvx,im*(km+1),mpi_real8, &
1052  mpi_sum,mpi_comm_world,ierr)
1053 
1054  end function zx_inty_zvx
1055 
1056  function vx_intz_zvx(zvx) ! Z積分
1057  !
1058  ! 3 次元格子点データのZ方向積分.
1059  !
1060  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1061  !(in) 3 次元格子点データ
1062 
1063  real(8), dimension(jc(ip),0:im-1) :: vx_IntZ_zvx
1064  !(out) Z積分された 2 次元VX(水平)格子点データ
1065 
1066  integer :: k
1067 
1068  vx_intz_zvx = 0.0d0
1069  do k=0,km
1070  vx_intz_zvx(:,:) = vx_intz_zvx(:,:) &
1071  + zvx(k,:,:) * z_z_weight(k)
1072  enddo
1073  end function vx_intz_zvx
1074 
1075  function x_intzy_zvx(zvx)
1076  !
1077  ! 3 次元格子点データのZV積分
1078  !
1079  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1080  !(in) 3 次元ZVX格子点データ
1081 
1082  real(8), dimension(0:im-1) :: x_IntZY_zvx
1083  !(out) ZV(子午面)積分された 1 次元X格子点データ
1084 
1085  real(8), dimension(0:im-1) :: x_IntZYTmp
1086  integer :: j, k, ierr
1087 
1088  x_intzytmp = 0.0d0
1089  do j=1,jc(ip)
1090  do k=0,km
1091  x_intzytmp = x_intzytmp &
1092  + zvx(k,j,:) * v_y_weight(j) * z_z_weight(k)
1093  enddo
1094  enddo
1095 
1096  CALL mpi_allreduce(x_intzytmp,x_intzy_zvx,im,mpi_real8, &
1097  mpi_sum,mpi_comm_world,ierr)
1098 
1099  end function x_intzy_zvx
1100 
1101  function v_intzx_zvx(zvx)
1102  !
1103  ! 3 次元格子点データのZX積分.
1104  !
1105  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1106  !(in) 3 次元ZVX格子点データ
1107 
1108  real(8), dimension(jc(ip)) :: v_IntZX_zvx
1109  !(out) ZX積分された 1 次元Y格子点データ
1110 
1111  integer :: i, k
1112 
1113  v_intzx_zvx = 0
1114  do i=0,im-1
1115  do k=0,km
1116  v_intzx_zvx = v_intzx_zvx &
1117  + zvx(k,:,i) * x_x_weight(i) * z_z_weight(k)
1118  enddo
1119  enddo
1120  end function v_intzx_zvx
1121 
1122  function z_intyx_zvx(zvx) ! VX(水平)積分
1123  !
1124  ! 3 次元格子点データのVX(水平)積分
1125  !
1126  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1127  !(in) 3 次元ZVX格子点データ
1128 
1129  real(8), dimension(0:km) :: z_IntYX_zvx
1130  !(out) VX(水平)積分された 1 次元Z格子点データ
1131 
1132  real(8), dimension(0:km) :: z_IntYXTmp
1133  integer :: i, j, ierr
1134 
1135  z_intyxtmp = 0
1136  do j=1,jc(ip)
1137  do i=0,im-1
1138  z_intyxtmp = z_intyxtmp &
1139  + zvx(:,j,i) * x_x_weight(i) * v_y_weight(j)
1140  enddo
1141  enddo
1142 
1143  CALL mpi_allreduce(z_intyxtmp,z_intyx_zvx,km+1,mpi_real8, &
1144  mpi_sum,mpi_comm_world,ierr)
1145 
1146  end function z_intyx_zvx
1147 
1148  function intzyx_zvx(zvx) ! ZVX(全領域)積分
1149  !
1150  ! 3 次元格子点データのZVX(全領域)積分
1151  !
1152  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1153  !(in) 3 次元格子点データ
1154 
1155  real(8) :: IntZYX_zvx
1156  !(out) 全領域積分値
1157 
1158  real(8) :: IntZYXTmp
1159  integer :: i, j, k, ierr
1160 
1161  intzyxtmp = 0
1162  do i=0,im-1
1163  do j=1,jc(ip)
1164  do k=0,km
1165  intzyxtmp = intzyxtmp &
1166  + zvx(k,j,i) * x_x_weight(i) &
1167  * v_y_weight(j) * z_z_weight(k)
1168  enddo
1169  enddo
1170  enddo
1171 
1172  CALL mpi_allreduce(intzyxtmp,intzyx_zvx,1,mpi_real8, &
1173  mpi_sum,mpi_comm_world,ierr)
1174 
1175  end function intzyx_zvx
1176 
1177  !----(入力データ zv)---
1178  function z_inty_zv(zv) ! Y積分
1179  !
1180  ! 2 次元(ZV)格子点データのY方向積分.
1181  !
1182  real(8), dimension(0:km,jc(ip)), intent(in) :: zv
1183  !(in) 2 次元ZV(子午面)格子点データ
1184 
1185  real(8), dimension(0:km) :: z_IntY_zv
1186  !(out) Y積分された 1 次元Z格子点データ
1187 
1188  real(8), dimension(0:km) :: z_IntYTmp
1189  integer :: j, ierr
1190 
1191  z_intytmp = 0.0d0
1192  do j=1,jc(ip)
1193  z_intytmp(:) = z_intytmp(:) + zv(:,j) * v_y_weight(j)
1194  enddo
1195 
1196  CALL mpi_allreduce(z_intytmp,z_inty_zv,km+1,mpi_real8, &
1197  mpi_sum,mpi_comm_world,ierr)
1198 
1199  end function z_inty_zv
1200 
1201  function v_intz_zv(zv) ! Z積分
1202  !
1203  ! 2 次元(ZV)格子点データのZ方向域積分.
1204  !
1205  real(8), dimension(0:km,jc(ip)), intent(in) :: zv
1206  !(in) 2 次元ZV格子点データ
1207 
1208  real(8), dimension(jc(ip)) :: v_IntZ_zv
1209  !(out) Z積分された 1 次元Y格子点データ
1210 
1211  integer :: k
1212 
1213  v_intz_zv = 0.0d0
1214  do k=0,km
1215  v_intz_zv(:) = v_intz_zv(:) &
1216  + zv(k,:) * z_z_weight(k)
1217  enddo
1218  end function v_intz_zv
1219 
1220  function intzy_zv(zv)
1221  !
1222  ! 2 次元(ZV)格子点データのZV積分
1223  !
1224  real(8), dimension(0:km,jc(ip)), intent(in) :: zv
1225  !(in) 2 次元ZV(子午面)格子点データ
1226 
1227  real(8) :: IntZY_zv
1228  !(out) 積分値
1229 
1230  real(8) :: IntZYTmp
1231  integer :: j, k, ierr
1232 
1233  intzytmp = 0
1234  do j=1,jc(ip)
1235  do k=0,km
1236  intzytmp = intzytmp &
1237  + zv(k,j) * v_y_weight(j) * z_z_weight(k)
1238  enddo
1239  enddo
1240 
1241  CALL mpi_allreduce(intzytmp,intzy_zv,1,mpi_real8, &
1242  mpi_sum,mpi_comm_world,ierr)
1243 
1244  end function intzy_zv
1245 
1246  !----(入力データ zx)---
1247  function z_intx_zx(zx)
1248  !
1249  ! 2 次元(ZX)格子点データのX方向積分.
1250  !
1251  real(8), dimension(0:km,0:im-1), intent(in) :: zx
1252  !(in) 2 次元ZV格子点データ
1253 
1254  real(8), dimension(0:km) :: z_IntX_zx
1255  !(out) X積分された 1 次元Z格子点データ
1256 
1257  integer :: i
1258 
1259  z_intx_zx = 0.0d0
1260  do i=0,im-1
1261  z_intx_zx(:) = z_intx_zx(:) + zx(:,i) * x_x_weight(i)
1262  enddo
1263 
1264  end function z_intx_zx
1265 
1266  function x_intz_zx(zx)
1267  !
1268  ! 2 次元(ZX)格子点データのZ方向積分.
1269  !
1270  real(8), dimension(0:km,0:im-1), intent(in) :: zx
1271  !(in) 2 次元ZV格子点データ
1272 
1273  real(8), dimension(0:im-1) :: x_IntZ_zx
1274  !(out) Z積分された 1 次元X格子点データ
1275 
1276  integer :: k
1277 
1278  x_intz_zx = 0.0d0
1279  do k=0,km
1280  x_intz_zx(:) = x_intz_zx(:) &
1281  + zx(k,:) * z_z_weight(k)
1282  enddo
1283 
1284  end function x_intz_zx
1285 
1286  function intzx_zx(zx) ! ZX積分
1287  !
1288  ! 2 次元(ZX)格子点データのZX積分
1289  !
1290  real(8), dimension(0:km,0:im-1), intent(in) :: zx
1291  !(in) 2 次元ZV格子点データ
1292 
1293  real(8) :: IntZX_zx
1294  !(out) 積分値
1295 
1296  integer :: i, k
1297 
1298  intzx_zx = 0
1299  do i=0,im-1
1300  do k=0,km
1301  intzx_zx = intzx_zx &
1302  + zx(k,i) * x_x_weight(i) * z_z_weight(k)
1303  enddo
1304  enddo
1305  end function intzx_zx
1306 
1307  !----(入力データ z)---
1308  function intz_z(z) ! Z積分
1309  !
1310  ! 1 次元(Z)格子点データのZ方向積分.
1311  !
1312  real(8), dimension(0:km), intent(in) :: z
1313  !(in) 1 次元Z格子点データ
1314 
1315  real(8) :: IntZ_z
1316  !(out) 積分値
1317 
1318  integer :: k
1319 
1320  intz_z = 0.0d0
1321  do k=0,km
1322  intz_z = intz_z + z(k) * z_z_weight(k)
1323  enddo
1324  end function intz_z
1325 
1326  !--------------- 平均計算 -----------------
1327  !----(入力データ zvx)---
1328  function zv_avrx_zvx(zvx) ! X積分
1329  !
1330  ! 3 次元格子点データのX方向平均.
1331  !
1332  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1333  !(in) 3 次元ZVX格子点データ
1334 
1335  real(8), dimension(0:km,jc(ip)) :: zv_AvrX_zvx
1336  !(out) X方向平均された 2 次元子午面格子点データ
1337 
1338  zv_avrx_zvx = zv_intx_zvx(zvx)/sum(x_x_weight)
1339 
1340  end function zv_avrx_zvx
1341 
1342  function zx_avry_zvx(zvx) ! Y平均
1343  !
1344  ! 3 次元格子点データのY方向域平均.
1345  !
1346  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1347  !(in) 3 次元格子点データ
1348 
1349  real(8), dimension(0:km,0:im-1) :: zx_AvrY_zvx
1350  !(out) Y平均された 2 次元ZX格子点データ
1351 
1352  zx_avry_zvx = zx_inty_zvx(zvx)/yl
1353 
1354  end function zx_avry_zvx
1355 
1356  function vx_avrz_zvx(zvx)
1357  !
1358  ! 3 次元格子点データのZ方向平均.
1359  !
1360  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1361  !(in) 3 次元格子点データ
1362 
1363  real(8), dimension(jc(ip),0:im-1) :: vx_AvrZ_zvx
1364  !(out) Z平均された 2 次元VX(水平)格子点データ
1365 
1366  vx_avrz_zvx = vx_intz_zvx(zvx)/sum(z_z_weight)
1367 
1368  end function vx_avrz_zvx
1369 
1370  function x_avrzy_zvx(zvx) ! ZV積分
1371  !
1372  ! 3 次元格子点データのZV平均
1373  !
1374  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1375  !(in) 3 次元ZVX格子点データ
1376 
1377  real(8), dimension(0:im-1) :: x_AvrZY_zvx
1378  !(out) ZV平均された 1 次元X格子点データ
1379 
1380  x_avrzy_zvx = x_intzy_zvx(zvx)/(yl*zl)
1381 
1382  end function x_avrzy_zvx
1383 
1384  function v_avrzx_zvx(zvx) ! ZX積分
1385  !
1386  ! 3 次元格子点データのZX平均.
1387  !
1388  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1389  !(in) 3 次元格子点データ
1390 
1391  real(8), dimension(jc(ip)) :: v_AvrZX_zvx
1392  !(out) ZX平均された 1 次元Y格子点データ
1393 
1394  v_avrzx_zvx = v_intzx_zvx(zvx) &
1395  /(sum(x_x_weight)*sum(z_z_weight))
1396 
1397  end function v_avrzx_zvx
1398 
1399  function z_avryx_zvx(zvx) ! VX(水平)積分
1400  !
1401  ! 3 次元格子点データのVX(水平)積分
1402  !
1403  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1404  !(in) 3 次元格子点データ
1405 
1406  real(8), dimension(0:km) :: z_AvrYX_zvx
1407  !(out) VX(水平)平均された 1 次元Z格子点データ
1408 
1409  z_avryx_zvx = z_intyx_zvx(zvx)/(xl*yl)
1410 
1411  end function z_avryx_zvx
1412 
1413  function avrzyx_zvx(zvx) ! ZVX(全領域)積分
1414  !
1415  ! 3 次元格子点データのZVX(全領域)積分
1416  !
1417  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
1418  !(in) 3 次元ZVX格子点データ
1419 
1420  real(8) :: AvrZYX_zvx
1421  !(out) 全領域平均値
1422 
1423  avrzyx_zvx = intzyx_zvx(zvx)/(xl*yl*zl)
1424 
1425  end function avrzyx_zvx
1426 
1427  !----(入力データ zv)---
1428  function z_avry_zv(zv)
1429  !
1430  ! 2 次元(ZV)格子点データのY方向平均.
1431  !
1432  real(8), dimension(0:km,jc(ip)), intent(in) :: zv
1433  !(in) 2 次元ZV格子点データ
1434 
1435  real(8), dimension(0:km) :: z_AvrY_zv
1436  !(out) Y平均された 1 次元Z格子点データ
1437 
1438  z_avry_zv = z_inty_zv(zv)/yl
1439 
1440  end function z_avry_zv
1441 
1442  function v_avrz_zv(zv)
1443  !
1444  ! 2 次元(ZV)格子点データのZ方向平均.
1445  !
1446  real(8), dimension(0:km,jc(ip)), intent(in) :: zv
1447  !(in) 2 次元ZV格子点データ
1448 
1449  real(8), dimension(jc(ip)) :: v_AvrZ_zv
1450  !(out) Z平均された 1 次元Y格子点データ
1451 
1452  v_avrz_zv = v_intz_zv(zv)/sum(z_z_weight)
1453 
1454  end function v_avrz_zv
1455 
1456  function avrzy_zv(zv) ! ZV平均
1457  !
1458  ! 2 次元(ZV)格子点データのZV平均
1459  !
1460  real(8), dimension(0:km,jc(ip)), intent(in) :: zv
1461  !(in) 2 次元ZV(子午面)格子点データ
1462 
1463  real(8) :: AvrZY_zv
1464  !(out) 平均値
1465 
1466  avrzy_zv = intzy_zv(zv)/(yl*zl)
1467 
1468  end function avrzy_zv
1469 
1470  !----(入力データ zx)---
1471  function z_avrx_zx(zx) ! X積分
1472  !
1473  ! 2 次元(ZX)格子点データのX方向平均.
1474  !
1475  real(8), dimension(0:km,0:im-1), intent(in) :: zx
1476  !(in) 2 次元ZV格子点データ
1477 
1478  real(8), dimension(0:km) :: z_AvrX_zx
1479  !(out) X平均された 1 次元Z格子点データ
1480 
1481  z_avrx_zx = z_intx_zx(zx)/sum(x_x_weight)
1482 
1483  end function z_avrx_zx
1484 
1485  function x_avrz_zx(zx) ! Z平均
1486  !
1487  ! 2 次元(ZX)格子点データのZ方向平均.
1488  !
1489  real(8), dimension(0:km,0:im-1), intent(in) :: zx
1490  !(in) 2 次元ZV格子点データ
1491 
1492  real(8), dimension(0:im-1) :: x_AvrZ_zx
1493  !(out) Z平均された 1 次元X格子点データ
1494 
1495  x_avrz_zx = x_intz_zx(zx)/sum(z_z_weight)
1496 
1497  end function x_avrz_zx
1498 
1499  function avrzx_zx(zx) ! ZX積分
1500  !
1501  ! 2 次元(ZX)格子点データのZX平均
1502  !
1503  real(8), dimension(0:km,0:im-1), intent(in) :: zx
1504  ! (in) 2 次元格子点データ
1505  real(8) :: AvrZX_zx
1506  ! 平均値
1507 
1508  avrzx_zx = intzx_zx(zx)/(sum(x_x_weight)*sum(z_z_weight))
1509 
1510  end function avrzx_zx
1511 
1512  !----(入力データ z)---
1513  function avrz_z(z)
1514  !
1515  ! 1 次元(Z)格子点データのZ方向平均.
1516  !
1517  real(8), dimension(0:km), intent(in) :: z
1518  !(in) 1 次元Z格子点データ
1519  real(8) :: AvrZ_z
1520  !(out) 平均値
1521 
1522  avrz_z = intz_z(z)/sum(z_z_weight)
1523 
1524  end function avrz_z
1525 
1526  !--------------- ポロイダル/トロイダルモデル用微分 -----------------
1527 
1528 !!$ function wt_KxRGrad_wt(wt)
1529 !!$ !
1530 !!$ ! 入力スペクトルデータにX微分 k×r・▽ = ∂/∂λを作用する.
1531 !!$ !
1532 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
1533 !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
1534 !!$
1535 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_KxRGrad_wt
1536 !!$ !(out) X微分を作用された 2 次元スペクトルデータ
1537 !!$
1538 !!$ wt_KxRGrad_wt = wa_Dlon_wa(wt)
1539 !!$
1540 !!$ end function wt_KxRGrad_wt
1541 !!$
1542 !!$ function zvx_KGrad_wt(wt) ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
1543 !!$ !
1544 !!$ ! 入力スペクトルデータに対応する格子データに軸方向微分
1545 !!$ !
1546 !!$ ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
1547 !!$ !
1548 !!$ ! を作用させた格子データが返される.
1549 !!$ ! ここでベクトル k は球の中心から北極向きの単位ベクトルである.
1550 !!$ !
1551 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
1552 !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
1553 !!$
1554 !!$ real(8), dimension(0:km,0:jm-1,0:im-1) :: yzx_KGrad_wt
1555 !!$ !(out) 軸方向微分を作用された 2 次元スペクトルデータ
1556 !!$
1557 !!$ xzv_KGrad_wt = cos(xyz_Y)*xyz_GradY_wt(wt) &
1558 !!$ + sin(xyz_Y)*xyz_wt(wt_Drad_wt(wt))
1559 !!$
1560 !!$ end function xyz_KGrad_wt
1561 !!$
1562 !!$ function wt_L2_wt(wt)
1563 !!$ !
1564 !!$ ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する.
1565 !!$ !
1566 !!$ ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる.
1567 !!$ ! 入力スペクトルデ ータに対応する格子点データに演算子
1568 !!$ !
1569 !!$ ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ)
1570 !!$ !
1571 !!$ ! を作用させたデータのスペクトル変換が返される.
1572 !!$ !
1573 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
1574 !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
1575 !!$
1576 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_L2_wt
1577 !!$ !(out) L^2 演算子を作用された 2 次元スペクトルデータ
1578 !!$
1579 !!$ wt_L2_wt = -wa_Lapla_wa(wt)
1580 !!$
1581 !!$ end function wt_L2_wt
1582 !!$
1583 !!$ function wt_L2Inv_wt(wt)
1584 !!$ !
1585 !!$ ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を
1586 !!$ ! 作用する.
1587 !!$ !
1588 !!$ ! スペクトルデータに L^2 演算子を作用させる関数 wt_L2_wt の逆計算を
1589 !!$ ! 行う関数である.
1590 !!$ !
1591 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
1592 !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
1593 !!$
1594 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_L2Inv_wt
1595 !!$ !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ
1596 !!$
1597 !!$ wt_L2Inv_wt = -wa_LaplaInv_wa(wt)
1598 !!$
1599 !!$ end function wt_L2Inv_wt
1600 !!$
1601 !!$ function wt_QOperator_wt(wt)
1602 !!$ !
1603 !!$ ! 入力スペクトルデータに対応する格子点データに演算子
1604 !!$ !
1605 !!$ ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2))
1606 !!$ !
1607 !!$ ! を作用させたデータのスペクトル変換が返される.
1608 !!$ !
1609 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
1610 !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
1611 !!$
1612 !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_QOperator_wt
1613 !!$ !(out) Q 演算子を作用された 2 次元スペクトルデータ
1614 !!$
1615 !!$ wt_QOperator_wt = &
1616 !!$ wt_xyz(xyz_KGrad_wt(wt) - xyz_KGrad_wt(wt_L2_wt(wt))/2) &
1617 !!$ - wt_L2_wt(wt_xyz(xyz_KGrad_wt(wt)))/2
1618 !!$
1619 !!$ end function wt_QOperator_wt
1620 
1621  function tef_zrot_zvx_zvx(zvx_VX,zvx_VY) ! z・(▽×v)
1622  !
1623  ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する.
1624  !
1625  ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す.
1626  !
1627  ! z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y
1628  !
1629  ! のスペクトル データが返される.
1630  !
1631  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX
1632  !(in) ベクトルのX成分
1633 
1634  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY
1635  !(in) ベクトルのY成分
1636 
1637  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_ZRot_zvx_zvx
1638  !(out) ベクトルの渦度とZベクトルの内積
1639 
1640  tef_zrot_zvx_zvx = tef_dx_tef(tef_zvx(zvx_vy)) &
1641  - tef_dy_tef(tef_zvx(zvx_vx))
1642 
1643  end function tef_zrot_zvx_zvx
1644 
1645  function tef_zrotrot_zvx_zvx_zvx(zvx_VX,zvx_VY,zvx_VZ)
1646  !
1647  ! ベクトル v に対して z・(▽×▽×v) を計算する.
1648  !
1649  ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分,
1650  ! Y成分, Z成分を表す.
1651  !
1652  ! z・(▽×▽×v) = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) )
1653  ! - ▽_H^2 v[z]
1654  !
1655  ! のスペクトルデータが返される.
1656  !
1657  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX
1658  !(in) ベクトルのX成分
1659 
1660  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY
1661  !(in) ベクトルのY成分
1662 
1663  real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VZ
1664  !(in) ベクトルのZ成分
1665 
1666  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_ZRotRot_zvx_zvx_zvx
1667  !(out) ベクトル v の z・(▽×▽×v)
1668 
1669  tef_zrotrot_zvx_zvx_zvx = &
1670  tef_dz_tef( tef_dx_tef(tef_zvx(zvx_vx)) &
1671  + tef_dy_tef(tef_zvx(zvx_vy)) ) &
1672  - tef_laplah_tef(tef_zvx(zvx_vz))
1673 
1674  end function tef_zrotrot_zvx_zvx_zvx
1675 
1676  subroutine tef_potential2vector(&
1677  zvx_VX,zvx_VY,zvx_VZ,tef_Torvel,tef_Polvel)
1678  !
1679  ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
1680  !
1681  ! v = ▽x(Ψz) + ▽x▽x(Φz)
1682  !
1683  ! の各成分を計算する
1684  !
1685  real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_VX
1686  !(out) ベクトル場のX成分
1687 
1688  real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_VY
1689  !(out) ベクトル場のY成分
1690 
1691  real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_VZ
1692  !(out) ベクトル場のZ成分
1693 
1694  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Torvel
1695  !(in) トロイダルポテンシャル
1696 
1697  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Polvel
1698  !(in) ポロイダルポテンシャル
1699 
1700  zvx_vx = zvx_tef( tef_dy_tef(tef_torvel) &
1701  + tef_dx_tef(tef_dz_tef(tef_polvel)) )
1702  zvx_vy = zvx_tef( - tef_dx_tef(tef_torvel) &
1703  + tef_dy_tef(tef_dz_tef(tef_polvel)) )
1704  zvx_vz = -zvx_tef(tef_laplah_tef(tef_polvel))
1705 
1706  end subroutine tef_potential2vector
1707 
1708  subroutine tef_potential2rotation(&
1709  zvx_RotVX,zvx_RotVY,zvx_RotVZ,tef_Torvel,tef_Polvel)
1710  !
1711  ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
1712  !
1713  ! v = ▽x(Ψz) + ▽x▽x(Φz)
1714  !
1715  ! に対して, その回転
1716  !
1717  ! ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z)
1718  !
1719  ! を計算する.
1720 
1721  ! ベクトル場の回転
1722  real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVX
1723  !(out) 回転のX成分
1724 
1725  real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVY
1726  !(out) 回転のY成分
1727 
1728  real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVZ
1729  !(out) 回転のZ成分
1730 
1731  ! 入力ベクトル場を表すポテンシャル
1732  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Torvel
1733  !(in) トロイダルポテンシャル
1734 
1735  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Polvel
1736  !(in) ポロイダルポテンシャル
1737 
1738  call tef_potential2vector( &
1739  zvx_rotvx,zvx_rotvy,zvx_rotvz, &
1740  -tef_lapla_tef(tef_polvel), tef_torvel)
1741 
1742  end subroutine tef_potential2rotation
1743 
1744  !--------------- 補間計算 -----------------
1745  function interpolate_tef(tef_data,x,y,z)
1746  !
1747  ! (x,y,z) における関数値を
1748  ! そのスペクトル係数 tef_data から補間計算する
1749  !
1750  real(8), intent(IN) :: tef_data(0:nm,-mm:mm,2*lc(ip))! スペクトルデータ
1751  real(8), intent(IN) :: x ! 補間する位置(X)
1752  real(8), intent(IN) :: y ! 補間する位置(Y)
1753  real(8), intent(IN) :: z ! 補間する位置(Z)
1754  real(8) :: Interpolate_tef ! 補間した値
1755 
1756  interpolate_tef = &
1757  interpolate_ef(ef_e2(a_interpolate_at(e2a_aef(tef_data),z)),x,y)
1758 
1759  end function interpolate_tef
1760 
1761  !--------------- ポロイダル/トロイダルモデル用スペクトル解析 ----------------
1762 
1763  function zef_toroidalenergyspectrum_tef(tef_TORPOT)
1764  !
1765  ! トロイダルポテンシャルから, トロイダルエネルギーの
1766  ! 水平 X 波数 l, Y 波数 m の各成分を計算する
1767  !
1768  ! * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分
1769  ! ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー
1770  ! スペクトルは (1/2)(l**2+m**2)|ψ(l,m,z)|^2 と計算される.
1771  !
1772  ! * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での
1773  ! 水平平均エネルギーに等しい.
1774  !
1775  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_TORPOT
1776  !(in) トロイダルポテンシャル
1777 
1778  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_ToroidalEnergySpectrum_tef
1779  !(out) エネルギースペクトルトロイダル成分
1780 
1781  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Work
1782  integer :: l, m
1783 
1784  zef_toroidalenergyspectrum_tef = 0.0d0
1785  zef_work = zef_tef(tef_torpot)
1786  do l=-lm,lm
1787  do m=-mm,mm
1788  if ( lf_l(l) > 0 ) then
1789  zef_toroidalenergyspectrum_tef(:,m,lf_l(l)) &
1790  = 0.5 * ((2*pi*l/xl)**2+(2*pi*m/yl)**2) &
1791  * ( zef_work(:,m,lf_l(l))**2 + zef_work(:,-m,lf_l(-l))**2 )
1792  endif
1793  enddo
1794  enddo
1795 
1796  end function zef_toroidalenergyspectrum_tef
1797 
1798 !!$ function nz_ToroidalEnergySpectrum_wt(wt_TORPOT)
1799 !!$ !
1800 !!$ ! トロイダルポテンシャルから, トロイダルエネルギーの
1801 !!$ ! 球面調和函数全波数の各成分を計算する.
1802 !!$ !
1803 !!$ ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分
1804 !!$ ! ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは
1805 !!$ ! Σ[m=-n]^n(1/2)n(n+1)4πr^2ψ(n,m,r)^2 と計算される.
1806 !!$ !
1807 !!$ ! * 全てのエネルギースペクトル成分の和をZ積分したもの(r^2の重み無し)
1808 !!$ ! が球殻内での全エネルギーに等しい.
1809 !!$ !
1810 !!$ real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_TORPOT
1811 !!$ !(in) トロイダルポテンシャル
1812 !!$
1813 !!$ real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wt
1814 !!$ !(out) エネルギースペクトルトロイダル成分
1815 !!$
1816 !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA ! 作業領域
1817 !!$ integer :: n, m
1818 !!$
1819 !!$ wz_DATA = wz_wt(wt_TORPOT)
1820 !!$ do n=0,nm
1821 !!$ nz_ToroidalEnergySpectrum_wt(n,:) &
1822 !!$ = 0.5 * n*(n+1)* (4*pi) * z_Z**2 &
1823 !!$ * sum(wz_Data(l_nm(n,(/(m,m=-n,n)/)),:)**2,1)
1824 !!$ enddo
1825 !!$
1826 !!$ end function nz_ToroidalEnergySpectrum_wt
1827 
1828  function zef_poloidalenergyspectrum_tef(tef_POLPOT)
1829  !
1830  ! ポロイダルポテンシャルから, ポロイダルエネルギーの
1831  ! 水平 X 波数 l, Y 波数 m の各成分を計算する.
1832  !
1833  ! * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分
1834  ! φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー
1835  ! スペクトルは
1836  !
1837  ! (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2}
1838  !
1839  ! と計算される.
1840  !
1841  ! * 全てのエネルギースペクトル成分の和をZ積分したもの
1842  ! が水平平均エネルギーに等しい.
1843  !
1844  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_POLPOT
1845  !(in) ポロイダルポテンシャル
1846 
1847  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_PoloidalEnergySpectrum_tef
1848  !(out) エネルギースペクトルポロイダル成分
1849 
1850 
1851  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Data ! 作業領域
1852  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_DData ! 作業領域
1853  integer :: l, m
1854 
1855  zef_poloidalenergyspectrum_tef = 0.0d0
1856  zef_data = zef_tef(tef_polpot)
1857  zef_ddata = zef_tef(tef_dz_tef(tef_polpot))
1858 
1859  do l=-lm,lm
1860  do m=-mm,mm
1861  if ( lf_l(l) > 0 ) then
1862  zef_poloidalenergyspectrum_tef(:,m,lf_l(l)) = &
1863  + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2) &
1864  *( zef_ddata(:,m,lf_l(l))**2 + zef_ddata(:,-m,lf_l(-l))**2 &
1865  + ((2*pi*l/xl)**2+(2*pi*m/yl)**2) &
1866  *( zef_data(:,m,lf_l(l))**2 + zef_data(:,-m,lf_l(-l))**2))
1867  endif
1868  enddo
1869  enddo
1870 
1871  end function zef_poloidalenergyspectrum_tef
1872 
1873 !!$ function nz_PoloidalEnergySpectrum_wt(wt_POLPOT)
1874 !!$ !
1875 !!$ ! ポロイダルポテンシャルから, ポロイダルエネルギーの
1876 !!$ ! 球面調和函数全波数の各成分を計算する
1877 !!$ !
1878 !!$ ! * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分
1879 !!$ ! φ(n,m,r)から全波数 n 成分のポロイダルエネルギースペクトルは
1880 !!$ !
1881 !!$ ! Σ[m=-n]^n ((1/2)n(n+1)4πr^2{[d(rφ(n,m,r))/dr]^2
1882 !!$ ! + n(n+1)φ(n,m,r)^2}
1883 !!$ !
1884 !!$ ! と計算される.
1885 !!$ !
1886 !!$ ! * 全ての全波数に対してのエネルギースペクトル成分の和をZ積分したもの
1887 !!$ ! (r^2の重み無し)が球殻内での全エネルギーに等しい.
1888 !!$ !
1889 !!$ real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_POLPOT
1890 !!$ !(in) ポロイダルポテンシャル
1891 !!$
1892 !!$ real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_wt
1893 !!$ !(out) エネルギースペクトルポロイダル成分
1894 !!$
1895 !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1 ! 作業領域
1896 !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2 ! 作業領域
1897 !!$ integer :: n, m
1898 !!$
1899 !!$ wz_Data1 = wz_wt(wt_POLPOT)
1900 !!$ wz_Data2 = wz_Z*wz_wt(wt_DZ_wt(wt_POLPOT)) & ! d(rφ)/dr
1901 !!$ + wz_wt(wt_POLPOT) ! = rdφ/dr+φ
1902 !!$
1903 !!$ do n=0,nm
1904 !!$ nz_PoloidalEnergySpectrum_wt(n,:) = &
1905 !!$ + 0.5* n*(n+1)* (4*pi) &
1906 !!$ *( sum(wz_Data2(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) &
1907 !!$ + n*(n+1)*sum(wz_Data1(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) )
1908 !!$ enddo
1909 !!$
1910 !!$ end function nz_PoloidalEnergySpectrum_wt
1911 
1912  function zef_toroidalenstrophyspectrum_tef(tef_TORPOT)
1913  !
1914  ! トロイダルポテンシャルから, トロイダルエンストロフィーの
1915  ! 水平 X 波数 l, Y 波数 m の各成分を計算する
1916  !
1917  ! * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分
1918  ! ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー
1919  ! スペクトルは
1920  !
1921  ! (1/2)(l**2+m**2){[dψ(n,m,z)/dz]^2 + (l**2+m**2)ψ(n,m,z)^2}
1922  !
1923  ! * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での
1924  ! 水平平均エネルギーに等しい.
1925  !
1926  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_TORPOT
1927  !(in) トロイダルポテンシャル
1928 
1929  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_ToroidalEnstrophySpectrum_tef
1930  !(out) エンストロフィースペクトルトロイダル成分
1931 
1932  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Data ! 作業領域
1933  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_DData ! 作業領域
1934  integer :: l, m
1935 
1936 
1937  zef_toroidalenstrophyspectrum_tef = 0.0d0
1938  zef_data = zef_tef(tef_torpot)
1939  zef_ddata = zef_tef(tef_dz_tef(tef_torpot))
1940 
1941  do l=-lm,lm
1942  do m=-mm,mm
1943  if ( lf_l(l) > 0 ) then
1944  zef_toroidalenstrophyspectrum_tef(:,m,lf_l(l)) = &
1945  + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2) &
1946  *( zef_ddata(:,m,lf_l(l))**2 + zef_ddata(:,-m,lf_l(-l))**2 &
1947  + ((2*pi*l/xl)**2+(2*pi*m/yl)**2) &
1948  *( zef_data(:,m,lf_l(l))**2 + zef_data(:,-m,lf_l(-l))**2))
1949  endif
1950  enddo
1951  enddo
1952 
1953  end function zef_toroidalenstrophyspectrum_tef
1954 
1955  function zef_poloidalenstrophyspectrum_tef(tef_POLPOT)
1956  !
1957  ! ポロイダルポテンシャルから, ポロイダルエンストロフィーの
1958  ! 水平 X 波数 l, Y 波数 m の各成分を計算する.
1959  !
1960  ! * X 波数 l, Y 波数 m のポロイダルポテンシャルエンストロフィーの
1961  ! スペクトル成分 φ(l,m,z)から X 波数 l, Y 波数 m 成分の
1962  ! ポロイダルエンストロフィースペクトルは
1963  !
1964  ! (1/2)(l**2+m**2){[d^φ(n,m,z)/dz^2] - (l**2+m**2)φ(n,m,z)]^2}
1965  !
1966  ! と計算される.
1967  !
1968  ! * 全てのエネルギースペクトル成分の和をZ積分したもの
1969  ! が水平平均エネルギーに等しい.
1970  !
1971  real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_POLPOT
1972  !(in) ポロイダルポテンシャル
1973 
1974  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_PoloidalEnstrophySpectrum_tef
1975  !(out) エネルギースペクトルポロイダル成分
1976 
1977  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Work
1978  integer :: l, m
1979 
1980  zef_poloidalenstrophyspectrum_tef = 0.0d0
1981  zef_work = zef_tef(tef_lapla_tef(tef_polpot))
1982  do l=-lm,lm
1983  do m=-mm,mm
1984  if ( lf_l(l) > 0 ) then
1985  zef_poloidalenstrophyspectrum_tef(:,m,lf_l(l)) &
1986  = 0.5 * ((2*pi*l/xl)**2+(2*pi*m/yl)**2) &
1987  * ( zef_work(:,m,lf_l(l))**2 + zef_work(:,-m,lf_l(-l))**2 )
1988  endif
1989  enddo
1990  enddo
1991 
1992  end function zef_poloidalenstrophyspectrum_tef
1993 
1994  !--------------- 境界値問題 -----------------
1995 
1996  subroutine tef_boundariestau(tef,values,cond)
1997  !
1998  ! スペクトルデータにディリクレ・ノイマン境界条件を適用する
1999  ! Chebyshev 空間での境界条件適用(タウ法)
2000  !
2001  ! チェビシェフ空間において境界条件を満たすべく高次の係数を
2002  ! 定める方法をとっている(タウ法).
2003  !
2004  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef
2005  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2006 
2007  real(8), dimension(2,-mm:mm,2*lc(ip)), intent(in), optional :: values
2008  !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える.
2009  ! 省略時は値/勾配 0 となる.
2010 
2011  character(len=2), intent(in), optional :: cond
2012  !(in) 境界条件. 省略時は 'DD'
2013  ! DD : 両端ディリクレ条件
2014  ! DN : 上端ディリクレ, 下端ノイマン条件
2015  ! ND : 上端ノイマン, 下端ディリクレ条件
2016  ! NN : 両端ノイマン条件
2017 
2018  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t
2019 
2020  real(8), dimension(2*lc(ip)*(2*mm+1),2) :: e22_values
2021  character(len=2) :: Bcond
2022 
2023  e2t = e2a_aef(tef)
2024 
2025  if (present(values)) then
2026  e22_values = e2a_aef(values)
2027  endif
2028 
2029  if (.not. present(cond)) then
2030  bcond='DD'
2031  else
2032  bcond = cond
2033  endif
2034 
2035  select case(bcond)
2036  case ('NN')
2037  if (present(values)) then
2038  call at_boundariestau_nn(e2t,e22_values)
2039  else
2040  call at_boundariestau_nn(e2t)
2041  endif
2042  case ('DN')
2043  if (present(values)) then
2044  call at_boundariestau_dn(e2t,e22_values)
2045  else
2046  call at_boundariestau_dn(e2t)
2047  endif
2048  case ('ND')
2049  if (present(values)) then
2050  call at_boundariestau_nd(e2t,e22_values)
2051  else
2052  call at_boundariestau_nd(e2t)
2053  endif
2054  case ('DD')
2055  if (present(values)) then
2056  call at_boundariestau_dd(e2t,e22_values)
2057  else
2058  call at_boundariestau_dd(e2t)
2059  endif
2060  case default
2061  call messagenotify('E','tef_BoundariesTau','B.C. not supported')
2062  end select
2063 
2064  tef = aef_e2a(e2t)
2065 
2066  end subroutine tef_boundariestau
2067 
2068  subroutine tef_boundariesgrid(tef,values,cond)
2069  !
2070  ! スペクトルデータにディリクレ・ノイマン境界条件を適用する
2071  ! 実空間での境界条件適用
2072  !
2073  ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
2074  ! 条件を課している(選点法). このルーチンを用いるためには
2075  ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を
2076  ! 等しくしておく必要がある.
2077  !
2078  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef
2079  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2080 
2081  real(8), dimension(2,-mm:mm,2*lc(ip)), intent(in), optional :: values
2082  !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える.
2083  ! 省略時は値/勾配 0 となる.
2084 
2085  character(len=2), intent(in), optional :: cond
2086  !(in) 境界条件. 省略時は 'DD'
2087  ! DD : 両端ディリクレ条件
2088  ! DN : 上端ディリクレ, 下端ノイマン条件
2089  ! ND : 上端ノイマン, 下端ディリクレ条件
2090  ! NN : 両端ノイマン条件
2091 
2092  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t
2093 
2094  real(8), dimension(2*lc(ip)*(2*mm+1),2) :: e22_values
2095  character(len=2) :: Bcond
2096 
2097  e2t = e2a_aef(tef)
2098 
2099  if (present(values)) then
2100  e22_values = e2a_aef(values)
2101  endif
2102 
2103  if (.not. present(cond)) then
2104  bcond = 'DD'
2105  else
2106  bcond = cond
2107  endif
2108 
2109  select case(bcond)
2110  case ('NN')
2111  if (present(values)) then
2112  call at_boundariesgrid_nn(e2t,e22_values)
2113  else
2114  call at_boundariesgrid_nn(e2t)
2115  endif
2116  case ('DN')
2117  if (present(values)) then
2118  call at_boundariesgrid_dn(e2t,e22_values)
2119  else
2120  call at_boundariesgrid_dn(e2t)
2121  endif
2122  case ('ND')
2123  if (present(values)) then
2124  call at_boundariesgrid_nd(e2t,e22_values)
2125  else
2126  call at_boundariesgrid_nd(e2t)
2127  endif
2128  case ('DD')
2129  if (present(values)) then
2130  call at_boundariesgrid_dd(e2t,e22_values)
2131  else
2132  call at_boundariesgrid_dd(e2t)
2133  endif
2134  case default
2135  call messagenotify('E','tef_BoundariesGrid','B.C. not supported')
2136  end select
2137 
2138  tef = aef_e2a(e2t)
2139 
2140  end subroutine tef_boundariesgrid
2141 
2142  subroutine tef_torboundariestau(tef_TOR,cond,new)
2143 
2144  ! トロイダル速度ポテンシャルに対して境界条件を適用する.
2145  ! Chebyshev 空間での境界条件適用
2146  !
2147  ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
2148  ! とっている(タウ法). トロイダルポテンシャルの境界条件は
2149  !
2150  ! ψ = 0 at boundaries (粘着条件)
2151  ! ∂ψ/∂z = 0 at boundaries (応力なし条件)
2152  !
2153  ! であるから tef_Boundaries で対応可能だが, 将来のため別途作成しておく.
2154  !
2155  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2156  !
2157  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR
2158  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2159 
2160  character(len=2), intent(in), optional :: cond
2161  !(in) 境界条件スイッチ. 省略時は 'RR'
2162  ! RR : 両端粘着条件
2163  ! RF : 上端粘着, 下端応力なし条件
2164  ! FR : 上端応力なし, 下端粘着条件
2165  ! FF : 両端応力なし条件
2166 
2167  logical, intent(IN), optional :: new
2168  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2169  ! default は false.
2170 
2171  real(8), dimension(:,:,:), allocatable :: alu
2172  integer, dimension(:,:), allocatable :: kp
2173 
2174  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2175  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2176 
2177  logical :: rigid1 = .true.
2178  logical :: rigid2 = .true.
2179  logical :: first = .true.
2180  logical :: new_matrix = .false.
2181  integer :: n
2182  save :: alu, kp, first
2183 
2184  select case (cond)
2185  case ('RR')
2186  rigid1 = .true. ; rigid2 = .true.
2187  case ('RF')
2188  rigid1 = .true. ; rigid2 = .false.
2189  case ('FR')
2190  rigid1 = .false. ; rigid2 = .true.
2191  case ('FF')
2192  rigid1 = .false. ; rigid2 = .false.
2193  case default
2194  call messagenotify('E','tef_TorBoundariesTau','B.C. not supported')
2195  end select
2196 
2197  if (.not. present(new)) then
2198  new_matrix=.false.
2199  else
2200  new_matrix=new
2201  endif
2202 
2203  if ( first .OR. new_matrix ) then
2204  first = .false.
2205 
2206  if ( allocated(alu) ) deallocate(alu)
2207  if ( allocated(kp) ) deallocate(kp)
2208  allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2209 
2210  do n=0,nm
2211  e2t_work= 0.0d0 ; e2t_work(:,n)= 1.0d0
2212  alu(:,:,n) = e2t_work
2213 
2214  ! 粘着条件
2215  ! 力学的条件粘着条件
2216  if ( rigid1 ) then
2217  e2z_work = az_at(e2t_work)
2218  else
2219  e2z_work = az_at(at_dz_at(e2t_work))
2220  endif
2221  alu(:,nm-1,n) = e2z_work(:,0)
2222 
2223  ! 力学的条件粘着条件
2224  if ( rigid2 ) then
2225  e2z_work = az_at(e2t_work)
2226  else
2227  e2z_work = az_at(at_dz_at(e2t_work))
2228  endif
2229  alu(:,nm,n) = e2z_work(:,km)
2230  enddo
2231 
2232  call ludecomp(alu,kp)
2233  call messagenotify('M','tef_TorBoundariesTau',&
2234  'Matrix to apply b.c. newly produced.')
2235  endif
2236 
2237  e2t_work = e2a_aef(tef_tor)
2238 
2239  e2t_work(:,nm-1) = 0.0d0
2240  e2t_work(:,nm) = 0.0d0
2241 
2242  tef_tor = aef_e2a(lusolve(alu,kp,e2t_work))
2243 
2244  end subroutine tef_torboundariestau
2245 
2246  subroutine tef_torboundariesgrid(tef_TOR,cond,new)
2247  !
2248  ! トロイダル速度ポテンシャルに対して境界条件を適用する.
2249  ! 鉛直実空間での境界条件適用.
2250  !
2251  ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
2252  ! 条件を課している(選点法). このルーチンを用いるためには
2253  ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を
2254  ! 等しくしておく必要がある.
2255  !
2256  ! トロイダル速度ポテンシャルの境界条件は
2257  !
2258  ! ψ = 0 at boundaries (粘着条件)
2259  ! ∂ψ/∂z = 0 at boundaries (応力なし条件)
2260  !
2261  ! であるので tef_Boundaries で対応可能だが, 将来のため別途作成しておく
2262  !
2263  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2264  !
2265  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR
2266  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2267 
2268  character(len=2), intent(in), optional :: cond
2269  !(in) 境界条件スイッチ. 省略時は 'RR'
2270  ! RR : 両端粘着条件
2271  ! RF : 上端粘着, 下端応力なし条件
2272  ! FR : 上端応力なし, 下端粘着条件
2273  ! FF : 両端応力なし条件
2274 
2275  logical, intent(IN), optional :: new
2276  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2277  ! default は false.
2278 
2279  real(8), dimension(:,:,:), allocatable :: alu
2280  integer, dimension(:,:), allocatable :: kp
2281 
2282  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2283  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2284 
2285  logical :: rigid1 = .true.
2286  logical :: rigid2 = .true.
2287  logical :: first = .true.
2288  logical :: new_matrix = .false.
2289  integer :: n
2290  save :: alu, kp, first
2291 
2292  select case (cond)
2293  case ('RR')
2294  rigid1 = .true. ; rigid2 = .true.
2295  case ('RF')
2296  rigid1 = .true. ; rigid2 = .false.
2297  case ('FR')
2298  rigid1 = .false. ; rigid2 = .true.
2299  case ('FF')
2300  rigid1 = .false. ; rigid2 = .false.
2301  case default
2302  call messagenotify('E','tef_TorBoundariesGrid','B.C. not supported')
2303  end select
2304 
2305  if (.not. present(new)) then
2306  new_matrix=.false.
2307  else
2308  new_matrix=new
2309  endif
2310 
2311  if ( first .OR. new_matrix ) then
2312  first = .false.
2313 
2314  if ( nm /= km ) then
2315  call messagenotify('E','tef_TorBoundariesGrid', &
2316  'Chebyshev truncation and number of grid points should be same.')
2317  endif
2318 
2319  if ( allocated(alu) ) deallocate(alu)
2320  if ( allocated(kp) ) deallocate(kp)
2321  allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2322 
2323  do n=0,nm
2324  e2t_work = 0.0d0 ; e2t_work(:,n)=1.0d0
2325  e2z_work = az_at(e2t_work)
2326 
2327  alu(:,:,n) = e2z_work ! 内部領域は値そのまま.
2328 
2329 
2330  ! 粘着条件
2331  ! 力学的条件粘着条件
2332  if ( rigid1 ) then
2333  e2z_work = az_at(e2t_work)
2334  else
2335  e2z_work=az_at(at_dz_at(e2t_work))
2336  endif
2337  alu(:,0,n) = e2z_work(:,0)
2338 
2339  ! 力学的条件粘着条件
2340  if ( rigid2 ) then
2341  e2z_work = az_at(e2t_work)
2342  else
2343  e2z_work=az_at(at_dz_at(e2t_work))
2344  endif
2345  alu(:,km,n) = e2z_work(:,km)
2346 
2347  enddo
2348  call ludecomp(alu,kp)
2349  call messagenotify('M','TorBoundariesGrid',&
2350  'Matrix to apply b.c. newly produced.')
2351  endif
2352 
2353  e2z_work = az_at(e2a_aef(tef_tor))
2354  e2z_work(:,0) = 0.0d0
2355  e2z_work(:,km) = 0.0d0
2356  tef_tor = aef_e2a(lusolve(alu,kp,e2z_work))
2357 
2358  end subroutine tef_torboundariesgrid
2359 
2360  function tef_laplapol2poltau_tef(tef,cond,new)
2361  !
2362  ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
2363  ! チェビシェフ空間でタウ法による境界条件を適用している.
2364  !
2365  ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
2366  !
2367  ! ▽^2Φ = f
2368  ! Φ = const. at boundaries.
2369  ! ∂Φ/∂z = 0 at boundaries (粘着条件)
2370  ! or ∂^2Φ/∂z^2 = 0 at boundaries (応力なし条件)
2371  !
2372  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2373  !
2374  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(in) :: tef
2375  !(in) 入力▽^2φ分布
2376 
2377  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaPol2PolTau_tef
2378  !(out) 出力ポロイダルポテンシャル分布
2379 
2380  character(len=2), intent(in), optional :: cond
2381  !(in) 境界条件スイッチ. 省略時は 'RR'
2382  ! RR : 両端粘着条件
2383  ! RF : 上端粘着, 下端応力なし条件
2384  ! FR : 上端応力なし, 下端粘着条件
2385  ! FF : 両端応力なし条件
2386 
2387  logical, intent(IN), optional :: new
2388  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2389  ! default は false.
2390 
2391  real(8), dimension(:,:,:), allocatable :: alu
2392  integer, dimension(:,:), allocatable :: kp
2393 
2394  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2395  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2396  logical :: rigid1 = .true.
2397  logical :: rigid2 = .true.
2398 
2399  logical :: first = .true.
2400  logical :: new_matrix = .false.
2401  integer :: n
2402  save :: alu, kp, first
2403 
2404  if ( present(cond) ) then
2405  select case (cond)
2406  case ('RR')
2407  rigid1 = .true. ; rigid2 = .true.
2408  case ('RF')
2409  rigid1 = .true. ; rigid2 = .false.
2410  case ('FR')
2411  rigid1 = .false. ; rigid2 = .true.
2412  case ('FF')
2413  rigid1 = .false. ; rigid2 = .false.
2414  case default
2415  call messagenotify('E','tef_LaplaPol2PolTau_tef','B.C. not supported')
2416  end select
2417  else
2418  rigid1 = .true. ; rigid2 = .true.
2419  endif
2420 
2421  if (.not. present(new)) then
2422  new_matrix=.false.
2423  else
2424  new_matrix=new
2425  endif
2426 
2427  if ( first .OR. new_matrix ) then
2428  first = .false.
2429 
2430  if ( allocated(alu) ) deallocate(alu)
2431  if ( allocated(kp) ) deallocate(kp)
2432  allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2433 
2434  alu = 0.0d0
2435  do n=0,nm-2
2436  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2437 
2438  ! 各水平波数に関して独立の式
2439  alu(:,:,n) = e2a_aef(tef_lapla_tef(aef_e2a(e2t_work)))
2440  enddo
2441 
2442  do n=0,nm
2443  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2444  e2z_work = az_at(e2t_work)
2445 
2446  ! 運動学的条件. 流線は境界で一定
2447  alu(:,nm,n) = e2z_work(:,0)
2448  alu(:,nm-1,n) = e2z_work(:,km)
2449 
2450  ! 力学的条件粘着条件
2451  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2452  if ( rigid1 ) then
2453  e2z_work=az_at(at_dz_at(e2t_work))
2454  else
2455  e2z_work=az_at(at_dz_at(at_dz_at(e2t_work)))
2456  endif
2457  alu(:,nm-2,n) = e2z_work(:,0)
2458 
2459  ! 力学的条件粘着条件
2460  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2461  if ( rigid2 ) then
2462  e2z_work=az_at(at_dz_at(e2t_work))
2463  else
2464  e2z_work=az_at(at_dz_at(at_dz_at(e2t_work)))
2465  endif
2466  alu(:,nm-3,n) = e2z_work(:,km)
2467  enddo
2468 
2469  call ludecomp(alu,kp)
2470 
2471  call messagenotify('M','tef_LaplaPol2PolTau_tef',&
2472  'Matrix to apply b.c. newly produced.')
2473  endif
2474 
2475  e2t_work = e2a_aef(tef)
2476  e2t_work(:,nm) = 0.0d0 ! 運動学的条件
2477  e2t_work(:,nm-1) = 0.0d0 ! 運動学的条件
2478  e2t_work(:,nm-2) = 0.0d0 ! 力学的条件
2479  e2t_work(:,nm-3) = 0.0d0 ! 力学的条件
2480 
2481  tef_laplapol2poltau_tef = aef_e2a(lusolve(alu,kp,e2t_work))
2482 
2483  end function tef_laplapol2poltau_tef
2484 
2485  function zef_laplapol2pol_zef(zef,cond,new)
2486  !
2487  ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
2488  !
2489  ! チェビシェフ格子点空間で境界条件を適用している.
2490  ! この関数を用いるためには wt_Initial にて設定する
2491  ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく
2492  ! しておく必要がある.
2493  !
2494  ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
2495  !
2496  ! ▽^2Φ = f
2497  ! Φ = const. at boundaries.
2498  ! ∂Φ/∂z = 0 at boundaries (粘着条件)
2499  ! or ∂^2Φ/∂z^2 = 0 at boundaries (応力なし条件)
2500  !
2501  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2502  !
2503  real(8), dimension(0:km,-mm:mm,2*lc(ip)),intent(in) :: zef
2504  !(in) 入力▽^2φ分布
2505 
2506  real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_LaplaPol2Pol_zef
2507  !(out) 出力ポロイダルポテンシャル分布
2508 
2509  character(len=2), intent(in), optional :: cond
2510  !(in) 境界条件スイッチ. 省略時は 'RR'
2511  ! RR : 両端粘着条件
2512  ! RF : 上端粘着, 下端応力なし条件
2513  ! FR : 上端応力なし, 下端粘着条件
2514  ! FF : 両端応力なし条件
2515 
2516  logical, intent(IN), optional :: new
2517  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2518  ! default は false.
2519 
2520  real(8), dimension(:,:,:), allocatable :: alu
2521  integer, dimension(:,:), allocatable :: kp
2522 
2523  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2524  logical :: rigid1 = .true.
2525  logical :: rigid2 = .true.
2526 
2527  logical :: first = .true.
2528  logical :: new_matrix = .false.
2529  integer :: k
2530  save :: alu, kp, first
2531 
2532  select case (cond)
2533  case ('RR')
2534  rigid1 = .true. ; rigid2 = .true.
2535  case ('RF')
2536  rigid1 = .true. ; rigid2 = .false.
2537  case ('FR')
2538  rigid1 = .false. ; rigid2 = .true.
2539  case ('FF')
2540  rigid1 = .false. ; rigid2 = .false.
2541  case default
2542  call messagenotify('E','zef_laplapol2pol_zef','B.C. not supported')
2543  end select
2544 
2545  if (.not. present(new)) then
2546  new_matrix=.false.
2547  else
2548  new_matrix=new
2549  endif
2550 
2551  if ( first .OR. new_matrix ) then
2552  first = .false.
2553 
2554  if ( nm /= km ) then
2555  call messagenotify('E','zef_LaplaPol2Pol_zef', &
2556  'Chebyshev truncation and number of grid points should be same.')
2557  endif
2558 
2559  if ( allocated(alu) ) deallocate(alu)
2560  if ( allocated(kp) ) deallocate(kp)
2561  allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:km),kp(2*lc(ip)*(2*mm+1),0:km))
2562 
2563  do k=0,km
2564  e2z_work = 0.0d0 ; e2z_work(:,k) = 1.0d0
2565 
2566  ! 各水平波数に関して独立の式
2567  alu(:,:,k) &
2568  = e2a_aef(zef_tef(tef_lapla_tef(tef_zef(aef_e2a(e2z_work)))))
2569  enddo
2570 
2571  do k=0,km
2572  e2z_work = 0.0d0 ; e2z_work(:,k) = 1.0d0
2573 
2574  ! 運動学的条件. 流線は境界で一定
2575  alu(:,0,k) = e2z_work(:,0)
2576  alu(:,km,k) = e2z_work(:,km)
2577 
2578  ! 力学的条件粘着条件
2579  e2z_work = 0.0d0 ; e2z_work(:,k) = 1.0d0
2580  if ( rigid1 ) then
2581  e2z_work=az_at(at_dz_at(at_az(e2z_work)))
2582  else
2583  e2z_work=az_at(at_dz_at(at_dz_at(at_az(e2z_work))))
2584  endif
2585  alu(:,1,k) = e2z_work(:,0)
2586 
2587  ! 力学的条件粘着条件
2588  e2z_work = 0.0d0 ; e2z_work(:,k) = 1.0d0
2589  if ( rigid2 ) then
2590  e2z_work=az_at(at_dz_at(at_az(e2z_work)))
2591  else
2592  e2z_work=az_at(at_dz_at(at_dz_at(at_az(e2z_work))))
2593  endif
2594  alu(:,km-1,k) = e2z_work(:,km)
2595  enddo
2596 
2597  call ludecomp(alu,kp)
2598 
2599  call messagenotify('M','zef_LaplaPol2Pol_zef',&
2600  'Matrix to apply b.c. newly produced.')
2601  endif
2602 
2603  e2z_work = e2a_aef(zef)
2604  e2z_work(:,1) = 0.0d0 ! 力学的条件
2605  e2z_work(:,km-1) = 0.0d0 ! 力学的条件
2606  e2z_work(:,0) = 0.0d0 ! 運動学的条件
2607  e2z_work(:,km) = 0.0d0 ! 運動学的条件
2608 
2609  e2z_work = lusolve(alu,kp,e2z_work)
2610 
2611  zef_laplapol2pol_zef = aef_e2a(e2z_work)
2612 
2613  end function zef_laplapol2pol_zef
2614 
2615  function tef_laplapol2polgrid_tef(tef,cond,new)
2616  !
2617  ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
2618  ! チェビシェフ格子点空間で境界条件を適用している.
2619  !
2620  ! この関数を用いるためには tef_Initial にて設定する
2621  ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく
2622  ! しておく必要がある.
2623  !
2624  ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
2625  !
2626  ! ▽^2Φ = f
2627  ! Φ = const. at boundaries.
2628  ! ∂Φ/∂z = 0 at boundaries (粘着条件)
2629  ! or ∂^2Φ/∂z^2 = 0 at boundaries (応力なし条件)
2630  !
2631  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2632  !
2633  ! 最終的にチェビシェフ係数の解が欲しい場合には, tef_LaplaPol2Pol_tef に
2634  ! 比べてチェビシェフ -- 格子点変換が 1 回分少なくて済む.
2635  !
2636  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(in) :: tef
2637  !(in) 入力▽^2φ分布
2638 
2639  real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaPol2PolGrid_tef
2640  !(out) 出力ポロイダルポテンシャル分布
2641 
2642  character(len=2), intent(in), optional :: cond
2643  !(in) 境界条件スイッチ. 省略時は 'RR'
2644  ! RR : 両端粘着条件
2645  ! RF : 上端粘着, 下端応力なし条件
2646  ! FR : 上端応力なし, 下端粘着条件
2647  ! FF : 両端応力なし条件
2648 
2649  logical, intent(IN), optional :: new
2650  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2651  ! default は false.
2652 
2653  real(8), dimension(:,:,:), allocatable :: alu
2654  integer, dimension(:,:), allocatable :: kp
2655 
2656  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2657  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2658  logical :: rigid1 = .true.
2659  logical :: rigid2 = .true.
2660 
2661  logical :: first = .true.
2662  logical :: new_matrix = .false.
2663  integer :: n
2664  save :: alu, kp, first
2665 
2666  select case (cond)
2667  case ('RR')
2668  rigid1 = .true. ; rigid2 = .true.
2669  case ('RF')
2670  rigid1 = .true. ; rigid2 = .false.
2671  case ('FR')
2672  rigid1 = .false. ; rigid2 = .true.
2673  case ('FF')
2674  rigid1 = .false. ; rigid2 = .false.
2675  case default
2676  call messagenotify('E','tef_LaplaPol2PolGrid_tef','B.C. not supported')
2677  end select
2678 
2679  if (.not. present(new)) then
2680  new_matrix=.false.
2681  else
2682  new_matrix=new
2683  endif
2684 
2685  if ( first .OR. new_matrix ) then
2686  first = .false.
2687 
2688  if ( nm /= km ) then
2689  call messagenotify('E','tef_LaplaPol2PolGrid_tef', &
2690  'Chebyshev truncation and number of grid points should be same.')
2691  endif
2692 
2693  if ( allocated(alu) ) deallocate(alu)
2694  if ( allocated(kp) ) deallocate(kp)
2695  allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2696 
2697  do n=0,nm
2698  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2699 
2700  ! 各水平波数に関して独立の式
2701  alu(:,:,n) = e2a_aef(zef_tef(tef_lapla_tef(aef_e2a(e2t_work))))
2702  enddo
2703 
2704  do n=0,nm
2705  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2706  e2z_work = az_at(e2t_work)
2707 
2708  ! 運動学的条件. 流線は境界で一定
2709  alu(:,0,n) = e2z_work(:,0)
2710  alu(:,km,n) = e2z_work(:,km)
2711 
2712  ! 力学的条件粘着条件
2713  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2714  if ( rigid1 ) then
2715  e2z_work=az_at(at_dz_at(e2t_work))
2716  else
2717  e2z_work=az_at(at_dz_at(at_dz_at(e2t_work)))
2718  endif
2719  alu(:,1,n) = e2z_work(:,0)
2720 
2721  ! 力学的条件粘着条件
2722  e2t_work = 0.0d0 ; e2t_work(:,n) = 1.0d0
2723  if ( rigid2 ) then
2724  e2z_work=az_at(at_dz_at(e2t_work))
2725  else
2726  e2z_work=az_at(at_dz_at(at_dz_at(e2t_work)))
2727  endif
2728  alu(:,km-1,n) = e2z_work(:,km)
2729  enddo
2730 
2731  call ludecomp(alu,kp)
2732 
2733  call messagenotify('M','tef_LaplaPol2PolGrid_tef',&
2734  'Matrix to apply b.c. newly produced.')
2735  endif
2736 
2737  e2z_work = az_at(e2a_aef(tef))
2738  e2z_work(:,1) = 0.0d0 ! 力学的条件
2739  e2z_work(:,km-1) = 0.0d0 ! 力学的条件
2740  e2z_work(:,0) = 0.0d0 ! 運動学的条件
2741  e2z_work(:,km) = 0.0d0 ! 運動学的条件
2742 
2743  tef_laplapol2polgrid_tef = aef_e2a(lusolve(alu,kp,e2z_work))
2744 
2745  end function tef_laplapol2polgrid_tef
2746 
2747  subroutine tef_tormagboundariestau(tef_TOR,new)
2748 
2749  ! 磁場トロイダルポテンシャルに対して境界条件を適用する.
2750  ! Chebyshev 空間での境界条件適用
2751  !
2752  ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
2753  ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ
2754  ! 対応している. その場合, 磁場トロイダルポテンシャルの境界条件は
2755  !
2756  ! 上側
2757  ! tef_psi = 0 at the outer boundary
2758  ! 下側
2759  ! tef_psi = 0 at the inner boundary
2760  !
2761  ! であるから tef_Boundaries で対応可能だが, 将来のため別途作成しておく.
2762  !
2763  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2764  !
2765  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR
2766  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2767 
2768  logical, intent(IN), optional :: new
2769  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2770  ! default は false.
2771 
2772  real(8), dimension(:,:,:), allocatable :: alu
2773  integer, dimension(:,:), allocatable :: kp
2774 
2775  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2776  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2777 
2778  logical :: first = .true.
2779  logical :: new_matrix = .false.
2780  integer :: n
2781  save :: alu, kp, first
2782 
2783  if (.not. present(new)) then
2784  new_matrix=.false.
2785  else
2786  new_matrix=new
2787  endif
2788 
2789  if ( first .OR. new_matrix ) then
2790  first = .false.
2791 
2792  if ( allocated(alu) ) deallocate(alu)
2793  if ( allocated(kp) ) deallocate(kp)
2794  allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2795 
2796  do n=0,nm
2797  e2t_work= 0.0d0 ; e2t_work(:,n)= 1.0d0
2798  alu(:,:,n) = e2t_work
2799 
2800  ! 非電気伝導体
2801  e2z_work = az_at(e2t_work)
2802  alu(:,nm-1,n) = e2z_work(:,0)
2803  alu(:,nm,n) = e2z_work(:,km)
2804 
2805  enddo
2806 
2807  call ludecomp(alu,kp)
2808  call messagenotify('M','tef_TormagBoundariesTau',&
2809  'Matrix to apply b.c. newly produced.')
2810  endif
2811 
2812  e2t_work = e2a_aef(tef_tor)
2813 
2814  e2t_work(:,nm-1) = 0.0d0
2815  e2t_work(:,nm) = 0.0d0
2816 
2817 
2818  tef_tor = aef_e2a(lusolve(alu,kp,e2t_work))
2819 
2820  end subroutine tef_tormagboundariestau
2821 
2822  subroutine tef_tormagboundariesgrid(tef_TOR,new)
2823  !
2824  ! 磁場トロイダルポテンシャルに対して境界条件を適用する.
2825  ! 鉛直実空間での境界条件適用.
2826  !
2827  ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
2828  ! 条件を課している(選点法). このルーチンを用いるためには
2829  ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を
2830  ! 等しくしておく必要がある.
2831  !
2832  ! 現在のところ境界物質が非電気伝導体の場合のみ対応している.
2833  ! その場合, 磁場トロイダルポテンシャルの境界条件は
2834  !
2835  ! 上側
2836  ! tef_psi = 0 at the outer boundary
2837  ! 下側
2838  ! tef_psi = 0 at the inner boundary
2839  !
2840  ! であるので tef_Boundaries で対応可能だが, 将来のため別途作成しておく
2841  !
2842  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2843  !
2844  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR
2845  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2846 
2847  logical, intent(IN), optional :: new
2848  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2849  ! default は false.
2850 
2851  real(8), dimension(:,:,:), allocatable :: alu
2852  integer, dimension(:,:), allocatable :: kp
2853 
2854  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2855  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2856 
2857  logical :: first = .true.
2858  logical :: new_matrix = .false.
2859  integer :: n
2860  save :: alu, kp, first
2861 
2862  if (.not. present(new)) then
2863  new_matrix=.false.
2864  else
2865  new_matrix=new
2866  endif
2867 
2868  if ( first .OR. new_matrix ) then
2869  first = .false.
2870 
2871  if ( nm /= km ) then
2872  call messagenotify('E','tef_TorMagBoundariesGrid', &
2873  'Chebyshev truncation and number of grid points should be same.')
2874  endif
2875 
2876  if ( allocated(alu) ) deallocate(alu)
2877  if ( allocated(kp) ) deallocate(kp)
2878  allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2879 
2880  do n=0,nm
2881  e2t_work = 0.0d0 ; e2t_work(:,n)=1.0d0
2882  e2z_work = az_at(e2t_work)
2883 
2884  alu(:,:,n) = e2z_work ! 内部領域は値そのまま.
2885 
2886  ! 非電気伝導体
2887  alu(:,0,n) = e2z_work(:,0)
2888  alu(:,km,n) = e2z_work(:,km)
2889 
2890  enddo
2891  call ludecomp(alu,kp)
2892  call messagenotify('M','TormagBoundariesGrid',&
2893  'Matrix to apply b.c. newly produced.')
2894  endif
2895 
2896  e2z_work = az_at(e2a_aef(tef_tor))
2897  e2z_work(:,0) = 0.0d0
2898  e2z_work(:,km) = 0.0d0
2899  tef_tor = aef_e2a(lusolve(alu,kp,e2z_work))
2900 
2901  end subroutine tef_tormagboundariesgrid
2902 
2903  subroutine tef_polmagboundariestau(tef_POL,new)
2904  !
2905  ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
2906  ! Chebyshev 空間での境界条件適用
2907  !
2908  ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
2909  ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ
2910  ! 対応している. その場合, 磁場ポロイダルポテンシャルの各水平スペクトル
2911  ! 成分 h にたいして境界条件が与えられ,
2912  !
2913  ! * 上側境界 : dh/dz + K h = 0
2914  ! * 下側境界 : dh/dz - K h = 0
2915  !
2916  ! である. ここで K=sqrt(l^2+m^2) は h の水平全波数である.
2917  !
2918  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
2919  !
2920  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_POL
2921  !(inout) 境界条件を適用するデータ. 修正された値を返す.
2922 
2923  logical, intent(IN), optional :: new
2924  !(in) true だと境界条件計算用行列を強制的に新たに作る.
2925  ! default は false.
2926 
2927  real(8), dimension(:,:,:), allocatable :: alu
2928  integer, dimension(:,:), allocatable :: kp
2929 
2930  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
2931  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
2932 
2933  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_kh
2934 
2935  logical :: first = .true.
2936  logical :: new_matrix = .false.
2937  real(8) :: rl, rm
2938  integer :: l, m, n, e2index
2939  save :: alu, kp, first
2940 
2941  if (.not. present(new)) then
2942  new_matrix=.false.
2943  else
2944  new_matrix=new
2945  endif
2946 
2947  if ( first .OR. new_matrix ) then
2948  first = .false.
2949 
2950  if ( allocated(alu) ) deallocate(alu)
2951  if ( allocated(kp) ) deallocate(kp)
2952 
2953  allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
2954 
2955  do n=0,nm
2956  e2t_work = 0.0d0 ; e2t_work(:,n)=1.0d0
2957 
2958  alu(:,:,n) = e2t_work(:,:) ! 内部領域は同じ
2959  enddo
2960 
2961  ! 非電気伝導体
2962  !e2t_kh = 0.0D0
2963  do m=-mm,mm
2964  do l=-lm,lm
2965  if ( lf_l(l) > 0 ) then
2966  e2index = (2*mm+1)*(lf_l(l)-1) + (m+mm+1)
2967  rl = 2*pi/xl*l ; rm = 2*pi/yl*m
2968  e2t_kh(e2index,:) = sqrt(rl**2+rm**2)
2969  endif
2970  enddo
2971  if ( lf_l(0) > 0 ) then
2972  ! e2index1 = (2*mm+1)*(lf_l(0)-1) + (m+mm+1)
2973  ! e2index2 = (2*mm+1)*(2*lc(ip))-1) + (m+mm+1)
2974  ! e2t_kh(e2index2,:) = e2t_kh(e2index1,:)
2975  e2t_kh((2*mm+1)*(2*lc(ip)-1)+(m+mm+1),:) = e2t_kh(m+mm+1,:)
2976  end if
2977  enddo
2978  do n=0,nm
2979  e2t_work = 0.0d0 ; e2t_work(:,n)=1.0d0
2980 
2981  e2z_work = az_at(at_dz_at(e2t_work) + e2t_kh * e2t_work)
2982  alu(:,nm-1,n) = e2z_work(:,0)
2983 
2984  e2z_work = az_at(at_dz_at(e2t_work) - e2t_kh * e2t_work)
2985  alu(:,nm,n) = e2z_work(:,km)
2986  enddo
2987 
2988  call ludecomp(alu,kp)
2989 
2990  call messagenotify('M','tef_PolmagBoundariesTau',&
2991  'Matrix to apply b.c. newly produced.')
2992  endif
2993 
2994  e2t_work = e2a_aef(tef_pol)
2995  e2t_work(:,nm-1) = 0.0d0
2996  e2t_work(:,nm) = 0.0d0
2997  tef_pol = aef_e2a(lusolve(alu,kp,e2t_work))
2998 
2999  end subroutine tef_polmagboundariestau
3000 
3001  subroutine tef_polmagboundariesgrid(tef_POL,new)
3002  !
3003  ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
3004  ! 鉛直実空間での境界条件適用.
3005  !
3006  ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
3007  ! 条件を課している(選点法). このルーチンを用いるためには
3008  ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を
3009  ! 等しくしておく必要がある.
3010  !
3011  ! 現在のところ境界物質が非電気伝導体の場合のみ対応している.
3012  ! その場合, 磁場ポロイダルポテンシャルの各水平スペクトル成分 h に
3013  ! たいして境界条件が与えられ,
3014  !
3015  ! * 上側境界 : dh/dz + K h = 0
3016  ! * 下側境界 : dh/dz - K h = 0
3017  !
3018  ! である. ここで K=sqrt(l^2+m^2) は h の水平全波数である.
3019  !
3020  ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
3021  !
3022  real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_POL
3023  !(inout) 境界条件を適用するデータ. 修正された値を返す.
3024 
3025  logical, intent(IN), optional :: new
3026  !(in) true だと境界条件計算用行列を強制的に新たに作る.
3027  ! default は false.
3028 
3029  real(8), dimension(:,:,:), allocatable :: alu
3030  integer, dimension(:,:), allocatable :: kp
3031 
3032  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
3033  real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work
3034 
3035  real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_kh
3036 
3037  logical :: first = .true.
3038  logical :: new_matrix = .false.
3039  real(8) :: rl, rm
3040  integer :: l, m, n, e2index
3041  save :: alu, kp, first
3042 
3043  if (.not. present(new)) then
3044  new_matrix=.false.
3045  else
3046  new_matrix=new
3047  endif
3048 
3049  if ( first .OR. new_matrix ) then
3050  first = .false.
3051 
3052  if ( nm /= km ) then
3053  call messagenotify('E','tef_PolMagBoundariesGrid', &
3054  'Chebyshev truncation and number of grid points should be same.')
3055  endif
3056 
3057  if ( allocated(alu) ) deallocate(alu)
3058  if ( allocated(kp) ) deallocate(kp)
3059  allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))
3060 
3061  do n=0,nm
3062  e2t_work = 0.0d0 ; e2t_work(:,n)=1.0d0
3063  e2z_work = az_at(e2t_work)
3064 
3065  alu(:,:,n) = e2z_work ! 内部領域は同じ
3066  enddo
3067 
3068  ! 非電気伝導体
3069  do m=-mm,mm
3070  do l=-lm,lm
3071  if ( lf_l(l) > 0 ) then
3072  e2index = (2*mm+1)*(lf_l(l)-1) + (m+mm+1)
3073  rl = 2*pi/xl*l ; rm = 2*pi/yl*m
3074  e2t_kh(e2index,:) = sqrt(rl**2+rm**2)
3075  endif
3076  enddo
3077  if ( lf_l(0) > 0 ) then
3078  ! e2index1 = (2*mm+1)*(lf_l(0)-1) + (m+mm+1)
3079  ! e2index2 = (2*mm+1)*(2*lc(ip))-1) + (m+mm+1)
3080  ! e2t_kh(e2index2,:) = e2t_kh(e2index1,:)
3081  e2t_kh((2*mm+1)*(2*lc(ip)-1)+(m+mm+1),:) = e2t_kh(m+mm+1,:)
3082  end if
3083  enddo
3084  do n=0,nm
3085  e2t_work = 0.0d0 ; e2t_work(:,n)=1.0d0
3086 
3087  e2z_work = az_at(at_dz_at(e2t_work) + e2t_kh * e2t_work)
3088  alu(:,0,n) = e2z_work(:,0)
3089 
3090  e2z_work = az_at(at_dz_at(e2t_work) - e2t_kh * e2t_work)
3091  alu(:,km,n) = e2z_work(:,km)
3092  enddo
3093 
3094  call ludecomp(alu,kp)
3095 
3096  call messagenotify('M','tef_PolmagBoundariesGrid',&
3097  'Matrix to apply b.c. newly produced.')
3098  endif
3099 
3100  e2z_work = az_at(e2a_aef(tef_pol))
3101  e2z_work(:,0) = 0.0d0
3102  e2z_work(:,km) = 0.0d0
3103  tef_pol = aef_e2a(lusolve(alu,kp,e2z_work))
3104 
3105  end subroutine tef_polmagboundariesgrid
3106 
3107 end module tee_mpi_module
at_module::at_dx_at
real(dp) function, dimension(size(at_data, 1), 0:size(at_data, 2) -1), public at_dx_at(at_data)
入力チェビシェフデータに X 微分を作用する(2 次元配列用).
Definition: at_module.f90:420
tee_mpi_module::zef_laplapol2pol_zef
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_laplapol2pol_zef(zef, cond, new)
Definition: tee_mpi_module.f90:2486
tee_mpi_module::tef_zef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_zef(zef)
Definition: tee_mpi_module.f90:533
ee_mpi_module::kc
integer, dimension(:), allocatable, save, public kc
Definition: ee_mpi_module.f90:198
ee_mpi_module::v_intx_vx
real(8) function, dimension(jc(ip)), public v_intx_vx(vx)
Definition: ee_mpi_module.f90:773
ee_mpi_module::je
integer, dimension(:), allocatable, save, public je
Definition: ee_mpi_module.f90:199
ee_mpi_module::inty_v
real(8) function, public inty_v(v)
Definition: ee_mpi_module.f90:832
at_module::at_ag
real(dp) function, dimension(size(ag_data, 1), 0:km), public at_ag(ag_data)
格子データからチェビシェフデータへ変換する(2 次元配列用).
Definition: at_module.f90:374
tee_mpi_module::tee_vmiss
real(8), public tee_vmiss
Definition: tee_mpi_module.f90:373
tee_mpi_module::zvx_x
real(8), dimension(:,:,:), allocatable, public zvx_x
Definition: tee_mpi_module.f90:368
tee_mpi_module::e2t_e2z
real(8) function, dimension(2 *lc(ip) *(2 *mm+1), 0:nm), public e2t_e2z(e2z)
Definition: tee_mpi_module.f90:559
tee_mpi_module::zvx_y
real(8), dimension(:,:,:), allocatable, public zvx_y
Definition: tee_mpi_module.f90:368
ee_mpi_module::kf_k
integer function, public kf_k(k)
Definition: ee_mpi_module.f90:337
tee_mpi_module::x_avrz_zx
real(8) function, dimension(0:im-1), public x_avrz_zx(zx)
Definition: tee_mpi_module.f90:1486
tee_mpi_module::zef_toroidalenergyspectrum_tef
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_toroidalenergyspectrum_tef(tef_TORPOT)
Definition: tee_mpi_module.f90:1764
tee_mpi_module::intzy_zv
real(8) function, public intzy_zv(zv)
Definition: tee_mpi_module.f90:1221
tee_mpi_module::tef_dx_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_dx_tef(tef)
Definition: tee_mpi_module.f90:692
tee_mpi_module::zx_avry_zvx
real(8) function, dimension(0:km, 0:im-1), public zx_avry_zvx(zvx)
Definition: tee_mpi_module.f90:1343
tee_mpi_module::avrzyx_zvx
real(8) function, public avrzyx_zvx(zvx)
Definition: tee_mpi_module.f90:1414
at_module
1 次元有限領域: チェビシェフ変換
Definition: at_module.f90:78
tee_mpi_module::zv_avrx_zvx
real(8) function, dimension(0:km, jc(ip)), public zv_avrx_zvx(zvx)
Definition: tee_mpi_module.f90:1329
ee_mpi_module::x_inty_vx
real(8) function, dimension(0:im-1), public x_inty_vx(vx)
Definition: ee_mpi_module.f90:794
tee_mpi_module::zef_toroidalenstrophyspectrum_tef
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_toroidalenstrophyspectrum_tef(tef_TORPOT)
Definition: tee_mpi_module.f90:1913
tee_mpi_module::tef_laplainvtaunn_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_laplainvtaunn_tef(tef, new)
Definition: tee_mpi_module.f90:922
tee_mpi_module::z_avryx_zvx
real(8) function, dimension(0:km), public z_avryx_zvx(zvx)
Definition: tee_mpi_module.f90:1400
tee_mpi_module::x_intz_zx
real(8) function, dimension(0:im-1), public x_intz_zx(zx)
Definition: tee_mpi_module.f90:1267
at_module::t_boundaries
subroutine, public t_boundaries(t_data, cond, values)
境界条件の適用(タウ法, 1 次元配列用) 境界条件と両境界での値をオプションで与える. デフォルトは DD...
Definition: at_module.f90:664
tee_mpi_module::aef_e2a
real(8) function, dimension(size(e2a, 2),-mm:mm, 2 *lc(ip)), public aef_e2a(e2a)
Definition: tee_mpi_module.f90:606
tee_mpi_module::tef_boundariestau
subroutine, public tef_boundariestau(tef, values, cond)
Definition: tee_mpi_module.f90:1997
ee_mpi_module::v_y_weight
real(8), dimension(:), allocatable, save, public v_y_weight
Definition: ee_mpi_module.f90:214
tee_mpi_module::tef_zrotrot_zvx_zvx_zvx
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_zrotrot_zvx_zvx_zvx(zvx_VX, zvx_VY, zvx_VZ)
Definition: tee_mpi_module.f90:1646
tee_mpi_module::v_avrzx_zvx
real(8) function, dimension(jc(ip)), public v_avrzx_zvx(zvx)
Definition: tee_mpi_module.f90:1385
tee_mpi_module::zvx_tef
real(8) function, dimension(0:km, jc(ip), 0:im-1), public zvx_tef(tef)
Definition: tee_mpi_module.f90:444
tee_mpi_module::zef_poloidalenergyspectrum_tef
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_poloidalenergyspectrum_tef(tef_POLPOT)
Definition: tee_mpi_module.f90:1829
tee_mpi_module::e2_ef
real(8) function, dimension(2 *lc(ip) *(2 *mm+1)), public e2_ef(ef)
Definition: tee_mpi_module.f90:638
tee_mpi_module::x_avrzy_zvx
real(8) function, dimension(0:im-1), public x_avrzy_zvx(zvx)
Definition: tee_mpi_module.f90:1371
ee_mpi_module::ks
integer, dimension(:), allocatable, save, public ks
Definition: ee_mpi_module.f90:197
ee_mpi_module::js
integer, dimension(:), allocatable, save, public js
Definition: ee_mpi_module.f90:200
tee_mpi_module::tef_tormagboundariesgrid
subroutine, public tef_tormagboundariesgrid(tef_TOR, new)
Definition: tee_mpi_module.f90:2823
ee_mpi_module::avry_v
real(8) function, public avry_v(v)
Definition: ee_mpi_module.f90:917
tee_mpi_module::zef_tef
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_tef(tef)
Definition: tee_mpi_module.f90:520
tee_mpi_module::e2z_e2t
real(8) function, dimension(2 *lc(ip) *(2 *mm+1), 0:km), public e2z_e2t(e2t)
Definition: tee_mpi_module.f90:546
tee_mpi_module::zef_poloidalenstrophyspectrum_tef
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_poloidalenstrophyspectrum_tef(tef_POLPOT)
Definition: tee_mpi_module.f90:1956
ee_mpi_module::v_y
real(8), dimension(:), allocatable, save, public v_y
Definition: ee_mpi_module.f90:211
tee_mpi_module::tef_boundariesgrid
subroutine, public tef_boundariesgrid(tef, values, cond)
Definition: tee_mpi_module.f90:2069
ee_mpi_module::v_avrx_vx
real(8) function, dimension(jc(ip)), public v_avrx_vx(vx)
Definition: ee_mpi_module.f90:869
tee_mpi_module::v_intzx_zvx
real(8) function, dimension(jc(ip)), public v_intzx_zvx(zvx)
Definition: tee_mpi_module.f90:1102
tee_mpi_module::tef_div_zvx_zvx_zvx
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_div_zvx_zvx_zvx(zvx_VX, zvx_VY, zvx_VZ)
Definition: tee_mpi_module.f90:822
tee_mpi_module::ef_e2
real(8) function, dimension(-mm:mm, 2 *lc(ip)), public ef_e2(e2)
Definition: tee_mpi_module.f90:664
ee_mpi_module::x_x
real(8), dimension(:), allocatable, save, public x_x
Definition: ee_mpi_module.f90:210
tee_mpi_module::vx_intz_zvx
real(8) function, dimension(jc(ip), 0:im-1), public vx_intz_zvx(zvx)
Definition: tee_mpi_module.f90:1057
at_module::at_boundaries
subroutine, public at_boundaries(at_data, cond, values)
境界条件の適用(タウ法, 2 次元配列用)
Definition: at_module.f90:619
tee_mpi_module::tee_mpi_initial
subroutine, public tee_mpi_initial(i, j, k, l, m, n, xmin, xmax, ymin, ymax, zmin, zmax)
Definition: tee_mpi_module.f90:383
ee_mpi_module::intx_x
real(8) function, public intx_x(x)
Definition: ee_mpi_module.f90:820
tee_mpi_module::zv_z
real(8), dimension(:,:), allocatable, public zv_z
Definition: tee_mpi_module.f90:369
ee_mpi_module::interpolate_ef
real(8) function, public interpolate_ef(ef_Data, x, y)
Definition: ee_mpi_module.f90:998
tee_mpi_module::zvx_zef
real(8) function, dimension(0:km, jc(ip), 0:im-1), public zvx_zef(zef)
Definition: tee_mpi_module.f90:470
tee_mpi_module
Definition: tee_mpi_module.f90:53
ee_mpi_module
Definition: ee_mpi_module.f90:21
tee_mpi_module::zx_z
real(8), dimension(:,:), allocatable, public zx_z
Definition: tee_mpi_module.f90:370
at_module::a_interpolate_at
real(dp) function, dimension(size(at_data, 1)), public a_interpolate_at(at_data, xval)
1 次元データの並んだ 2次元配列の補間
Definition: at_module.f90:576
tee_mpi_module::z_intyx_zvx
real(8) function, dimension(0:km), public z_intyx_zvx(zvx)
Definition: tee_mpi_module.f90:1123
tee_mpi_module::tef_potential2rotation
subroutine, public tef_potential2rotation(zvx_RotVX, zvx_RotVY, zvx_RotVZ, tef_Torvel, tef_Polvel)
Definition: tee_mpi_module.f90:1710
tee_mpi_module::v_avrz_zv
real(8) function, dimension(jc(ip)), public v_avrz_zv(zv)
Definition: tee_mpi_module.f90:1443
tee_mpi_module::zx_x
real(8), dimension(:,:), allocatable, public zx_x
Definition: tee_mpi_module.f90:370
at_module::t_g
real(dp) function, dimension(0:km), public t_g(g_data)
格子データからチェビシェフデータへ変換する(1 次元配列用).
Definition: at_module.f90:399
ee_mpi_module::ef_vx
real(8) function, dimension(-lm:lm, 2 *kc(ip)), public ef_vx(vx)
Definition: ee_mpi_module.f90:408
at_module::g_x
real(dp), dimension(:), allocatable, save, public g_x
格子点座標(X)を格納した 1 次元配列
Definition: at_module.f90:218
tee_mpi_module::tef_laplah_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_laplah_tef(tef)
Definition: tee_mpi_module.f90:772
tee_mpi_module::tef_laplapol2polgrid_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_laplapol2polgrid_tef(tef, cond, new)
Definition: tee_mpi_module.f90:2616
tee_mpi_module::avrzx_zx
real(8) function, public avrzx_zx(zx)
Definition: tee_mpi_module.f90:1500
tee_mpi_module::x_intzy_zvx
real(8) function, dimension(0:im-1), public x_intzy_zvx(zvx)
Definition: tee_mpi_module.f90:1076
ee_mpi_module::intyx_vx
real(8) function, public intyx_vx(vx)
Definition: ee_mpi_module.f90:744
tee_mpi_module::tef_dz_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_dz_tef(tef)
Definition: tee_mpi_module.f90:734
tee_mpi_module::zv_y
real(8), dimension(:,:), allocatable, public zv_y
Definition: tee_mpi_module.f90:369
tee_mpi_module::zv_intx_zvx
real(8) function, dimension(0:km, jc(ip)), public zv_intx_zvx(zvx)
Definition: tee_mpi_module.f90:1014
tee_mpi_module::vx_avrz_zvx
real(8) function, dimension(jc(ip), 0:im-1), public vx_avrz_zvx(zvx)
Definition: tee_mpi_module.f90:1357
at_module::ag_at
real(dp) function, dimension(size(at_data, 1), 0:im), public ag_at(at_data)
チェビシェフデータから格子データへ変換する(2 次元配列用).
Definition: at_module.f90:328
at_module::at_initial
subroutine, public at_initial(i, k, xmin, xmax)
チェビシェフ変換の格子点数, 波数, 領域の大きさを設定する.
Definition: at_module.f90:268
at_module::t_dx_t
real(dp) function, dimension(size(t_data)), public t_dx_t(t_data)
入力チェビシェフデータに X 微分を作用する(1 次元配列用).
Definition: at_module.f90:467
tee_mpi_module::e2a_aef
real(8) function, dimension(2 *lc(ip) *(2 *mm+1), size(aef, 1)), public e2a_aef(aef)
Definition: tee_mpi_module.f90:572
tee_mpi_module::tef_dy_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_dy_tef(tef)
Definition: tee_mpi_module.f90:713
ee_mpi_module::x_avry_vx
real(8) function, dimension(0:im-1), public x_avry_vx(vx)
Definition: ee_mpi_module.f90:886
tee_mpi_module::tef_tormagboundariestau
subroutine, public tef_tormagboundariestau(tef_TOR, new)
Definition: tee_mpi_module.f90:2748
at_module::g_x_weight
real(dp), dimension(:), allocatable, save, public g_x_weight
重み座標を格納した 1 次元配列
Definition: at_module.f90:220
tee_mpi_module::z_avrx_zx
real(8) function, dimension(0:km), public z_avrx_zx(zx)
Definition: tee_mpi_module.f90:1472
at_module::g_t
real(dp) function, dimension(0:im), public g_t(t_data)
チェビシェフデータから格子データへ変換する(1 次元配列用).
Definition: at_module.f90:357
ee_mpi_module::avryx_vx
real(8) function, public avryx_vx(vx)
Definition: ee_mpi_module.f90:852
ee_mpi_module::vx_ef
real(8) function, dimension(js(ip):je(ip), 0:im-1), public vx_ef(ef)
Definition: ee_mpi_module.f90:359
tee_mpi_module::zef_z
real(8), dimension(:,:,:), allocatable, public zef_z
Definition: tee_mpi_module.f90:371
tee_mpi_module::tef_zvx
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_zvx(zvx)
Definition: tee_mpi_module.f90:457
tee_mpi_module::tef_torboundariesgrid
subroutine, public tef_torboundariesgrid(tef_TOR, cond, new)
Definition: tee_mpi_module.f90:2247
tee_mpi_module::zvx_z
real(8), dimension(:,:,:), allocatable, public zvx_z
Definition: tee_mpi_module.f90:368
tee_mpi_module::tef_polmagboundariestau
subroutine, public tef_polmagboundariestau(tef_POL, new)
Definition: tee_mpi_module.f90:2904
tee_mpi_module::z_avry_zv
real(8) function, dimension(0:km), public z_avry_zv(zv)
Definition: tee_mpi_module.f90:1429
ee_mpi_module::vx_x
real(8), dimension(:,:), allocatable, save, public vx_x
Definition: ee_mpi_module.f90:216
tee_mpi_module::tef_polmagboundariesgrid
subroutine, public tef_polmagboundariesgrid(tef_POL, new)
Definition: tee_mpi_module.f90:3002
ee_mpi_module::vx_y
real(8), dimension(:,:), allocatable, save, public vx_y
Definition: ee_mpi_module.f90:217
tee_mpi_module::tef_laplahinv_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_laplahinv_tef(tef)
Definition: tee_mpi_module.f90:797
ee_mpi_module::ef_laplainv_ef
real(8) function, dimension(-lm:lm, 2 *kc(ip)), public ef_laplainv_ef(ef)
Definition: ee_mpi_module.f90:594
ee_mpi_module::ef_lapla_ef
real(8) function, dimension(-lm:lm, 2 *kc(ip)), public ef_lapla_ef(ef)
Definition: ee_mpi_module.f90:564
ee_mpi_module::avrx_x
real(8) function, public avrx_x(x)
Definition: ee_mpi_module.f90:903
ee_mpi_module::x_x_weight
real(8), dimension(:), allocatable, save, public x_x_weight
Definition: ee_mpi_module.f90:213
tee_mpi_module::interpolate_tef
real(8) function, public interpolate_tef(tef_data, x, y, z)
Definition: tee_mpi_module.f90:1746
tee_mpi_module::intz_z
real(8) function, public intz_z(z)
Definition: tee_mpi_module.f90:1309
tee_mpi_module::tef_zrot_zvx_zvx
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_zrot_zvx_zvx(zvx_VX, zvx_VY)
Definition: tee_mpi_module.f90:1622
tee_mpi_module::z_inty_zv
real(8) function, dimension(0:km), public z_inty_zv(zv)
Definition: tee_mpi_module.f90:1179
tee_mpi_module::tef_lapla_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_lapla_tef(tef)
Definition: tee_mpi_module.f90:751
tee_mpi_module::tef_torboundariestau
subroutine, public tef_torboundariestau(tef_TOR, cond, new)
Definition: tee_mpi_module.f90:2143
tee_mpi_module::z_intx_zx
real(8) function, dimension(0:km), public z_intx_zx(zx)
Definition: tee_mpi_module.f90:1248
tee_mpi_module::zx_inty_zvx
real(8) function, dimension(0:km, 0:im-1), public zx_inty_zvx(zvx)
Definition: tee_mpi_module.f90:1033
tee_mpi_module::intzx_zx
real(8) function, public intzx_zx(zx)
Definition: tee_mpi_module.f90:1287
tee_mpi_module::avrz_z
real(8) function, public avrz_z(z)
Definition: tee_mpi_module.f90:1514
ee_mpi_module::ee_mpi_initial
subroutine, public ee_mpi_initial(i_in, j_in, k_in, l_in, xmin_in, xmax_in, ymin_in, ymax_in)
Definition: ee_mpi_module.f90:233
ee_mpi_module::jc
integer, dimension(:), allocatable, save, public jc
Definition: ee_mpi_module.f90:201
tee_mpi_module::avrzy_zv
real(8) function, public avrzy_zv(zv)
Definition: tee_mpi_module.f90:1457
ee_mpi_module::ke
integer, dimension(:), allocatable, save, public ke
Definition: ee_mpi_module.f90:196
tee_mpi_module::tef_potential2vector
subroutine, public tef_potential2vector(zvx_VX, zvx_VY, zvx_VZ, tef_Torvel, tef_Polvel)
Definition: tee_mpi_module.f90:1678
tee_mpi_module::tef_laplainvtaudd_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_laplainvtaudd_tef(tef, new)
Definition: tee_mpi_module.f90:844
tee_mpi_module::intzyx_zvx
real(8) function, public intzyx_zvx(zvx)
Definition: tee_mpi_module.f90:1149
tee_mpi_module::v_intz_zv
real(8) function, dimension(jc(ip)), public v_intz_zv(zv)
Definition: tee_mpi_module.f90:1202
ee_mpi_module::ef_dy_ef
real(8) function, dimension(-lm:lm, 2 *kc(ip)), public ef_dy_ef(ef)
Definition: ee_mpi_module.f90:655
ee_mpi_module::ef_dx_ef
real(8) function, dimension(-lm:lm, 2 *kc(ip)), public ef_dx_ef(ef)
Definition: ee_mpi_module.f90:627
tee_mpi_module::tef_laplapol2poltau_tef
real(8) function, dimension(0:nm,-mm:mm, 2 *lc(ip)), public tef_laplapol2poltau_tef(tef, cond, new)
Definition: tee_mpi_module.f90:2361
tee_mpi_module::zef_zvx
real(8) function, dimension(0:km,-mm:mm, 2 *lc(ip)), public zef_zvx(zvx)
Definition: tee_mpi_module.f90:495