!-- !---------------------------------------------------------------------- ! Copyright (c) 2016 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 ut_mpi_module_svpack ! ! spml/ut_mpi_module_svpack モジュールは球面上および球殻内での流体運動を ! スペクトル法と MPI 並列化によって数値計算するための Fortran90 ! 関数を提供するものである. ! ! 水平方向に球面調和函数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 関数, サブルーチンの名前と機能は wt_mpi_module のものと同じである. ! したがって use 文を wt_mpi_module から wt_mpi_module_supack に ! 変更するだけで SUPACK-MPI の機能が使えるようになる. ! ! ただし l_nm, nm_l の使い方には注意されたい. wt_mpi_module の l_nm, nm_l は ! wt_Initial で初期化しなくとも用いることができる(結果が切断波数に依らない)が, ! wt_mpi_module_supack のものは初期化したのちにしか使うことができない. ! ! 内部で ua_mpi_module_svpack を用いている. ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして ! ISPACK2 の Fortran77 サブルーチンを用いている. ! ! !履歴 2016/07/29 竹広真一 wt_mpi_module_supack より svpack 用に改造 ! !凡例 ! データ種類と index ! x : 経度 y : 緯度 z : 動径 h : 動径(分散格子) ! w : 球面調和関数スペクトル ! u : 球面調和関数スペクトル(分散データ) ! n : 球面調和関数スペクトル(水平全波数) ! m : 球面調和関数スペクトル(水平東西波数) ! t : チェビシェフ関数スペクトル ! a : 任意の次元 ! ! xyz : 3 次元格子点データ ! xvz : 3 次元分散格子点データ ! xy : 水平 2 次元格子点データ ! yz : 子午面 2 次元格子点データ ! xz : 緯度面 2 次元格子点データ ! ! wz : 水平スペクトル動径格子点データ ! wt : スペクトルデータ ! !++ module ut_mpi_module_svpack ! != ut_mpi_module_svpack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id: wt_mpi_module_supack.f90 598 2013-08-20 03:23:44Z takepiro $ ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== 概要 ! ! spml/ut_mpi_module_svpack モジュールは球面上および球殻内での流体運動を ! スペクトル法と MPI 並列化によって数値計算するための Fortran90 ! 関数を提供するものである. ! ! 水平方向に球面調和函数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ua_mpi_module_svpack を用いている. ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして ! ISPACK2 の Fortran77 サブルーチンを用いている. ! ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭(wt_, nmz_, nz_, xyz_, xvz_ wz_, w_, xy_, x_, y_, z_, a_)は, ! 返す値の形を示している. ! wt_ :: スペクトルデータ(球面調和函数・チェビシェフ変換) ! nmz_:: 水平スペクトルデータ(全波数 n, 帯状波数各成分, 動径) ! nz_ :: 水平スペクトルデータ(全波数 n, 動径) ! xyz_ :: 3 次元格子点データ(経度・緯度・動径) ! xvz_ :: 3 次元分散格子点データ(経度・緯度・動径) ! wz_ :: 水平スペクトル, 動径格子点データ ! ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 (_wt, _xyz, _xvz, _wz_, _w, _xy, _xv, _x, _y, _z, _a) は, ! 入力変数の形がスペクトルデータおよび格子点データであることを示している. ! _wt :: スペクトルデータ ! _xyz :: 3 次元格子点データ ! _xyz_xyz :: 2 つの3 次元格子点データ, ... ! _xvz :: 3 次元分散格子点データ ! _xvz_xvz :: 2 つの3 次元分散格子点データ, ... ! !=== 各データの種類の説明 ! ! * xyh : 3 次元格子点データ(経度・緯度・分散動径) ! * 変数の種類と次元は real(8), dimension(0:im-1,jm,kc). ! * im, jm, km はそれぞれ経度, 緯度, 動径座標の格子点数であり, ! サブルーチン wt_mpi_Initial にてあらかじめ設定しておく. ! ! * ut : スペクトルデータ ! * 変数の種類と次元は real(8), dimension(nc,0:lm). ! * nm は球面調和函数の最大全波数, lm はチェビシェフ多項式の最大次数 ! であり, サブルーチン ut_mpi_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * nmz : 水平スペクトルデータの並んだ 3 次元配列. ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,0:km). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wt_Initial にて ! あらかじめ設定しておく. ! ! * nz : スペクトルデータの並んだ 2 次元配列. ! * 変数の種類と次元は real(8), dimension(0:nm,0:km). ! 第 1 次元が水平全波数を表す. ! * nm は球面調和函数の最大全波数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wz : 水平スペクトル, 動径格子点データ. ! * 変数の種類と次元は real(8), dimension(nc,0:km). ! ! * wt_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyz_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * xvz_ で始まる関数が返す値は 3 次元分散格子点データに同じ. ! ! * wz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! wt_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_Lon, x_Lon_Weight, :: 格子点座標(経度)と重みを格納した 1 次元配列 ! y_Lat, y_Lat_Weight :: 格子点座標(緯度)と重みを格納した 1 次元配列 ! h_Rad, h_Rad_Weight :: 格子点座標(同形)と重みを格納した 1 次元配列 ! xyz_Lon, xyz_Lat, xyz_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (格子点データ型 3 次元配列) ! xyh_Lon, xyh_Lat, xyh_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (分散格子点データ型 3 次元配列) ! !==== 基本変換 ! ! xyz_wt, wt_xyz :: スペクトルデータと 3 次元格子データの間の変換 (球面調和函数, チェビシェフ変換) ! ! xvz_wt, wt_xvz :: スペクトルデータと 3 次元分散格子データの間の変換 (球面調和函数, チェビシェフ変換) ! ! xyz_wz, wz_xyz :: 3 次元格子データと水平スペクトル・動径格子データとの間の変換 (球面調和函数) ! ! xvz_wz, wz_xvz :: 3 次元分散格子データと水平スペクトル・動径格子データとの間の変換 (球面調和函数) ! ! wz_wt, wt_wz :: スペクトルデータと水平スペクトル・動径格子データとの間の変換 (チェビシェフ変換) ! ! w_xy, xy_w :: スペクトルデータと 2 次元水平格子データの間の変換(球面調和函数変換) ! w_xv, xv_w :: スペクトルデータと 2 次元水平分散格子データの間の変換(球面調和函数変換) ! ! az_at, at_az :: 同時に複数個行う (チェビシェフ変換)格子データとチェビシェフデータの間の変換を ! ! 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φ・∂/∂λを作用させる ! xvz_GradLon_wt :: スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λを作用させる ! ! xyz_GradLat_wt :: スペクトルデータに勾配型緯度微分 1/r・∂/∂φを作用させる ! xvz_GradLat_wt :: スペクトルデータに勾配型緯度微分 1/r・∂/∂φを作用させる(分散格子) ! wt_DivLon_xyz :: 格子データに発散型経度微分 1/rcosφ・∂/∂λを作用させる ! wt_DivLon_xvz :: 分散格子データに発散型経度微分 1/rcosφ・∂/∂λを作用させるa ! wt_DivLat_xyz :: 格子データに発散型緯度微分 1/rcosφ・∂(g cosφ)/∂φを作用させる ! wt_DivLat_xvz :: 分散格子データに発散型緯度微分 1/rcosφ・∂(g cosφ)/∂φを作用させる(分散格子) ! ! wt_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに発散を作用させる ! wt_Div_xvz_xvz_xvz :: ベクトル成分である 3 つの分散格子データに発散を作用させる ! xyz_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに発散を作用させる ! xvz_Div_xvz_xvz_xvz :: ベクトル成分である 3 つの分散格子データに発散を作用させる ! xyz_RotLon_wt_wt :: ベクトル場の回転の経度成分を計算する ! xvz_RotLon_wt_wt :: ベクトル場の回転の経度成分を計算する(分散格子) ! xyz_RotLat_wt_wt :: ベクトル場の回転の緯度成分を計算する ! xvz_RotLat_wt_wt :: ベクトル場の回転の緯度成分を計算する(分散格子) ! wt_RotRad_xyz_xyz :: ベクトル場の回転の動径成分を計算する ! wt_RotRad_xvz_xvz :: ベクトル場の回転の動径成分を計算する(分散格子) ! !==== トロイダルポロイダル計算用微分 ! ! wt_KxRGrad_wt :: スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用させる ! ! xyz_KGrad_wt :: スペクトルデータに軸方向微分 k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r を作用させる ! xvz_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_RadRot_xvz_xvz :: ベクトル v の渦度と動径ベクトル r の内積 r・(▽×v) を計算する(分散格子) ! wt_RadRotRot_xyz_xyz_xyz :: ベクトルの v の r・(▽×▽×v) を計算する ! wt_RadRotRot_xvz_xvz_xvz :: ベクトルの v の r・(▽×▽×v) を計算する(分散格子) ! wt_Potential2Vector :: トロイダルポロイダルポテンシャルからベクトル場を計算する ! wt_Potential2VectorMPI :: トロイダルポロイダルポテンシャルからベクトル場を計算する(分散格子) ! ! wt_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される非発散ベクトル場の回転の各成分を計算する ! wt_Potential2RotationMPI :: トロイダルポロイダルポテンシャルで表される非発散ベクトル場の回転の各成分を計算する(分散格子) ! ! !==== 非線形計算 ! ! wt_VGradV :: ベクトル v から v・▽v を計算する ! wt_VGradVMPI :: ベクトル v から v・▽v を計算する(分散格子) ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! nmz_ToroidalEnergySpectrum_wt, :: トロイダルポテンシャルからエネルギーの ! nz_ToroidalEnergySpectrum_wt :: 球面調和函数各成分を計算する ! ! nmz_PoloidalEnergySpectrum_wt, :: ポロイダルポテンシャルからエネルギーの ! nz_PoloidalEnergySpectrum_wt :: 球面調和函数各成分を計算する ! !==== 境界値問題 ! ! wt_BoundariesTau, :: ディリクレ, ノイマン境界条件を適用する ! wt_BoundariesGrid, ::(タウ法, 選点法) ! wt_Boundaries :: ! ! wt_TorBoundariesTau, :: 速度トロイダルポテンシャルの境界条件を ! wt_TorBoundariesGrid, :: 適用する(タウ法,選点法) │ ! wt_TorBoundaries :: ! ! wz_LaplaPol2Pol_wz, :: 速度ポロイダルポテンシャルΦを▽^2Φから ! wt_LaplaPol2Pol_wt :: 求める (入出力がそれぞれチェビシェフ格子点, ! :: チェビシェフ係数) ! ! wt_TorMagBoundariesTau, :: 磁場トロイダルポテンシャルの境界条件を ! wt_TorMagBoundariesGrid, :: 適用する(タウ法, 選点法) ! wt_TorMagBoundaries :: ! ! wt_PolMagBoundariesTau, :: 磁場トロイダルポテンシャル境界の境界条件を ! wt_PolMagBoundariesGrid, :: 適用する(タウ法, 選点法) ! wt_PolMagBoundaries :: ! !==== 積分・平均(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_xvz, AvrLonLatRad_xvz :: 3 次元格子点データの全領域積分および平均 ! ! z_IntLonLat_xvz, z_AvrLonLat_xvz :: 3 次元格子点データの緯度経度(水平・球面)積分および平均 ! ! v_IntLonRad_xvz, v_AvrLonRad_xvz :: 3 次元格子点データの緯度動径積分および平均 ! ! z_IntLatRad_xvz, z_AvrLatRad_xvz :: 3 次元格子点データの経度動径(子午面)積分および平均 ! ! vz_IntLon_xvz, vz_AvrLon_xvz :: 3 次元格子点データの経度方向積分および平均 ! xz_IntLat_xvz, xz_AvrLat_xvz :: 3 次元格子点データの緯度方向積分および平均 ! xy_IntRad_xvz, xy_AvrRad_xvz :: 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)格子点データの動径方向積分および ! :: 平均 ! ! IntLonLat_xv, AvrLonLat_xv :: 2 次元格子点データの水平(球面)積分および平均 ! IntLonRad_xz, AvrLonRad_xz :: 2 次元(XZ)格子点データの経度動径積分 ! :: および平均 ! IntLatRad_vz, AvrLatRad_vz :: 2 次元(YZ)格子点データの緯度動径(子午面) ! :: 積分および平均 ! v_IntLon_xv, v_AvrLon_xv :: 水平 2 次元(球面)格子点データの経度方向 ! :: 積分および平均 ! v_IntLat_xv, x_AvrLat_xv :: 水平2 次元(球面)格子点データの緯度方向積分 ! :: および平均 ! z_IntLon_xz, z_AvrLon_xz :: 2 次元(XZ)格子点データの経度方向積分および ! :: 平均 ! x_IntRad_xz, x_AvrRad_xz :: 2 次元(XZ)格子点データの動径方向積分および ! :: 平均 ! z_IntLat_vz, z_AvrLat_vz :: 2 次元(YZ)格子点データの緯度方向積分および ! :: 平均 ! v_IntRad_vz, v_AvrRad_vz :: 2 次元(YZ)格子点データの動径方向積分および ! :: 平均 ! !==== 積分・平均(1 次元データ) ! ! IntLon_x, AvrLon_x :: 1 次元(X)格子点データの経度方向積分および平均 ! IntLat_y, AvrLat_y :: 1 次元(Y)格子点データの緯度方向積分および平均 ! IntRad_z, AvrRad_z :: 1 次元(Z)格子点データの動径方向積分および平均 ! ! IntLat_v, AvrLat_v :: 1 次元(Y)格子点データの緯度方向積分および平均 ! !==== 補間計算 ! ! Interpolate_ut :: スペクトルデータから任意の点の値を補間する. ! use dc_message use lumatrix use mpi use ua_base_mpi_module_svpack, only : ua_base_mpi_Initial, kc, k1, k2, nc, & l_nm, nm_l, kp_kl, xyb_ua, ua_xyb use w_module_svpack, only : w_Initial, & x_Lon, y_Lat, x_Lon_Weight, y_Lat_Weight, xy_Lon, xy_Lat, & IntLonLat_xy, y_IntLon_xy, x_IntLat_xy, IntLon_x, IntLat_y, & AvrLonLat_xy, y_AvrLon_xy, x_AvrLat_xy, AvrLon_x, AvrLat_y use w_interpolate_mpi_module_supack, only : Interpolate_u => Interpolate_w use ua_deriv_mpi_module_svpack, only : & ua_deriv_mpi_Initial, ua_deriv_mpi_Finalize, & ua_Lapla_ua, ua_LaplaInv_ua, ua_DLon_ua, & xyb_GradLon_ua, xyb_GradLat_ua, & ua_DivLon_xyb, ua_DivLat_xyb, & ua_Div_xyb_xyb !!$ use ua_interpolate_mpi_module_svpack !!$ use ua_spectrum_mpi_module_svpack !!$ use ua_mpi_module_svpack use at_module, z_RAD => g_X, z_RAD_WEIGHT => g_X_WEIGHT, & at_az => at_ag, az_at => ag_at, & t_z => t_g, z_t => g_t, & t_Dr_t => t_Dx_t, at_Dr_at => at_Dx_at implicit none private public ut_mpi_Initial public kc public nc public x_Lon, x_Lon_Weight public y_Lat, y_Lat_Weight public z_Rad, z_Rad_Weight public h_Rad, h_Rad_Weight public l_nm, nm_l public kh_kz public xy_Lon, xy_Lat public xyh_Lon, xyh_Lat, xyh_Rad public uz_Rad public ut_VMiss !!$ public u_xy, xy_u !!$ public at_Dr_at, t_Dr_t, az_at, at_az public uz_ut, ut_uz public xyh_ut, ut_xyh public xyh_uz, uz_xyh public ut_DRad_ut, ut_DivRad_ut, ut_RotRad_ut, ut_Lapla_ut public xyh_GradLon_ut, xyh_Gradlat_ut public ut_DivLon_xyh, ut_DivLat_xyh public ut_Div_xyh_xyh_xyh, xyh_Div_xyh_xyh_xyh public xyh_RotLon_ut_ut, xyh_RotLat_ut_ut, ut_RotRad_xyh_xyh public ut_KxRGrad_ut, xyh_KGrad_ut public ut_QOperatorMPI_ut public ut_L2_ut, ut_L2Inv_ut public ut_RadRot_xyh_xyh, ut_RadRotRot_xyh_xyh_xyh public ut_Potential2vectorMPI, ut_Potential2RotationMPI public yh_IntLon_xyh, xh_IntLat_xyh, xy_IntRad_xyh public x_IntLatRad_xyh, y_IntLonRad_xyh, h_IntLonLat_xyh public IntLonLatRad_xyh public x_IntLat_xy, y_IntLon_xy, IntLonLat_xy public h_IntLat_yh, y_IntRad_yh, IntLatRad_yh public h_IntLon_xh, x_IntRad_xh, IntLonRad_xh public IntLon_x, IntLat_y, IntRad_h public yh_AvrLon_xyh, xh_AvrLat_xyh, xy_AvrRad_xyh public x_AvrLatRad_xyh, y_AvrLonRad_xyh, h_AvrLonLat_xyh public AvrLonLatRad_xyh public x_AvrLat_xy, y_AvrLon_xy, AvrLonLat_xy public h_AvrLat_yh, y_AvrRad_yh, AvrLatRad_yh public h_AvrLon_xh, x_AvrRad_xh, AvrLonRad_xh public AvrLon_x, AvrLat_y, AvrRad_h !!$ public wt_VGradVMPI public Interpolate_ut public nmz_ToroidalEnergySpectrum_ut, nz_ToroidalEnergySpectrum_ut public nmz_PoloidalEnergySpectrum_ut, nz_PoloidalEnergySpectrum_ut public ut_Boundaries, ut_TorBoundaries public ut_LaplaPol2Pol_ut, uz_LaplaPol2Pol_uz public ut_TormagBoundaries, ut_PolmagBoundaries public ut_BoundariesTau, ut_TorBoundariesTau, ut_LaplaPol2PolTau_ut public ut_TormagBoundariesTau, ut_PolmagBoundariesTau public ut_BoundariesGrid, ut_TorBoundariesGrid, ut_LaplaPol2PolGrid_ut public ut_TormagBoundariesGrid, ut_PolmagBoundariesGrid interface ut_Boundaries module procedure ut_BoundariesTau end interface interface ut_TorBoundaries module procedure ut_TorBoundariesTau end interface interface ut_LaplaPol2Pol_ut module procedure ut_LaplaPol2PolTau_ut end interface interface ut_TorMagBoundaries module procedure ut_TorMagBoundariesTau end interface interface ut_PolMagBoundaries module procedure ut_PolMagBoundariesTau end interface integer :: im=64, jm=32, km=16 ! 格子点の設定(経度, 緯度, 動径) integer :: nm=21, lm=16 ! 切断波数の設定(水平, 動径) real(8) :: ri=0.0, ro=1.0 ! 球殻内外半径 real(8), parameter :: pi=3.1415926535897932385D0 real(8), dimension(:), allocatable :: h_Rad real(8), dimension(:), allocatable :: h_Rad_Weight real(8), dimension(:,:,:), allocatable :: xyh_LON, xyh_LAT, xyh_RAD ! 座標 real(8), dimension(:,:), allocatable :: uz_RAD ! 座標 real(8) :: ut_VMiss = -999.0 ! 欠損値 save im, jm, km, nm, lm, ri, ro, ut_VMiss contains !--------------- 初期化 ----------------- subroutine ut_mpi_Initial(i,j,k,n,l,r_in,r_out,omp,ua_init) ! ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! integer,intent(in) :: i ! 格子点数(経度λ) integer,intent(in) :: j ! 格子点数(緯度φ) integer,intent(in) :: k ! 格子点数(動径 r) integer,intent(in) :: n ! 切断波数(水平全波数) integer,intent(in) :: l ! 切断波数(動径波数) real(8),intent(in) :: r_in ! 球殻内半径 real(8),intent(in) :: r_out ! 球殻外半径 integer,intent(in), optional :: omp ! OPENMP での最大スレッド数 logical,intent(in), optional :: ua_init ! wa_initial スイッチ logical :: ua_initialize=.true. ! wa_initial スイッチ im = i ; jm = j ; km = k nm = n ; lm = l ri = r_in ; ro = r_out if ( present(ua_init) ) then ua_initialize = ua_init else ua_initialize = .true. endif if ( present(omp) ) then call MessageNotify('M','ut_mpi_Initial', & 'Optional argument "omp" is dummy in ut_mpi_module_svpack') endif if ( ua_initialize ) then ! call ua_mpi_Initial(nm,im,jm,km) call w_Initial(nm,im,jm) call ua_base_mpi_Initial(nm,im,jm,km+1) call ua_deriv_mpi_Initial endif call at_Initial(km,lm,r_in,r_out) allocate(h_Rad(kc)) allocate(h_Rad_Weight(kc)) allocate(xyh_Lon(0:im-1,jm,kc)) allocate(xyh_Lat(0:im-1,jm,kc)) allocate(xyh_Rad(0:im-1,jm,kc)) allocate(uz_Rad(nc,0:km)) uz_Rad = spread(z_Rad,1,nc) h_Rad = z_Rad(k1-1:k2-1) z_Rad_Weight = z_Rad_Weight * z_Rad**2 ! r^2 dr の積分重み h_Rad_Weight = z_Rad_Weight(k1-1:k2-1) xyh_Lon = spread(xy_Lon,3,kc) xyh_Lat = spread(xy_Lat,3,kc) xyh_Rad = spread(spread(z_Rad(k1-1:k2-1),1,jm),1,im) call MessageNotify('M','ut_mpi_initial', & 'ut_mpi_module_svpack (2016/08/04) is initialized') end subroutine ut_mpi_Initial function kh_kz(kz) integer, intent(IN) :: kz ! チェビシェフ格子点番号(全層) integer :: kh_kz ! 層番号(分割プロセス中) kh_kz = kp_kl(kz+1) end function kh_kz !--------------- 基本変換 ----------------- function xyh_ut(ut) ! ! スペクトルデータから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(0:im-1,jm,kc) :: xyh_ut !(out) 3 次元経度緯度動径格子点データ real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ xyh_ut = xyb_ua(uz_ut(ut)) end function xyh_ut function ut_xyh(xyh) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(nc,0:lm) :: ut_xyh !(out) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ ut_xyh = ut_uz(ua_xyb(xyh)) end function ut_xyh function xyh_uz(uz) ! ! 水平スペクトル・動径格子点データから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(0:im-1,jm,kc) :: xyh_uz !(out) 3 次元経度緯度動径格子点データ real(8), dimension(nc,0:km), intent(in) :: uz !(in) 2 次元球面調和函数スペクトル・動径格子点データ xyh_uz = xyb_ua(uz) end function xyh_uz function uz_xyh(xyh) ! ! 3 次元格子データから水平スペクトル・動径格子点データへ(正)変換する. ! real(8), dimension(nc,0:km) :: uz_xyh !(out) 2 次元球面調和函数スペクトル・動径格子点データ real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ uz_xyh = ua_xyb(xyh) end function uz_xyh function uz_ut(ut) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:km) :: uz_ut !(out) 2 次元球面調和函数スペクトル・動径格子点データ uz_ut = az_at(ut) end function uz_ut function ut_uz(uz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(nc,0:km), intent(in) :: uz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(8), dimension(nc,0:lm) :: ut_uz !(out) 2 次元球面調和函数チェビシェフスペクトルデータ ut_uz = at_az(uz) end function ut_uz !--------------- 微分計算 ----------------- function ut_DRad_ut(ut) ! ! 入力スペクトルデータに動径微分 ∂/∂r を作用する. ! ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_DRad_ut !(in) 動径微分された2 次元球面調和函数チェビシェフスペクトルデータ ut_DRad_ut = at_Dr_at(ut) end function ut_DRad_ut function ut_DivRad_ut(ut) ! ! 入力スペクトルデータに発散型動径微分 ! ! 1/r^2 ∂/∂r (r^2 .)= ∂/∂r + 2/r ! ! を作用する. ! ! スペクトルデータの発散型動径微分とは, 対応する格子点データに ! 発散型動径微分を作用させたデータのスペクトル変換のことである. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_DivRad_ut !(out) 発散型動径微分を作用された 2 次元スペクトルデータ ut_DivRad_ut = ut_Drad_ut(ut) + ut_uz(2/uz_rad*uz_ut(ut)) end function ut_DivRad_ut function ut_RotRad_ut(ut) ! ! 入力スペクトルデータに回転型動径微分 ! ! 1/r ∂(r.)/∂r = ∂(.)/∂r + (.)/r ! ! を作用する. ! ! スペクトルデータの回転型動径微分とは, 対応する格子点データに ! 回転型動径微分を作用させたデータのスペクトル変換のことである. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_RotRad_ut !(out) 回転型動径微分を作用された 2 次元スペクトルデータ ut_RotRad_ut = ut_Drad_ut(ut) + ut_uz(1/uz_Rad*uz_ut(ut)) end function ut_RotRad_ut function ut_Lapla_ut(ut) ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = 1/r^2 cos^2φ・∂^2/∂λ^2 ! + 1/r^2 cosφ・∂/∂φ(cosφ∂/∂φ) ! + 1/r^2 ∂/∂r (r^2 ∂/∂r) ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_Lapla_ut !(out) ラプラシアンを作用された 2 次元スペクトルデータ ut_Lapla_ut = ut_DivRad_ut(ut_Drad_ut(ut)) & + ut_uz(uz_ut(ua_Lapla_ua(ut))/uz_Rad**2) end function ut_Lapla_ut function xyh_GradLon_ut(ut) ! ! スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λ ! を作用させる. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,jm,kc) :: xyh_GradLon_ut !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ xyh_GradLon_ut = xyb_GradLon_ua(uz_ut(ut))/xyh_Rad end function xyh_GradLon_ut function xyh_GradLat_ut(ut) ! ! スペクトルデータに勾配型経度微分 1/r ∂/∂φ を作用させる. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,jm,kc) :: xyh_GradLat_ut !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ xyh_GradLat_ut = xyb_GradLat_ua(uz_ut(ut))/xyh_Rad end function xyh_GradLat_ut function ut_DivLon_xyh(xyh) ! ! 格子点データに発散型経度微分 1/rcosφ・∂/∂λ を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(nc,0:lm) :: ut_DivLon_xyh !(out) 発散型経度微分を作用された 2 次元スペクトルデータ ut_DivLon_xyh = ut_uz(ua_DivLon_xyb(xyh/xyh_Rad)) end function ut_DivLon_xyh function ut_DivLat_xyh(xyh) ! ! 格子データに発散型緯度微分 1/rcosφ・∂(f cosφ)/∂φ を ! 作用させたスペクトルデータを返す. ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(nc,0:lm) :: ut_DivLat_xyh !(out) 発散型緯度微分を作用された 2 次元スペクトルデータ ut_DivLat_xyh = ut_uz(ua_divlat_xyb(xyh/xyh_Rad)) end function ut_DivLat_xyh function ut_Div_xyh_xyh_xyh(xyh_Vlon,xyh_Vlat,xyh_Vrad) ! ! べクトル成分である 3 つの格子データに発散を作用させた ! スペクトルデータを返す. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表し, 発散は ! ! 1/rcosφ・∂u/∂λ + 1/rcosφ・∂(v cosφ)/∂φ ! + 1/r^2 ∂/∂r (r^2 w) ! ! と計算される. ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vrad !(in) ベクトル場の動径成分 real(8), dimension(nc,0:lm) :: ut_Div_xyh_xyh_xyh !(out) ベクトル場の発散 ut_Div_xyh_xyh_xyh = ut_DivLon_xyh(xyh_Vlon) & + ut_DivLat_xyh(xyh_Vlat) & + ut_DivRad_ut(ut_xyh(xyh_Vrad)) end function ut_Div_xyh_xyh_xyh function xyh_Div_xyh_xyh_xyh(xyh_Vlon,xyh_Vlat,xyh_Vrad) ! ! ベクトル成分である 3 つの格子データに発散を作用させる. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表す. ! ! 極の特異性を回避するためにベクトル場に cosφ/r の重みをかけて ! 計算している. ! ! div V = (r/cosφ)・div (Vcosφ/r) + V_φtanφ/r + V_r/r ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vrad !(in) ベクトル場の動径成分 real(8), dimension(0:im-1,jm,kc) :: xyh_Div_xyh_xyh_xyh !(out) ベクトル場の発散 xyh_Div_xyh_xyh_xyh & = xyh_Rad/cos(xyh_Lat) & * xyh_ut(ut_Div_xyh_xyh_xyh(xyh_VLon*cos(xyh_Lat)/xyh_Rad, & xyh_VLat*cos(xyh_Lat)/xyh_Rad, & xyh_VRad*cos(xyh_Lat)/xyh_Rad ))& + xyh_VLat*tan(xyh_Lat)/xyh_Rad & + xyh_VRad/xyh_Rad end function xyh_Div_xyh_xyh_xyh function xyh_RotLon_ut_ut(ut_Vrad,ut_Vlat) ! ! ベクトル場の動径成分, 緯度成分である第 1, 2 引数 Vrad, Vlat から ! 回転の経度成分 ! ! 1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r を計算する. ! ! を計算する ! real(8), dimension(nc,0:lm), intent(in) :: ut_Vrad !(in) ベクトル場の動径成分 real(8), dimension(nc,0:lm), intent(in) :: ut_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,jm,kc) :: xyh_RotLon_ut_ut !(out) ベクトル場の回転の経度成分 xyh_RotLon_ut_ut = xyh_GradLat_ut(ut_Vrad) & - xyh_ut(ut_RotRad_ut(ut_Vlat)) end function xyh_RotLon_ut_ut function xyh_RotLat_ut_ut(ut_Vlon,ut_Vrad) ! ! ベクトル場の経度成分, 動径成分である第 1, 2 引数 Vlon, Vrad から ! 回転の緯度成分 ! ! 1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ ! ! を計算する. ! real(8), dimension(nc,0:lm), intent(in) :: ut_Vlon !(in) ベクトル場の経度成分 real(8), dimension(nc,0:lm), intent(in) :: ut_Vrad !(in) ベクトル場の動径成分 real(8), dimension(0:im-1,jm,kc) :: xyh_RotLat_ut_ut !(out) ベクトル場の回転の緯度成分 xyh_RotLat_ut_ut = xyh_ut(ut_RotRad_ut(ut_Vlon)) & - xyh_GradLon_ut(ut_Vrad) end function xyh_RotLat_ut_ut function ut_RotRad_xyh_xyh(xyh_Vlat,xyh_Vlon) ! ! ベクトルの緯度成分, 経度成分である第 1, 2 引数 Vlat, Vlon に対して ! ベクトル場の回転の動径成分 ! ! 1/rcosφ・∂Vlat/∂λ - 1/rcosφ・∂(Vlon cosφ)/∂φ ! ! を計算する. ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_Vlon !(in) ベクトル場の経度成分 real(8), dimension(nc,0:lm) :: ut_RotRad_xyh_xyh !(out) ベクトル場の回転の動径成分 ut_RotRad_xyh_xyh = ut_DivLon_xyh(xyh_Vlat) & - ut_DivLat_xyh(xyh_Vlon) end function ut_RotRad_xyh_xyh !--------------- ポロイダル/トロイダルモデル用微分 ----------------- function ut_KxRGrad_ut(ut) ! ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_KxRGrad_ut !(out) 経度微分を作用された 2 次元スペクトルデータ ut_KxRGrad_ut = ua_Dlon_ua(ut) end function ut_KxRGrad_ut function ut_L2_ut(ut) ! ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する. ! ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. ! 入力スペクトルデ ータに対応する格子点データに演算子 ! ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用させたデータのスペクトル変換が返される. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_L2_ut !(out) L^2 演算子を作用された 2 次元スペクトルデータ ut_L2_ut = -ua_Lapla_ua(ut) end function ut_L2_ut function ut_L2Inv_ut(ut) ! ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を ! 作用する. ! ! スペクトルデータに L^2 演算子を作用させる関数 ut_L2_ut の逆計算を ! 行う関数である. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_L2Inv_ut !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ ut_L2Inv_ut = -ua_LaplaInv_ua(ut) end function ut_L2Inv_ut function xyh_KGrad_ut(ut) ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r ! ! 入力スペクトルデータに対応する格子データに軸方向微分 ! ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r ! ! を作用させた格子データが返される. ! ここでベクトル k は球の中心から北極向きの単位ベクトルである. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,jm,kc) :: xyh_KGrad_ut !(out) 軸方向微分を作用された 2 次元スペクトルデータ xyh_KGrad_ut = cos(xyh_Lat)*xyh_GradLat_ut(ut) & + sin(xyh_Lat)*xyh_ut(ut_Drad_ut(ut)) end function xyh_KGrad_ut function ut_QOperatorMPI_ut(ut) ! ! 入力スペクトルデータに対応する格子点データに演算子 ! ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) ! ! を作用させたデータのスペクトル変換が返される. ! real(8), dimension(nc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nc,0:lm) :: ut_QOperatorMPI_ut !(out) Q 演算子を作用された 2 次元スペクトルデータ ut_QOperatorMPI_ut = & ut_xyh(xyh_KGrad_ut(ut) - xyh_KGrad_ut(ut_L2_ut(ut))/2) & - ut_L2_ut(ut_xyh(xyh_KGrad_ut(ut)))/2 end function ut_QOperatorMPI_ut function ut_RadRot_xyh_xyh(xyh_VLON,xyh_VLAT) ! r・(▽×v) ! ! ベクトルの渦度と動径ベクトルの内積 r・(▽×v) を計算する. ! ! 第 1, 2 引数(v[λ], v[φ])がそれぞれベクトルの経度成分, 緯度成分を表す. ! ! r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ ! ! のスペクトル データが返される. ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_VLAT !(in) ベクトルの緯度成分 real(8), dimension(nc,0:lm) :: ut_RadRot_xyh_xyh !(out) ベクトルの渦度と動径ベクトルの内積 ut_RadRot_xyh_xyh = ut_uz(ua_DivLon_xyb(xyh_VLAT) & - ua_DivLat_xyb(xyh_VLON)) end function ut_RadRot_xyh_xyh function ut_RadRotRot_xyh_xyh_xyh(xyh_VLON,xyh_VLAT,xyh_VRAD) ! ! ベクトル v に対して r・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分, ! 緯度成分, 動径成分を表す. ! ! r・(▽×▽×v) = 1/r ∂/∂r (r・( 1/cosφ・∂v[λ]/∂λ ! + 1/cosφ・∂(v[φ] cosφ)/∂φ ) ) ! + L^2 v[r]/r ! ! のスペクトルデータが返される. ! real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_VLAT !(in) ベクトルの緯度成分 real(8), dimension(0:im-1,jm,kc), intent(in) :: xyh_VRAD !(in) ベクトルの動径成分 real(8), dimension(nc,0:lm) :: ut_RadRotRot_xyh_xyh_xyh !(out) ベクトル v の r・(▽×▽×v) ut_RadRotRot_xyh_xyh_xyh = & ut_RotRad_ut(ut_uz( & (ua_DivLon_xyb(xyh_VLON)+ ua_DivLat_xyb(xyh_VLAT)))) & + ut_L2_ut(ut_xyh(xyh_VRAD/xyh_RAD)) end function ut_RadRotRot_xyh_xyh_xyh subroutine ut_Potential2VectorMPI(& xyh_VLON,xyh_VLAT,xyh_VRAD,ut_TORPOT,ut_POLPOT) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψr) + ▽x▽x(Φr) ! ! の各成分を計算する ! real(8), dimension(0:im-1,jm,kc) :: xyh_VLON !(out) ベクトル場の経度成分 real(8), dimension(0:im-1,jm,kc) :: xyh_VLAT !(out) ベクトル場の緯度成分 real(8), dimension(0:im-1,jm,kc) :: xyh_VRAD !(out) ベクトル場の動径成分 real(8), dimension(nc,0:lm), intent(in) :: ut_TORPOT !(in) トロイダルポテンシャル real(8), dimension(nc,0:lm), intent(in) :: ut_POLPOT !(in) ポロイダルポテンシャル xyh_VLON = xyh_RAD * xyh_GradLat_ut(ut_TORPOT) & + xyb_GradLon_ua(uz_ut(ut_RotRad_ut(ut_POLPOT))) xyh_VLAT = - xyh_RAD * xyh_GradLon_ut(ut_TORPOT) & + xyb_GradLat_ua(uz_ut(ut_RotRad_ut(ut_POLPOT))) xyh_VRAD = xyh_ut(ut_L2_ut(ut_POLPOT))/xyh_RAD end subroutine ut_Potential2VectorMPI subroutine ut_Potential2RotationMPI(& xyh_RotVLON,xyh_RotVLAT,xyh_RotVRAD,ut_TORPOT,ut_POLPOT) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψr) + ▽x▽x(Φr) ! ! に対して, その回転 ! ! ▽xv = ▽x▽x(Ψr) + ▽x▽x▽x(Φr) = ▽x▽x(Ψr) - ▽x((▽^2Φ)r) ! ! を計算する. ! ベクトル場の回転 real(8), dimension(0:im-1,jm,kc), intent(OUT) :: xyh_RotVLON !(out) 回転の経度成分 real(8), dimension(0:im-1,jm,kc), intent(OUT) :: xyh_RotVLAT !(out) 回転の緯度成分 real(8), dimension(0:im-1,jm,kc), intent(OUT) :: xyh_RotVRAD !(out) 回転の動径成分 ! 入力ベクトル場を表すポテンシャル real(8), dimension(nc,0:lm), intent(in) :: ut_TORPOT !(in) トロイダルポテンシャル real(8), dimension(nc,0:lm), intent(in) :: ut_POLPOT !(in) ポロイダルポテンシャル call ut_Potential2VectorMPI( & xyh_RotVLON,xyh_RotVLAT,xyh_RotVRAD, & -ut_Lapla_ut(ut_POLPOT), ut_TORPOT) end subroutine ut_Potential2RotationMPI !!$ !------------------- 非線形項計算 ---------------------- !!$ subroutine wt_VGradVMPI(xvz_VGRADV_LON,xvz_VGRADV_LAT,xvz_VGRADV_RAD, & !!$ xvz_VLON,xvz_VLAT,xvz_VRAD ) !!$ ! !!$ ! ベクトル場の v・▽v を計算する. !!$ ! !!$ ! ベクトル場 v=(v[λ],v[φ],v[r]) に対するv・▽vの各成分は !!$ ! 次のように計算される. !!$ ! !!$ ! (v・▽v)[λ] = ▽・(v[λ]v) + v[λ]v[r]/r - v[λ]v[φ]tan(φ)/r !!$ ! (v・▽v)[φ] = ▽・(v[φ]v) + v[φ]v[r]/r - v[λ]^2tan(φ)/r !!$ ! (v・▽v)[r] = ▽・(v[r]v) + (v[λ]^2+v[φ]^2)/r !!$ ! !!$ ! 非発散速度場に対してはポテンシャルから wt_Potential2Rotation を !!$ ! 用いて回転を計算し, 恒等式 v・▽v = ▽(v[2^/2) - vx▽xv を !!$ ! 用いる方がよいだろう. !!$ ! !!$ real(8), dimension(0:im-1,jc,0:km),intent(out) :: xvz_VGRADV_LON !!$ !(out) (v・▽v) 経度成分 !!$ !!$ real(8), dimension(0:im-1,jc,0:km),intent(out) :: xvz_VGRADV_LAT !!$ !(out) (v・▽v) 緯度成分 !!$ !!$ real(8), dimension(0:im-1,jc,0:km),intent(out) :: xvz_VGRADV_RAD !!$ !(out) (v・▽v) 動径成分 !!$ !!$ real(8), dimension(0:im-1,jc,0:km),intent(in) :: xvz_VLON !!$ !(in) ベクトル場 v の経度成分 !!$ !!$ real(8), dimension(0:im-1,jc,0:km),intent(in) :: xvz_VLAT !!$ !(in) ベクトル場 v の緯度成分 !!$ !!$ real(8), dimension(0:im-1,jc,0:km),intent(in) :: xvz_VRAD !!$ !(in) ベクトル場 v の動径成分 !!$ !!$ xvz_VGRADV_LON = & !!$ xvz_Div_xvz_xvz_xvz( & !!$ xvz_VLON * xvz_VLON, xvz_VLON*xvz_VLAT, xvz_VLON*xvz_VRAD ) & !!$ + xvz_VLON*xvz_VRAD/xvz_RAD & !!$ - xvz_VLON*xvz_VLAT*tan(xvz_LAT)/xvz_RAD !!$ !!$ xvz_VGRADV_LAT = & !!$ xvz_Div_xvz_xvz_xvz( & !!$ xvz_VLAT*xvz_VLON, xvz_VLAT*xvz_VLAT, xvz_VLAT*xvz_VRAD ) & !!$ + xvz_VLAT*xvz_VRAD/xvz_RAD & !!$ + xvz_VLON**2*tan(xvz_LAT)/xvz_RAD !!$ !!$ xvz_VGRADV_RAD = & !!$ xvz_Div_xvz_xvz_xvz( & !!$ xvz_VRAD*xvz_VLON, xvz_VRAD*xvz_VLAT, xvz_VRAD*xvz_VRAD ) & !!$ - (xvz_VLON**2 + xvz_VLAT**2)/xvz_RAD !!$ !!$ end subroutine wt_VGradVMPI !--------------- 積分計算 ----------------- !----(入力データ xyh)--- function yh_IntLon_xyh(xyh) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)積分. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm,kc) :: yh_IntLon_xyh !(out) 経度方向(帯状)積分された 2 次元子午面格子点データ integer :: i yh_IntLon_xyh = 0.0d0 do i=0,im-1 yh_IntLon_xyh(:,:) = yh_IntLon_xyh(:,:) & + xyh(i,:,:) * x_Lon_Weight(i) enddo end function yh_IntLon_xyh function xh_IntLat_xyh(xyh) ! ! 3 次元格子点データの緯度方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) cosφ dφ を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,kc) :: xh_IntLat_xyh !(out) 緯度積分された 2 次元緯度動径格子点データ. ! 緯度円格子点データ integer :: j xh_IntLat_xyh = 0.0d0 do j=1,jm xh_IntLat_xyh(:,:) = xh_IntLat_xyh(:,:) & + xyh(:,j,:) * y_Lat_Weight(j) enddo end function xh_IntLat_xyh function xy_IntRad_xyh(xyh) ! 動径積分 ! ! 3 次元格子点データの動径方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dr を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,1:jm) :: xy_IntRad_xyh !(out) 動径積分された 2 次元経度緯度(水平, 球面)格子点データ real(8), dimension(0:im-1,1:jm) :: xy_IntRadTmp integer :: k, ierr xy_IntRadTmp = 0.0d0 do k=1,kc xy_IntRadTmp(:,:) = xy_IntRadTmp(:,:) & + xyh(:,:,k) * h_Rad_Weight(k) enddo CALL MPI_ALLREDUCE(xy_IntRadTMP,xy_IntRad_xyh,im*jm,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function xy_IntRad_xyh function x_IntLatRad_xyh(xyh) ! ! 3 次元格子点データの緯度動径(子午面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dφdr ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1) :: x_IntLatRad_xyh !(out) 緯度動径(子午面)積分された 1 次元経度格子点データ real(8), dimension(0:im-1) :: x_IntLatRadTmp integer :: j, k, ierr x_IntLatRadTmp = 0.0D0 do k=1,kc do j=1,jm x_IntLatRadTmp = x_IntLatRadTmp & + xyh(:,j,k) * y_Lat_Weight(j) * h_Rad_Weight(k) enddo enddo CALL MPI_ALLREDUCE(x_IntLatRadTMP,x_IntLatRad_xyh,im,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function x_IntLatRad_xyh function y_IntLonRad_xyh(xyh) ! ! 3 次元格子点データの経度動径(緯度円)積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dλdr を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm) :: y_IntLonRad_xyh !(out) 経度動径(緯度円)積分された 1 次元緯度格子点データ real(8), dimension(1:jm) :: y_IntLonRadTmp integer :: i, k, ierr y_IntLonRadTmp = 0 do k=1,kc do i=0,im-1 y_IntLonRadTmp = y_IntLonRadTmp & + xyh(i,:,k) * x_Lon_Weight(i) * h_Rad_Weight(k) enddo enddo CALL MPI_ALLREDUCE(y_IntLonRadTMP,y_IntLonRad_xyh,jm,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function y_IntLonRad_xyh function h_IntLonLat_xyh(xyh) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(kc) :: h_IntLonLat_xyh !(out) 緯度経度(水平, 球面)積分された 1 次元動径格子点データ integer :: i, j h_IntLonLat_xyh = 0 do j=1,jm do i=0,im-1 h_IntLonLat_xyh = h_IntLonLat_xyh & + xyh(i,j,:) * x_Lon_Weight(i) * y_Lat_Weight(j) enddo enddo end function h_IntLonLat_xyh function IntLonLatRad_xyh(xyh) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8) :: IntLonLatRad_xyh !(out) 全球積分値 real(8) :: IntLonLatRadTmp integer :: i, j, k, ierr IntLonLatRadTmp = 0 do k=1,kc do j=1,jm do i=0,im-1 IntLonLatRadTmp = IntLonLatRadTmp & + xyh(i,j,k) * x_Lon_Weight(i) & * y_Lat_Weight(j) * h_Rad_Weight(k) enddo enddo enddo CALL MPI_ALLREDUCE(IntLonLatRadTMP,IntLonLatRad_xyh,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntLonLatRad_xyh !----(入力データ yh)--- function h_IntLat_yh(yh) ! 緯度積分 ! ! 2 次元(YH)格子点データの緯度方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) cosφ dφ を計算する. ! real(8), dimension(jm,kc), intent(in) :: yh !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(kc) :: h_IntLat_yh !(out) 緯度積分された 1 次元動径格子点データ integer :: j h_IntLat_yh = 0.0d0 do j=1,jm h_IntLat_yh(:) = h_IntLat_yh(:) + yh(j,:) * y_Lat_Weight(j) enddo end function h_IntLat_yh function y_IntRad_yh(yh) ! 動径積分 ! ! 2 次元(YH)格子点データの動径方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) r^2dr を計算する. ! real(8), dimension(1:jm,kc), intent(in) :: yh !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(1:jm) :: y_IntRad_yh !(out) 動径積分された 1 次元緯度格子点データ real(8), dimension(1:jm) :: y_IntRadTmp integer :: k, ierr y_IntRadTmp = 0.0d0 do k=1,kc y_IntRadTmp(:) = y_IntRadTmp(:) & + yh(:,k) * h_Rad_Weight(k) enddo CALL MPI_ALLREDUCE(y_IntRadTMP,y_IntRad_yh,jm,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function y_IntRad_yh function IntLatRad_yh(yh) ! ! 2 次元(YH)格子点データの緯度動径積分(子午面)および平均 ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2cosφ dφdr を計算する. ! real(8), dimension(1:jm,kc), intent(in) :: yh !(in) 2 次元緯度動径(子午面)格子点データ real(8) :: IntLatRad_yh !(out) 積分値 real(8) :: IntLatRadTmp integer :: j, k, ierr IntLatRadTmp = 0 do k=1,kc do j=1,jm IntLatRadTmp = IntLatRadTmp & + yh(j,k) * y_Lat_Weight(j) * h_Rad_Weight(k) enddo enddo CALL MPI_ALLREDUCE(IntLatRadTMP,IntLatRad_yh,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntLatRad_yh !----(入力データ xh)--- function h_IntLon_xh(xh) ! ! 2 次元(XH)格子点データの経度方向積分. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ を計算する. ! real(8), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(8), dimension(kc) :: h_IntLon_xh !(out) 経度積分された 1 次元動径格子点データ integer :: i h_IntLon_xh = 0.0d0 do i=0,im-1 h_IntLon_xh(:) = h_IntLon_xh(:) + xh(i,:) * x_Lon_Weight(i) enddo end function h_IntLon_xh function x_IntRad_xh(xh) ! ! 2 次元(XH)格子点データの動径方向域積分. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r) r^2dr を計算する. ! real(8), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(8), dimension(0:im-1) :: x_IntRad_xh !(out) 動径積分された 1 次元経度格子点データ real(8), dimension(0:im-1) :: x_IntRadTmp integer :: k, ierr x_IntRadTmp = 0.0d0 do k=1,kc x_IntRadTmp(:) = x_IntRadTmp(:) & + xh(:,k) * h_Rad_Weight(k) enddo CALL MPI_ALLREDUCE(x_IntRadTMP,x_IntRad_xh,im,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function x_IntRad_xh function IntLonRad_xh(xh) ! 経度動径(緯度円)積分 ! ! 2 次元(XH)格子点データの経度動径積分 ! ! 2 次元データ f(λ,r) に対して∫f(λ,r) r^2dλdr を計算する. ! real(8), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(8) :: IntLonRad_xh !(out) 積分値 real(8) :: IntLonRadTmp integer :: i, k, ierr IntLonRadTmp = 0 do k=1,kc do i=0,im-1 IntLonRadTmp = IntLonRadTmp & + xh(i,k) * x_Lon_Weight(i) * h_Rad_Weight(k) enddo enddo CALL MPI_ALLREDUCE(IntLonRadTMP,IntLonRad_xh,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntLonRad_xh !----(入力データ h)--- function IntRad_h(h) ! 動径積分 ! ! 1 次元(Z)格子点データの動径方向域積分. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr を計算する. ! real(8), dimension(kc), intent(in) :: h !(in) 1 次元動径格子点データ real(8) :: IntRad_h !(out) 積分値 real(8) :: IntRadTMP integer :: k, ierr IntRadTmp = 0.0d0 do k=1,kc IntRadTmp = IntRadTmp + h(k) * h_Rad_Weight(k) enddo CALL MPI_ALLREDUCE(IntRadTMP,IntRad_h,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntRad_h !--------------- 平均計算 ----------------- !----(入力データ xyh)--- function yh_AvrLon_xyh(xyh) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ/2π を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm,kc) :: yh_AvrLon_xyh !(out) 経度方向(帯状)平均された 2 次元子午面格子点データ yh_AvrLon_xyh = yh_IntLon_xyh(xyh)/sum(x_Lon_Weight) end function yh_AvrLon_xyh function xh_AvrLat_xyh(xyh) ! 緯度積分 ! ! 3 次元格子点データの緯度方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)cosφ dφ/2 を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,kc) :: xh_AvrLat_xyh !(out) 緯度平均された 2 次元緯度動径格子点データ xh_AvrLat_xyh = xh_IntLat_xyh(xyh)/sum(y_Lat_Weight) end function xh_AvrLat_xyh function xy_AvrRad_xyh(xyh) ! ! 3 次元格子点データの動径方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dr/((r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,1:jm) :: xy_AvrRad_xyh !(out) 動径平均された 2 次元経度緯度(水平, 球面)格子点データ ! 水平格子点データ xy_AvrRad_xyh = xy_IntRad_xyh(xyh)/sum(z_Rad_Weight) end function xy_AvrRad_xyh function x_AvrLatRad_xyh(xyh) ! 緯度動径(子午面)積分 ! ! 3 次元格子点データの緯度動径(子午面)平均 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1) :: x_AvrLatRad_xyh !(out) 緯度動径(子午面)平均された 1 次元経度格子点データ x_AvrLatRad_xyh = x_IntLatRad_xyh(xyh) & /( sum(y_Lat_Weight)*sum(z_Rad_Weight) ) end function x_AvrLatRad_xyh function y_AvrLonRad_xyh(xyh) ! 経度動径(緯度円)積分 ! ! 3 次元格子点データの経度動径(緯度円)平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm) :: y_AvrLonRad_xyh !(out) 経度動径(緯度円)平均された 1 次元緯度格子点データ y_AvrLonRad_xyh = y_IntLonRad_xyh(xyh) & /(sum(x_Lon_Weight)*sum(z_Rad_Weight)) end function y_AvrLonRad_xyh function h_AvrLonLat_xyh(xyh) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ /4π ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8), dimension(kc) :: h_AvrLonLat_xyh !(out) 緯度経度(水平, 球面)平均された 1 次元動径格子点データ h_AvrLonLat_xyh = h_IntLonLat_xyh(xyh) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight)) end function h_AvrLonLat_xyh function AvrLonLatRad_xyh(xyh) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr /(4π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,kc), intent(in) :: xyh !(in) 3 次元経度緯度動径格子点データ real(8) :: AvrLonLatRad_xyh !(out) 全球平均値 AvrLonLatRad_xyh = IntLonLatRad_xyh(xyh) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(z_Rad_Weight)) end function AvrLonLatRad_xyh !----(入力データ yh)--- function h_AvrLat_yh(yh) ! ! 2 次元(YH)格子点データの緯度方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) cosφ dφ/2 を計算する. ! real(8), dimension(1:jm,kc), intent(in) :: yh !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(kc) :: h_AvrLat_yh !(out) 緯度平均された 1 次元動径格子点データ h_AvrLat_yh = h_IntLat_yh(yh)/sum(y_Lat_Weight) end function h_AvrLat_yh function y_AvrRad_yh(yh) ! ! 2 次元(YH)格子点データの動径方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2dr /((r[o]^3-r[i]^3)/3) ! を計算する. ! real(8), dimension(1:jm,kc), intent(in) :: yh !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(1:jm) :: y_AvrRad_yh !(out) 動径平均された 1 次元緯度格子点データ y_AvrRad_yh = y_IntRad_yh(yh)/sum(z_Rad_Weight) end function y_AvrRad_yh function AvrLatRad_yh(yh) ! 緯度動径(子午面)積分 ! ! 2 次元(YH)格子点データの緯度動径(子午面)平均 ! ! 2 次元データ f(φ,r) に対して ! ! ∫f(φ,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(1:jm,kc), intent(in) :: yh !(in) 2 次元緯度動径(子午面)格子点データ real(8) :: AvrLatRad_yh !(out) 平均値 AvrLatRad_yh = IntLatRad_yh(yh)/(sum(y_Lat_Weight)*sum(z_Rad_Weight)) end function AvrLatRad_yh !----(入力データ xh)--- function h_AvrLon_xh(xh) ! 経度(帯状)積分 ! ! 2 次元(XH)格子点データの経度方向平均. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ/2π を計算する. ! real(8), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(8), dimension(kc) :: h_AvrLon_xh !(out) 経度平均された 1 次元動径格子点データ h_AvrLon_xh = h_IntLon_xh(xh)/sum(x_Lon_Weight) end function h_AvrLon_xh function x_AvrRad_xh(xh) ! 動径積分 ! ! 2 次元(XH)格子点データの動径方向域平均. ! ! 2 次元データ f(λ,r) に対して ! ! ∫f(λ,r) r^2dr /((r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(8), dimension(0:im-1) :: x_AvrRad_xh !(out) 動径平均された 1 次元経度格子点データ x_AvrRad_xh = x_IntRad_xh(xh)/sum(z_Rad_Weight) end function x_AvrRad_xh function AvrLonRad_xh(xh) ! 経度動径(緯度円)積分 ! ! 2 次元(XH)格子点データの経度動径平均 ! ! 2 次元データ f(λ,r) に対して ! ! ∫f(λ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,kc), intent(in) :: xh ! (in) 2 次元格子点データ real(8) :: AvrLonRad_xh ! 積分値 AvrLonRad_xh = IntLonRad_xh(xh)/(sum(x_Lon_Weight)*sum(z_Rad_Weight)) end function AvrLonRad_xh !----(入力データ z)--- function AvrRad_h(z) ! ! 1 次元(Z)格子点データの動径方向域平均. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr /((r[o]^3-r[i]^3)/3) を ! 計算する. ! real(8), dimension(kc), intent(in) :: z !(in) 1 次元動径格子点データ real(8) :: AvrRad_h !(out) 平均値 AvrRad_h = IntRad_h(z)/sum(z_Rad_Weight) end function AvrRad_h !--------------- 補間計算 ----------------- function Interpolate_ut(ut_data,alon,alat,arad) ! ! 緯度 alon, 経度 alat 動径 arad における関数値を ! その球面調和変換係数 wa_data から補間計算する ! real(8), intent(IN) :: ut_data(nc,0:lm) ! スペクトルデータ real(8), intent(IN) :: alon ! 補間する位置(経度) real(8), intent(IN) :: alat ! 補間する位置(緯度) real(8), intent(IN) :: arad ! 補間する位置(動径) real(8) :: Interpolate_ut ! 補間した値 Interpolate_ut = & Interpolate_u(a_Interpolate_at(ut_data,arad),alon,alat) end function Interpolate_ut !--------------- ポロイダル/トロイダルモデル用スペクトル解析 ---------------- function nmz_ToroidalEnergySpectrum_ut(ut_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する ! ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(n,m,r)から全波数 n, 帯状波数 m 成分のトロイダルエネルギー ! スペクトルは (1/2)n(n+1)4πr^2|ψ(n,m,r)|^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和を動径積分したもの(r^2の重み無し) ! が球殻内での全エネルギーに等しい. ! ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される. ! ut_VMiss によって設定できる (初期値は -999.0) ! real(8), intent(in) :: ut_TORPOT(nc,0:lm) !(in) 流線関数(スペクトルデータ) real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_ToroidalEnergySpectrum_ut !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間) real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_Work ! 作業領域 real(8), dimension(nc,0:km) ::uz_DATA ! 作業領域 integer :: l, n, m, nm_ary(2) integer :: ierr uz_DATA = uz_ut(ut_TORPOT) nmz_work = 0.0D0 do l=1,nc nm_ary = nm_l(l) n = nm_ary(1) ; m = abs(nm_ary(2)) if ( n .ge. 0 ) then if ( m == 0 ) then nmz_Work(n,0,:) = nmz_Work(n,0,:) & + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * uz_Data(l,:)**2 else nmz_Work(n,m,:) = nmz_Work(n,m,:) & + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * uz_Data(l,:)**2 nmz_Work(n,-m,:) = nmz_Work(n,-m,:) & + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * uz_Data(l,:)**2 endif endif enddo CALL MPI_ALLREDUCE(nmz_Work,nmz_ToroidalEnergySpectrum_ut, & (nm+1)*(2*nm+1)*(km+1), & MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR) do m=1,nm do n=0,m-1 nmz_ToroidalEnergySpectrum_ut(n,m,:) = ut_Vmiss nmz_ToroidalEnergySpectrum_ut(n,-m,:) = ut_Vmiss enddo enddo end function nmz_ToroidalEnergySpectrum_ut function nz_ToroidalEnergySpectrum_ut(ut_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 球面調和函数全波数の各成分を計算する. ! ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは ! Σ[m=-n]^n(1/2)n(n+1)4πr^2ψ(n,m,r)^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和を動径積分したもの(r^2の重み無し) ! が球殻内での全エネルギーに等しい. ! real(8), dimension(nc,0:lm), intent(in) :: ut_TORPOT !(in) トロイダルポテンシャル real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_ut !(out) エネルギースペクトルトロイダル成分 real(8), dimension(nc,0:km) ::uz_DATA ! 作業領域 real(8), dimension(0:nm,0:km) ::nz_Work integer :: l, n, m, nm_ary(2) real(8) :: factor integer :: ierr nz_Work = 0.0D0 uz_DATA = uz_ut(ut_TORPOT) do l=1,nc nm_ary = nm_l(l) n = nm_ary(1) ; m = abs(nm_ary(2)) if ( n .ge. 0 ) then if ( m == 0 ) then factor = 1.0D0 else factor = 2.0D0 endif nz_Work(n,:) = nz_Work(n,:) & + factor * 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * uz_Data(l,:)**2 ! + 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * uz_Data(l,:)**2 endif enddo CALL MPI_ALLREDUCE(nz_Work,nz_ToroidalEnergySpectrum_ut,(nm+1)*(km+1), & MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR) end function nz_ToroidalEnergySpectrum_ut function nmz_PoloidalEnergySpectrum_ut(ut_POLPOT) ! ! ポロイダルポテンシャルから, ポロイダルエネルギーの ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する. ! ! * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分 ! φ(n,m,r)から全波数 n, 帯状波数 m 成分のポロイダルエネルギー ! スペクトルは ! ! (1/2)n(n+1)4π{|d(rφ(n,m,r))/dr|^2 + n(n+1)|φ(n,m,r)|^2} ! ! と計算される. ! ! * 全てのエネルギースペクトル成分の和を動径積分したもの(r^2の重み無し) ! が球殻内での全エネルギーに等しい. ! ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される. ! 欠損値の値はモジュール変数 ut_VMiss によって設定できる ! (初期値は -999.0) ! real(8), intent(in) :: ut_POLPOT(nc,0:lm) !(in) 流線関数(スペクトルデータ) real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_PoloidalEnergySpectrum_ut !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間) real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_Work ! 作業領域 real(8), dimension(nc,0:km) ::uz_DATA1 ! 作業領域 real(8), dimension(nc,0:km) ::uz_DATA2 ! 作業領域 integer :: l, n, m, nm_ary(2) integer :: ierr uz_Data1 = uz_ut(ut_POLPOT) uz_Data2 = uz_Rad*uz_ut(ut_DRad_ut(ut_POLPOT)) & ! d(rφ)/dr + uz_ut(ut_POLPOT) ! = rdφ/dr+φ nmz_work = 0.0D0 do l=1,nc nm_ary = nm_l(l) n = nm_ary(1) ; m = abs(nm_ary(2)) if ( n .ge. 0 ) then if ( m == 0 ) then nmz_Work(n,0,:) = nmz_Work(n,0,:) & + 0.5 * n*(n+1) * (4*pi) & * ( uz_Data2(l,:)**2 + n*(n+1)*uz_Data1(l,:)**2 ) else nmz_Work(n,m,:) = nmz_Work(n,m,:) & + 0.5 * n*(n+1) * (4*pi) & * ( uz_Data2(l,:)**2 + n*(n+1)*uz_Data1(l,:)**2 ) nmz_Work(n,-m,:) = nmz_Work(n,-m,:) & + 0.5 * n*(n+1) * (4*pi) & * ( uz_Data2(l,:)**2 + n*(n+1)*uz_Data1(l,:)**2 ) endif endif enddo CALL MPI_ALLREDUCE(nmz_Work,nmz_PoloidalEnergySpectrum_ut, & (nm+1)*(2*nm+1)*(km+1), & MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR) do m=1,nm do n=0,m-1 nmz_PoloidalEnergySpectrum_ut(n,m,:) = ut_Vmiss nmz_PoloidalEnergySpectrum_ut(n,-m,:) = ut_Vmiss enddo enddo end function nmz_PoloidalEnergySpectrum_ut function nz_PoloidalEnergySpectrum_ut(ut_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} ! ! と計算される. ! ! * 全ての全波数に対してのエネルギースペクトル成分の和を動径積分したもの ! (r^2の重み無し)が球殻内での全エネルギーに等しい. ! real(8), dimension(nc,0:lm), intent(in) :: ut_POLPOT !(in) ポロイダルポテンシャル real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_ut !(out) エネルギースペクトルポロイダル成分 real(8), dimension(nc,0:km) ::uz_DATA1 ! 作業領域 real(8), dimension(nc,0:km) ::uz_DATA2 ! 作業領域 real(8), dimension(0:nm,0:km) ::nz_Work integer :: l, n, m, nm_ary(2) real(8) :: factor integer :: ierr nz_Work = 0.0D0 uz_Data1 = uz_ut(ut_POLPOT) uz_Data2 = uz_Rad*uz_ut(ut_DRad_ut(ut_POLPOT)) & ! d(rφ)/dr + uz_ut(ut_POLPOT) ! = rdφ/dr+φ do l=1,nc nm_ary = nm_l(l) n = nm_ary(1) ; m = abs(nm_ary(2)) if ( n .ge. 0 ) then if ( m == 0 ) then factor = 1.0D0 else factor = 2.0D0 endif nz_Work(n,:) = nz_Work(n,:) & + factor * 0.5* n*(n+1)* (4*pi) & *( uz_Data2(l,:)**2 + n*(n+1)*uz_Data1(l,:)**2 ) endif enddo CALL MPI_ALLREDUCE(nz_Work,nz_PoloidalEnergySpectrum_ut,(nm+1)*(km+1), & MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR) end function nz_PoloidalEnergySpectrum_ut !--------------- 境界値問題 ----------------- subroutine ut_BoundariesTau(ut,values,cond) ! ! スペクトルデータにディリクレ・ノイマン境界条件を適用する ! Chebyshev 空間での境界条件適用(タウ法) ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を ! 定める方法をとっている(タウ法). ! real(8), dimension(nc,0:lm),intent(inout) :: ut !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,2), intent(in), optional :: values !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える. ! 省略時は値/勾配 0 となる. character(len=2), intent(in), optional :: cond !(in) 境界条件. 省略時は 'DD' ! DD : 両端ディリクレ条件 ! DN : 上端ディリクレ, 下端ノイマン条件 ! ND : 上端ノイマン, 下端ディリクレ条件 ! NN : 両端ノイマン条件 if (.not. present(cond)) then if (present(values)) then call at_BoundariesTau_DD(ut,values) else call at_BoundariesTau_DD(ut) endif return endif select case(cond) case ('NN') if (present(values)) then call at_BoundariesTau_NN(ut,values) else call at_BoundariesTau_NN(ut) endif case ('DN') if (present(values)) then call at_BoundariesTau_DN(ut,values) else call at_BoundariesTau_DN(ut) endif case ('ND') if (present(values)) then call at_BoundariesTau_ND(ut,values) else call at_BoundariesTau_ND(ut) endif case ('DD') if (present(values)) then call at_BoundariesTau_DD(ut,values) else call at_BoundariesTau_DD(ut) endif case default call MessageNotify('E','ut_BoundariesTau','B.C. not supported') end select end subroutine ut_BoundariesTau subroutine ut_BoundariesGrid(ut,values,cond) ! ! スペクトルデータにディリクレ・ノイマン境界条件を適用する ! 実空間での境界条件適用 ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! real(8), dimension(nc,0:lm),intent(inout) :: ut !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,2), intent(in), optional :: values !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える. ! 省略時は値/勾配 0 となる. character(len=2), intent(in), optional :: cond !(in) 境界条件. 省略時は 'DD' ! DD : 両端ディリクレ条件 ! DN : 上端ディリクレ, 下端ノイマン条件 ! ND : 上端ノイマン, 下端ディリクレ条件 ! NN : 両端ノイマン条件 if (.not. present(cond)) then if (present(values)) then call at_boundariesGrid_DD(ut,values) else call at_boundariesGrid_DD(ut) endif return endif select case(cond) case ('NN') if (present(values)) then call at_BoundariesGrid_NN(ut,values) else call at_BoundariesGrid_NN(ut) endif case ('DN') if (present(values)) then call at_BoundariesGrid_DN(ut,values) else call at_BoundariesGrid_DN(ut) endif case ('ND') if (present(values)) then call at_BoundariesGrid_ND(ut,values) else call at_BoundariesGrid_ND(ut) endif case ('DD') if (present(values)) then call at_BoundariesGrid_DD(ut,values) else call at_BoundariesGrid_DD(ut) endif case default call MessageNotify('E','ut_BoundariesGrid','B.C. not supported') end select end subroutine ut_BoundariesGrid subroutine ut_TorBoundariesTau(ut_TORPOT,values,cond,new) ! ! 速度トロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用. ! ! 速度トロイダルポテンシャルΨに対して与えられる境界条件は ! ! * 粘着条件 : Ψ = Ψb(lon,lat). Ψb は境界球面での速度分布. ! default は 0(静止状態). ! ! * 応力なし条件 : ∂(Ψ/r)/∂r = 0. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(inout) :: ut_TORPOT !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,2), intent(in), optional :: values !(in) 両端境界でのトロイダルポテンシャル ! 粘着条件の時のみ有効 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(0:lm) :: t_data real(8), dimension(0:km) :: z_data logical :: rigid1, rigid2 ! 境界条件 logical :: first = .true. logical :: new_matrix = .false. integer :: l save :: alu, kp, first if (.not. present(cond)) then rigid1=.TRUE. ; rigid2=.TRUE. else 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','ut_TorBoundariesTau','B.C. not supported') end select 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(0:lm,0:lm),kp(0:lm)) do l=0,lm t_data = 0.0D0 ; t_data(l) = 1.0D0 alu(:,l) = t_data ! 力学的条件粘着条件 if ( rigid1 ) then z_data = z_t(t_data) else z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad))) endif alu(lm-1,l) = z_data(0) ! 境界 k=0 での条件式代入 if ( rigid2 ) then z_data = z_t(t_data) else z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad))) endif alu(lm,l) = z_data(km) ! 境界 k=km での条件式代入 enddo call ludecomp(alu,kp) if ( rigid1 .AND. present(values) ) then call MessageNotify('M','ut_TorBoundariesTau',& 'Toroidal potential at k=0 was given by the optional variable.') else if ( rigid1 .AND. (.NOT.present(values)) ) then call MessageNotify('M','ut_TorBoundariesTau',& 'Toroidal potential at k=0 was set to zero.') else if ( (.NOT. rigid1) .AND. present(values) ) then call MessageNotify('W','ut_TorBoundariesTau',& 'Boundary value k=0 cannot be set under stress-free condition.') endif if ( rigid2 .AND. present(values) ) then call MessageNotify('M','ut_TorBoundariesTau',& 'Toroidal potential at k=0 was given by the optional variable.') else if ( rigid2 .AND. (.NOT.present(values)) ) then call MessageNotify('M','ut_TorBoundariesTau',& 'Toroidal potential at k=0 was set to zero.') else if ( (.NOT. rigid2) .AND. present(values) ) then call MessageNotify('W','ut_TorBoundariesTau',& 'Boundary value k=0 cannot be set under stress-free condition.') endif call MessageNotify('M','ut_TorBoundariesTau',& 'Matrix to apply b.c. newly produced.') endif if ( rigid1 .AND. present(values) ) then ut_torpot(:,lm-1) = values(:,1) else ut_torpot(:,lm-1) = 0.0D0 endif if ( rigid2 .AND. present(values) ) then ut_torpot(:,lm) = values(:,2) else ut_torpot(:,lm) = 0.0D0 endif ut_torpot = lusolve(alu,kp,ut_TORPOT) end subroutine ut_TorBoundariesTau subroutine ut_TorBoundariesGrid(ut_TORPOT,values,cond,new) ! ! 速度トロイダルポテンシャルに対して境界条件を適用する. ! 実空間での境界条件適用 ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! ! 速度トロイダルポテンシャルΨに対して与えられる境界条件は ! ! * 粘着条件 : Ψ = Ψb(lon,lat). Ψb は境界球面での速度分布. ! default は 0 (静止状態). ! ! * 応力なし条件 : ∂(Ψ/r)/∂r = 0. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(inout) :: ut_TORPOT !(inout) 境界条件を適用するデータ. 修正された値を返す. real(8), dimension(nc,2), intent(in), optional :: values !(in) 両端境界でのトロイダルポテンシャル ! 粘着条件の時のみ有効 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(nc,0:km):: uz_TORPOT real(8), dimension(:,:), allocatable :: alu integer, dimension(:), allocatable :: kp real(8), dimension(0:lm) :: t_data real(8), dimension(0:km) :: z_data logical :: rigid1, rigid2 ! 境界条件 logical :: first = .true. logical :: new_matrix = .false. integer :: l save :: alu, kp, first if (.not. present(cond)) then rigid1=.TRUE. ; rigid2=.TRUE. else 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','ut_TorBoundariesGrid','B.C. not supported') end select endif if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( lm /= km ) then call MessageNotify('E','TorBoundariesGrid', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(0:km,0:lm),kp(0:lm)) do l=0,lm t_data = 0.0D0 ; t_data(l)=1.0D0 alu(:,l) = z_t(t_data) if ( rigid1 ) then z_data = z_t(t_data) else z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad))) endif alu(0,l) = z_data(0) ! 境界 k=0 での条件式代入 if ( rigid2 ) then z_data = z_t(t_data) else z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad))) endif alu(km,l) = z_data(km) ! 境界 k=km での条件式代入 end do call ludecomp(alu,kp) if ( rigid1 .AND. present(values) ) then call MessageNotify('M','ut_TorBoundariesGrid',& 'Toroidal potential at k=0 was given by the optional variable.') else if ( rigid1 .AND. (.NOT.present(values)) ) then call MessageNotify('M','ut_TorBoundariesGrid',& 'Toroidal potential at k=0 was set to zero.') else if ( (.NOT. rigid1) .AND. present(values) ) then call MessageNotify('W','ut_TorBoundariesGrid',& 'Boundary value at k=0 cannot be set under stress-free condition.') endif if ( rigid2 .AND. present(values) ) then call MessageNotify('M','ut_TorBoundariesGrid',& 'Toroidal potential at k=km was given by the optional variable.') else if ( rigid2 .AND. (.NOT.present(values)) ) then call MessageNotify('M','ut_TorBoundariesGrid',& 'Toroidal potential at k=km was set to zero.') else if ( (.NOT. rigid2) .AND. present(values) ) then call MessageNotify('W','ut_TorBoundariesGrid',& 'Boundary value at k=km cannot be set under stress-free condition.') endif call MessageNotify('M','ut_TorBoundariesGrid',& 'Matrix to apply b.c. newly produced.') endif uz_TorPot = uz_ut(ut_TorPot) if ( rigid1 .AND. present(values) ) then uz_TorPot(:,0) = values(:,1) else uz_TorPot(:,0) = 0.0D0 endif if ( rigid2 .AND. present(values) ) then uz_TorPot(:,km) = values(:,2) else uz_TorPot(:,km) = 0.0D0 endif ut_torpot = lusolve(alu,kp,uz_TorPot) end subroutine ut_TorBoundariesGrid function uz_LaplaPol2Pol_uz(uz,cond,new) ! ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する. ! ! チェビシェフ格子点空間で境界条件を適用している. ! この関数を用いるためには ut_Initial にて設定する ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく ! しておく必要がある. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂r = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂r^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:km),intent(in) :: uz !(in) 入力▽^2φ分布 real(8), dimension(nc,0:km) :: uz_LaplaPol2Pol_uz !(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(nc,0:km) :: uz_work real(8), dimension(0:km,0:km) :: gg real(8), dimension(0:km,0:km) :: gg_work logical :: rigid1, rigid2 ! 境界条件 logical :: first = .true. logical :: new_matrix = .false. integer :: k,n save :: alu, kp, first if (.not. present(cond)) then rigid1=.TRUE. ; rigid2=.TRUE. else 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','uz_LaplaPol2Pol_ut','B.C. not supported') end select endif if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( lm /= km ) then call MessageNotify('E','uz_LaplaPol2Pol_uz', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(nc,0:km,0:km),kp(nc,0:km)) do k=0,km uz_work = 0.0D0 ; uz_work(:,k) = 1.0D0 ! 各水平波数に関して独立の式 alu(:,:,k) = uz_ut(ut_lapla_ut(ut_uz(uz_work))) enddo ! 運動学的条件. 流線は境界で一定 gg = 0.0D0 do k=0,km gg(k,k)=1.0D0 enddo do n=1,nc alu(n,0,:) = gg(:,0) alu(n,km,:) = gg(:,km) enddo ! 力学的条件粘着条件 if ( rigid1 ) then gg_work=az_at(at_dr_at(at_az(gg))) else gg_work=az_at(at_dr_at(at_dr_at(at_az(gg)))) endif do n=1,nc alu(n,1,:) = gg_work(:,0) enddo ! 力学的条件粘着条件 if ( rigid2 ) then gg_work=az_at(at_dr_at(at_az(gg))) else gg_work=az_at(at_dr_at(at_dr_at(at_az(gg)))) endif do n=1,nc alu(n,km-1,:) = gg_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','uz_LaplaPol2Pol_uz',& 'Matrix to apply b.c. newly produced.') endif uz_work = uz uz_work(:,1) = 0.0D0 ! 力学的条件 uz_work(:,km-1) = 0.0D0 ! 力学的条件 uz_work(:,0) = 0.0D0 ! 運動学的条件 uz_work(:,km) = 0.0D0 ! 運動学的条件 uz_laplapol2pol_uz = lusolve(alu,kp,uz_work) end function uz_LaplaPol2Pol_uz function ut_LaplaPol2PolTau_ut(ut,cond,new) ! ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する. ! チェビシェフタウ法で境界条件を適用している. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂r = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂r^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(in) :: ut !(in) 入力▽^2φ分布 real(8), dimension(nc,0:lm) :: ut_LaplaPol2PolTau_ut !(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(nc,0:km) :: uz_work real(8), dimension(nc,0:lm) :: ut_work logical :: rigid1, rigid2 ! 境界条件 logical :: first = .true. logical :: new_matrix = .false. integer :: l save :: alu, kp, first if (.not. present(cond)) then rigid1=.TRUE. ; rigid2=.TRUE. else 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','ut_LaplaPol2PolTau_ut','B.C. not supported') end select 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(nc,0:lm,0:lm),kp(nc,0:lm)) alu = 0.0D0 do l=0,lm-2 ut_work = 0.0D0 ; ut_work(:,l) = 1.0D0 ! 各水平波数に関して独立の式 alu(:,:,l) = ut_Lapla_ut(ut_work) enddo ! 運動学的条件. 流線は境界で一定 do l=0,lm ut_work=0.0D0 ut_work(:,l) = 1.0D0 uz_work = uz_ut(ut_work) alu(:,lm-1,l) = uz_work(:,0) alu(:,lm,l) = uz_work(:,km) ! 力学的条件粘着条件 if ( rigid1 ) then uz_work=uz_ut(ut_DRad_ut(ut_work)) else uz_work=uz_ut(ut_DRad_ut(ut_DRad_ut(ut_work))) endif alu(:,lm-2,l) = uz_work(:,0) ! 力学的条件粘着条件 if ( rigid2 ) then uz_work=uz_ut(ut_DRad_ut(ut_work)) else uz_work=uz_ut(ut_DRad_ut(ut_DRad_ut(ut_work))) endif alu(:,lm-3,l) = uz_work(:,km) end do call ludecomp(alu,kp) call MessageNotify('M','ut_LaplaPol2PolTau_ut',& 'Matrix to apply b.c. newly produced.') endif ut_work = ut ut_work(:,lm) = 0.0D0 ! 運動学的条件 ut_work(:,lm-1) = 0.0D0 ! 運動学的条件 ut_work(:,lm-2) = 0.0D0 ! 力学的条件 ut_work(:,lm-3) = 0.0D0 ! 力学的条件 ut_LaplaPol2PolTau_ut = lusolve(alu,kp,ut_work) end function ut_LaplaPol2PolTau_ut function ut_LaplaPol2PolGrid_ut(ut,cond,new) ! ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する. ! チェビシェフ格子点空間で境界条件を適用している. ! ! この関数を用いるためには ut_Initial にて設定する ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく ! しておく必要がある. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂r = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂r^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! ! 最終的にチェビシェフ係数の解が欲しい場合には, uz_LaplaPol2Pol_uz に ! 比べてチェビシェフ -- 格子点変換が 1 回分少なくて済む. ! real(8), dimension(nc,0:lm),intent(in) :: ut !(in) 入力▽^2φ分布 real(8), dimension(nc,0:lm) :: ut_LaplaPol2PolGrid_ut !(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(nc,0:km) :: uz_work real(8), dimension(nc,0:lm) :: ut_work real(8), dimension(0:lm,0:lm) :: tt_I real(8), dimension(0:lm,0:km) :: tz_work logical :: rigid1, rigid2 ! 境界条件 logical :: first = .true. logical :: new_matrix = .false. integer :: l,n save :: alu, kp, first if (.not. present(cond)) then rigid1=.TRUE. ; rigid2=.TRUE. else 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','ut_LaplaPol2PolGrid_ut','B.C. not supported') end select endif if (.not. present(new)) then new_matrix=.false. else new_matrix=new endif if ( first .OR. new_matrix ) then first = .false. if ( lm /= km ) then call MessageNotify('E','ut_LaplaPol2PolGrid_ut', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) allocate(alu(nc,0:km,0:lm),kp(nc,0:lm)) do l=0,lm ut_work = 0.0D0 ; ut_work(:,l) = 1.0D0 ! 各水平波数に関して独立の式 alu(:,:,l) = uz_ut(ut_Lapla_ut(ut_work)) enddo ! 運動学的条件. 流線は境界で一定 tt_I = 0.0D0 do l=0,lm tt_I(l,l)=1.0D0 enddo ! 非電気伝導体 tz_work = az_at(tt_I) do n=1,nc alu(n,0,:) = tz_work(:,0) alu(n,km,:) = tz_work(:,km) enddo ! 力学的条件粘着条件 if ( rigid1 ) then tz_work=az_at(at_Dr_at(tt_I)) else tz_work=az_at(at_Dr_at(at_Dr_at(tt_I))) endif do n=1,nc alu(n,1,:) = tz_work(:,0) enddo ! 力学的条件粘着条件 if ( rigid2 ) then tz_work=az_at(at_Dr_at(tt_I)) else tz_work=az_at(at_Dr_at(at_Dr_at(tt_I))) endif do n=1,nc alu(n,km-1,:) = tz_work(:,km) enddo call ludecomp(alu,kp) call MessageNotify('M','ut_LaplaPol2PolGrid_ut',& 'Matrix to apply b.c. newly produced.') endif uz_work = uz_ut(ut) uz_work(:,1) = 0.0D0 ! 力学的条件 uz_work(:,km-1) = 0.0D0 ! 力学的条件 uz_work(:,0) = 0.0D0 ! 運動学的条件 uz_work(:,km) = 0.0D0 ! 運動学的条件 ut_LaplaPol2PolGrid_ut = lusolve(alu,kp,uz_work) end function ut_LaplaPol2PolGrid_ut subroutine ut_TormagBoundariesTau(ut_TOR,new) ! 磁場トロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ ! 対応している. その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 外側 ! ut_psi = 0 at the outer boundary ! 内側 ! ut_psi = 0 at the inner boundary ! ! であるから ut_Boundaries で対応可能だが, 将来のため別途作成しておく. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(inout) :: ut_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: ut_I real(8), dimension(:,:), allocatable :: uz_PSI logical :: first = .true. logical :: new_matrix = .false. integer :: l 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(ut_I) ) deallocate(ut_I) if ( allocated(uz_PSI) ) deallocate(uz_PSI) allocate(alu(nc,0:lm,0:lm),kp(nc,0:lm)) allocate(ut_I(nc,0:lm),uz_PSI(nc,0:km)) do l=0,lm ut_I = 0.0D0 ; ut_I(:,l)=1.0D0 alu(:,:,l) = ut_I ! 非電気伝導体 uz_PSI = uz_ut(ut_I) alu(:,lm-1,l) = uz_PSI(:,0) alu(:,lm,l) = uz_PSI(:,km) enddo call ludecomp(alu,kp) deallocate(ut_I,uz_PSI) call MessageNotify('M','TormagBoundariesTau',& 'Matrix to apply b.c. newly produced.') endif ut_TOR(:,lm-1) = 0.0D0 ut_TOR(:,lm) = 0.0D0 ut_TOR = lusolve(alu,kp,ut_TOR) end subroutine ut_TormagBoundariesTau subroutine ut_TormagBoundariesGrid(ut_TOR,new) ! ! 磁場トロイダルポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場トロイダルポテンシャルの境界条件は ! ! 外側 ! ut_psi = 0 at the outer boundary ! 内側 ! ut_psi = 0 at the inner boundary ! ! であるので ut_Boundaries で対応可能だが, 将来のため別途作成しておく ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(inout) :: ut_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: ut_I real(8), dimension(:,:), allocatable :: uz_PSI real(8), dimension(nc,0:km) :: uz_TOR logical :: first = .true. logical :: new_matrix = .false. integer :: l 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 ( lm /= km ) then call MessageNotify('E','TorMagBoundariesGrid', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) if ( allocated(ut_I) ) deallocate(ut_I) if ( allocated(uz_PSI) ) deallocate(uz_PSI) allocate(alu(nc,0:km,0:lm),kp(nc,0:lm)) allocate(ut_I(nc,0:lm),uz_PSI(nc,0:km)) do l=0,lm ut_I = 0.0D0 ; ut_I(:,l)=1.0D0 uz_PSI = uz_ut(ut_I) alu(:,:,l) = uz_PSI ! 非電気伝導体 alu(:,0,l) = uz_PSI(:,0) alu(:,km,l) = uz_PSI(:,km) enddo call ludecomp(alu,kp) deallocate(ut_I,uz_PSI) call MessageNotify('M','TormagBoundariesGrid',& 'Matrix to apply b.c. newly produced.') endif uz_TOR = uz_ut(ut_TOR) uz_TOR(:,0) = 0.0D0 uz_TOR(:,km) = 0.0D0 ut_TOR = lusolve(alu,kp,uz_TOR) end subroutine ut_TormagBoundariesGrid subroutine ut_PolmagBoundariesTau(ut_POL,new) ! ! 磁場ポロイダルポテンシャルに対して境界条件を適用する. ! Chebyshev 空間での境界条件適用 ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ ! 対応している. その場合, 磁場ポロイダルポテンシャルの各水平スペクトル ! 成分 h にたいして境界条件が与えられ, ! ! * 外側境界 : dh/dr + (n+1)h/r = 0 ! * 内側境界 : dh/dr - nh/r = 0 ! ! である. ここで n は h の水平全波数である. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(inout) :: ut_POL !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: ut_I real(8), dimension(:,:), allocatable :: uz_PSI real(8), dimension(:,:), allocatable :: uz_DPSIDR real(8), dimension(:), allocatable :: u_n logical :: first = .true. logical :: new_matrix = .false. integer :: l, n, 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(ut_I) ) deallocate(ut_I) if ( allocated(uz_PSI) ) deallocate(uz_PSI) if ( allocated(uz_DPSIDR) ) deallocate(uz_DPSIDR) if ( allocated(u_n) ) deallocate(u_n) allocate(alu(nc,0:lm,0:lm),kp(nc,0:lm)) allocate(ut_I(nc,0:lm),u_n(nc)) allocate(uz_PSI(nc,0:km),uz_DPSIDR(nc,0:km)) do n=1,nc nn=nm_l(n) u_n(n) = nn(1) enddo do l=0,lm ut_I = 0.0D0 ; ut_I(:,l) = 1.0D0 alu(:,:,l) = ut_I uz_PSI = uz_ut(ut_I) uz_DPSIDR = uz_ut(ut_DRad_ut(ut_I)) alu(:,lm-1,l) = uz_DPSIDR(:,0) + (u_n+1) * uz_PSI(:,0)/z_RAD(0) alu(:,lm,l) = uz_DPSIDR(:,km) - u_n * uz_PSI(:,km)/z_RAD(km) enddo call ludecomp(alu,kp) deallocate(ut_I,uz_PSI,uz_DPSIDR,u_n) call MessageNotify('M','PolmagBoundariesTau',& 'Matrix to apply b.c. newly produced.') endif ut_POL(:,lm-1) = 0.0D0 ut_POL(:,lm) = 0.0D0 ut_POL = lusolve(alu,kp,ut_POL) end subroutine ut_PolmagBoundariesTau subroutine ut_PolmagBoundariesGrid(ut_POL,new) ! ! 磁場ポロイダルポテンシャルに対して境界条件を適用する. ! 鉛直実空間での境界条件適用. ! ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように ! 条件を課している(選点法). このルーチンを用いるためには ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を ! 等しくしておく必要がある. ! ! 現在のところ境界物質が非電気伝導体の場合のみ対応している. ! その場合, 磁場ポロイダルポテンシャルの各水平スペクトル成分 h に ! たいして境界条件が与えられ, ! ! * 外側境界 : dh/dr + (n+1)h/r = 0 ! * 内側境界 : dh/dr - nh/r = 0 ! ! である. ここで n は h の水平全波数である. ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(8), dimension(nc,0:lm),intent(inout) :: ut_POL !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(8), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(8), dimension(:,:), allocatable :: ut_I real(8), dimension(:,:), allocatable :: uz_PSI real(8), dimension(:,:), allocatable :: uz_DPSIDR real(8), dimension(:), allocatable :: u_n real(8), dimension(nc,0:km) :: uz_POL logical :: first = .true. logical :: new_matrix = .false. integer :: l, n, 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 ( lm /= km ) then call MessageNotify('E','PolMagBoundariesGrid', & 'Chebyshev truncation and number of grid points should be same.') endif if ( allocated(alu) ) deallocate(alu) if ( allocated(kp) ) deallocate(kp) if ( allocated(ut_I) ) deallocate(ut_I) if ( allocated(uz_PSI) ) deallocate(uz_PSI) if ( allocated(uz_DPSIDR) ) deallocate(uz_DPSIDR) if ( allocated(u_n) ) deallocate(u_n) allocate(alu(nc,0:lm,0:lm),kp(nc,0:lm)) allocate(ut_I(nc,0:lm),u_n(nc)) allocate(uz_PSI(nc,0:km),uz_DPSIDR(nc,0:km)) do n=1,nc nn=nm_l(n) u_n(n) = nn(1) enddo do l=0,lm ut_I = 0.0D0 ; ut_I(:,l) = 1.0D0 uz_PSI = uz_ut(ut_I) alu(:,:,l) = uz_PSI uz_DPSIDR = uz_ut(ut_DRad_ut(ut_I)) alu(:,0,l) = uz_DPSIDR(:,0) + (u_n+1) * uz_PSI(:,0)/z_RAD(0) alu(:,km,l) = uz_DPSIDR(:,km) - u_n * uz_PSI(:,km)/z_RAD(km) enddo call ludecomp(alu,kp) deallocate(ut_I,uz_PSI,uz_DPSIDR,u_n) call MessageNotify('M','PolmagBoundariesGrid',& 'Matrix to apply b.c. newly produced.') endif uz_POL = uz_ut(ut_POL) uz_POL(:,0) = 0.0D0 uz_POL(:,km) = 0.0D0 ut_POL = lusolve(alu,kp,uz_POL) end subroutine ut_PolmagBoundariesGrid end module ut_mpi_module_svpack