!-- !---------------------------------------------------------------------- ! Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- !表題 wtq_mpi_module_supack ! ! spml/wtq_mpi_module_supack モジュールは球とその外側の球殻領域での ! 流体運動をスペクトル法と MPI 並列化によってによって数値計算するための ! Fortran90 関数を提供するものである. ! ! 水平方向に球面調和函数変換および動径方向にチェビシェフ変換ならびに ! Matsushima and Marcus (1994) で提唱された多項式を用いる場合の ! スペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wt_mpi_module_supack, wq_mpi_module_supack を用いている. ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして ISPACK2 の ! Fortran77 サブルーチンを用いている. ! ! !履歴 2016/03/31 竹広真一 wtq_module_sjpack より改変 ! !凡例 ! データ種類と index ! x : 経度 y : 緯度 v : 緯度(分散格子) ! z : 動径(球殻内) r : 動径(球内) ! w : 球面調和関数スペクトル ! n : 球面調和関数スペクトル(水平全波数) ! m : 球面調和関数スペクトル(帯状波数) ! t : チェビシェフ関数スペクトル(球殻内) ! q : 動径スペクトル関数スペクトル(球内) ! a : 任意の次元 ! ! xyz : 球殻内 3 次元格子点データ ! xvr : 球殻内の 3 次元分散格子点データ ! xyr : 球内の 3 次元格子点データ ! xvr : 球内の 3 次元分散格子点データ ! xy : 水平 2 次元格子点データ ! xv : 水平 2 次元分散格子点データ ! yz : 球殻内子午面 2 次元格子点データ ! vz : 球殻内子午面 2 次元分散格子点データ ! yr : 球内子午面 2 次元格子点データ ! vr : 球内子午面 2 次元分散格子点データ ! xz : 球殻内緯度面 2 次元格子点データ ! xr : 球内緯度面 2 次元格子点データ ! ! wz : 水平スペクトル動径格子点データ(球殻内) ! wr : 水平スペクトル動径格子点データ(球内) ! wt : スペクトルデータ(球殻内) ! wq : スペクトルデータ(球内) ! !++ module wtq_mpi_module_supack ! != wtq_module_supack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/wtq_mpi_module_supack モジュールは球とその外側の球殻領域での ! 流体運動をスペクトル法とMPI 並列化によってによって数値計算するための ! Fortran90 関数を提供するものである. ! ! 水平方向に球面調和函数変換および動径方向にチェビシェフ変換ならびに ! Matsushima and Marcus (1994) で提唱された多項式を用いる場合の ! スペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wt_mpi_module_supack, wq_mpi_module_supack を用いている. ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして 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 次元格子点データ, ... ! !=== 各データの種類の説明 ! ! * xvz : 3 次元分散格子点データ(経度・緯度・動径-球殻) ! * 変数の種類と次元は real(8), dimension(0:im-1,jc,0:kmo). ! * im, kmo はそれぞれ経度, 動径座標の格子点数であり, ! * jc はこのプロセスで保有する緯度格子点数である. ! サブルーチン wtq_mpi_Initial を呼ぶと jc が設定される. ! ! * wt : スペクトルデータ-球殻 ! * 変数の種類と次元は real(8), dimension(nc,0:lmo). ! * nc はこのプロセスで保有するスペクトルデータ数である. ! サブルーチン wtq_mpi_Initial を呼ぶと nc が設定される. ! * lmo はチェビシェフ多項式の最大次数 ! であり, サブルーチン wtq_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * nmz : 水平スペクトルデータの並んだ 3 次元配列(球殻). ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,0:kmo). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wt_Initial にて ! あらかじめ設定しておく. ! ! * nz : スペクトルデータの並んだ 2 次元配列(球殻). ! * 変数の種類と次元は real(8), dimension(0:nm,0:kmo). ! * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wz : 水平スペクトル, 動径格子点データ(球殻). ! * 変数の種類と次元は real(8), dimension(nc,0:kmo). ! ! * xvr : 3 次元分散格子点データ(経度・緯度・動径-球) ! * 変数の種類と次元は real(8), dimension(0:im-1,jc,kmi). ! * im, kmi はそれぞれ経度, 動径座標の格子点数であり, ! * jc はこのプロセスで保有する緯度格子点数である. ! サブルーチン wtq_mpi_Initial を呼ぶと jc が設定される. ! ! * wq : スペクトルデータ(球) ! * 変数の種類と次元は real(8), dimension(nc,0:lmi). ! * nc はこのプロセスで保有するスペクトルデータ数である. ! サブルーチン wtq_mpi_Initial を呼ぶと nc が設定される. ! * lmi は Matsushima-Marcus 多項式の最大次数であり, ! サブルーチン wtq_mpi_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * nmr : 水平スペクトルデータの並んだ 3 次元配列(球). ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,kmi). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wtq_mpi_Initial にて ! あらかじめ設定しておく. ! ! * nr : スペクトルデータの並んだ 2 次元配列(球). ! * 変数の種類と次元は real(8), dimension(0:nm,kmi). ! * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり, ! サブルーチン wqt_mpi_Initial にてあらかじめ設定しておく. ! ! * wr : 水平スペクトル, 動径格子点データ(球). ! * 変数の種類と次元は real(8), dimension(nc,kmi). ! ! * wt_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyz_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * wz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * wq_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyr_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * wr_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! wtq_mpi_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_QOperatorMPI_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_QOperatorMPI_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_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_mpi_module_supack use wq_mpi_module_supack, only : wq_mpi_initial, & r_Rad, r_Rad_Weight, wr_Rad, wq_VMiss, & xvr_Lon, xvr_Lat, xvr_Rad, & wq_RadDRad_wq, q_RadDRad_q, & wq_Rad2_wq, q_Rad2_q, wq_Rad2Inv_wq, q_Rad2Inv_q, & wr_wq, wq_wr, xvr_wq, wq_xvr, xvr_wr, wr_xvr, & wr_DivRad_wq, wr_DivRad_xvr, & wr_RotDRad_wq, wr_RotDRad_wr, wr_RotDRad_wr, & wq_RotDRad_wr, wq_Lapla_wq, & xvr_GradLon_wq, xvr_GradLat_wq, wr_DivLon_xvr, wr_DivLat_xvr, & wr_Div_xvr_xvr_xvr, xvr_Div_xvr_xvr_xvr, & xvr_RotLon_wq_wq, xvr_RotLat_wq_wq, wr_RotRad_xvr_xvr, & vr_IntLon_xvr, xr_IntLat_xvr, xv_IntRad_xvr, & x_IntLatRad_xvr, v_IntLonRad_xvr, r_IntLonLat_xvr, & IntLonLatRad_xvr, & r_IntLat_vr, v_IntRad_vr, IntLatRad_vr, & r_IntLon_xr, x_IntRad_xr, IntLonRad_xr, & IntRad_r, & vr_AvrLon_xvr, xr_AvrLat_xvr, xv_AvrRad_xvr, & x_AvrLatRad_xvr, v_AvrLonRad_xvr, r_AvrLonLat_xvr, & AvrLonLatRad_xvr, & r_AvrLat_vr, v_AvrRad_vr, AvrLatRad_vr, & r_AvrLon_xr, x_AvrRad_xr, AvrLonRad_xr, & AvrLon_x, AvrLat_v, AvrRad_r, & wq_KxRGrad_wq, xvr_KGrad_wq, & wq_L2_wq, wq_L2Inv_wq, wq_QOperatorMPI_wq, & wr_RadRot_xvr_xvr, wr_RadRotRot_xvr_xvr_xvr, & wq_RadRotRot_xvr_xvr_xvr, & wq_Potential2vectorMPI, wq_Potential2RotationMPI, & 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_mpi_Initial public jc public nc public x_Lon, x_Lon_Weight public v_Lat, v_Lat_Weight public z_Rad, z_Rad_Weight public l_nm, nm_l public xv_Lon, xv_Lat public xvz_Lon, xvz_Lat, xvz_Rad public wz_Rad public wt_VMiss public w_xv, xv_w public at_Dr_at, t_Dr_t, az_at, at_az public wz_wt, wt_wz public xvz_wt, wt_xvz, xvz_wz, wz_xvz public wt_DRad_wt, wt_DivRad_wt, wt_RotRad_wt, wt_Lapla_wt public xvz_GradLon_wt, xvz_GradLat_wt public wt_DivLon_xvz, wt_DivLat_xvz public wt_Div_xvz_xvz_xvz, xvz_Div_xvz_xvz_xvz public xvz_RotLon_wt_wt, xvz_RotLat_wt_wt, wt_RotRad_xvz_xvz public vz_IntLon_xvz, xz_IntLat_xvz, xv_IntRad_xvz public x_IntLatRad_xvz, v_IntLonRad_xvz, z_IntLonLat_xvz public IntLonLatRad_xvz public x_IntLat_xv, v_IntLon_xv, IntLonLat_xv public z_IntLat_vz, v_IntRad_vz, IntLatRad_vz public z_IntLon_xz, x_IntRad_xz, IntLonRad_xz public IntLon_x, IntLat_v, IntRad_z public vz_AvrLon_xvz, xz_AvrLat_xvz, xv_AvrRad_xvz public x_AvrLatRad_xvz, v_AvrLonRad_xvz, z_AvrLonLat_xvz public AvrLonLatRad_xvz public x_AvrLat_xv, v_AvrLon_xv, AvrLonLat_xv public z_AvrLat_vz, v_AvrRad_vz, AvrLatRad_vz public z_AvrLon_xz, x_AvrRad_xz, AvrLonRad_xz public AvrLon_x, AvrLat_v, AvrRad_z public wt_KxRGrad_wt, xvz_KGrad_wt public wt_L2_wt, wt_L2Inv_wt, wt_QOperatorMPI_wt public wt_RadRot_xvz_xvz, wt_RadRotRot_xvz_xvz_xvz public wt_Potential2vectorMPI, wt_Potential2RotationMPI 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 xvr_Lon, xvr_Lat, xvr_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 xvr_wq, wq_xvr, xvr_wr, wr_xvr public wr_DivRad_wq, wr_DivRad_xvr public wr_RotDRad_wq, wr_RotDRad_wr, wq_Lapla_wq public wq_RotDRad_wr public xvr_GradLon_wq, xvr_GradLat_wq public wr_DivLon_xvr, wr_DivLat_xvr public wr_Div_xvr_xvr_xvr, xvr_Div_xvr_xvr_xvr public xvr_RotLon_wq_wq, xvr_RotLat_wq_wq, wr_RotRad_xvr_xvr public vr_IntLon_xvr, xr_IntLat_xvr, xv_IntRad_xvr public x_IntLatRad_xvr, v_IntLonRad_xvr, r_IntLonLat_xvr public IntLonLatRad_xvr public r_IntLat_vr, v_IntRad_vr, IntLatRad_vr public r_IntLon_xr, x_IntRad_xr, IntLonRad_xr public IntRad_r public vr_AvrLon_xvr, xr_AvrLat_xvr, xv_AvrRad_xvr public x_AvrLatRad_xvr, v_AvrLonRad_xvr, r_AvrLonLat_xvr public AvrLonLatRad_xvr public r_AvrLat_vr, v_AvrRad_vr, AvrLatRad_vr public r_AvrLon_xr, x_AvrRad_xr, AvrLonRad_xr public AvrRad_r public wq_KxRGrad_wq, xvr_KGrad_wq public wq_L2_wq, wq_L2Inv_wq, wq_QOperatorMPI_wq public wr_RadRot_xvr_xvr, wr_RadRotRot_xvr_xvr_xvr public wq_RadRotRot_xvr_xvr_xvr public wq_Potential2vectorMPI, wq_Potential2RotationMPI 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_mpi_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 :: nmary(2), nl im = i ; jm = j ; kmo = ko ; kmi = ki nm = n ; lmo = lo ; lmi=li ri = r_in ; ro = r_out if ( present(np) ) then call MessageNotify('M','wtq_mpi_Initial', & 'Optional argument "NP_IN" is dummy in wq_mpi_module_supack') endif if ( kmo+1 > kmi ) then call wt_mpi_Initial(im,jm,kmo,nm,lmo,ri,ro) call wq_mpi_Initial(im,jm,kmi,nm,lmi,ri,wa_init=.false.) else call wq_mpi_Initial(im,jm,kmi,nm,lmi,ri) call wt_mpi_Initial(im,jm,kmo,nm,lmo,ri,ro,wa_init=.false.) endif allocate(nd(nc)) do nl=1,nc nmary = nm_l(nl) if ( nmary(1) .ge. 0 ) then nd(nl) = nmary(1) ! 水平 n 次は r^n から else nd(nl) = 0 endif enddo call MessageNotify('M','wtq_mpi_initial', & 'wtq_mpi_module_supack (2016/04/02) is initialized') end subroutine wtq_mpi_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(nc,0:lmo),intent(inout) :: wt_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,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(nc,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(nc,0:lmo+lmi+1,0:lmo+lmi+1)) allocate(kp(nc,0:lmo+lmi+1)) allocate(wt_I(nc,0:lmo),wz_PSI(nc,0:kmo)) allocate(wq_I(nc,0:lmi),wr_PSI(nc,kmi)) allocate(wz_DPSIDR(nc,0:kmo)) allocate(wr_DPSIDR(nc,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,nc 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,nc 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,nc 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,nc 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(nc,0:kmo),intent(inout) :: wz_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,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(nc,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(nc,0:kmo+kmi,0:kmo+kmi)) allocate(kp(nc,0:kmo+kmi)) allocate(wz_I(nc,0:kmo),wz_PSI(nc,0:kmo)) allocate(wr_I(nc,kmi),wr_PSI(nc,kmi)) allocate(wz_DPSIDR(nc,0:kmo)) allocate(wr_DPSIDR(nc,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(nc,0:lmo),intent(inout) :: wt_Pol !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,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(nc,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(nc,0:lmo+lmi+1,0:lmo+lmi+1)) allocate(kp(nc,0:lmo+lmi+1)) allocate(wt_I(nc,0:lmo),wz_PSI(nc,0:kmo)) allocate(wq_I(nc,0:lmi),wr_PSI(nc,kmi)) allocate(wz_DPSIDR(nc,0:kmo)) allocate(wr_DPSIDR(nc,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,nc 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,nc 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,nc 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,nc 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(nc,0:kmo),intent(inout) :: wz_Pol !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,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(nc,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(nc,0:kmo+kmi,0:kmo+kmi)) allocate(kp(nc,0:kmo+kmi)) allocate(wz_I(nc,0:kmo)) allocate(wr_I(nc,kmi)) allocate(wz_DPSIDR(nc,0:kmo)) allocate(wr_DPSIDR(nc,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))) do n=1,nc 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_mpi_module_supack