!-- !---------------------------------------------------------------------- ! Copyright(c) 2009-2016 SPMDODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 tee_mpi_module ! ! spml/tee_module モジュールは平行平板間での 3 次元流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向にフーリエ数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ee_module, at_module を用いている. 最下部ではフーリエ ! およびチェビシェフ変換のエンジンとして ISPACK の Fortran77 ! サブルーチンを用いている. ! ! !履歴 2009/12/19 竹広真一 新規作成 ! 2010/03/10 佐々木洋平 threadprivate 削除(コンパイラ依存) ! 2010/03/27 竹広真一 goto 文を削除 ! 2011/12/05 竹広真一 tee_Div_zyx_zyx_zyx を public 変数に ! 2012/07/08 竹広真一 tee_LaplaPol2PolTau_tee を追加 ! 2012/08/24 竹広真一 tee_LaplaPol2Pol_tee public に追加 ! 2013/08/20 竹広真一 gnu fortran 対応 ! 2016/09/22 竹広真一 tee_module を mpi 並列化 ! !凡例 ! データ種類と index ! x : 水平(X) y : 水平(Y) v : 水平(Y,MPI 分割) z : 鉛直 ! e : フーリエ変換スペクトル ! f : フーリエ変換スペクトル(MPI 分割 ! l : フーリエ関数スペクトル(X 方向波数) ! m : フーリエ関数スペクトル(Y 方向波数) ! t : チェビシェフ関数スペクトル ! a : 任意の次元 ! ! zvx : 3 次元格子点データ(MPI 分割) ! vx : 水平 2 次元格子点データ(MPI 分割) ! zv : 鉛直 2 次元格子点データ(MPI 分割) ! zx : 鉛直 2 次元格子点データ ! ! zef : 水平スペクトル鉛直格子点データ(MPI 分割) ! tef : スペクトルデータ鉛直チェビシェフデータ(MPI 分割) ! !++ module tee_mpi_module ! != tee_mpi_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/tee_mpi_module モジュールは平行平板間での 3 次元流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向にフーリエ数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ee_mpi_module, at_module を用いている. 最下部ではフーリエ ! およびチェビシェフ変換のエンジンとして ISPACK の Fortran77 ! サブルーチンを用いている. ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (tef_, zvx_, zef_, ef_, vx_, x_, v_, z_, a_) は, ! 返す値の形を示している. ! tef_ :: スペクトルデータ(2 重フーリエ・チェビシェフ変換) ! zvx_ :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・) ! zef_ :: 水平スペクトル, 鉛直格子点データ ! e2a_ :: 1 次元化した水平スペクトル, 任意座標データ ! aef_ :: 任意座標データ, 水平スペクトルデータ ! ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 (tef_, zvx_, zef_, ef_, vx_, x_, v_, z_, a_) は, 入力変数の ! 形がスペクトルデータおよび格子点データであることを示している. ! _tef :: スペクトルデータ(2 重フーリエ・チェビシェフ変換) ! _zvx :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・) ! _zef :: 水平スペクトル, 鉛直格子点データ ! _e2a :: 1 次元化した水平スペクトル, 任意座標データ ! _aef :: 任意座標データ, 水平スペクトルデータ ! !=== 各データの種類の説明 ! ! * zvx : 3 次元格子点データ(鉛直, 水平 2 次元, MPI 分割) ! * 変数の種類と次元は real(8), dimension(0:km,jc(ip),0:im-1). ! * im, km はそれぞれ水平 X, 鉛直 Z 座標の格子点数である. ! サブルーチン tee_mpi_Initial にてあらかじめ設定しておく. ! * jc(ip) は tee_mpi_initial にて全 Y 格子点数 jm を設定後に使用できる. ! ! * tef : スペクトルデータ ! * 変数の種類と次元は real(8), dimension(0:nm,-mm:mm,2*lc(ip)). ! * mm は Y 方向最大波数, nm はチェビシェフ多項式の最大次数 ! であり, サブルーチン tee_Initial にてあらかじめ設定しておく. ! * lc(ip) は tee_mpi_initial にて全 Y 方向最大波数 lm を設定後に使用できる. ! ! * zef : 水平スペクトル, 鉛直格子点データ. ! * 変数の種類と次元は real(8), dimension(0:km,-mm:mm,2*lc(ip)). ! ! * e2a : 1 次元化した水平スペクトル, 任意座標データ ! * 変数の種類と次元は real(8), dimension((2*mm+1)*2*lc(ip),:) ! ! * aef : 任意座標データ, 水平スペクトルデータ ! * 変数の種類と次元は real(8), dimension(:,-mm:mm,2*lc(ip)). ! ! * tef_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * zvx_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * zef_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! tee_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_X, v_Y, z_Z :: 格子点座標(水平 X,Y, 鉛直 Z 座標)を ! 格納した1 次元配列 ! x_X_Weight, v_X_Weight, z_Z_Weight :: 重み座標を格納した 1 次元配列 ! zvx_X, zvx_Y, zyx_Z :: 格子点データの水平鉛直座標(X,Y,Z) ! (格子点データ型 3 次元配列) ! vx_X, vx_Y :: 格子点データの水平座標(X,Y) ! zv_Z, zv_Y :: 格子点データの鉛直水平座標(Z,Y) ! zx_Z, zx_X :: 格子点データの鉛直水平座標(Z,X) ! (格子点データ型 2 次元配列) ! !==== 基本変換 ! ! zvx_tee, tee_zvx :: スペクトルデータと 3 次元格子データの間の変換 ! (2 重フーリエ, チェビシェフ変換) ! zvx_zee, zee_zvx :: 3 次元格子データと水平スペクトル・鉛直格子データとの間 ! の変換 (2 重フーリエ変換) ! zef_tee, tee_zef :: スペクトルデータと水平スペクトル・鉛直格子データとの間 ! の変換 (チェビシェフ変換) ! ee_vx, vx_ee :: スペクトルデータと 2 次元水平格子データの間の変換 ! (2 重フーリエ変換) ! az_at, at_az :: 同時に複数個行う (チェビシェフ変換)格子データと ! チェビシェフデータの間の変換を ! e2a_aef, aef_e2a :: 水平スペクトル軸を 1 次元化し転置, 転置 ２次元化する. ! !==== 微分 ! ! tef_Dx_tef :: スペクトルデータに動径微分∂/∂x を作用させる ! tef_Dy_tef :: スペクトルデータに動径微分∂/∂y を作用させる ! tef_Dz_tef :: スペクトルデータに動径微分∂/∂z を作用させる ! ! tef_Lapla_tef :: スペクトルデータにラプラシアンを作用させる ! tef_LaplaH_tef :: スペクトルデータに水平ラプラシアンを作用させる ! tef_LaplaHInv_tef :: スペクトルデータに逆水平ラプラシアンを作用させる ! ! tef_Div_zvx_zvx_zvx :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! !==== トロイダルポロイダル計算用微分 ! ! tef_ZRot_zvx_zvx :: ベクトル v の渦度と動径ベクトル r の内積 ! z・(▽×v) を計算する ! tef_ZRotRot_zvx_zvx_zvx :: ベクトルの v の z・(▽×▽×v) を計算する ! tef_Potential2Vector :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! tef_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! zef_ToroidalEnergySpectrum_tef, zk_ToroidalEnergySpectrum_tef :: ! トロイダルポテンシャルからエネルギーのフーリエ各成分を計算する ! zef_PoloidalEnergySpectrum_tef, zk_PoloidalEnergySpectrum_tef :: ! ポロイダルポテンシャルからエネルギーのフーリエ各成分を計算する ! !==== 境界値問題 ! ! tef_BoundariesTau, tef_BoundariesGrid, tef_Boundaries :: ! ディリクレ, ノイマン境界条件を適用する(タウ法, 選点法) ! ! tef_TorBoundariesTau, tef_TorBoundariesGrid, tef_TorBoundaries :: ! 速度トロイダルポテンシャルの境界条件を適用する(タウ法,選点法) ! ! tef_LaplaPol2PolTau_tef, zee_LaplaPol2Pol_zee, tef_LaplaPol2PolGrid_tef :: ! 速度ポロイダルポテンシャルΦを▽^2Φから求める ! (入出力がそれぞれチェビシェフ格子点,チェビシェフ係数) ! ! tef_TorMagBoundariesTau, tef_TorMagBoundariesGrid, tef_TorMagBoundaries :: ! 磁場トロイダルポテンシャルの境界条件を適用する(タウ法, 選点法) ! ! tef_PolMagBoundariesTau, tef_PolMagBoundariesGrid, tef_PolMagBoundaries :: ! 磁場トロイダルポテンシャル境界の境界条件を適用する(タウ法, 選点法) ! !==== 積分・平均(3 次元データ) ! ! IntZYX_zvx, AvrZYX_zvx :: 3 次元格子点データの全領域積分および平均 ! z_IntYX_zvx, z_AvrYX_zvx :: 3 次元格子点データの水平積分および平均 ! v_IntZX_zvx, v_AvrZX_zvx :: 3 次元格子点データのZX積分および平均 ! z_IntZY_zvx, z_AvrZY_zvx :: 3 次元格子点データのZY積分および平均 ! zv_IntX_zvx, zv_AvrX_zvx :: 3 次元格子点データの水平X方向積分および平均 ! zx_IntY_zvx, zx_AvrY_zvx :: 3 次元格子点データの水平Y方向積分および平均 ! zx_IntZ_zvx, zx_AvrZ_zvx :: 3 次元格子点データの鉛直方向積分および平均 ! !==== 積分・平均(2 次元データ) ! ! IntYX_vx, AvrYX_vx :: 2 次元格子点データの水平積分および平均 ! IntZX_zx, AvrZX_zx :: 2 次元(ZX)格子点データのZX積分および平均 ! IntZY_zv, AvrZY_zv :: 2 次元(ZY)格子点データのZY積分および平均 ! y_IntX_vx, y_AvrX_vx :: 水平 2 次元格子点データのX方向積分および平均 ! x_IntY_vx, x_AvrY_vx :: 水平2 次元格子点データのY方向積分および平均 ! z_IntX_zx, z_AvrX_zx :: 2 次元(ZX)格子点データのX方向積分および平均 ! x_IntZ_zx, x_AvrZ_zx :: 2 次元(ZX)格子点データのZ方向積分および平均 ! z_IntY_zv, z_AvrY_zv :: 2 次元(ZY)格子点データのY方向積分および平均 ! y_IntZ_zv, y_AvrZ_zv :: 2 次元(ZY)格子点データのZ方向積分および平均 ! !==== 積分・平均(1 次元データ) ! ! IntX_x, AvrX_x :: 1 次元(X)格子点データのX方向積分および平均 ! IntY_v, AvrY_v :: 1 次元(Y)格子点データのY方向積分および平均 ! IntZ_z, AvrZ_z :: 1 次元(Z)格子点データのZ方向積分および平均 ! !==== 補間計算 ! ! Interpolate_tee :: スペクトルデータから任意の点の値を補間する. ! use dc_message use mpi use lumatrix use ee_mpi_module, ls=>ks, le=>ke, lc=>kc, lf_l=>kf_k use at_module, only : z_Z => g_X, z_Z_WEIGHT => g_X_WEIGHT, & at_Initial, at_az => at_ag, az_at => ag_at, & t_z => t_g, z_t => g_t, & t_Dz_t => t_Dx_t, at_Dz_at => at_Dx_at, & a_Interpolate_at, & at_Boundaries, t_Boundaries, & at_BoundariesTau_DD, at_BoundariesTau_DN, & at_BoundariesTau_ND, at_BoundariesTau_NN, & at_BoundariesGrid_DD, at_BoundariesGrid_DN, & at_BoundariesGrid_ND, at_BoundariesGrid_NN implicit none private public tee_mpi_Initial public ls, le, lc public js, je, jc public lf_l public x_X, x_X_Weight public v_Y, v_Y_Weight public z_Z, z_Z_Weight public vx_X, vx_Y, zv_Z, zv_Y, zx_X, zx_Z public zvx_X, zvx_Y, zvx_Z public zef_Z public tee_VMiss public ef_vx, vx_ef public at_Dz_at, t_Dz_t, az_at, at_az, z_t, t_z public zvx_tef, tef_zvx, zvx_zef, zef_zvx, zef_tef, tef_zef public e2z_e2t, e2t_e2z public aef_e2a, e2a_aef, ef_e2, e2_ef public tef_Dx_tef, tef_Dy_tef, tef_Dz_tef public tef_Lapla_tef, tef_LaplaH_tef, tef_LaplaHInv_tef public tef_Div_zvx_zvx_zvx public zv_IntX_zvx, zx_IntY_zvx, vx_IntZ_zvx public x_IntZY_zvx, v_IntZX_zvx, z_IntYX_zvx public IntZYX_zvx public x_IntY_vx, v_IntX_vx, IntYX_vx public z_IntY_zv, v_IntZ_zv, IntZY_zv public z_IntX_zx, x_IntZ_zx, IntZX_zx public IntX_x, IntY_v, IntZ_z public zv_AvrX_zvx, zx_AvrY_zvx, vx_AvrZ_zvx public x_AvrZY_zvx, v_AvrZX_zvx, z_AvrYX_zvx public AvrZYX_zvx public x_AvrY_vx, v_AvrX_vx, AvrYX_vx public z_AvrY_zv, v_AvrZ_zv, AvrZY_zv public z_AvrX_zx, x_AvrZ_zx, AvrZX_zx public AvrX_x, AvrY_v, AvrZ_z public tef_ZRot_zvx_zvx, tef_ZRotRot_zvx_zvx_zvx public tef_Potential2vector, tef_Potential2Rotation public Interpolate_tef public zef_ToroidalEnergySpectrum_tef ! nz_ToroidalEnergySpectrum_wt public zef_PoloidalEnergySpectrum_tef ! nz_PoloidalEnergySpectrum_wt public tef_Boundaries, zef_LaplaPol2Pol_zef public tef_TorBoundaries, tef_LaplaPol2Pol_tef public tef_TormagBoundaries, tef_PolmagBoundaries public tef_BoundariesTau, tef_LaplaPol2PolTau_tef, tef_TorBoundariesTau public tef_TormagBoundariesTau, tef_PolmagBoundariesTau public tef_BoundariesGrid, tef_LaplaPol2PolGrid_tef, tef_TorBoundariesGrid public tef_TormagBoundariesGrid, tef_PolmagBoundariesGrid public at_Boundaries, t_Boundaries interface tef_Boundaries module procedure tef_BoundariesTau end interface interface tef_TorBoundaries module procedure tef_TorBoundariesTau end interface interface tef_LaplaPol2Pol_tef module procedure tef_LaplaPol2PolTau_tef end interface interface tef_TorMagBoundaries module procedure tef_TorMagBoundariesTau end interface interface tef_PolMagBoundaries module procedure tef_PolMagBoundariesTau end interface integer :: im=32, jm=32, km=16 ! 格子点の設定(X,Y,Z) integer :: lm=10, mm=10, nm=16 ! 切断波数の設定(水平X,Y, 鉛直Z) real(8) :: xl, yl, zl ! 領域の大きさ(水平X,Y, 鉛直Z) real(8), parameter :: pi=3.1415926535897932385D0 real(8), dimension(:,:,:), allocatable :: zvx_X, zvx_Y, zvx_Z ! 座標 real(8), dimension(:,:), allocatable :: zv_Z, zv_Y ! 座標 real(8), dimension(:,:), allocatable :: zx_Z, zx_X ! 座標 real(8), dimension(:,:,:), allocatable :: zef_Z ! 座標 real(8) :: tee_VMiss = -999.0 ! 欠損値 ! integer :: np ! MPI プロセス数 integer :: ip ! MPI プロセス ID save im, jm, km, lm, mm, nm, xl, yl, zl contains !--------------- 初期化 ----------------- subroutine tee_mpi_Initial(i,j,k,l,m,n,xmin,xmax,ymin,ymax,zmin,zmax) ! ! スペクトル変換の格子点数, 波数, 各座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! integer,intent(in) :: i ! 格子点数(水平X) integer,intent(in) :: j ! 格子点数(水平Y) integer,intent(in) :: k ! 格子点数(鉛直Z) integer,intent(in) :: l ! 切断波数(水平X波数) integer,intent(in) :: m ! 切断波数(水平Y波数) integer,intent(in) :: n ! 切断波数(鉛直Z波数) real(8),intent(in) :: xmin, xmax ! X 方向領域 real(8),intent(in) :: ymin, ymax ! Y 方向領域 real(8),intent(in) :: zmin, zmax ! Z 方向領域 integer :: ierr im = i ; jm = j ; km = k lm = l ; mm = m ; nm = n xl = xmax - xmin yl = ymax - ymin zl = zmax - zmin CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR) call ee_mpi_Initial(im,jm,lm,mm,xmin,xmax,ymin,ymax) call at_Initial(km,nm,zmin,zmax) allocate(zvx_X(0:km,jc(ip),0:im-1)) allocate(zvx_Y(0:km,jc(ip),0:im-1)) allocate(zvx_Z(0:km,jc(ip),0:im-1)) allocate(zv_Z(0:km,jc(ip))) allocate(zv_Y(0:km,jc(ip))) allocate(zx_Z(0:km,0:im-1)) allocate(zx_X(0:km,0:im-1)) allocate(zef_Z(0:km,-mm:mm,2*lc(ip))) zvx_X = spread(vx_X,1,km+1) zvx_Y = spread(vx_Y,1,km+1) zvx_Z = spread(spread(z_Z,2,jc(ip)),3,im) zv_Z = spread(z_Z,2,jc(ip)) zv_Y = spread(v_Y,1,km+1) zx_Z = spread(z_Z,2,im) zx_X = spread(x_X,1,km+1) zef_Z = spread(spread(z_Z,2,2*mm+1),3,2*lc(ip)) call MessageNotify('M','tee_mpi_initial', & 'tee_mpi_module (2016/09/22) is initialized') end subroutine tee_mpi_Initial !--------------- 基本変換 ----------------- function zvx_tef(tef) ! ! スペクトルデータから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 2 重フーリエチェビシェフスペクトルデータ real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_tef !(out) 3 次元格子点データ zvx_tef = zvx_zef(zef_tef(tef)) end function zvx_tef function tef_zvx(zvx) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_zvx !(out) 2 重フーリエチェビシェフスペクトルデータ tef_zvx = tef_zef(zef_zvx(zvx)) end function tef_zvx function zvx_zef(zef) ! ! 水平スペクトル・鉛直格子点データから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(0:km,-mm:mm,2*lc(ip)), intent(in) :: zef !(in) 2 次元水平スペクトル・鉛直格子点データ real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_zef !(out) 3 次元格子点データ real(8), dimension(-mm:mm,2*lc(ip)) :: ef real(8), dimension(jc(ip),0:im-1) :: vx integer :: k do k = 0, km !zyx_zef(k,:,:) = yx_ee(zef(k,:,:)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! ef = zef(k,-mm:mm,1:2*lc(ip)) vx = vx_ef(ef) zvx_zef(k,:,:) = vx enddo end function zvx_zef function zef_zvx(zvx) ! ! 3 次元格子データから水平スペクトル・鉛直格子点データへ(正)変換する. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_zvx !(out) 2 次元スペクトル・鉛直格子点データ real(8), dimension(-mm:mm,2*lc(ip)) :: ef real(8), dimension(jc(ip),0:im-1) :: vx integer :: k do k = 0, km !zef_zyx(k,:,:) = ee_yx(zyx(k,:,:)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! vx = zvx(k,1:jc(ip),0:im-1) ef = ef_vx(vx) zef_zvx(k,:,:) = ef enddo end function zef_zvx function zef_tef(tef) ! ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 2 次元水平スペクトル鉛直チェビシェフスペクトルデータ real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_tef !(out) 2 次元水平スペクトル・鉛直格子点データ zef_tef = aef_e2a(e2z_e2t(e2a_aef(tef))) end function zef_tef function tef_zef(zef) ! ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:km,-mm:mm,2*lc(ip)), intent(in) :: zef !(in) 2 次元水平スペクトル・鉛直格子点データ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_zef !(out) 2 次元水平鉛直チェビシェフスペクトルデータ tef_zef = aef_e2a(e2t_e2z(e2a_aef(zef))) end function tef_zef function e2z_e2t(e2t) ! ! 鉛直スペクトルを格子点に変換する ! real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t !(in) 1 次元化された水平スペクトル・鉛直スペクトル real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_e2t !(out) 1 次元化された水平スペクトル・鉛直格子座標 e2z_e2t = az_at(e2t) end function e2z_e2t function e2t_e2z(e2z) ! ! 鉛直格子点をスペクトルに変換する ! real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z !(in) 1 次元化された水平スペクトル・鉛直格子座標 real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_e2z !(out) 1 次元化された水平スペクトル・鉛直スペクトル e2t_e2z = at_az(e2z) end function e2t_e2z function e2a_aef(aef) ! ! 水平スペクトルを 1 次元化転置する. ! real(8), dimension(:,-mm:,:), intent(in) :: aef !(in) 任意座標・水平 2 次元スペクトルデータ dimension(:,-mm:mm,2*lc(ip)) real(8), dimension(2*lc(ip)*(2*mm+1),size(aef,1)) :: e2a_aef !(out) 1 次元化された水平スペクトル・任意座標データ integer :: j,k,l,m if ( size(aef,2) /= 2*mm+1 ) & call MessageNotify('E','e2a_aef',& '2nd dimension of input data invalid') if ( size(aef,3) /= 2*lc(ip) ) & call MessageNotify('E','e2a_aef',& '3rd dimension of input data invalid') !e2a_aee = transpose(reshape(aee,(/size(aee,1),2*lc(ip)*(2*mm+1)/))) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do k=1,size(aef,1) do l=1,2*lc(ip) do m=-mm,mm j = (m+mm+1) + (2*mm+1)*(l-1) e2a_aef(j,k) = aef(k,m,l) enddo enddo enddo end function e2a_aef function aef_e2a(e2a) ! ! 水平スペクトルを転置展開する. ! real(8), dimension(:,:),intent(in) :: e2a !(in) 1 次元化された水平スペクトル・任意座標データ ! dimmension((2*mm*1)*(2*lm*1),:) real(8), dimension(size(e2a,2),-mm:mm,2*lc(ip)) :: aef_e2a !(out) 任意座標・水平 2 次元スペクトルデータ integer :: j,k,l,m if ( size(e2a,1) /= (2*mm+1)*2*lc(ip) ) & call MessageNotify('E','aef_e2a',& '1st dimension of input data invalid') !aee_e2a = reshape(transpose(e2a),(/size(e2a,2),2*mm+1,2*lc(ip)/)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do k=1,size(e2a,2) do l=1,2*lc(ip) do m=-mm,mm j = (m+mm+1) + (2*mm+1)*(l-1) aef_e2a(k,m,l) = e2a(j,k) enddo enddo enddo end function aef_e2a function e2_ef(ef) ! ! 水平スペクトルを 1 次元化転置する. ! real(8), dimension(-mm:,:), intent(in) :: ef !(in) 任意座標・水平 2 次元スペクトルデータ dimension(-mm:mm,2*lc(ip)) real(8), dimension(2*lc(ip)*(2*mm+1)) :: e2_ef !(out) 1 次元化された水平スペクトル・任意座標データ integer :: j,l,m !e2_ee = reshape(e2a_aee(reshape(ee,(/1,2*mm+1,2*lc(ip)/))), & ! (/lm*(2*mm+1)/)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do l=1,2*lc(ip) do m=-mm,mm j = (m+mm+1)+ (2*mm+1)*(l-1) e2_ef(j) = ef(m,l) enddo enddo end function e2_ef function ef_e2(e2) ! ! 水平スペクトルを転置展開する. ! real(8), dimension(:),intent(in) :: e2 !(in) 1 次元化された水平スペクトル・任意座標データ ! dimmension((2*mm*1)*2*lc(ip),:) real(8), dimension(-mm:mm,2*lc(ip)) :: ef_e2 !(out) 任意座標・水平 2 次元スペクトルデータ integer :: j,l,m !ee_e2 = reshape(aee_e2a(reshape(e2,(/(2*mm+1)*2*lc(ip),1/))), & ! (/2*mm+1,2*lc(ip)/)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do l=1,2*lc(ip) do m=-mm,mm j = (m+mm+1) + (2*mm+1)*(l-1) ef_e2(m,l) = e2(j) enddo enddo end function ef_e2 !--------------- 微分計算 ----------------- function tef_Dx_tef(tef) ! ! 入力スペクトルデータに水平微分 ∂/∂x を作用する. ! ! スペクトルデータのX微分とは, 対応する格子点データにX微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Dx_tef !(in) X微分されたスペクトルデータ integer :: n do n=0,nm tef_Dx_tef(n,:,:) = ef_Dx_ef(tef(n,:,:)) enddo end function tef_Dx_tef function tef_Dy_tef(tef) ! ! 入力スペクトルデータに水平微分 ∂/∂y を作用する. ! ! スペクトルデータのY微分とは, 対応する格子点データにY微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Dy_tef !(in) X微分されたスペクトルデータ integer :: n do n=0,nm tef_Dy_tef(n,:,:) = ef_Dy_ef(tef(n,:,:)) enddo end function tef_Dy_tef function tef_Dz_tef(tef) ! ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する. ! ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Dz_tef !(in) 鉛直微分されたスペクトルデータ tef_Dz_tef = aef_e2a(at_Dz_at(e2a_aef(tef))) end function tef_Dz_tef function tef_Lapla_tef(tef) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2 ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Lapla_tef !(out) ラプラシアンを作用された 2 次元スペクトルデータ tef_Lapla_tef = tef_LaplaH_tef(tef) + tef_Dz_tef(tef_Dz_tef(tef)) end function tef_Lapla_tef function tef_LaplaH_tef(tef) ! ! 入力スペクトルデータに水平ラプラシアン ! ! ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2 ! ! を作用する. ! ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaH_tef !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ integer :: n do n=0,nm tef_LaplaH_tef(n,:,:) = ef_Lapla_ef(tef(n,:,:)) enddo end function tef_LaplaH_tef function tef_LaplaHInv_tef(tef) ! ! 入力スペクトルデータに逆水平ラプラシアン ! ! ▽^-2_H = (∂^2/∂X^2 + ∂^2/∂Y^2)^-1 ! ! を作用する. ! ! スペクトルデータの逆水平ラプラシアンとは, 対応する格子点データに ! 逆水平ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaHInv_tef !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ integer :: n do n=0,nm tef_LaplaHInv_tef(n,:,:) = ef_LaplaInv_ef(tef(n,:,:)) enddo end function tef_LaplaHInv_tef function tef_Div_zvx_zvx_zvx(zvx_VX,zvx_VY,zvx_VZ) ! ! べクトル成分である 3 つの格子データに発散を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX !(in) ベクトル場のX成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY !(in) ベクトル場のY成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VZ !(in) ベクトル場のZ成分 real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_Div_zvx_zvx_zvx !(out) ベクトル場の発散 tef_Div_zvx_zvx_zvx = tef_Dx_tef(tef_zvx(zvx_VX)) & + tef_Dy_tef(tef_zvx(zvx_VY)) & + tef_Dz_tef(tef_zvx(zvx_VZ)) end function tef_Div_zvx_zvx_zvx !--------------- 積分計算 ----------------- !----(入力データ zvx)--- function zv_IntX_zvx(zvx) ! X積分 ! ! 3 次元格子点データのX方向積分. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(0:km,jc(ip)) :: zv_IntX_zvx !(out) X方向積分された 2 次元ZV格子点データ integer :: i zv_IntX_zvx = 0.0d0 do i=0,im-1 zv_IntX_zvx(:,:) = zv_IntX_zvx(:,:) & + zvx(:,:,i) * x_X_Weight(i) enddo end function zv_IntX_zvx function zx_IntY_zvx(zvx) ! ! 3 次元格子点データのY方向域積分. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(0:km,0:im-1) :: zx_IntY_zvx !(out) Y積分された 2 次元ZV格子点データ. real(8), dimension(0:km,0:im-1) :: zx_IntYTmp integer :: j, ierr zx_IntYTmp = 0.0d0 do j=1,jc(ip) zx_IntYTmp(:,:) = zx_IntYTmp(:,:) & + zvx(:,j,:) * v_Y_Weight(j) enddo CALL MPI_ALLREDUCE(zx_IntYTmp,zx_IntY_zvx,im*(km+1),MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function zx_IntY_zvx function vx_IntZ_zvx(zvx) ! Z積分 ! ! 3 次元格子点データのZ方向積分. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(jc(ip),0:im-1) :: vx_IntZ_zvx !(out) Z積分された 2 次元VX(水平)格子点データ integer :: k vx_IntZ_zvx = 0.0d0 do k=0,km vx_IntZ_zvx(:,:) = vx_IntZ_zvx(:,:) & + zvx(k,:,:) * z_Z_Weight(k) enddo end function vx_IntZ_zvx function x_IntZY_zvx(zvx) ! ! 3 次元格子点データのZV積分 ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元ZVX格子点データ real(8), dimension(0:im-1) :: x_IntZY_zvx !(out) ZV(子午面)積分された 1 次元X格子点データ real(8), dimension(0:im-1) :: x_IntZYTmp integer :: j, k, ierr x_IntZYtmp = 0.0D0 do j=1,jc(ip) do k=0,km x_IntZYTmp = x_IntZYTmp & + zvx(k,j,:) * v_Y_Weight(j) * z_Z_Weight(k) enddo enddo CALL MPI_ALLREDUCE(x_IntZYTmp,x_IntZY_zvx,im,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function x_IntZY_zvx function v_IntZX_zvx(zvx) ! ! 3 次元格子点データのZX積分. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元ZVX格子点データ real(8), dimension(jc(ip)) :: v_IntZX_zvx !(out) ZX積分された 1 次元Y格子点データ integer :: i, k v_IntZX_zvx = 0 do i=0,im-1 do k=0,km v_IntZX_zvx = v_IntZX_zvx & + zvx(k,:,i) * x_X_Weight(i) * z_Z_Weight(k) enddo enddo end function v_IntZX_zvx function z_IntYX_zvx(zvx) ! VX(水平)積分 ! ! 3 次元格子点データのVX(水平)積分 ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元ZVX格子点データ real(8), dimension(0:km) :: z_IntYX_zvx !(out) VX(水平)積分された 1 次元Z格子点データ real(8), dimension(0:km) :: z_IntYXTmp integer :: i, j, ierr z_IntYXTmp = 0 do j=1,jc(ip) do i=0,im-1 z_IntYXTmp = z_IntYXTmp & + zvx(:,j,i) * x_X_Weight(i) * v_Y_Weight(j) enddo enddo CALL MPI_ALLREDUCE(z_IntYXTmp,z_IntYX_zvx,km+1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function z_IntYX_zvx function IntZYX_zvx(zvx) ! ZVX(全領域)積分 ! ! 3 次元格子点データのZVX(全領域)積分 ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8) :: IntZYX_zvx !(out) 全領域積分値 real(8) :: IntZYXTmp integer :: i, j, k, ierr IntZYXTmp = 0 do i=0,im-1 do j=1,jc(ip) do k=0,km IntZYXTmp = IntZYXTmp & + zvx(k,j,i) * x_X_Weight(i) & * v_Y_Weight(j) * z_Z_Weight(k) enddo enddo enddo CALL MPI_ALLREDUCE(IntZYXTmp,IntZYX_zvx,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntZYX_zvx !----(入力データ zv)--- function z_IntY_zv(zv) ! Y積分 ! ! 2 次元(ZV)格子点データのY方向積分. ! real(8), dimension(0:km,jc(ip)), intent(in) :: zv !(in) 2 次元ZV(子午面)格子点データ real(8), dimension(0:km) :: z_IntY_zv !(out) Y積分された 1 次元Z格子点データ real(8), dimension(0:km) :: z_IntYTmp integer :: j, ierr z_IntYTmp = 0.0d0 do j=1,jc(ip) z_IntYTmp(:) = z_IntYTmp(:) + zv(:,j) * v_Y_Weight(j) enddo CALL MPI_ALLREDUCE(z_IntYTmp,z_IntY_zv,km+1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function z_IntY_zv function v_IntZ_zv(zv) ! Z積分 ! ! 2 次元(ZV)格子点データのZ方向域積分. ! real(8), dimension(0:km,jc(ip)), intent(in) :: zv !(in) 2 次元ZV格子点データ real(8), dimension(jc(ip)) :: v_IntZ_zv !(out) Z積分された 1 次元Y格子点データ integer :: k v_IntZ_zv = 0.0d0 do k=0,km v_IntZ_zv(:) = v_IntZ_zv(:) & + zv(k,:) * z_Z_Weight(k) enddo end function v_IntZ_zv function IntZY_zv(zv) ! ! 2 次元(ZV)格子点データのZV積分 ! real(8), dimension(0:km,jc(ip)), intent(in) :: zv !(in) 2 次元ZV(子午面)格子点データ real(8) :: IntZY_zv !(out) 積分値 real(8) :: IntZYTmp integer :: j, k, ierr IntZYTmp = 0 do j=1,jc(ip) do k=0,km IntZYTmp = IntZYTmp & + zv(k,j) * v_Y_Weight(j) * z_Z_Weight(k) enddo enddo CALL MPI_ALLREDUCE(IntZYTmp,IntZY_zv,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntZY_zv !----(入力データ zx)--- function z_IntX_zx(zx) ! ! 2 次元(ZX)格子点データのX方向積分. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZV格子点データ real(8), dimension(0:km) :: z_IntX_zx !(out) X積分された 1 次元Z格子点データ integer :: i z_IntX_zx = 0.0d0 do i=0,im-1 z_IntX_zx(:) = z_IntX_zx(:) + zx(:,i) * x_X_Weight(i) enddo end function z_IntX_zx function x_IntZ_zx(zx) ! ! 2 次元(ZX)格子点データのZ方向積分. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZV格子点データ real(8), dimension(0:im-1) :: x_IntZ_zx !(out) Z積分された 1 次元X格子点データ integer :: k x_IntZ_zx = 0.0d0 do k=0,km x_IntZ_zx(:) = x_IntZ_zx(:) & + zx(k,:) * z_Z_Weight(k) enddo end function x_IntZ_zx function IntZX_zx(zx) ! ZX積分 ! ! 2 次元(ZX)格子点データのZX積分 ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZV格子点データ real(8) :: IntZX_zx !(out) 積分値 integer :: i, k IntZX_zx = 0 do i=0,im-1 do k=0,km IntZX_zx = IntZX_zx & + zx(k,i) * x_X_Weight(i) * z_Z_Weight(k) enddo enddo end function IntZX_zx !----(入力データ z)--- function IntZ_z(z) ! Z積分 ! ! 1 次元(Z)格子点データのZ方向積分. ! real(8), dimension(0:km), intent(in) :: z !(in) 1 次元Z格子点データ real(8) :: IntZ_z !(out) 積分値 integer :: k IntZ_z = 0.0d0 do k=0,km IntZ_z = IntZ_z + z(k) * z_Z_Weight(k) enddo end function IntZ_z !--------------- 平均計算 ----------------- !----(入力データ zvx)--- function zv_AvrX_zvx(zvx) ! X積分 ! ! 3 次元格子点データのX方向平均. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元ZVX格子点データ real(8), dimension(0:km,jc(ip)) :: zv_AvrX_zvx !(out) X方向平均された 2 次元子午面格子点データ zv_AvrX_zvx = zv_IntX_zvx(zvx)/sum(x_X_Weight) end function zv_AvrX_zvx function zx_AvrY_zvx(zvx) ! Y平均 ! ! 3 次元格子点データのY方向域平均. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(0:km,0:im-1) :: zx_AvrY_zvx !(out) Y平均された 2 次元ZX格子点データ zx_AvrY_zvx = zx_IntY_zvx(zvx)/yl end function zx_AvrY_zvx function vx_AvrZ_zvx(zvx) ! ! 3 次元格子点データのZ方向平均. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(jc(ip),0:im-1) :: vx_AvrZ_zvx !(out) Z平均された 2 次元VX(水平)格子点データ vx_AvrZ_zvx = vx_IntZ_zvx(zvx)/sum(z_Z_Weight) end function vx_AvrZ_zvx function x_AvrZY_zvx(zvx) ! ZV積分 ! ! 3 次元格子点データのZV平均 ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元ZVX格子点データ real(8), dimension(0:im-1) :: x_AvrZY_zvx !(out) ZV平均された 1 次元X格子点データ x_AvrZY_zvx = x_IntZY_zvx(zvx)/(yl*zl) end function x_AvrZY_zvx function v_AvrZX_zvx(zvx) ! ZX積分 ! ! 3 次元格子点データのZX平均. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(jc(ip)) :: v_AvrZX_zvx !(out) ZX平均された 1 次元Y格子点データ v_AvrZX_zvx = v_IntZX_zvx(zvx) & /(sum(x_X_Weight)*sum(z_Z_Weight)) end function v_AvrZX_zvx function z_AvrYX_zvx(zvx) ! VX(水平)積分 ! ! 3 次元格子点データのVX(水平)積分 ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元格子点データ real(8), dimension(0:km) :: z_AvrYX_zvx !(out) VX(水平)平均された 1 次元Z格子点データ z_AvrYX_zvx = z_IntYX_zvx(zvx)/(xl*yl) end function z_AvrYX_zvx function AvrZYX_zvx(zvx) ! ZVX(全領域)積分 ! ! 3 次元格子点データのZVX(全領域)積分 ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx !(in) 3 次元ZVX格子点データ real(8) :: AvrZYX_zvx !(out) 全領域平均値 AvrZYX_zvx = IntZYX_zvx(zvx)/(xl*yl*zl) end function AvrZYX_zvx !----(入力データ zv)--- function z_AvrY_zv(zv) ! ! 2 次元(ZV)格子点データのY方向平均. ! real(8), dimension(0:km,jc(ip)), intent(in) :: zv !(in) 2 次元ZV格子点データ real(8), dimension(0:km) :: z_AvrY_zv !(out) Y平均された 1 次元Z格子点データ z_AvrY_zv = z_IntY_zv(zv)/yl end function z_AvrY_zv function v_AvrZ_zv(zv) ! ! 2 次元(ZV)格子点データのZ方向平均. ! real(8), dimension(0:km,jc(ip)), intent(in) :: zv !(in) 2 次元ZV格子点データ real(8), dimension(jc(ip)) :: v_AvrZ_zv !(out) Z平均された 1 次元Y格子点データ v_AvrZ_zv = v_IntZ_zv(zv)/sum(z_Z_Weight) end function v_AvrZ_zv function AvrZY_zv(zv) ! ZV平均 ! ! 2 次元(ZV)格子点データのZV平均 ! real(8), dimension(0:km,jc(ip)), intent(in) :: zv !(in) 2 次元ZV(子午面)格子点データ real(8) :: AvrZY_zv !(out) 平均値 AvrZY_zv = IntZY_zv(zv)/(yl*zl) end function AvrZY_zv !----(入力データ zx)--- function z_AvrX_zx(zx) ! X積分 ! ! 2 次元(ZX)格子点データのX方向平均. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZV格子点データ real(8), dimension(0:km) :: z_AvrX_zx !(out) X平均された 1 次元Z格子点データ z_AvrX_zx = z_IntX_zx(zx)/sum(x_X_Weight) end function z_AvrX_zx function x_AvrZ_zx(zx) ! Z平均 ! ! 2 次元(ZX)格子点データのZ方向平均. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZV格子点データ real(8), dimension(0:im-1) :: x_AvrZ_zx !(out) Z平均された 1 次元X格子点データ x_AvrZ_zx = x_IntZ_zx(zx)/sum(z_Z_Weight) end function x_AvrZ_zx function AvrZX_zx(zx) ! ZX積分 ! ! 2 次元(ZX)格子点データのZX平均 ! real(8), dimension(0:km,0:im-1), intent(in) :: zx ! (in) 2 次元格子点データ real(8) :: AvrZX_zx ! 平均値 AvrZX_zx = IntZX_zx(zx)/(sum(x_X_Weight)*sum(z_Z_Weight)) end function AvrZX_zx !----(入力データ z)--- function AvrZ_z(z) ! ! 1 次元(Z)格子点データのZ方向平均. ! real(8), dimension(0:km), intent(in) :: z !(in) 1 次元Z格子点データ real(8) :: AvrZ_z !(out) 平均値 AvrZ_z = IntZ_z(z)/sum(z_Z_Weight) end function AvrZ_z !--------------- ポロイダル/トロイダルモデル用微分 ----------------- !!$ function wt_KxRGrad_wt(wt) !!$ ! !!$ ! 入力スペクトルデータにX微分 k×r・▽ = ∂/∂λを作用する. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_KxRGrad_wt !!$ !(out) X微分を作用された 2 次元スペクトルデータ !!$ !!$ wt_KxRGrad_wt = wa_Dlon_wa(wt) !!$ !!$ end function wt_KxRGrad_wt !!$ !!$ function zvx_KGrad_wt(wt) ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r !!$ ! !!$ ! 入力スペクトルデータに対応する格子データに軸方向微分 !!$ ! !!$ ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r !!$ ! !!$ ! を作用させた格子データが返される. !!$ ! ここでベクトル k は球の中心から北極向きの単位ベクトルである. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(0:km,0:jm-1,0:im-1) :: yzx_KGrad_wt !!$ !(out) 軸方向微分を作用された 2 次元スペクトルデータ !!$ !!$ xzv_KGrad_wt = cos(xyz_Y)*xyz_GradY_wt(wt) & !!$ + sin(xyz_Y)*xyz_wt(wt_Drad_wt(wt)) !!$ !!$ end function xyz_KGrad_wt !!$ !!$ function wt_L2_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する. !!$ ! !!$ ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. !!$ ! 入力スペクトルデ ータに対応する格子点データに演算子 !!$ ! !!$ ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ) !!$ ! !!$ ! を作用させたデータのスペクトル変換が返される. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_L2_wt !!$ !(out) L^2 演算子を作用された 2 次元スペクトルデータ !!$ !!$ wt_L2_wt = -wa_Lapla_wa(wt) !!$ !!$ end function wt_L2_wt !!$ !!$ function wt_L2Inv_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を !!$ ! 作用する. !!$ ! !!$ ! スペクトルデータに L^2 演算子を作用させる関数 wt_L2_wt の逆計算を !!$ ! 行う関数である. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_L2Inv_wt !!$ !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ !!$ !!$ wt_L2Inv_wt = -wa_LaplaInv_wa(wt) !!$ !!$ end function wt_L2Inv_wt !!$ !!$ function wt_QOperator_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに対応する格子点データに演算子 !!$ ! !!$ ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) !!$ ! !!$ ! を作用させたデータのスペクトル変換が返される. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_QOperator_wt !!$ !(out) Q 演算子を作用された 2 次元スペクトルデータ !!$ !!$ wt_QOperator_wt = & !!$ wt_xyz(xyz_KGrad_wt(wt) - xyz_KGrad_wt(wt_L2_wt(wt))/2) & !!$ - wt_L2_wt(wt_xyz(xyz_KGrad_wt(wt)))/2 !!$ !!$ end function wt_QOperator_wt function tef_ZRot_zvx_zvx(zvx_VX,zvx_VY) ! z・(▽×v) ! ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する. ! ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す. ! ! z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y ! ! のスペクトル データが返される. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX !(in) ベクトルのX成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY !(in) ベクトルのY成分 real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_ZRot_zvx_zvx !(out) ベクトルの渦度とZベクトルの内積 tef_ZRot_zvx_zvx = tef_DX_tef(tef_zvx(zvx_VY)) & - tef_DY_tef(tef_zvx(zvx_VX)) end function tef_ZRot_zvx_zvx function tef_ZRotRot_zvx_zvx_zvx(zvx_VX,zvx_VY,zvx_VZ) ! ! ベクトル v に対して z・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分, ! Y成分, Z成分を表す. ! ! z・(▽×▽×v) = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) ) ! - ▽_H^2 v[z] ! ! のスペクトルデータが返される. ! real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX !(in) ベクトルのX成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY !(in) ベクトルのY成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VZ !(in) ベクトルのZ成分 real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_ZRotRot_zvx_zvx_zvx !(out) ベクトル v の z・(▽×▽×v) tef_ZRotRot_zvx_zvx_zvx = & tef_DZ_tef( tef_DX_tef(tef_zvx(zvx_VX)) & + tef_DY_tef(tef_zvx(zvx_VY)) ) & - tef_LaplaH_tef(tef_zvx(zvx_VZ)) end function tef_ZRotRot_zvx_zvx_zvx subroutine tef_Potential2Vector(& zvx_VX,zvx_VY,zvx_VZ,tef_Torvel,tef_Polvel) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψz) + ▽x▽x(Φz) ! ! の各成分を計算する ! real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_VX !(out) ベクトル場のX成分 real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_VY !(out) ベクトル場のY成分 real(8), dimension(0:km,jc(ip),0:im-1) :: zvx_VZ !(out) ベクトル場のZ成分 real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Torvel !(in) トロイダルポテンシャル real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Polvel !(in) ポロイダルポテンシャル zvx_VX = zvx_tef( tef_DY_tef(tef_Torvel) & + tef_DX_tef(tef_DZ_tef(tef_Polvel)) ) zvx_VY = zvx_tef( - tef_DX_tef(tef_Torvel) & + tef_DY_tef(tef_DZ_tef(tef_Polvel)) ) zvx_VZ = -zvx_tef(tef_LaplaH_tef(tef_Polvel)) end subroutine tef_Potential2Vector subroutine tef_Potential2Rotation(& zvx_RotVX,zvx_RotVY,zvx_RotVZ,tef_Torvel,tef_Polvel) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψz) + ▽x▽x(Φz) ! ! に対して, その回転 ! ! ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z) ! ! を計算する. ! ベクトル場の回転 real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVX !(out) 回転のX成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVY !(out) 回転のY成分 real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVZ !(out) 回転のZ成分 ! 入力ベクトル場を表すポテンシャル real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Torvel !(in) トロイダルポテンシャル real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Polvel !(in) ポロイダルポテンシャル call tef_Potential2Vector( & zvx_RotVX,zvx_RotVY,zvx_RotVZ, & -tef_Lapla_tef(tef_Polvel), tef_Torvel) end subroutine tef_Potential2Rotation !--------------- 補間計算 ----------------- function Interpolate_tef(tef_data,x,y,z) ! ! (x,y,z) における関数値を ! そのスペクトル係数 tef_data から補間計算する ! real(8), intent(IN) :: tef_data(0:nm,-mm:mm,2*lc(ip))! スペクトルデータ real(8), intent(IN) :: x ! 補間する位置(X) real(8), intent(IN) :: y ! 補間する位置(Y) real(8), intent(IN) :: z ! 補間する位置(Z) real(8) :: Interpolate_tef ! 補間した値 Interpolate_tef = & Interpolate_ef(ef_e2(a_Interpolate_at(e2a_aef(tef_data),z)),x,y) end function Interpolate_tef !--------------- ポロイダル/トロイダルモデル用スペクトル解析 ---------------- function zef_ToroidalEnergySpectrum_tef(tef_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 水平 X 波数 l, Y 波数 m の各成分を計算する ! ! * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー ! スペクトルは (1/2)(l**2+m**2)|ψ(l,m,z)|^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での ! 水平平均エネルギーに等しい. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_TORPOT !(in) トロイダルポテンシャル real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_ToroidalEnergySpectrum_tef !(out) エネルギースペクトルトロイダル成分 real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Work integer :: l, m zef_ToroidalEnergySpectrum_tef = 0.0D0 zef_Work = zef_tef(tef_Torpot) do l=-lm,lm do m=-mm,mm if ( lf_l(l) > 0 ) then zef_ToroidalEnergySpectrum_tef(:,m,lf_l(l)) & = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) & * ( zef_Work(:,m,lf_l(l))**2 + zef_Work(:,-m,lf_l(-l))**2 ) endif enddo enddo end function zef_ToroidalEnergySpectrum_tef !!$ function nz_ToroidalEnergySpectrum_wt(wt_TORPOT) !!$ ! !!$ ! トロイダルポテンシャルから, トロイダルエネルギーの !!$ ! 球面調和函数全波数の各成分を計算する. !!$ ! !!$ ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分 !!$ ! ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは !!$ ! Σ[m=-n]^n(1/2)n(n+1)4πr^2ψ(n,m,r)^2 と計算される. !!$ ! !!$ ! * 全てのエネルギースペクトル成分の和をZ積分したもの(r^2の重み無し) !!$ ! が球殻内での全エネルギーに等しい. !!$ ! !!$ real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_TORPOT !!$ !(in) トロイダルポテンシャル !!$ !!$ real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wt !!$ !(out) エネルギースペクトルトロイダル成分 !!$ !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA ! 作業領域 !!$ integer :: n, m !!$ !!$ wz_DATA = wz_wt(wt_TORPOT) !!$ do n=0,nm !!$ nz_ToroidalEnergySpectrum_wt(n,:) & !!$ = 0.5 * n*(n+1)* (4*pi) * z_Z**2 & !!$ * sum(wz_Data(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) !!$ enddo !!$ !!$ end function nz_ToroidalEnergySpectrum_wt function zef_PoloidalEnergySpectrum_tef(tef_POLPOT) ! ! ポロイダルポテンシャルから, ポロイダルエネルギーの ! 水平 X 波数 l, Y 波数 m の各成分を計算する. ! ! * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分 ! φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー ! スペクトルは ! ! (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2} ! ! と計算される. ! ! * 全てのエネルギースペクトル成分の和をZ積分したもの ! が水平平均エネルギーに等しい. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_POLPOT !(in) ポロイダルポテンシャル real(8), dimension(0:km,-nm:nm,2*lc(ip)) :: zef_PoloidalEnergySpectrum_tef !(out) エネルギースペクトルポロイダル成分 real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Data ! 作業領域 real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_DData ! 作業領域 integer :: l, m zef_PoloidalEnergySpectrum_tef = 0.0D0 zef_Data = zef_tef(tef_POLPOT) zef_DData = zef_tef(tef_DZ_tef(tef_POLPOT)) do l=-lm,lm do m=-mm,mm if ( lf_l(l) > 0 ) then zef_PoloidalEnergySpectrum_tef(:,m,lf_l(l)) = & + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2) & *( zef_DData(:,m,lf_l(l))**2 + zef_DData(:,-m,lf_l(-l))**2 & + ((2*pi*l/xl)**2+(2*pi*m/yl)**2) & *( zef_Data(:,m,lf_l(l))**2 + zef_Data(:,-m,lf_l(-l))**2)) endif enddo enddo end function zef_PoloidalEnergySpectrum_tef !!$ function nz_PoloidalEnergySpectrum_wt(wt_POLPOT) !!$ ! !!$ ! ポロイダルポテンシャルから, ポロイダルエネルギーの !!$ ! 球面調和函数全波数の各成分を計算する !!$ ! !!$ ! * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分 !!$ ! φ(n,m,r)から全波数 n 成分のポロイダルエネルギースペクトルは !!$ ! !!$ ! Σ[m=-n]^n ((1/2)n(n+1)4πr^2{[d(rφ(n,m,r))/dr]^2 !!$ ! + n(n+1)φ(n,m,r)^2} !!$ ! !!$ ! と計算される. !!$ ! !!$ ! * 全ての全波数に対してのエネルギースペクトル成分の和をZ積分したもの !!$ ! (r^2の重み無し)が球殻内での全エネルギーに等しい. !!$ ! !!$ real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_POLPOT !!$ !(in) ポロイダルポテンシャル !!$ !!$ real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_wt !!$ !(out) エネルギースペクトルポロイダル成分 !!$ !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1 ! 作業領域 !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2 ! 作業領域 !!$ integer :: n, m !!$ !!$ wz_Data1 = wz_wt(wt_POLPOT) !!$ wz_Data2 = wz_Z*wz_wt(wt_DZ_wt(wt_POLPOT)) & ! d(rφ)/dr !!$ + wz_wt(wt_POLPOT) ! = rdφ/dr+φ !!$ !!$ do n=0,nm !!$ nz_PoloidalEnergySpectrum_wt(n,:) = & !!$ + 0.5* n*(n+1)* (4*pi) & !!$ *( sum(wz_Data2(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) & !!$ + n*(n+1)*sum(wz_Data1(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) ) !!$ enddo !!$ !!$ end function nz_PoloidalEnergySpectrum_wt !--------------- 境界値問題 ----------------- subroutine tef_BoundariesTau(tef,values,cond) ! ! スペクトルデータにディリクレ・ノイマン境界条件を適用する ! Chebyshev 空間での境界条件適用(タウ法) ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を ! 定める方法をとっている(タウ法). ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(2,-mm:mm,2*lc(ip)), intent(in), optional :: values !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える. ! 省略時は値/勾配 0 となる. character(len=2), intent(in), optional :: cond !(in) 境界条件. 省略時は 'DD' ! DD : 両端ディリクレ条件 ! DN : 上端ディリクレ, 下端ノイマン条件 ! ND : 上端ノイマン, 下端ディリクレ条件 ! NN : 両端ノイマン条件 real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t real(8), dimension(2*lc(ip)*(2*mm+1),2) :: e22_values character(len=2) :: Bcond e2t = e2a_aef(tef) if (present(values)) then e22_values = e2a_aef(values) endif if (.not. present(cond)) then Bcond='DD' else Bcond = cond endif select case(Bcond) case ('NN') if (present(values)) then call at_BoundariesTau_NN(e2t,e22_values) else call at_BoundariesTau_NN(e2t) endif case ('DN') if (present(values)) then call at_BoundariesTau_DN(e2t,e22_values) else call at_BoundariesTau_DN(e2t) endif case ('ND') if (present(values)) then call at_BoundariesTau_ND(e2t,e22_values) else call at_BoundariesTau_ND(e2t) endif case ('DD') if (present(values)) then call at_BoundariesTau_DD(e2t,e22_values) else call at_BoundariesTau_DD(e2t) endif case default call MessageNotify('E','tef_BoundariesTau','B.C. not supported') end select tef = aef_e2a(e2t) end subroutine tef_BoundariesTau subroutine tef_BoundariesGrid(tef,values,cond) ! ! スペクトルデータにディリクレ・ノイマン境界条件を適用する ! 実空間での境界条件適用 ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(2,-mm:mm,2*lc(ip)), intent(in), optional :: values !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える. ! 省略時は値/勾配 0 となる. character(len=2), intent(in), optional :: cond !(in) 境界条件. 省略時は 'DD' ! DD : 両端ディリクレ条件 ! DN : 上端ディリクレ, 下端ノイマン条件 ! ND : 上端ノイマン, 下端ディリクレ条件 ! NN : 両端ノイマン条件 real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t real(8), dimension(2*lc(ip)*(2*mm+1),2) :: e22_values character(len=2) :: Bcond e2t = e2a_aef(tef) if (present(values)) then e22_values = e2a_aef(values) endif if (.not. present(cond)) then BCond = 'DD' else BCond = cond endif select case(Bcond) case ('NN') if (present(values)) then call at_BoundariesGrid_NN(e2t,e22_values) else call at_BoundariesGrid_NN(e2t) endif case ('DN') if (present(values)) then call at_BoundariesGrid_DN(e2t,e22_values) else call at_BoundariesGrid_DN(e2t) endif case ('ND') if (present(values)) then call at_BoundariesGrid_ND(e2t,e22_values) else call at_BoundariesGrid_ND(e2t) endif case ('DD') if (present(values)) then call at_BoundariesGrid_DD(e2t,e22_values) else call at_BoundariesGrid_DD(e2t) endif case default call MessageNotify('E','tef_BoundariesGrid','B.C. not supported') end select tef = aef_e2a(e2t) end subroutine tef_BoundariesGrid subroutine tef_TorBoundariesTau(tef_TOR,cond,new) ! トロイダル速度ポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). トロイダルポテンシャルの境界条件は ! ! ψ = 0 at boundaries (粘着条件) ! ∂ψ/∂z = 0 at boundaries (応力なし条件) ! ! であるから tef_Boundaries で対応可能だが, 将来のため別途作成しておく. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. character(len=2), intent(in), optional :: cond !(in) 境界条件スイッチ. 省略時は 'RR' ! RR : 両端粘着条件 ! RF : 上端粘着, 下端応力なし条件 ! FR : 上端応力なし, 下端粘着条件 ! FF : 両端応力なし条件 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: rigid1 = .true. logical :: rigid2 = .true. logical :: first = .true. logical :: new_matrix = .false. integer :: n save :: alu, kp, first select case (cond) case ('RR') rigid1 = .TRUE. ; rigid2 = .TRUE. case ('RF') rigid1 = .TRUE. ; rigid2 = .FALSE. case ('FR') rigid1 = .FALSE. ; rigid2 = .TRUE. case ('FF') rigid1 = .FALSE. ; rigid2 = .FALSE. case default call MessageNotify('E','tef_TorBoundariesTau','B.C. not supported') end select 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) allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work= 0.0D0 ; e2t_work(:,n)= 1.0D0 alu(:,:,n) = e2t_work ! 粘着条件 ! 力学的条件粘着条件 if ( rigid1 ) then e2z_work = az_at(e2t_work) else e2z_work = az_at(at_Dz_at(e2t_work)) endif alu(:,nm-1,n) = e2z_work(:,0) ! 力学的条件粘着条件 if ( rigid2 ) then e2z_work = az_at(e2t_work) else e2z_work = az_at(at_Dz_at(e2t_work)) endif alu(:,nm,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','tef_TorBoundariesTau',& 'Matrix to apply b.c. newly produced.') endif e2t_work = e2a_aef(tef_Tor) e2t_work(:,nm-1) = 0.0D0 e2t_work(:,nm) = 0.0D0 tef_Tor = aef_e2a(lusolve(alu,kp,e2t_work)) end subroutine tef_TorBoundariesTau subroutine tef_TorBoundariesGrid(tef_TOR,cond,new) ! ! トロイダル速度ポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! ! トロイダル速度ポテンシャルの境界条件は ! ! ψ = 0 at boundaries (粘着条件) ! ∂ψ/∂z = 0 at boundaries (応力なし条件) ! ! であるので tef_Boundaries で対応可能だが, 将来のため別途作成しておく ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. character(len=2), intent(in), optional :: cond !(in) 境界条件スイッチ. 省略時は 'RR' ! RR : 両端粘着条件 ! RF : 上端粘着, 下端応力なし条件 ! FR : 上端応力なし, 下端粘着条件 ! FF : 両端応力なし条件 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: rigid1 = .true. logical :: rigid2 = .true. logical :: first = .true. logical :: new_matrix = .false. integer :: n save :: alu, kp, first select case (cond) case ('RR') rigid1 = .TRUE. ; rigid2 = .TRUE. case ('RF') rigid1 = .TRUE. ; rigid2 = .FALSE. case ('FR') rigid1 = .FALSE. ; rigid2 = .TRUE. case ('FF') rigid1 = .FALSE. ; rigid2 = .FALSE. case default call MessageNotify('E','tef_TorBoundariesGrid','B.C. not supported') end select if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( nm /= km ) then call MessageNotify('E','tef_TorBoundariesGrid', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0 e2z_work = az_at(e2t_work) alu(:,:,n) = e2z_work ! 内部領域は値そのまま. ! 粘着条件 ! 力学的条件粘着条件 if ( rigid1 ) then e2z_work = az_at(e2t_work) else e2z_work=az_at(at_Dz_at(e2t_work)) endif alu(:,0,n) = e2z_work(:,0) ! 力学的条件粘着条件 if ( rigid2 ) then e2z_work = az_at(e2t_work) else e2z_work=az_at(at_Dz_at(e2t_work)) endif alu(:,km,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','TorBoundariesGrid',& 'Matrix to apply b.c. newly produced.') endif e2z_work = az_at(e2a_aef(tef_Tor)) e2z_work(:,0) = 0.0D0 e2z_work(:,km) = 0.0D0 tef_TOR = aef_e2a(lusolve(alu,kp,e2z_work)) end subroutine tef_TorBoundariesGrid function tef_LaplaPol2PolTau_tef(tef,cond,new) ! ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する. ! チェビシェフ空間でタウ法による境界条件を適用している. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂z = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂z^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(in) :: tef !(in) 入力▽^2φ分布 real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaPol2PolTau_tef !(out) 出力ポロイダルポテンシャル分布 character(len=2), intent(in), optional :: cond !(in) 境界条件スイッチ. 省略時は 'RR' ! RR : 両端粘着条件 ! RF : 上端粘着, 下端応力なし条件 ! FR : 上端応力なし, 下端粘着条件 ! FF : 両端応力なし条件 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: rigid1 = .true. logical :: rigid2 = .true. logical :: first = .true. logical :: new_matrix = .false. integer :: n save :: alu, kp, first if ( present(cond) ) then select case (cond) case ('RR') rigid1 = .TRUE. ; rigid2 = .TRUE. case ('RF') rigid1 = .TRUE. ; rigid2 = .FALSE. case ('FR') rigid1 = .FALSE. ; rigid2 = .TRUE. case ('FF') rigid1 = .FALSE. ; rigid2 = .FALSE. case default call MessageNotify('E','tef_LaplaPol2PolTau_tef','B.C. not supported') end select else rigid1 = .TRUE. ; rigid2 = .TRUE. endif 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) allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) alu = 0.0d0 do n=0,nm-2 e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 ! 各水平波数に関して独立の式 alu(:,:,n) = e2a_aef(tef_lapla_tef(aef_e2a(e2t_work))) enddo do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 e2z_work = az_at(e2t_work) ! 運動学的条件. 流線は境界で一定 alu(:,nm,n) = e2z_work(:,0) alu(:,nm-1,n) = e2z_work(:,km) ! 力学的条件粘着条件 e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 if ( rigid1 ) then e2z_work=az_at(at_Dz_at(e2t_work)) else e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work))) endif alu(:,nm-2,n) = e2z_work(:,0) ! 力学的条件粘着条件 e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 if ( rigid2 ) then e2z_work=az_at(at_Dz_at(e2t_work)) else e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work))) endif alu(:,nm-3,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','tef_LaplaPol2PolTau_tef',& 'Matrix to apply b.c. newly produced.') endif e2t_work = e2a_aef(tef) e2t_work(:,nm) = 0.0D0 ! 運動学的条件 e2t_work(:,nm-1) = 0.0D0 ! 運動学的条件 e2t_work(:,nm-2) = 0.0D0 ! 力学的条件 e2t_work(:,nm-3) = 0.0D0 ! 力学的条件 tef_LaplaPol2PolTau_tef = aef_e2a(lusolve(alu,kp,e2t_work)) end function tef_LaplaPol2PolTau_tef function zef_LaplaPol2Pol_zef(zef,cond,new) ! ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する. ! ! チェビシェフ格子点空間で境界条件を適用している. ! この関数を用いるためには wt_Initial にて設定する ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく ! しておく必要がある. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂z = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂z^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:km,-mm:mm,2*lc(ip)),intent(in) :: zef !(in) 入力▽^2φ分布 real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_LaplaPol2Pol_zef !(out) 出力ポロイダルポテンシャル分布 character(len=2), intent(in), optional :: cond !(in) 境界条件スイッチ. 省略時は 'RR' ! RR : 両端粘着条件 ! RF : 上端粘着, 下端応力なし条件 ! FR : 上端応力なし, 下端粘着条件 ! FF : 両端応力なし条件 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: rigid1 = .true. logical :: rigid2 = .true. logical :: first = .true. logical :: new_matrix = .false. integer :: k save :: alu, kp, first select case (cond) case ('RR') rigid1 = .TRUE. ; rigid2 = .TRUE. case ('RF') rigid1 = .TRUE. ; rigid2 = .FALSE. case ('FR') rigid1 = .FALSE. ; rigid2 = .TRUE. case ('FF') rigid1 = .FALSE. ; rigid2 = .FALSE. case default call MessageNotify('E','zef_laplapol2pol_zef','B.C. not supported') end select if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( nm /= km ) then call MessageNotify('E','zef_LaplaPol2Pol_zef', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:km),kp(2*lc(ip)*(2*mm+1),0:km)) do k=0,km e2z_work = 0.0D0 ; e2z_work(:,k) = 1.0D0 ! 各水平波数に関して独立の式 alu(:,:,k) & = e2a_aef(zef_tef(tef_lapla_tef(tef_zef(aef_e2a(e2z_work))))) enddo do k=0,km e2z_work = 0.0D0 ; e2z_work(:,k) = 1.0D0 ! 運動学的条件. 流線は境界で一定 alu(:,0,k) = e2z_work(:,0) alu(:,km,k) = e2z_work(:,km) ! 力学的条件粘着条件 e2z_work = 0.0D0 ; e2z_work(:,k) = 1.0D0 if ( rigid1 ) then e2z_work=az_at(at_Dz_at(at_az(e2z_work))) else e2z_work=az_at(at_Dz_at(at_Dz_at(at_az(e2z_work)))) endif alu(:,1,k) = e2z_work(:,0) ! 力学的条件粘着条件 e2z_work = 0.0D0 ; e2z_work(:,k) = 1.0D0 if ( rigid2 ) then e2z_work=az_at(at_Dz_at(at_az(e2z_work))) else e2z_work=az_at(at_Dz_at(at_Dz_at(at_az(e2z_work)))) endif alu(:,km-1,k) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','zef_LaplaPol2Pol_zef',& 'Matrix to apply b.c. newly produced.') endif e2z_work = e2a_aef(zef) e2z_work(:,1) = 0.0D0 ! 力学的条件 e2z_work(:,km-1) = 0.0D0 ! 力学的条件 e2z_work(:,0) = 0.0D0 ! 運動学的条件 e2z_work(:,km) = 0.0D0 ! 運動学的条件 e2z_work = lusolve(alu,kp,e2z_work) zef_laplapol2pol_zef = aef_e2a(e2z_work) end function zef_LaplaPol2Pol_zef function tef_LaplaPol2PolGrid_tef(tef,cond,new) ! ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する. ! チェビシェフ格子点空間で境界条件を適用している. ! ! この関数を用いるためには tef_Initial にて設定する ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく ! しておく必要がある. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂z = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂z^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! ! 最終的にチェビシェフ係数の解が欲しい場合には, tef_LaplaPol2Pol_tef に ! 比べてチェビシェフ -- 格子点変換が 1 回分少なくて済む. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(in) :: tef !(in) 入力▽^2φ分布 real(8), dimension(0:nm,-mm:mm,2*lc(ip)) :: tef_LaplaPol2PolGrid_tef !(out) 出力ポロイダルポテンシャル分布 character(len=2), intent(in), optional :: cond !(in) 境界条件スイッチ. 省略時は 'RR' ! RR : 両端粘着条件 ! RF : 上端粘着, 下端応力なし条件 ! FR : 上端応力なし, 下端粘着条件 ! FF : 両端応力なし条件 logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: rigid1 = .true. logical :: rigid2 = .true. logical :: first = .true. logical :: new_matrix = .false. integer :: n save :: alu, kp, first select case (cond) case ('RR') rigid1 = .TRUE. ; rigid2 = .TRUE. case ('RF') rigid1 = .TRUE. ; rigid2 = .FALSE. case ('FR') rigid1 = .FALSE. ; rigid2 = .TRUE. case ('FF') rigid1 = .FALSE. ; rigid2 = .FALSE. case default call MessageNotify('E','tef_LaplaPol2PolGrid_tef','B.C. not supported') end select if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( nm /= km ) then call MessageNotify('E','tef_LaplaPol2PolGrid_tef', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 ! 各水平波数に関して独立の式 alu(:,:,n) = e2a_aef(zef_tef(tef_lapla_tef(aef_e2a(e2t_work)))) enddo do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 e2z_work = az_at(e2t_work) ! 運動学的条件. 流線は境界で一定 alu(:,0,n) = e2z_work(:,0) alu(:,km,n) = e2z_work(:,km) ! 力学的条件粘着条件 e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 if ( rigid1 ) then e2z_work=az_at(at_Dz_at(e2t_work)) else e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work))) endif alu(:,1,n) = e2z_work(:,0) ! 力学的条件粘着条件 e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0 if ( rigid2 ) then e2z_work=az_at(at_Dz_at(e2t_work)) else e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work))) endif alu(:,km-1,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','tef_LaplaPol2PolGrid_tef',& 'Matrix to apply b.c. newly produced.') endif e2z_work = az_at(e2a_aef(tef)) e2z_work(:,1) = 0.0D0 ! 力学的条件 e2z_work(:,km-1) = 0.0D0 ! 力学的条件 e2z_work(:,0) = 0.0D0 ! 運動学的条件 e2z_work(:,km) = 0.0D0 ! 運動学的条件 tef_LaplaPol2PolGrid_tef = aef_e2a(lusolve(alu,kp,e2z_work)) end function tef_LaplaPol2PolGrid_tef subroutine tef_TormagBoundariesTau(tef_TOR,new) ! 磁場トロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ ! 対応している. その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 上側 ! tef_psi = 0 at the outer boundary ! 下側 ! tef_psi = 0 at the inner boundary ! ! であるから tef_Boundaries で対応可能だが, 将来のため別途作成しておく. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: first = .true. logical :: new_matrix = .false. integer :: n 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) allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work= 0.0D0 ; e2t_work(:,n)= 1.0D0 alu(:,:,n) = e2t_work ! 非電気伝導体 e2z_work = az_at(e2t_work) alu(:,nm-1,n) = e2z_work(:,0) alu(:,nm,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','tef_TormagBoundariesTau',& 'Matrix to apply b.c. newly produced.') endif e2t_work = e2a_aef(tef_Tor) e2t_work(:,nm-1) = 0.0D0 e2t_work(:,nm) = 0.0D0 tef_Tor = aef_e2a(lusolve(alu,kp,e2t_work)) end subroutine tef_TormagBoundariesTau subroutine tef_TormagBoundariesGrid(tef_TOR,new) ! ! 磁場トロイダルポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 上側 ! tef_psi = 0 at the outer boundary ! 下側 ! tef_psi = 0 at the inner boundary ! ! であるので tef_Boundaries で対応可能だが, 将来のため別途作成しておく ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work logical :: first = .true. logical :: new_matrix = .false. integer :: n 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 ( nm /= km ) then call MessageNotify('E','tef_TorMagBoundariesGrid', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0 e2z_work = az_at(e2t_work) alu(:,:,n) = e2z_work ! 内部領域は値そのまま. ! 非電気伝導体 alu(:,0,n) = e2z_work(:,0) alu(:,km,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','TormagBoundariesGrid',& 'Matrix to apply b.c. newly produced.') endif e2z_work = az_at(e2a_aef(tef_Tor)) e2z_work(:,0) = 0.0D0 e2z_work(:,km) = 0.0D0 tef_TOR = aef_e2a(lusolve(alu,kp,e2z_work)) end subroutine tef_TormagBoundariesGrid subroutine tef_PolmagBoundariesTau(tef_POL,new) ! ! 磁場ポロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ ! 対応している. その場合, 磁場ポロイダルポテンシャルの各水平スペクトル ! 成分 h にたいして境界条件が与えられ, ! ! * 上側境界 : dh/dz + K h = 0 ! * 下側境界 : dh/dz - K h = 0 ! ! である. ここで K=sqrt(l^2+m^2) は h の水平全波数である. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_POL !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_kh logical :: first = .true. logical :: new_matrix = .false. real(8) :: rl, rm integer :: l, m, n, e2index 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) allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0 alu(:,:,n) = e2t_work(:,:) ! 内部領域は同じ enddo ! 非電気伝導体 do m=-mm,mm do l=-lm,lm if ( lf_l(l) > 0 ) then e2index = (2*mm+1)*(lf_l(l)-1) + (m+mm+1) rl = 2*PI/xl*l ; rm = 2*PI/yl*m e2t_kh(e2index,:) = sqrt(rl**2+rm**2) endif enddo enddo do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0 e2z_work = az_at(at_Dz_at(e2t_work) + e2t_kh * e2t_work) alu(:,nm-1,n) = e2z_work(:,0) e2z_work = az_at(at_Dz_at(e2t_work) - e2t_kh * e2t_work) alu(:,nm,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','tef_PolmagBoundariesTau',& 'Matrix to apply b.c. newly produced.') endif e2t_work = e2a_aef(tef_POL) e2t_work(:,nm-1) = 0.0D0 e2t_work(:,nm) = 0.0D0 tef_POL = aef_e2a(lusolve(alu,kp,e2t_work)) end subroutine tef_PolmagBoundariesTau subroutine tef_PolmagBoundariesGrid(tef_POL,new) ! ! 磁場ポロイダルポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場ポロイダルポテンシャルの各水平スペクトル成分 h に ! たいして境界条件が与えられ, ! ! * 上側境界 : dh/dz + K h = 0 ! * 下側境界 : dh/dz - K h = 0 ! ! である. ここで K=sqrt(l^2+m^2) は h の水平全波数である. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout) :: tef_POL !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_kh logical :: first = .true. logical :: new_matrix = .false. real(8) :: rl, rm integer :: l, m, n, e2index 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 ( nm /= km ) then call MessageNotify('E','tef_PolMagBoundariesGrid', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm)) do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0 e2z_work = az_at(e2t_work) alu(:,:,n) = e2z_work ! 内部領域は同じ enddo ! 非電気伝導体 do m=-mm,mm do l=-lm,lm if ( lf_l(l) > 0 ) then e2index = (2*mm+1)*(lf_l(l)-1) + (m+mm+1) rl = 2*PI/xl*l ; rm = 2*PI/yl*m e2t_kh(e2index,:) = sqrt(rl**2+rm**2) endif enddo enddo do n=0,nm e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0 e2z_work = az_at(at_Dz_at(e2t_work) + e2t_kh * e2t_work) alu(:,0,n) = e2z_work(:,0) e2z_work = az_at(at_Dz_at(e2t_work) - e2t_kh * e2t_work) alu(:,km,n) = e2z_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','tef_PolmagBoundariesGrid',& 'Matrix to apply b.c. newly produced.') endif e2z_work = az_at(e2a_aef(tef_POL)) e2z_work(:,0) = 0.0D0 e2z_work(:,km) = 0.0D0 tef_POL = aef_e2a(lusolve(alu,kp,e2z_work)) end subroutine tef_PolmagBoundariesGrid end module tee_mpi_module