!-- !---------------------------------------------------------------------- ! Copyright (c) 2010--2012 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- !表題 wtq_module_sjpack ! ! spml/wtq_module モジュールは球とその外側の球殻領域での流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向に球面調和函数変換および動径方向にチェビシェフ変換ならびに ! Matsushima and Marcus (1994) で提唱された多項式を用いる場合の ! スペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wt_module_sjpack, wa_module_sjpack を用いている. ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして ! ISPACK の Fortran77 サブルーチンを用いている. ! ! !履歴 2010/04/17 竹広真一 wtu_module より改変 ! 2010/04/18 竹広真一 wtq_module より改変 ! 2010/04/20 竹広真一 wq_Boundary, wq_TorBoundary, ! wq_TormagBoundary, wq_PolmagBoundary を ! use 文から削除(for intel fortran) ! 2010/04/28 佐々木洋平 threadprivate 分削除 ! 2012/07/08 竹広真一 wt_LaplaPol2PolTau_wt 追加 ! 2012/08/24 竹広真一 wt_LaplaPol2Pol_wt 追加 ! !凡例 ! データ種類と index ! x : 経度 y : 緯度 z : 動径(球殻内) r : 動径(球内) ! w : 球面調和関数スペクトル ! n : 球面調和関数スペクトル(水平全波数) ! m : 球面調和関数スペクトル(帯状波数) ! t : チェビシェフ関数スペクトル(球殻内) ! q : 動径スペクトル関数スペクトル(球内) ! a : 任意の次元 ! ! xyz : 球殻内 3 次元格子点データ ! xyr : 球内の 3 次元格子点データ ! xy : 水平 2 次元格子点データ ! yr : 子午面 2 次元格子点データ ! xr : 緯度面 2 次元格子点データ ! ! wz : 水平スペクトル動径格子点データ(球殻内) ! wr : 水平スペクトル動径格子点データ(球内) ! wt : スペクトルデータ(球殻内) ! wq : スペクトルデータ(球内) ! !++ module wtq_module_sjpack ! != wtq_module_sjpack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/wtq_module モジュールは球とその外側の球殻領域での流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向に球面調和函数変換および動径方向にチェビシェフ変換ならびに ! Matsushima and Marcus (1994) で提唱された多項式を用いる場合の ! スペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wt_module_sjpack, wa_module_sjpack を用いている. ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして ! ISPACK の Fortran77 サブルーチンを用いている. ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 ! (wq_, nmr_, nr_, xyr_, wr_, w_, xy_, x_, y_, r_, a_, wt_, nmz_, nz_, ! xyz_, wz_, w_, xy_, x_, y_, z_, a_) は, 返す値の形を示している. ! wt_ :: スペクトルデータ(球面調和函数・チェビシェフ変換-球殻) ! nmz_ :: 水平スペクトルデータ(全波数 n, 帯状波数各成分, 動径-球殻) ! nz_ :: 水平スペクトルデータ(全波数 n, 動径-球殻) ! xyz_ :: 3 次元格子点データ(経度・緯度・動径-球殻) ! wz_ :: 水平スペクトル, 動径格子点データ, 球殻 ! wq_ :: スペクトルデータ(球面調和函数・Matsushima Marcus 多項式) ! nmr_ :: 水平スペクトルデータ(全波数 n, 帯状波数各成分, 動径-球) ! nr_ :: 水平スペクトルデータ(全波数 n, 動径-球) ! xyr_ :: 3 次元格子点データ(経度・緯度・動径-球) ! wr_ :: 水平スペクトル, 動径格子点データ, 球 ! ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 ! (wt_, xyz_, wz_, w_, xy_, x_, y_, z_, a_, ! wq_, xyz_, wr_, w_, xy_, x_, y_, r_, a_) は, 入力変数の形が ! スペクトルデータおよび格子点データであることを示している. ! _wt :: スペクトルデータ ! _xyz :: 3 次元格子点データ ! _xyz_xyz :: 2 つの3 次元格子点データ, ... ! _wq :: スペクトルデータ ! _xyr :: 3 次元格子点データ ! _xyr_xyr :: 2 つの3 次元格子点データ, ... ! !=== 各データの種類の説明 ! ! ! * xyz : 3 次元格子点データ(経度・緯度・動径-球殻) ! * 変数の種類と次元は real(8), dimension(0:im-1,1:jm,0:km). ! * im, jm, km はそれぞれ経度, 緯度, 動径座標の格子点数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wt : スペクトルデータ-球殻 ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),0:lm). ! * nm は球面調和函数の最大全波数, lm はチェビシェフ多項式の最大次数 ! であり, サブルーチン wt_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * nmz : 水平スペクトルデータの並んだ 3 次元配列(球殻). ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,0:km). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wt_Initial にて ! あらかじめ設定しておく. ! ! * nz : スペクトルデータの並んだ 2 次元配列(球殻). ! * 変数の種類と次元は real(8), dimension(0:nm,0:km). ! * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wz : 水平スペクトル, 動径格子点データ(球殻). ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),0:km). ! ! * xyr : 3 次元格子点データ(経度・緯度・動径-球) ! * 変数の種類と次元は real(8), dimension(0:im-1,1:jm,0:km). ! * im, jm, km はそれぞれ経度, 緯度, 動径座標の格子点数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wq : スペクトルデータ(球) ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),0:lm). ! * nm は球面調和函数の最大全波数, lm はチェビシェフ多項式の最大次数 ! であり, サブルーチン wtq_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * nmr : 水平スペクトルデータの並んだ 3 次元配列(球). ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,0:km). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wt_Initial にて ! あらかじめ設定しておく. ! ! * nr : スペクトルデータの並んだ 2 次元配列(球). ! * 変数の種類と次元は real(8), dimension(0:nm,0:km). ! * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wr : 水平スペクトル, 動径格子点データ(球). ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),0:km). ! ! * wt_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyz_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * wz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * wq_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyr_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * wr_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! wtq_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_Lon, y_Lat :: 格子点座標(緯度, 経度)を格納した1 次元配列 ! z_Rad :: 格子点座標(動径, 球殻)を格納した1 次元配列 ! r_Rad :: 格子点座標(動径, 球)を格納した1 次元配列 ! x_Lon_Weight, y_Lat_Weight :: 重み座標を格納した 1 次元配列 ! z_Rad_Weight :: 重み座標を格納した 1 次元配列 ! r_Rad_Weight :: 重み座標を格納した 1 次元配列 ! xyz_Lon, xyz_Lat :: 格子点データの経度・緯度(X,Y) ! (格子点データ型 3 次元配列) ! xyz_Rad :: 格子点データの動径(球殻)(Z) ! (格子点データ型 3 次元配列) ! xyr_Lon, xyr_Lat :: 格子点データの経度・緯度(X,Y) ! (格子点データ型 3 次元配列) ! xyr_Rad :: 格子点データの動径(球) ! (格子点データ型 3 次元配列) ! !==== 基本変換 ! ! xyz_wt, wt_xyz :: スペクトルデータと 3 次元格子データの間の変換 ! (球面調和函数, チェビシェフ変換) ! xyz_wz, wz_xyz :: 3 次元格子データと水平スペクトル・動径格子データとの ! 間の変換 (球面調和函数) ! wz_wt, wt_wz :: スペクトルデータと水平スペクトル・動径格子データとの ! 間の変換 (チェビシェフ変換) ! xyr_wq, wq_xyr :: スペクトルデータと 3 次元格子データの ! 間の変換 (球面調和函数, チェビシェフ変換) ! xyr_wr, wr_xyr :: 3 次元格子データと水平スペクトル・動径格子データとの ! 間の変換 (球面調和函数) ! wr_wq, wq_wr :: スペクトルデータと水平スペクトル・動径格子データとの ! 間の変換 (M.M. 多項式変換) ! w_xy, xy_w :: スペクトルデータと 2 次元水平格子データの ! 間の変換(球面調和函数変換) ! l_nm, nm_l :: スペクトルデータの格納位置と全波数・帯状波数の変換 ! !==== 微分 ! ! wt_DRad_wt :: スペクトルデータに動径微分 ! ∂/∂r を作用させる ! wt_DivRad_wt :: スペクトルデータに発散型動径微分 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! wt_RotRad_wt :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wt_Lapla_wt :: スペクトルデータにラプラシアンを作用させる ! xyz_GradLon_wt :: スペクトルデータに勾配型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! xyz_GradLat_wt :: スペクトルデータに勾配型緯度微分 ! 1/r・∂/∂φを作用させる ! wt_DivLon_xyz :: 格子データに発散型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! wt_DivLat_xyz :: 格子データに発散型緯度微分 ! 1/rcosφ・∂(g cosφ)/∂φを作用させる ! wt_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyz_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyz_RotLon_wt_wt :: ベクトル場の回転の経度成分を計算する ! xyz_RotLat_wt_wt :: ベクトル場の回転の緯度成分を計算する ! wt_RotRad_xyz_xyz :: ベクトル場の回転の動径成分を計算する ! ! wq_RadDRad_wq :: スペクトルデータに動径微分 r∂/∂r を作用させる ! wr_DivRad_wq :: スペクトルデータに発散型動径微分 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! wr_DivRad_xyr :: スペクトルデータに発散型動径微分 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! wr_RotDRad_wq :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wr_RotDRad_wr :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wq_RotDRad_wr :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wq_Lapla_wq :: スペクトルデータにラプラシアンを作用させる ! xyr_GradLon_wq :: スペクトルデータに勾配型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! xyr_GradLat_wq :: スペクトルデータに勾配型緯度微分 ! 1/r・∂/∂φを作用させる ! wr_DivLon_xyr :: 格子データに発散型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! wr_DivLat_xyr :: 格子データに発散型緯度微分 ! 1/rcosφ・∂(g cosφ)/∂φを作用させる ! wr_Div_xyr_xyr_xyr :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyr_Div_xyr_xyr_xyr :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyr_RotLon_wq_wq :: ベクトル場の回転の経度成分を計算する ! xyr_RotLat_wq_wq :: ベクトル場の回転の緯度成分を計算する ! wr_RotRad_xyr_xyr :: ベクトル場の回転の動径成分を計算する ! !==== トロイダルポロイダル計算用微分 ! ! wt_KxRGrad_wt :: スペクトルデータに経度微分 ! k×r・▽ = ∂/∂λを作用させる ! xyz_KGrad_wt :: スペクトルデータに軸方向微分 ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r を ! 作用させる ! wt_L2_wt :: スペクトルデータに ! L2 演算子 = -水平ラプラシアンを ! 作用させる ! wt_L2Inv_wt :: スペクトルデータに ! L2 演算子の逆 = -逆水平ラプラシアンを ! 作用させる ! wt_QOperator_wt :: スペクトルデータに演算子 ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) を ! 作用させる ! wt_RadRot_xyz_xyz :: ベクトル v の渦度と動径ベクトル r の内積 ! r・(▽×v) を計算する ! wt_RadRotRot_xyz_xyz_xyz :: ベクトルの v の r・(▽×▽×v) を計算する ! wt_Potential2Vector :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! wt_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! ! wq_KxRGrad_wq :: スペクトルデータに経度微分 ! k×r・▽ = ∂/∂λを作用させる ! xyr_KGrad_wq :: スペクトルデータに軸方向微分 ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r を作用させる ! wq_L2_wq :: スペクトルデータに L2 演算子 = -水平ラプラシアンを ! 作用させる ! wq_L2Inv_wq :: スペクトルデータに L2 演算子の逆 = -逆水平ラプラシアン ! を作用させる ! wq_QOperator_wq :: スペクトルデータに演算子 ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) を作用させる ! wq_RadRot_xyr_xyr :: ベクトル v の渦度と動径ベクトル r の内積 r・(▽×v) を ! 計算する ! wr_RadRotRot_xyr_xyr_xyr :: ベクトルの v の r・(▽×▽×v) を計算する ! wq_RadRotRot_xyr_xyr_xyr :: ベクトルの v の r・(▽×▽×v) を計算する ! wq_Potential2Vector :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! wq_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! nmz_ToroidalEnergySpectrum_wt, nz_ToroidalEnergySpectrum_wt :: ! トロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! nmz_PoloidalEnergySpectrum_wt, nz_PoloidalEnergySpectrum_wt :: ! ポロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! nmr_ToroidalEnergySpectrum_wq, nr_ToroidalEnergySpectrum_wq :: ! トロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! nmr_PoloidalEnergySpectrum_wq, nr_PoloidalEnergySpectrum_wq :: ! ポロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! !==== 境界値問題 ! ! wt_BoundariesTau, wt_BoundariesGrid, wt_Boundaries :: ! タウ法, 選点法 ! wt_TorBoundariesTau, wt_TorBoundariesGrid, wt_TorBoundaries :: ! 速度トロイダルポテンシャルの境界条件を球殻領域に適用する ! (タウ法,選点法) ! wz_LaplaPol2Pol_wz, wt_LaplaPol2PolTau_wt, wt_LaplaPol2Pol_wt :: ! 球殻領域の速度ポロイダルポテンシャルΦを▽^2Φから求める ! (入出力がそれぞれチェビシェフ格子点,チェビシェフ係数) ! wt_TorMagBoundariesTau, wt_TorMagBoundariesGrid, wt_TorMagBoundaries :: ! 磁場トロイダルポテンシャルの境界条件を球殻領域に適用する ! (タウ法, 選点法) ! wt_PolMagBoundariesTau, wt_PolMagBoundariesGrid, wt_PolMagBoundaries :: ! 磁場トロイダルポテンシャル境界の境界条件を球殻領域に適用する ! (タウ法, 選点法) ! wq_BoundaryTau, wr_BoundaryGrid, wq_Boundary :: ! ディリクレ, ノイマン境界条件を適用する(タウ法, 選点法) ! wq_TorBoundaryTau, wr_TorBoundaryGrid, wq_TorBoundary :: ! 速度トロイダルポテンシャルの境界条件を適用する(タウ法,選点法) ! wr_LaplaPol2Pol_wr, wq_LaplaPol2Pol_wq :: ! 速度ポロイダルポテンシャルΦを▽^2Φから求める ! (入出力がそれぞれ格子点およびスペクトル係数) ! wq_TorMagBoundaryTau, wr_TorMagBoundaryGrid, wq_TorMagBoundary :: ! 磁場トロイダルポテンシャルの境界条件を適用する(タウ法, 選点法) ! wq_PolMagBoundaryTau, wr_PolMagBoundaryGrid, wq_PolMagBoundary :: ! 磁場トロイダルポテンシャル境界の境界条件を適用する(タウ法, 選点法) ! wtq_TormagBoundariesTau, wtq_TormagBoundariesGrid :: ! 磁場トロイダルポテンシャル境界の境界条件を球および球殻領域全体に ! 適用する(タウ法, 選点法) ! !==== 積分・平均(3 次元データ) ! ! IntLonLatRad_xyz, AvrLonLatRad_xyz :: 3 次元格子点データの ! 全領域積分および平均 ! z_IntLonLat_xyz, z_AvrLonLat_xyz :: 3 次元格子点データの ! 緯度経度(水平・球面)積分および平均 ! y_IntLonRad_xyz, y_AvrLonRad_xyz :: 3 次元格子点データの ! 緯度動径積分および平均 ! z_IntLatRad_xyz, z_AvrLatRad_xyz :: 3 次元格子点データの ! 経度動径(子午面)積分および平均 ! yz_IntLon_xyz, yz_AvrLon_xyz :: 3 次元格子点データの ! 経度方向積分および平均 ! xz_IntLat_xyz, xz_AvrLat_xyz :: 3 次元格子点データの ! 緯度方向積分および平均 ! xz_IntRad_xyz, xz_AvrRad_xyz :: 3 次元格子点データの ! 動径方向積分および平均 ! IntLonLatRad_xyr, AvrLonLatRad_xyr :: 3 次元格子点データの ! 全領域積分および平均 ! z_IntLonLat_xyr, z_AvrLonLat_xyr :: 3 次元格子点データの ! 緯度経度(水平・球面)積分および平均 ! y_IntLonRad_xyr, y_AvrLonRad_xyr :: 3 次元格子点データの ! 緯度動径積分および平均 ! z_IntLatRad_xyr, z_AvrLatRad_xyr :: 3 次元格子点データの ! 経度動径(子午面)積分および平均 ! yz_IntLon_xyr, yz_AvrLon_xyr :: 3 次元格子点データの ! 経度方向積分および平均 ! xz_IntLat_xyr, xz_AvrLat_xyr :: 3 次元格子点データの ! 緯度方向積分および平均 ! xz_IntRad_xyr, xz_AvrRad_xyr :: 3 次元格子点データの ! 動径方向積分および平均 ! !==== 積分・平均(2 次元データ) ! ! IntLonLat_xy, AvrLonLat_xy :: 2 次元格子点データの水平(球面)積分および平均 ! ! IntLonRad_xz, AvrLonRad_xz :: 2 次元(XZ)格子点データの経度動径積分 ! および平均 ! IntLatRad_yz, AvrLatRad_yz :: 2 次元(YZ)格子点データの緯度動径(子午面) ! 積分および平均 ! y_IntLon_xy, y_AvrLon_xy :: 水平 2 次元(球面)格子点データの経度方向 ! 積分および平均 ! x_IntLat_xy, x_AvrLat_xy :: 水平2 次元(球面)格子点データの緯度方向積分 ! および平均 ! z_IntLon_xz, z_AvrLon_xz :: 2 次元(XZ)格子点データの経度方向積分および ! 平均 ! x_IntRad_xz, x_AvrRad_xz :: 2 次元(XZ)格子点データの動径方向積分および ! 平均 ! z_IntLat_yz, z_AvrLat_yz :: 2 次元(YZ)格子点データの緯度方向積分および ! 平均 ! y_IntRad_yz, y_AvrRad_yz :: 2 次元(YZ)格子点データの動径方向積分および ! 平均 ! IntLonRad_xr, AvrLonRad_xr :: 2 次元(XZ)格子点データの経度動径積分 ! および平均 ! IntLatRad_yr, AvrLatRad_yr :: 2 次元(YZ)格子点データの緯度動径(子午面) ! 積分および平均 ! r_IntLon_xr, r_AvrLon_xr :: 2 次元(XZ)格子点データの経度方向積分および ! 平均 ! x_IntRad_xr, x_AvrRad_xr :: 2 次元(XZ)格子点データの動径方向積分および ! 平均 ! r_IntLat_yr, r_AvrLat_yr :: 2 次元(YZ)格子点データの緯度方向積分および ! 平均 ! y_IntRad_yr, y_AvrRad_yr :: 2 次元(YZ)格子点データの動径方向積分および ! 平均 ! !==== 積分・平均(1 次元データ) ! ! IntLon_x, AvrLon_x :: 1 次元(X)格子点データの経度方向積分および平均 ! IntLat_y, AvrLat_y :: 1 次元(Y)格子点データの緯度方向積分および平均 ! IntRad_z, AvrRad_z :: 1 次元(Z)格子点データの動径方向積分および平均 ! IntRad_r, AvrRad_r :: 1 次元(Z)格子点データの動径方向積分および平均 ! ! use dc_message use lumatrix use wt_module_sjpack use wq_module_sjpack, only : wq_initial, & r_Rad, r_Rad_Weight, xyr_Lon, xyr_Lat, xyr_Rad, wr_Rad, wq_VMiss, & wq_RadDRad_wq, q_RadDRad_q, & wq_Rad2_wq, q_Rad2_q, wq_Rad2Inv_wq, q_Rad2Inv_q, & wr_wq, wq_wr, xyr_wq, wq_xyr, xyr_wr, wr_xyr, & wr_DivRad_wq, wr_DivRad_xyr, wr_RotDRad_wq, wr_RotDRad_wr, wr_RotDRad_wr, & wq_RotDRad_wr, wq_Lapla_wq, & xyr_GradLon_wq, xyr_gradlat_wq, wr_DivLon_xyr, wr_DivLat_xyr, & wr_Div_xyr_xyr_xyr, xyr_Div_xyr_xyr_xyr, & xyr_RotLon_wq_wq, xyr_RotLat_wq_wq, wr_RotRad_xyr_xyr, & yr_IntLon_xyr, xr_IntLat_xyr, xy_IntRad_xyr, & x_IntLatRad_xyr, y_IntLonRad_xyr, r_IntLonLat_xyr, & IntLonLatRad_xyr, & r_IntLat_yr, y_IntRad_yr, IntLatRad_yr, & r_IntLon_xr, x_IntRad_xr, IntLonRad_xr, & IntRad_r, & yr_AvrLon_xyr, xr_AvrLat_xyr, xy_AvrRad_xyr, & x_AvrLatRad_xyr, y_AvrLonRad_xyr, r_AvrLonLat_xyr, & AvrLonLatRad_xyr, & r_AvrLat_yr, y_AvrRad_yr, AvrLatRad_yr, & r_AvrLon_xr, x_AvrRad_xr, AvrLonRad_xr, & AvrLon_x, AvrLat_y, AvrRad_r, & wq_KxRGrad_wq, xyr_KGrad_wq, wq_L2_wq, wq_L2Inv_wq, wq_QOperator_wq, & wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr, & wq_RadRotRot_xyr_xyr_xyr, & wq_Potential2vector, wq_Potential2Rotation, & nmr_ToroidalEnergySpectrum_wq, nr_ToroidalEnergySpectrum_wq,& nmr_PoloidalEnergySpectrum_wq, nr_PoloidalEnergySpectrum_wq, & wq_BoundaryTau, wq_TorBoundaryTau, wq_LaplaPol2PolTau_wq, & wq_TormagBoundaryTau, wq_PolmagBoundaryTau, & wr_BoundaryGrid, wr_TorBoundaryGrid, & wr_TormagBoundaryGrid, wr_PolmagBoundaryGrid implicit none private public wtq_Initial public x_Lon, x_Lon_Weight public y_Lat, y_Lat_Weight public z_Rad, z_Rad_Weight public l_nm, nm_l public xy_Lon, xy_Lat public xyz_Lon, xyz_Lat, xyz_Rad public wz_Rad public wt_VMiss public w_xy, xy_w public at_Dr_at, t_Dr_t, az_at, at_az public xyz_wt, wt_xyz, xyz_wz, wz_xyz, wz_wt, wt_wz public wt_DRad_wt, wt_DivRad_wt, wt_RotRad_wt, wt_Lapla_wt public xyz_GradLon_wt, xyz_gradlat_wt public wt_DivLon_xyz, wt_DivLat_xyz public wt_Div_xyz_xyz_xyz, xyz_Div_xyz_xyz_xyz public xyz_RotLon_wt_wt, xyz_RotLat_wt_wt, wt_RotRad_xyz_xyz public yz_IntLon_xyz, xz_IntLat_xyz, xy_IntRad_xyz public x_IntLatRad_xyz, y_IntLonRad_xyz, z_IntLonLat_xyz public IntLonLatRad_xyz public x_IntLat_xy, y_IntLon_xy, IntLonLat_xy public z_IntLat_yz, y_IntRad_yz, IntLatRad_yz public z_IntLon_xz, x_IntRad_xz, IntLonRad_xz public IntLon_x, IntLat_y, IntRad_z public yz_AvrLon_xyz, xz_AvrLat_xyz, xy_AvrRad_xyz public x_AvrLatRad_xyz, y_AvrLonRad_xyz, z_AvrLonLat_xyz public AvrLonLatRad_xyz public x_AvrLat_xy, y_AvrLon_xy, AvrLonLat_xy public z_AvrLat_yz, y_AvrRad_yz, AvrLatRad_yz public z_AvrLon_xz, x_AvrRad_xz, AvrLonRad_xz public AvrLon_x, AvrLat_y, AvrRad_z public wt_KxRGrad_wt, xyz_KGrad_wt, wt_L2_wt, wt_L2Inv_wt, wt_QOperator_wt public wt_RadRot_xyz_xyz, wt_RadRotRot_xyz_xyz_xyz public wt_Potential2vector, wt_Potential2Rotation public wt_VGradV public Interpolate_wt public nmz_ToroidalEnergySpectrum_wt, nz_ToroidalEnergySpectrum_wt public nmz_PoloidalEnergySpectrum_wt, nz_PoloidalEnergySpectrum_wt public wt_Boundaries, wt_TorBoundaries public wt_LaplaPol2Pol_wt, wz_LaplaPol2Pol_wz public wt_TormagBoundaries, wt_PolmagBoundaries public wt_BoundariesTau, wt_TorBoundariesTau, wt_LaplaPol2PolTau_wt public wt_TormagBoundariesTau, wt_PolmagBoundariesTau public wt_BoundariesGrid, wt_TorBoundariesGrid, wt_LaplaPol2PolGrid_wt public wt_TormagBoundariesGrid, wt_PolmagBoundariesGrid public r_Rad, r_Rad_Weight public xyr_Lon, xyr_Lat, xyr_Rad public wr_Rad public wq_VMiss public wq_RadDRad_wq, q_RadDRad_q, wr_wq, wq_wr public wq_Rad2_wq, q_Rad2_q, wq_Rad2Inv_wq, q_Rad2Inv_q public xyr_wq, wq_xyr, xyr_wr, wr_xyr public wr_DivRad_wq, wr_DivRad_xyr, wr_RotDRad_wq, wr_RotDRad_wr, wq_Lapla_wq public wq_RotDRad_wr public xyr_GradLon_wq, xyr_GradLat_wq public wr_DivLon_xyr, wr_DivLat_xyr public wr_Div_xyr_xyr_xyr, xyr_Div_xyr_xyr_xyr public xyr_RotLon_wq_wq, xyr_RotLat_wq_wq, wr_RotRad_xyr_xyr public yr_IntLon_xyr, xr_IntLat_xyr, xy_IntRad_xyr public x_IntLatRad_xyr, y_IntLonRad_xyr, r_IntLonLat_xyr public IntLonLatRad_xyr public r_IntLat_yr, y_IntRad_yr, IntLatRad_yr public r_IntLon_xr, x_IntRad_xr, IntLonRad_xr public IntRad_r public yr_AvrLon_xyr, xr_AvrLat_xyr, xy_AvrRad_xyr public x_AvrLatRad_xyr, y_AvrLonRad_xyr, r_AvrLonLat_xyr public AvrLonLatRad_xyr public r_AvrLat_yr, y_AvrRad_yr, AvrLatRad_yr public r_AvrLon_xr, x_AvrRad_xr, AvrLonRad_xr public AvrRad_r public wq_KxRGrad_wq, xyr_KGrad_wq, wq_L2_wq, wq_L2Inv_wq, wq_QOperator_wq public wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr public wq_RadRotRot_xyr_xyr_xyr public wq_Potential2vector, wq_Potential2Rotation public nmr_ToroidalEnergySpectrum_wq, nr_ToroidalEnergySpectrum_wq public nmr_PoloidalEnergySpectrum_wq, nr_PoloidalEnergySpectrum_wq public wq_Boundary, wq_TorBoundary, wq_LaplaPol2Pol_wq ! wr_LaplaPol2Pol_wr public wq_TormagBoundary, wq_PolmagBoundary public wq_BoundaryTau, wq_TorBoundaryTau, wq_LaplaPol2PolTau_wq public wq_TormagBoundaryTau, wq_PolmagBoundaryTau public wr_BoundaryGrid, wr_TorBoundaryGrid public wr_TormagBoundaryGrid, wr_PolmagBoundaryGrid public wtq_TormagBoundaries, wtq_TormagBoundariesTau, wtq_TormagBoundariesGrid public wtq_PolmagBoundaries, wtq_PolmagBoundariesTau, wtq_PolmagBoundariesGrid interface wq_Boundary module procedure wq_BoundaryTau end interface interface wq_TorBoundary module procedure wq_TorBoundaryTau end interface interface wq_LaplaPol2Pol_wq module procedure wq_LaplaPol2PolTau_wq end interface interface wq_TorMagBoundary module procedure wq_TorMagBoundaryTau end interface interface wq_PolMagBoundary module procedure wq_PolMagBoundaryTau end interface interface wtq_TorMagBoundaries module procedure wtq_TorMagBoundariesTau end interface interface wtq_PolMagBoundaries module procedure wtq_PolMagBoundariesTau end interface integer :: im=64, jm=32 ! 水平格子点の設定(経度, 緯度) integer :: nm=21 ! 水平切断波数の設定 integer :: kmo=16, kmi=16 ! 鉛直格子点の設定(球殻, 球) integer :: lmo=16, lmi=31 ! 鉛直切断波数の設定(球殻, 球) integer, dimension(:), allocatable :: nd ! 重み r^n の指数 real(8) :: ri=1.0D0, ro=2.0D0 ! 内外球半径 real(8), parameter :: pi=3.1415926535897932385D0 save im, jm, kmo, kmi, nm, lmo, lmi, ri, ro contains !--------------- 初期化 ----------------- subroutine wtq_Initial(i,j,ki,ko,n,li,lo,r_in,r_out,np) ! ! スペクトル変換の格子点数, 波数, 球および球殻動径座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! ! np に 1 より大きな値を指定すれば ISPACK の球面調和函数変換 ! OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには, ! 実行時に環境変数 OMP_NUM_THREADS を np 以下の数字に設定する等の ! システムに応じた準備が必要となる. ! ! np に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない. ! ! integer,intent(in) :: i ! 格子点数(経度λ) integer,intent(in) :: j ! 格子点数(緯度φ) integer,intent(in) :: ki ! 格子点数(動径 r) integer,intent(in) :: ko ! 格子点数(動径 z) integer,intent(in) :: n ! 切断波数(水平全波数) integer,intent(in) :: li ! 切断波数(球動径波数) integer,intent(in) :: lo ! 切断波数(球殻動径波数) real(8),intent(in) :: r_in, r_out ! 球半径 integer,intent(in), optional :: np ! OPENMP での最大スレッド数 integer :: nn, mm im = i ; jm = j ; kmo = ko ; kmi = ki nm = n ; lmo = lo ; lmi=li ri = r_in ; ro = r_out if ( present(np) ) then if ( kmo+1 > kmi ) then call wt_Initial(im,jm,kmo,nm,lmo,ri,ro,np) call wq_Initial(im,jm,kmi,nm,lmi,ri,np,wa_init=.false.) else call wq_Initial(im,jm,kmi,nm,lmi,ri,np) call wt_Initial(im,jm,kmo,nm,lmo,ri,ro,np,wa_init=.false.) endif else if ( kmo+1 > kmi ) then call wt_Initial(im,jm,kmo,nm,lmo,ri,ro) call wq_Initial(im,jm,kmi,nm,lmi,ri,wa_init=.false.) else call wq_Initial(im,jm,kmi,nm,lmi,ri) call wt_Initial(im,jm,kmo,nm,lmo,ri,ro,wa_init=.false.) endif endif allocate(nd((nm+1)*(nm+1))) do nn=0,nm do mm=-nn,nn nd(l_nm(nn,mm)) = nn ! 水平 n 次は r^n から enddo enddo call MessageNotify('M','wtq_initial','wtq_module_sjpack (2012/08/24) is initialized') end subroutine wtq_Initial !--------------- 境界値問題 ----------------- subroutine wtq_TormagBoundariesTau(wt_TOR,wq_TOR,Pmo,Pmi,new) ! ! 磁場トロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). 現在のところ外側境界物質が非電気伝導体の場合のみ ! 対応している. その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 外側 ! Psi_o = 0 at the outer boundary ! 球--球殻境界 ! Psi_o = Psi_i, Pm_o DrDPsi_o = Pm_i DrDPsi_i at the boundary ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension((nm+1)*(nm+1),0:lmo),intent(inout) :: wt_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension((nm+1)*(nm+1),0:lmi),intent(inout) :: wq_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8),intent(in) :: Pmo !(in) 球殻の磁気プランドル数 real(8),intent(in) :: Pmi !(in) 内球の磁気プランドル数 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: wt_I real(8), dimension(:,:), allocatable :: wz_PSI real(8), dimension(:,:), allocatable :: wz_DPSIDR real(8), dimension(:,:), allocatable :: wq_I real(8), dimension(:,:), allocatable :: wr_PSI real(8), dimension(:,:), allocatable :: wr_DPSIDR real(8), dimension((nm+1)*(nm+1),0:lmo+lmi+1) :: wtq_PSI logical :: first = .true. logical :: new_matrix = .false. integer :: n, l, lend save :: alu, kp, first if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) if ( allocated(wt_I) ) deallocate(wt_I) if ( allocated(wz_PSI) ) deallocate(wz_PSI) if ( allocated(wq_I) ) deallocate(wq_I) if ( allocated(wr_PSI) ) deallocate(wr_PSI) allocate(alu((nm+1)*(nm+1),0:lmo+lmi+1,0:lmo+lmi+1)) allocate(kp((nm+1)*(nm+1),0:lmo+lmi+1)) allocate(wt_I((nm+1)*(nm+1),0:lmo),wz_PSI((nm+1)*(nm+1),0:kmo)) allocate(wq_I((nm+1)*(nm+1),0:lmi),wr_PSI((nm+1)*(nm+1),kmi)) allocate(wz_DPSIDR((nm+1)*(nm+1),0:kmo)) allocate(wr_DPSIDR((nm+1)*(nm+1),kmi)) alu = 0.0D0 do l=0,lmo+lmi+1 wt_I = 0.0D0 alu(:,l,l) = 1.0D0 ! 低波数はそのまま enddo do l=0,lmo wt_I = 0.0D0 wt_I(:,l) = 1.0D0 wz_PSI = wz_wt(wt_I) wz_DPSIDR = wz_wt(wt_DRad_wt(wt_I)) alu(:,lmo-1,l) = wz_PSI(:,0) ! 外側境界 alu(:,lmo,l) = Pmo * wz_DPSIDR(:,kmo) ! 球--球殻境界 do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then lend = lmi else lend = lmi-1 endif alu(n,lmo+lend+1,l) = wz_PSI(n,kmo) ! 球--球殻境界 enddo enddo do l=0,lmi wq_I = 0.0D0 wq_I(:,l) = 1.0D0 wr_PSI = wr_wq(wq_I) wr_DPSIDR = wr_wq(wq_RadDRad_wq(wq_I))/wr_Rad do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then lend = lmi else lend = lmi-1 endif alu(n,lmo,lmo+1+l) = -Pmi * wr_DPSIDR(n,kmi) ! 球--球殻境界 alu(n,lmo+lend+1,lmo+1+l) = - wr_PSI(n,kmi) ! 球--球殻境界 enddo enddo ! 関係ないところを 0 で埋める. do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then lend = lmi else lend = lmi-1 endif do l=0,nd(n)-1 alu(n,lmo+lend,l) = 0.0D0 enddo do l=nd(n)+1,lmi,2 alu(n,lmo+lend,l) = 0.0D0 enddo enddo call ludecomp(alu,kp) deallocate(wt_I,wz_PSI,wz_DPSIDR,wq_I,wr_PSI,wr_DPSIDR) call MessageNotify('M','wtq_TormagBoundaryTau',& 'Matrix to apply b.c. newly produced.') endif wtq_PSI(:,0:lmo) = wt_TOR wtq_PSI(:,lmo+1:lmo+lmi+1) = wq_TOR wtq_PSI(:,lmo-1) = 0.0D0 ! 境界条件に対応する行を 0 にする wtq_PSI(:,lmo) = 0.0D0 ! 境界条件に対応する行を 0 にする do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then wtq_PSI(n,lmo+lmi+1) = 0.0D0 ! 境界条件に対応する行を 0 にする else wtq_PSI(n,lmo+lmi) = 0.0D0 ! 境界条件に対応する行を 0 にする endif enddo wtq_PSI = lusolve(alu,kp,wtq_PSI) wt_TOR = wtq_PSI(:,0:lmo) wq_TOR = wtq_PSI(:,lmo+1:lmo+lmi+1) end subroutine wtq_TormagBoundariesTau subroutine wtq_TormagBoundariesGrid(wz_TOR,wr_TOR,Pmo,Pmi,new) ! ! 磁場トロイダルポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). ! ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 外側 ! Psi_o = 0 at the outer boundary ! 球--球殻境界 ! Psi_o = Psi_i, Pm_o DrDPsi_o = Pm_i DrDPsi_i at the boundary ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension((nm+1)*(nm+1),0:kmo),intent(inout) :: wz_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension((nm+1)*(nm+1),kmi),intent(inout) :: wr_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8),intent(in) :: Pmo !(in) 球殻の磁気プランドル数 real(8),intent(in) :: Pmi !(in) 内球の磁気プランドル数 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: wz_I real(8), dimension(:,:), allocatable :: wz_PSI real(8), dimension(:,:), allocatable :: wz_DPSIDR real(8), dimension(:,:), allocatable :: wr_I real(8), dimension(:,:), allocatable :: wr_PSI real(8), dimension(:,:), allocatable :: wr_DPSIDR real(8), dimension((nm+1)*(nm+1),0:kmo+kmi) :: wzr_PSI logical :: first = .true. logical :: new_matrix = .false. integer :: k save :: alu, kp, first if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) if ( allocated(wz_I) ) deallocate(wz_I) if ( allocated(wz_PSI) ) deallocate(wz_PSI) if ( allocated(wr_I) ) deallocate(wr_I) if ( allocated(wr_PSI) ) deallocate(wr_PSI) allocate(alu((nm+1)*(nm+1),0:kmo+kmi,0:kmo+kmi)) allocate(kp((nm+1)*(nm+1),0:kmo+kmi)) allocate(wz_I((nm+1)*(nm+1),0:kmo),wz_PSI((nm+1)*(nm+1),0:kmo)) allocate(wr_I((nm+1)*(nm+1),kmi),wr_PSI((nm+1)*(nm+1),kmi)) allocate(wz_DPSIDR((nm+1)*(nm+1),0:kmo)) allocate(wr_DPSIDR((nm+1)*(nm+1),kmi)) alu = 0.0D0 do k=0,kmo+kmi alu(:,k,k) = 1.0D0 ! 内部領域は値そのまま. enddo ! 境界条件 do k=0,kmo wz_I = 0.0D0 wz_I(:,k) = 1.0D0 wz_DPSIDR = wz_wt(wt_DRad_wt(wt_wz(wz_I))) alu(:,0,k) = wz_I(:,0) ! 外側境界 alu(:,kmo,k) = Pmo * wz_DPSIDR(:,kmo) ! 球--球殻境界 alu(:,kmo+kmi,k) = wz_I(:,kmo) ! 球--球殻境界 enddo do k=1,kmi wr_I = 0.0D0 wr_I(:,k) = 1.0D0 wr_DPSIDR = wr_wq(wq_RadDRad_wq(wq_wr(wr_I)))/wr_Rad alu(:,kmo,kmo+k) = -Pmi * wr_DPSIDR(:,kmi) alu(:,kmo+kmi,kmo+k) = - wr_I(:,kmi) enddo call ludecomp(alu,kp) deallocate(wz_I,wz_PSI,wz_DPSIDR,wr_I,wr_PSI,wr_DPSIDR) call MessageNotify('M','wtq_TormagBoundaryGrid',& 'Matrix to apply b.c. newly produced.') endif wzr_PSI(:,0:kmo) = wz_TOR wzr_PSI(:,kmo+1:kmo+kmi) = wr_TOR wzr_PSI(:,0) = 0.0D0 wzr_PSI(:,kmo) = 0.0D0 wzr_PSI(:,kmo+kmi) = 0.0D0 wzr_PSI = lusolve(alu,kp,wzr_PSI) wz_TOR = wzr_PSI(:,0:kmo) wr_TOR = wzr_PSI(:,kmo+1:kmo+kmi) end subroutine wtq_TormagBoundariesGrid subroutine wtq_PolmagBoundariesTau(wt_Pol,wq_Pol,new) ! ! 磁場ポロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). 現在のところ外側境界物質が非電気伝導体の場合のみ ! 対応している. その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 現在のところ外側境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 外側 ! dPol_o/dr +(n+1)/r Pol_o= 0 at the outer boundary ! 球--球殻境界 ! Pol_o = Phi_i, DrDPol_o = DrDPol_i at the boundary ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension((nm+1)*(nm+1),0:lmo),intent(inout) :: wt_Pol !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension((nm+1)*(nm+1),0:lmi),intent(inout) :: wq_Pol !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: wt_I real(8), dimension(:,:), allocatable :: wz_PSI real(8), dimension(:,:), allocatable :: wz_DPSIDR real(8), dimension(:,:), allocatable :: wq_I real(8), dimension(:,:), allocatable :: wr_PSI real(8), dimension(:,:), allocatable :: wr_DPSIDR real(8), dimension((nm+1)*(nm+1),0:lmo+lmi+1) :: wtq_Pol logical :: first = .true. logical :: new_matrix = .false. integer :: n,l,nn(2), lend save :: alu, kp, first if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) if ( allocated(wt_I) ) deallocate(wt_I) if ( allocated(wz_PSI) ) deallocate(wz_PSI) if ( allocated(wq_I) ) deallocate(wq_I) if ( allocated(wr_PSI) ) deallocate(wr_PSI) allocate(alu((nm+1)*(nm+1),0:lmo+lmi+1,0:lmo+lmi+1)) allocate(kp((nm+1)*(nm+1),0:lmo+lmi+1)) allocate(wt_I((nm+1)*(nm+1),0:lmo),wz_PSI((nm+1)*(nm+1),0:kmo)) allocate(wq_I((nm+1)*(nm+1),0:lmi),wr_PSI((nm+1)*(nm+1),kmi)) allocate(wz_DPSIDR((nm+1)*(nm+1),0:kmo)) allocate(wr_DPSIDR((nm+1)*(nm+1),kmi)) alu = 0.0D0 do l=0,lmo+lmi+1 alu(:,l,l) = 1.0d0 ! 低波数はそのまま. enddo ! 球殻外側境界 do l=0,lmo wt_I = 0.0D0 wt_I(:,l) = 1.0D0 wz_PSI = wz_wt(wt_I) wz_DPSIDR = wz_wt(wt_DRad_wt(wt_I)) do n=1,(nm+1)*(nm+1) nn=nm_l(n) alu(n,lmo-1,l) = wz_DPSIDR(n,0) + (nn(1)+1)*wz_PSI(n,0)/z_RAD(0) enddo enddo ! 球--球殻境界 do l=0,lmo wt_I = 0.0D0 wt_I(:,l) = 1.0D0 wz_PSI = wz_wt(wt_I) wz_DPSIDR = wz_wt(wt_DRad_wt(wt_I)) alu(:,lmo,l) = wz_DPSIDR(:,kmo) ! 球--球殻境界 do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then lend = lmi else lend = lmi-1 endif alu(n,lmo+lend+1,l) = wz_PSI(n,kmo) ! 球--球殻境界 enddo enddo do l=0,lmi wq_I = 0.0D0 wq_I(:,l) = 1.0D0 wr_PSI = wr_wq(wq_I) wr_DPSIDR = wr_wq(wq_RadDRad_wq(wq_I))/wr_Rad do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then lend = lmi else lend = lmi-1 endif alu(n,lmo,lmo+1+l) = - wr_DPSIDR(n,kmi) ! 球--球殻境界 alu(n,lmo+lend+1,lmo+1+l) = - wr_PSI(n,kmi) ! 球--球殻境界 enddo enddo call ludecomp(alu,kp) deallocate(wt_I,wz_PSI,wz_DPSIDR,wq_I,wr_PSI,wr_DPSIDR) call MessageNotify('M','PolmagBoundaryTau',& 'Matrix to apply b.c. newly produced.') endif wtq_Pol(:,0:lmo) = wt_Pol wtq_Pol(:,lmo+1:lmo+lmi+1) = wq_Pol wtq_Pol(:,lmo-1) = 0.0D0 wtq_Pol(:,lmo) = 0.0D0 do n=1,(nm+1)**2 if ( mod(nd(n),2) .eq. mod(lmi,2) ) then wtq_Pol(n,lmo+lmi+1) = 0.0D0 ! 境界条件に対応する行を 0 にする else wtq_Pol(n,lmo+lmi) = 0.0D0 ! 境界条件に対応する行を 0 にする endif enddo wtq_Pol = lusolve(alu,kp,wtq_Pol) wt_Pol = wtq_Pol(:,0:lmo) wq_Pol = wtq_Pol(:,lmo+1:lmo+lmi+1) end subroutine wtq_PolmagBoundariesTau subroutine wtq_PolmagBoundariesGrid(wz_Pol,wr_Pol,new) ! ! 磁場ポロイダルポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). ! ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 外側 ! dPol_o/dr +(n+1)/r Pol_o= 0 at the outer boundary ! 球--球殻境界 ! Pol_o = Phi_i, DrDPol_o = DrDPol_i at the boundary ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension((nm+1)*(nm+1),0:kmo),intent(inout) :: wz_Pol !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension((nm+1)*(nm+1),kmi),intent(inout) :: wr_Pol !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: wz_I real(8), dimension(:,:), allocatable :: wz_DPSIDR real(8), dimension(:,:), allocatable :: wr_I real(8), dimension(:,:), allocatable :: wr_DPSIDR real(8), dimension((nm+1)*(nm+1),0:kmo+kmi) :: wzr_Pol logical :: first = .true. logical :: new_matrix = .false. integer :: n,k,nn(2) save :: alu, kp, first if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) if ( allocated(wz_I) ) deallocate(wz_I) if ( allocated(wr_I) ) deallocate(wr_I) allocate(alu((nm+1)*(nm+1),0:kmo+kmi,0:kmo+kmi)) allocate(kp((nm+1)*(nm+1),0:kmo+kmi)) allocate(wz_I((nm+1)*(nm+1),0:kmo)) allocate(wr_I((nm+1)*(nm+1),kmi)) allocate(wz_DPSIDR((nm+1)*(nm+1),0:kmo)) allocate(wr_DPSIDR((nm+1)*(nm+1),kmi)) alu = 0.0D0 alu = 0.0D0 do k=0,kmo+kmi alu(:,k,k) = 1.0D0 ! 内部領域は値そのまま. enddo ! 球殻外側境界 do k=0,kmo wz_I = 0.0D0 wz_I(:,k) = 1.0D0 wz_DPSIDR = wz_wt(wt_DRad_wt(wt_wz(wz_I))) do n=1,(nm+1)*(nm+1) nn=nm_l(n) alu(n,0,k) = wz_DPSIDR(n,0) + (nn(1)+1) * wz_I(n,0)/z_RAD(0) enddo enddo ! 球--球殻境界 do k=0,kmo wz_I = 0.0D0 wz_I(:,k) = 1.0D0 wz_DPSIDR = wz_wt(wt_DRad_wt(wt_wz(wz_I))) alu(:,kmo,k) = wz_DPSIDR(:,kmo) alu(:,kmo+kmi,k) = wz_I(:,kmo) enddo do k=1,kmi wr_I = 0.0D0 wr_I(:,k) = 1.0D0 wr_DPSIDR = wr_wq(wq_RadDRad_wq(wq_wr(wr_I)))/wr_Rad alu(:,kmo,kmo+k) = -wr_DPSIDR(:,kmi) alu(:,kmo+kmi,kmo+k) = - wr_I(:,kmi) enddo call ludecomp(alu,kp) deallocate(wz_I,wz_DPSIDR,wr_I,wr_DPSIDR) call MessageNotify('M','wtq_PolmagBoundaryGrid',& 'Matrix to apply b.c. newly produced.') endif wzr_Pol(:,0:kmo) = wz_Pol wzr_Pol(:,kmo+1:kmo+kmi) = wr_Pol wzr_Pol(:,0) = 0.0D0 wzr_Pol(:,kmo) = 0.0D0 wzr_Pol(:,kmo+kmi) = 0.0D0 wzr_Pol = lusolve(alu,kp,wzr_Pol) wz_Pol = wzr_Pol(:,0:kmo) wr_Pol = wzr_Pol(:,kmo+1:kmo+kmi) end subroutine wtq_PolmagBoundariesGrid end module wtq_module_sjpack