!-- !---------------------------------------------------------------------- ! Copyright (c) 2019--2021 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 ut_mpi_module_mint ! ! spml/ut_mpi_module_mint モジュールは球面上および球殻内での流体運動を ! スペクトル法と MPI 並列化によって数値計算するための Fortran90 ! 関数を提供するものである. ! ! 水平方向に球面調和函数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ua_mpi_module_mint を用いている. 最下部では球面調和変換 ! として SYPACK/ISPACK3 を, チェビシェフ変換のエンジンとして ISPACK の ! サブルーチンを用いている. ! ! !履歴 2002/05/19 竹広真一 yt_module よりモジュール名変更. ! 格子点次元を区別する必要があるため. ! 2002/06/10 竹広真一 ポロイダル磁場境界条件ルーチン追加 ! 2002/11/10 竹広真一 実空間での境界条件ルーチン追加 ! 2002/11/24 竹広真一 チェビシェフ係数入出力の ! ラプラシアン逆解きルーチンを追加 ! 2002/11/28 竹広真一 VGradV の計算変更. ! 2005/01/09 竹広真一 msgdmp -> MessageNotify に変更 ! 2005/01/09 竹広真一 ベクトル場の回転の各成分を計算する関数を追加. ! 2005/02/19 竹広真一 トロイダル粘着境界条件にて非零速度場を ! 与えられるようにオプションを追加 ! 2005/02/19 竹広真一 変数 xy_Lon, xy_Lat を追加 ! 2005/03/18 竹広真一 積分・平均関数を追加 ! z_IntLon_xz, x_IntRad_xz, IntLonRad_xz ! z_InRad_xz, y_IntRad_yz, IntLatRad_yz ! IntRad_z ! z_IntLon_xz, x_IntRad_xz, IntLonRad_xz ! z_InRad_xz, y_IntRad_yz, IntLatRad_yz ! IntRad_z ! 2005/04/24 竹広真一 スペクトル解析計算ルーチンを追加 ! nmz_ToroidalEnergySpectrum_wt ! nz_ToroidalEnergySpectrum_wt ! nmz_PoloidalEnergySpectrum_wt ! nz_PoloidalEnergySpectrum_wt ! 内部変数欠損値 wt_VMiss を追加 ! 2005/07/09 竹広真一 OPENMP 版変換ルーチンに対応 ! 2006/03/03 竹広真一 配列サイズミス修正 ! wt_wz, nmz_ToroidalEnergySpectrum_wt, ! nmz_PoloidalEnergySpectrum_wt ! 2006/03/08 竹広真一 コメントを RDoc 用に修正 ! 2006/03/19 竹広真一 変数・手続き群の要約をコメントに追加 ! 2007/08/11 竹広真一 境界値ルーチンのコメントを増強 ! 2007/09/15 竹広真一 微分ルーチンのコメントを修正 ! 2007/11/02 竹広真一 補間計算ルーチン追加 ! 2007/11/11 竹広真一 境界条件関数サブルーチンに行列設定スイッチ導入 ! 2007/11/21 竹広真一 初期化サブルーチンメッセージ出力 ! 2008/01/07 竹広真一 xyz_RotLon_wt_wt, xyz_RotLat_wt_wt バグ修正 ! 2008/01/13 竹広真一 wa_initial スイッチ導入 ! 2008/05/29 竹広真一 MPI 並列化 ! 2008/06/11 竹広真一 OPENMP スイッチ導入 ! 2008/06/11 竹広真一 配列添え字変更 ! 2009/02/27 佐々木洋平 RDoc 用のコメントを修正 ! 2010/01/07 佐々木洋平 RDoc 用のドキュメント修正, ! include 'mpif.h' -> use mpi ! 2011/09/08 竹広真一 コメント修正 ! 2011/09/12 竹広真一 wt_mpi_Initial に wa_init スイッチ導入 ! 2012/05/25 竹広真一 wt_QoperatorMPI_wt を追加 ! 2012/08/24 竹広真一 wt_LaplaPol2Pol_wt を追加 ! 2013/08/20 竹広真一 gnu fortran 対応 ! 2014/03/10 竹広真一 wt_mpi_module_nm に改造 ! 2014/03/17 竹広真一 mint パラメターに対応 ! 2014/03/21 竹広真一 OpenMP 版 ! 2015/02/15 竹広真一 Interpolate 関数追加 ! 2016/03/02 竹広真一 nz_ToroidalEnergySpectrum_wt, ! nz_PoloidalEnergySpectrum_wt 追加 ! 2016/11/04 竹広真一 nmz_ToroidalEnergySpectrum_wt, ! nmz_PoloidalEnergySpectrum_wt 追加 ! 2017/02/07 竹広真一 wt_mpi_Finalize 追加 ! 2017/10/26 竹広真一 ut_mpi_module_nmomp に改造 ! 2017/11/07 竹広真一 ut_mpi_module_nmomp_lp に改造 ! 2017/11/11 竹広真一 ut_wh, wh_ut を追加 ! 2018/01/28 竹広真一 ut_mpi_module_nmomp に改名 ! 2018/06/14 竹広真一 ut_mpi_module_extomp に改名, 変数定義変更 ! 2019/05/28 竹広真一 ut_mpi_module に改造 ! 2020/01/25 竹広真一 u_xv, xv_u を ua_mpi_module に移動 ! 2020/11/11 竹広真一 セクター計算オプション導入 ! 2021/11/01 竹広真一 Potentail2VectorMPI を ua_Potentail2VectorMPI で計算 ! 2021/11/15 竹広真一 ut_Vector2VorDivMPI, ut_RadRot_RadRotRotRot を作成 ! !凡例 ! データ種類と index ! x : 経度 y : 緯度 v : 緯度(分散格子) z : 動径 h : 動径(分散格子) ! w : 球面調和関数スペクトル ! n : 球面調和関数スペクトル(水平全波数) ! m : 球面調和関数スペクトル(帯状波数) ! t : チェビシェフ関数スペクトル ! a : 任意の次元 ! ! xvh : 3 次元分散格子点データ ! xv : 水平 2 次元格子点データ ! vh : 子午面 2 次元格子点データ ! xh : 緯度面 2 次元格子点データ ! ! wz : 水平スペクトル動径格子点データ ! wt : スペクトルデータ ! !++ module ut_mpi_module_mint ! != ut_mpi_module_mint ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id: wt_mpi_module_mint.f90 598 2013-08-20 03:23:44Z takepiro $ ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== 概要 ! ! spml/wt_mpi_module_mint モジュールは球面上および球殻内での流体運動を ! スペクトル法と MPI 並列化によって数値計算するための Fortran90 ! 関数を提供するものである. ! ! 水平方向に球面調和函数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ua_mpi_module_mint を用いている. 最下部では球面調和変換 ! として SYPACK/ISPACK3 を, チェビシェフ変換のエンジンとして ISPACK の ! サブルーチンを用いている. ! ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭(ut_, nmz_, nz_, xyz_, xvh_ uz_, w_, xv_, x_, y_, z_, a_)は, ! 返す値の形を示している. ! ut_ :: スペクトルデータ(球面調和函数・チェビシェフ変換) ! xvh_ :: 3 次元分散格子点データ(経度・緯度・動径) ! uz_ :: 水平スペクトル, 動径格子点データ ! ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 (_ut, _xyz, _xvh, _uz_, _w, _xy, _xv, _x, _y, _z, _a) は, ! 入力変数の形がスペクトルデータおよび格子点データであることを示している. ! _ut :: スペクトルデータ ! _xvh :: 3 次元分散格子点データ ! _xvh_xvh :: 2 つの3 次元分散格子点データ, ... ! !=== 各データの種類の説明 ! ! * xvh : 3 次元分散格子点データ(経度・緯度・動径) ! * 変数の種類と次元は real(8), dimension(0:im-1,jc,kc). ! * im, は経度の格子点数である. ! * jc, kc はこのプロセスで保有する緯度・動径格子点数である. ! サブルーチン ut_mpi_Initial を呼ぶと jc が設定される. ! ! * wh : スペクトルデータ ! * 変数の種類と次元は real(8), dimension(nc,kc). ! * nc はこのプロセスで保有する球面調和函数成分の数, ! lm はチェビシェフ多項式の最大次数であり, ! サブルーチン ut_mpi_Initial を呼ぶと nc が設定される. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * ut : スペクトルデータ ! * 変数の種類と次元は real(8), dimension(nmc,0:lm). ! * nmc はこのプロセスで保有する球面調和函数成分の数, ! lm はチェビシェフ多項式の最大次数であり, ! サブルーチン ut_mpi_Initial を呼ぶと nc が設定される. ! * 水平スペクトルデータの格納のされ方は関数 lu_nm, nm_lu によって調べる ! ことができる. ! ! * uz : 水平スペクトル, 動径格子点データ. ! * 変数の種類と次元は real(8), dimension(nmc,0:km). ! ! * ut_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xvh_ で始まる関数が返す値は 3 次元分散格子点データに同じ. ! ! * uz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! ut_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_Lon, x_Lon_Weight, :: 格子点座標(経度)と重みを格納した 1 次元配列 ! y_Lat, y_Lat_Weight :: 格子点座標(緯度)と重みを格納した 1 次元配列 ! v_Lat, v_Lat_Weight :: 分散格子点座標(緯度)と重みを格納した 1 次元配列 ! z_Rad, z_Rad_Weight :: 格子点座標(同形)と重みを格納した 1 次元配列 ! xvh_Lon, xvh_Lat, xvh_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (分散格子点データ型 3 次元配列) ! !==== 基本変換 ! ! xvh_ut, ut_xvh :: スペクトルデータと 3 次元分散格子データの間の変換 (球面調和函数, チェビシェフ変換) ! ! xvh_uz, uz_xvh :: 3 次元分散格子データと水平スペクトル・動径格子データとの間の変換 (球面調和函数) ! ! uz_ut, ut_uz :: スペクトルデータと水平スペクトル・動径格子データとの間の変換 (チェビシェフ変換) ! ! w_xv, xv_w :: スペクトルデータと 2 次元水平分散格子データの間の変換(球面調和函数変換) ! ! az_at, at_az :: 同時に複数個行う (チェビシェフ変換)格子データとチェビシェフデータの間の変換を ! ! l_nm, nm_l :: スペクトルデータの格納位置と全波数・帯状波数の変換 ! !==== 微分 ! ! ut_DRad_ut :: スペクトルデータに動径微分∂/∂r を作用させる ! ut_DivRad_ut :: スペクトルデータに発散型動径微分 1 /r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! ! ut_RotRad_ut :: スペクトルデータに回転型動径微分 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! ! ut_Lapla_ut :: スペクトルデータにラプラシアンを作用させる ! ! xvh_GradLon_ut :: スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λを作用させる ! ! xvh_GradLat_ut :: スペクトルデータに勾配型緯度微分 1/r・∂/∂φを作用させる(分散格子) ! ut_DivLon_xvh :: 分散格子データに発散型経度微分 1/rcosφ・∂/∂λを作用させる ! ut_DivLat_xvh :: 分散格子データに発散型緯度微分 1/rcosφ・∂(g cosφ)/∂φを作用させる(分散格子) ! ! ut_Div_xvh_xvh_xvh :: ベクトル成分である 3 つの分散格子データに発散を作用させる ! xvh_Div_xvh_xvh_xvh :: ベクトル成分である 3 つの分散格子データに発散を作用させる ! xvh_RotLon_ut_ut :: ベクトル場の回転の経度成分を計算する(分散格子) ! xvh_RotLat_ut_ut :: ベクトル場の回転の緯度成分を計算する(分散格子) ! ut_RotRad_xyz_xyz :: ベクトル場の回転の動径成分を計算する ! ut_RotRad_xvh_xvh :: ベクトル場の回転の動径成分を計算する(分散格子) ! !==== トロイダルポロイダル計算用微分 ! ! ut_L2_ut :: スペクトルデータに L2 演算子 = -水平ラプラシアンを作用させる ! ! ut_L2Inv_ut :: スペクトルデータに L2 演算子の逆 = -逆水平ラプラシアンを作用させる ! ! ut_RadRot_xvh_xvh :: ベクトル v の渦度と動径ベクトル r の内積 r・(▽×v) を計算する(分散格子) ! ut_RadRotRot_xvh_xvh_xvh :: ベクトルの v の r・(▽×▽×v) を計算する(分散格子) ! ut_Potential2VectorMPI :: トロイダルポロイダルポテンシャルからベクトル場を計算する(分散格子) ! ! ut_Potential2RotationMPI :: トロイダルポロイダルポテンシャルで表される非発散ベクトル場の回転の各成分を計算する(分散格子) ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! nmz_ToroidalEnergySpectrum_ut, :: トロイダルポテンシャルからエネルギーの ! nz_ToroidalEnergySpectrum_ut :: 球面調和函数各成分を計算する ! ! nmz_PoloidalEnergySpectrum_ut, :: ポロイダルポテンシャルからエネルギーの ! nz_PoloidalEnergySpectrum_ut :: 球面調和函数各成分を計算する ! !==== 境界値問題 ! ! ut_BoundariesTau, :: ディリクレ, ノイマン境界条件を適用する ! ut_BoundariesGrid, ::(タウ法, 選点法) ! ut_Boundaries :: ! ! ut_TorBoundariesTau, :: 速度トロイダルポテンシャルの境界条件を ! ut_TorBoundariesGrid, :: 適用する(タウ法,選点法) ! ut_TorBoundaries :: ! ! uz_LaplaPol2Pol_uz, :: 速度ポロイダルポテンシャルΦを▽^2Φから ! ut_LaplaPol2Pol_ut :: 求める (入出力がそれぞれチェビシェフ格子点, ! :: チェビシェフ係数) ! ! ut_TorMagBoundariesTau, :: 磁場トロイダルポテンシャルの境界条件を ! ut_TorMagBoundariesGrid, :: 適用する(タウ法, 選点法) ! ut_TorMagBoundaries :: ! ! ut_PolMagBoundariesTau, :: 磁場トロイダルポテンシャル境界の境界条件を ! ut_PolMagBoundariesGrid, :: 適用する(タウ法, 選点法) ! ut_PolMagBoundaries :: ! !==== 積分・平均(3 次元データ) ! ! IntLonLatRad_xvh, AvrLonLatRad_xvh :: 3 次元格子点データの全領域積分および平均 ! ! h_IntLonLat_xvh, h_AvrLonLat_xvh :: 3 次元格子点データの緯度経度(水平・球面)積分および平均 ! ! v_IntLonRad_xvh, v_AvrLonRad_xvh :: 3 次元格子点データの緯度動径積分および平均 ! ! h_IntLatRad_xvh, h_AvrLatRad_xvh :: 3 次元格子点データの経度動径(子午面)積分および平均 ! ! vh_IntLon_xvh, vz_AvrLon_xvh :: 3 次元格子点データの経度方向積分および平均 ! xh_IntLat_xvh, xz_AvrLat_xvh :: 3 次元格子点データの緯度方向積分および平均 ! !==== 積分・平均(2 次元データ) ! ! IntLonLat_xv, AvrLonLat_xv :: 2 次元格子点データの水平(球面)積分および平均 ! IntLonRad_xh, AvrLonRad_xh :: 2 次元(XZ)格子点データの経度動径積分 ! :: および平均 ! IntLatRad_vh, AvrLatRad_vh :: 2 次元(YZ)格子点データの緯度動径(子午面) ! :: 積分および平均 ! v_IntLon_xv, v_AvrLon_xv :: 水平 2 次元(球面)格子点データの経度方向 ! :: 積分および平均 ! v_IntLat_xv, x_AvrLat_xv :: 水平2 次元(球面)格子点データの緯度方向積分 ! :: および平均 ! h_IntLon_xh, h_AvrLon_xh :: 2 次元(XZ)格子点データの経度方向積分および ! :: 平均 ! h_IntRad_xh, x_AvrRad_xh :: 2 次元(XZ)格子点データの動径方向積分および ! :: 平均 ! h_IntLat_vh, z_AvrLat_vh :: 2 次元(YZ)格子点データの緯度方向積分および ! :: 平均 ! v_IntRad_vh, v_AvrRad_vh :: 2 次元(YZ)格子点データの動径方向積分および ! :: 平均 ! !==== 積分・平均(1 次元データ) ! ! IntLon_x, AvrLon_x :: 1 次元(X)格子点データの経度方向積分および平均 ! IntLat_v, AvrLat_v :: 1 次元(Y)格子点データの緯度方向積分および平均 ! IntRad_h, AvrRad_h :: 1 次元(Z)格子点データの動径方向積分および平均 ! !==== 補間計算 ! ! Interpolate_ut :: スペクトルデータから任意の点の値を補間する. ! use dc_types, only: DP use dc_message use lumatrix use mpi use ua_mpi_module_mint, & jc_ua_mpi=>jc, nc_ua_mpi=>nc, & ks_ua_mpi=>ks, ke_ua_mpi=>ke, kc_ua_mpi=>kc, & nms_ua_mpi=>nms, nme_ua_mpi=>nme, nmc_ua_mpi=>nmc use ua_mpi_module_base_mint, only: & & MPI_COMM_H, MPI_COMM_V, iproc, iproc_h, iproc_v use w_mpi_module_mint, only : y_Lat, y_Lat_Weight, & w_Vector2VorDivMPI, & Interpolate_w 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, ut_mpi_Finalize public nproc_h, nproc_v public iproc_h, iproc_v public jc ! 緯度方向分散格子点数 public ks,ke,kc public nc ! MPI 分割水平波数(w_) public nms,nme,nmc ! MPI 分割水平波数(u_) public l_nm, nm_l public lu_nm, nm_lu public kp_kl public x_Lon, x_Lon_Weight public y_Lat, y_Lat_Weight public v_Lat, v_Lat_Weight public z_Rad, z_Rad_Weight public h_Rad, h_Rad_Weight public xv_Lon, xv_Lat public xvh_Lon, xvh_Lat, xvh_Rad public wh_Rad public uz_Rad public ut_VMiss public u_xv, xv_u public at_Dr_at, t_Dr_t, az_at, at_az !public xvh_wt, wt_xvh, xvh_wh, wh_xvh public xvh_ut, ut_xvh public ut_wh, wh_ut, ua_wa, wa_ua public wz_uz, uz_wz, wz_wh, wh_wz !public wz_wh, wh_wz, wz_uz, uz_wz, az_ah, ah_az public wz_wt, wt_wz, wt_ut, ut_wt, uz_ut, ut_uz public ut_DRad_ut, ut_DivRad_ut, ut_RotRad_ut public ut_Lapla_ut, ua_LaplaH_ua, ua_LaplaHInv_ua public xvh_GradLon_ut, xvh_Gradlat_ut public ut_DivLon_xvh, ut_DivLat_xvh public ut_Div_xvh_xvh_xvh, xvh_Div_xvh_xvh_xvh public xvh_RotLon_ut_ut, xvh_RotLat_ut_ut, ut_RotRad_xvh_xvh public vh_IntLon_xvh, xh_IntLat_xvh, xv_IntRad_xvh public x_IntLatRad_xvh, v_IntLonRad_xvh, h_IntLonLat_xvh public IntLonLatRad_xvh public h_IntLon_xh, x_IntRad_xh, IntLonRad_xh public IntLon_x, IntRad_h, IntRad_z public x_IntLat_xv, v_IntLon_xv, IntLonLat_xv public h_IntLat_vh, v_IntRad_vh, IntLatRad_vh public IntLat_v public vh_AvrLon_xvh, xh_AvrLat_xvh, xv_AvrRad_xvh public x_AvrLatRad_xvh, v_AvrLonRad_xvh, h_AvrLonLat_xvh public AvrLonLatRad_xvh public h_AvrLon_xh, x_AvrRad_xh, AvrLonRad_xh public AvrLon_x, AvrRad_h, AvrRad_z public x_AvrLat_xv, v_AvrLon_xv, AvrLonLat_xv public h_AvrLat_vh, v_AvrRad_vh, AvrLatRad_vh public AvrLat_v public ut_L2_ut, ut_L2Inv_ut public ut_RadRot_xvh_xvh, ut_RadRotRot_xvh_xvh_xvh public ut_RadRot_RadRotRot public ut_Potential2vectorMPI, ut_Potential2RotationMPI public ut_Vector2VorDivMPI 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 Interpolate_ut ! ! 緯度 alon, 経度 alat 動径 arad における関数値を ! そのスペクトルデータ wt_data から補間計算する ! ! 入力する緯度経度座標は, ! １点, 経度複数点緯度動径1点 ! の 2 種類 ! module procedure Interpolate_array000_ut module procedure a_Interpolate_array100_ut end interface 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 ! 切断波数の設定(水平, 動径) integer :: mi=1 ! 経度方向対称性 real(DP) :: ri=0.0, ro=1.0 ! 球殻内外半径 real(DP), parameter :: pi=3.1415926535897932385D0 integer :: jc, nc integer :: ks, ke, kc integer :: nms, nme, nmc real(DP), dimension(:,:,:), allocatable :: xvh_LON, xvh_LAT, xvh_RAD ! 座標 real(DP), dimension(:,:), allocatable :: wh_RAD ! 座標 real(DP), dimension(:,:), allocatable :: uz_RAD ! 座標 real(DP), dimension(:), allocatable :: h_rad ! 座標 real(DP), dimension(:), allocatable :: h_rad_weight ! 座標 real(DP) :: ut_VMiss = -999.0 ! 欠損値 !!$ integer :: ks, ke, kc !!$ integer, allocatable :: ksp(:), kep(:), kcp(:) !!$ !!$ integer :: nc !!$ integer :: nms, nme, nmc !!$ integer, allocatable :: nmsp(:), nmep(:), nmcp(:) save im, jm, km, nm, lm, mi, ri, ro, ut_VMiss save nc, jc, ks, ke, kc, nms, nme, nmc contains !--------------- 初期化 ----------------- subroutine ut_mpi_Initial(i,j,k,n,l,r_in,r_out,np,npv,mint) ! ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! integer,intent(in) :: i ! 格子点数(経度λ) integer,intent(in) :: j ! 格子点数(緯度φ) integer,intent(in) :: k ! 格子点数(動径 r) integer,intent(in) :: n ! 切断波数(水平全波数) integer,intent(in) :: l ! 切断波数(動径波数) real(DP),intent(in) :: r_in ! 球殻内半径 real(DP),intent(in) :: r_out ! 球殻外半径 integer,intent(in), optional :: np ! OPENMP での最大スレッド数(dummy) integer,intent(in), optional :: npv ! 鉛直並列プロセス数 integer,intent(in), optional :: mint ! 経度方向対称性 logical :: wa_initialize=.true. ! wa_initial スイッチ integer :: npv_in ! 鉛直並列プロセス数 im = i ; jm = j ; km = k nm = n ; lm = l ri = r_in ; ro = r_out if (present(npv)) then npv_in = npv else npv_in = 1 endif if (present(mint)) then mi = mint else mi = 1 endif if ( wa_initialize ) then if ( present(np) ) then call ua_mpi_Initial(nm,im,jm,km+1,np=np,npv=npv_in,mint=mi) else call ua_mpi_Initial(nm,im,jm,km+1,npv=npv_in,mint=mi) endif endif call at_Initial(km,lm,r_in,r_out) nc = nc_ua_mpi ; jc = jc_ua_mpi ks = ks_ua_mpi; ke = ke_ua_mpi; kc = kc_ua_mpi nms = nms_ua_mpi; nme = nme_ua_mpi; nmc = nmc_ua_mpi allocate(h_Rad(kc)) allocate(h_Rad_Weight(kc)) allocate(xvh_Lon(0:im-1,jc,kc)) allocate(xvh_Lat(0:im-1,jc,kc)) allocate(xvh_Rad(0:im-1,jc,kc)) allocate(wh_Rad(nc,kc)) allocate(uz_Rad(nmc,0:km)) z_Rad_Weight = z_Rad_Weight * z_Rad**2 ! r^2 dr の積分重み h_Rad = z_Rad(ks-1:ke-1) h_Rad_Weight = z_Rad_Weight(ks-1:ke-1) xvh_Lon = spread(xv_Lon,3,kc) xvh_Lat = spread(xv_Lat,3,kc) xvh_Rad = spread(spread(h_Rad,1,jc),1,im) wh_Rad = spread(h_Rad,1,nc) uz_Rad = spread(z_Rad,1,nmc) call MessageNotify('M','ut_mpi_initial', & 'ut_mpi_module_mint (2021/11/15) is initialized') end subroutine ut_mpi_Initial !--------------- 基本変換 ----------------- function ut_xvh(xvh) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(DP), dimension(nmc,0:lm) :: ut_xvh !(out) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ !ut_xvh = ut_uz(uz_wz(wz_wh(wh_xvh(xvh)))) ut_xvh = ut_uz(ua_xvb(xvh)) end function ut_xvh function xvh_ut(ut) ! ! 水平スペクトル・動径格子点データから 3 次元格子点データへ(逆)変換する. ! real(DP), dimension(0:im-1,jc,kc) :: xvh_ut !(out) 3 次元経度緯度動径格子点データ real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数スペクトル・動径格子点データ xvh_ut = xvb_ua(uz_ut(ut)) end function xvh_ut !!$ function wt_xvh(xvh) !!$ ! !!$ ! 3 次元格子点データからスペクトルデータへ(正)変換する. !!$ ! !!$ real(DP), dimension(nc,0:lm) :: wt_xvh !!$ !(out) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !!$ !(in) 3 次元経度緯度動径格子点データ !!$ !!$ wt_xvh = wt_uz(ua_xvb(xvh)) !!$ !!$ end function wt_xvh !!$ !!$ function xvh_wt(wt) !!$ ! !!$ ! 水平スペクトル・動径格子点データから 3 次元格子点データへ(逆)変換する. !!$ ! !!$ real(DP), dimension(0:im-1,jc,kc) :: xvh_wt !!$ !(out) 3 次元経度緯度動径格子点データ !!$ !!$ real(DP), dimension(nc,0:lm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数スペクトル・動径格子点データ !!$ !!$ xvh_wt = xvb_ua(uz_wt(wt)) !!$ !!$ end function xvh_wt !!$ !!$ function wh_xvh(xvh) !!$ ! !!$ ! 3 次元格子データから水平スペクトル・動径格子点データへ(正)変換する. !!$ ! !!$ real(DP), dimension(nc,kc) :: wh_xvh !!$ !(out) 2 次元球面調和函数スペクトル・動径格子点データ !!$ !!$ real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !!$ !(in) 3 次元経度緯度動径格子点データ !!$ !!$ wh_xvh = wb_xvb(xvh) !!$ !!$ end function wh_xvh !!$ !!$ function xvh_wh(wh) !!$ ! !!$ ! 水平スペクトル・動径格子点データから 3 次元格子データへ(逆)変換する. !!$ ! !!$ real(DP), dimension(0:im-1,jc,kc) :: xvh_wh !!$ !(out) 2 次元球面調和函数スペクトル・動径格子点データ !!$ !!$ real(DP), dimension(nc,kc), intent(in) :: wh !!$ !(in) 3 次元経度緯度動径格子点データ !!$ !!$ xvh_wh = xvb_wb(wh) !!$ !!$ end function xvh_wh !!$ !!$ function u_xv(xv) !!$ ! !!$ ! 2 次元格子データから水平スペクトル(正)変換する. !!$ ! !!$ real(DP), dimension(nmc) :: u_xv !!$ !(out) 2 次元球面調和函数スペクトル・動径格子点データ !!$ !!$ real(DP), dimension(0:im-1,jc), intent(in) :: xv !!$ !(in) 3 次元経度緯度動径格子点データ !!$ !!$ u_xv = reshape(ua_wa(reshape(w_xv(xv),(/nc,1/))),(/nmc/)) !!$ !!$ end function u_xv !!$ !!$ function xv_u(u) !!$ ! !!$ ! 水平スペクトルデータから 2 次元格子データへ(逆)変換する. !!$ ! !!$ real(DP), dimension(0:im-1,jc) :: xv_u !!$ !(out) 2 次元球面調和函数スペクトル・動径格子点データ !!$ !!$ real(DP), dimension(nmc), intent(in) :: u !!$ !(in) 3 次元経度緯度動径格子点データ !!$ !!$ xv_u = xv_w(reshape(wa_ua(reshape(u,(/nmc,1/))),(/nc/))) !!$ !!$ end function xv_u function ut_wh(wh) ! ! 水平スペクトル・動径格子点データを集める. ! real(DP), dimension(nc,kc), intent(in) :: wh !(in) 2 次元球面調和函数スペクトル・動径格子点データ(分割) real(DP), dimension(nmc,0:lm) :: ut_wh !(out) 2 次元球面調和函数スペクトル・動径格子点データ !ut_wh = ut_uz(uz_wz(wz_wh(wh))) ut_wh = ut_uz(uz_wh(wh)) end function ut_wh function wh_ut(ut) ! ! 水平スペクトル・動径格子点データを分割する. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(DP), dimension(nc,kc) :: wh_ut !(out) 2 次元球面調和函数スペクトル・動径格子点データ(分割) wh_ut = wh_uz(uz_ut(ut)) end function wh_ut function wz_wh(wh) ! ! 水平スペクトル・動径格子点データを集める. ! real(DP), dimension(nc,kc), intent(in) :: wh !(in) 2 次元球面調和函数スペクトル・動径格子点データ(分割) real(DP), dimension(nc,0:km) :: wz_wh !(out) 2 次元球面調和函数スペクトル・動径格子点データ wz_wh = aa_ab(wh) end function wz_wh function wh_wz(wz) ! ! 水平スペクトル・動径格子点データを分割する. ! real(DP), dimension(nc,0:km), intent(in) :: wz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(DP), dimension(nc,kc) :: wh_wz !(out) 2 次元球面調和函数スペクトル・動径格子点データ(分割) wh_wz = ab_aa(wz) end function wh_wz function wz_wt(wt) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(DP), dimension(nc,0:lm), intent(in) :: wt !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nc,0:km) :: wz_wt !(out) 2 次元球面調和函数スペクトル・動径格子点データ wz_wt = az_at(wt) end function wz_wt function wt_wz(wz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(DP), dimension(nc,0:km), intent(in) :: wz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(DP), dimension(nc,0:lm) :: wt_wz !(out) 2 次元球面調和函数チェビシェフスペクトルデータ wt_wz = at_az(wz) end function wt_wz function ut_wt(wt) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(DP), dimension(nc,0:lm), intent(in) :: wt !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:lm) :: ut_wt !(out) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割) ut_wt = ua_wa(wt) end function ut_wt function wt_ut(ut) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割) real(DP), dimension(nc,0:lm) :: wt_ut !(out) 2 次元球面調和函数チェビシェフスペクトルデータ wt_ut = wa_ua(ut) end function wt_ut function uz_wz(wz) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(DP), dimension(nc,0:km), intent(in) :: wz !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:km) :: uz_wz !(out) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割) uz_wz = ua_wa(wz) end function uz_wz function wz_uz(uz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(DP), dimension(nmc,0:km), intent(in) :: uz !(in) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割) real(DP), dimension(nc,0:km) :: wz_uz !(out) 2 次元球面調和函数チェビシェフスペクトルデータ wz_uz = wa_ua(uz) end function wz_uz function uz_ut(ut) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:km) :: uz_ut !(out) 2 次元球面調和函数スペクトル・動径格子点データ uz_ut = az_at(ut) end function uz_ut function ut_uz(uz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(DP), dimension(nmc,0:km), intent(in) :: uz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(DP), dimension(nmc,0:lm) :: ut_uz !(out) 2 次元球面調和函数チェビシェフスペクトルデータ ut_uz = at_az(uz) end function ut_uz !!$ function wa_ua(ua) !!$ ! !!$ ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. !!$ ! !!$ real(DP), dimension(:,:), intent(in) :: ua !!$ !(in) 2 次元球面調和函数(水平再分割) (nc,:) !!$ real(DP), dimension(nc,size(ua,2)) :: wa_ua !!$ !(out) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(DP), dimension(size(ua,2),nmc) :: au !!$ real(DP), dimension(size(ua,2),nc) :: aw !!$ !!$ integer :: a_size !!$ integer :: IRQSTS !!$ integer :: IST( MPI_STATUS_SIZE ) !!$ !!$ integer :: ierr, ip_v !!$ integer, parameter :: itag=110 !!$ ! integer, allocatable ::snd_rank(:) !!$ integer, allocatable ::snd_rank_v(:) !!$ ! integer, allocatable ::rcv_rank(:) !!$ integer, allocatable ::rcv_rank_v(:) !!$ integer, allocatable ::rcv_idx(:) !!$ logical :: first=.true. !!$ !!$ save first, snd_rank_v, rcv_rank_v, rcv_idx !, snd_rank, rcv_rank !!$ !!$ if ( first ) then !!$ first = .false. !!$ ! allocate(snd_rank(1:nproc_v-1)) !!$ allocate(snd_rank_v(1:nproc_v-1)) !!$ ! allocate(rcv_rank(1:nproc_v-1)) !!$ allocate(rcv_rank_v(1:nproc_v-1)) !!$ allocate(rcv_idx(0:nproc_v-1)) !!$ !!$ do ip_v=1,nproc_v-1 !!$ rcv_rank_v(ip_v) = mod(iproc_v+ip_v,nproc_v) !!$ ! rcv_rank(ip_v) = myid_h + nproc_h*rcv_rank_v(ip_v) !!$ snd_rank_v(ip_v) = mod(iproc_v-ip_v+nproc_v,nproc_v) !!$ ! snd_rank(ip_v) = myid_h + nproc_h*snd_rank_v(ip_v) !!$ enddo !!$ !!$ rcv_idx(0) = 1 !!$ do ip_v=1,nproc_v-1 !!$ rcv_idx(ip_v) = rcv_idx(ip_v-1)+nmcp(ip_v-1) !!$ enddo !!$ endif !!$ !!$ au = transpose(ua) !!$ a_size = size(ua,2) !!$ !!$ call MPI_Barrier(MPI_COMM_V,IERR) !!$ do ip_v=1,nproc_v-1 !!$ call MPI_ISEND(au, a_size*nmc, MPI_REAL8, rcv_rank_v(ip_v),& !!$ itag,MPI_COMM_V,IRQSTS,IERR) !!$ call MPI_RECV(aw(1,rcv_idx(snd_rank_v(ip_v))), & !!$ a_size*nmcp(snd_rank_v(ip_v)), MPI_REAL8, & !!$ snd_rank_v(ip_v), itag, MPI_COMM_V, IST, IERR) !!$ CALL MPI_WAIT(IRQSTS,IST,IERR) !!$ enddo !!$ !!$ aw(:,rcv_idx(iproc_v):rcv_idx(iproc_v)+nmcp(iproc_v)-1) = au !!$ !!$ wa_ua = transpose(aw) !!$ !!$ end function wa_ua !!$ !!$ function ua_wa(wa) !!$ ! !!$ ! 水平スペクトルデータを再分割する !!$ ! !!$ real(DP), dimension(:,:), intent(in) :: wa !!$ !(in) 2 次元球面調和函数データ (nc,:) !!$ real(DP), dimension(nmc,size(wa,2)) :: ua_wa !!$ !(out) 2 次元球面調和函数データ(水平再分割) !!$ !!$ ua_wa = wa(nms:nme,:) !!$ !!$ end function ua_wa !!$ function uz_wh(wh) ! ! 水平スペクトル・動径格子点データを集める. ! real(DP), dimension(nc,kc), intent(in) :: wh !(in) 動径格子点データ(分割) real(DP), dimension(nmc,0:km) :: uz_wh !(out) 動径格子点データ uz_wh = ua_wb(wh) end function uz_wh function wh_uz(uz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(DP), dimension(nmc,0:km), intent(in) :: uz !(in) 2 次元球面調和函数(水平再分割) (nc,:) real(DP), dimension(nc,kc) :: wh_uz !(out) 2 次元球面調和函数チェビシェフスペクトルデータ wh_uz = wb_ua(uz) end function wh_uz !--------------- 微分計算 ----------------- function xvh_GradLon_ut(ut) ! ! スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λ ! を作用させる. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(0:im-1,jc,kc) :: xvh_GradLon_ut !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ xvh_GradLon_ut = xvb_GradLon_ua(uz_ut(ut))/xvh_Rad end function xvh_GradLon_ut function xvh_GradLat_ut(ut) ! ! スペクトルデータに勾配型経度微分 1/r ∂/∂φ を作用させる. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(0:im-1,jc,kc) :: xvh_GradLat_ut !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ xvh_GradLat_ut = xvb_GradLat_ua(uz_ut(ut))/xvh_Rad end function xvh_GradLat_ut function ut_DivLon_xvh(xvh) ! ! 格子点データに発散型経度微分 1/rcosφ・∂/∂λ を作用させた ! スペクトルデータを返す. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(nmc,0:lm) :: ut_DivLon_xvh !(out) 発散型経度微分を作用された 2 次元スペクトルデータ ut_DivLon_xvh = ut_uz(ua_DivLon_xvb(xvh/xvh_Rad)) end function ut_DivLon_xvh function ut_DivLat_xvh(xvh) ! ! 格子データに発散型緯度微分 1/rcosφ・∂(f cosφ)/∂φ を ! 作用させたスペクトルデータを返す. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(nmc,0:lm) :: ut_DivLat_xvh !(out) 発散型緯度微分を作用された 2 次元スペクトルデータ ut_DivLat_xvh = ut_uz(ua_divlat_xvb(xvh/xvh_Rad)) end function ut_DivLat_xvh function ut_Div_xvh_xvh_xvh(xvh_Vlon,xvh_Vlat,xvh_Vrad) ! ! べクトル成分である 3 つの格子データに発散を作用させた ! スペクトルデータを返す. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表し, 発散は ! ! 1/rcosφ・∂u/∂λ + 1/rcosφ・∂(v cosφ)/∂φ ! + 1/r^2 ∂/∂r (r^2 w) ! ! と計算される. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlon !(in) ベクトル場の経度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlat !(in) ベクトル場の緯度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vrad !(in) ベクトル場の動径成分 real(DP), dimension(nmc,0:lm) :: ut_Div_xvh_xvh_xvh !(out) ベクトル場の発散 ut_Div_xvh_xvh_xvh = ut_DivLon_xvh(xvh_Vlon) & + ut_DivLat_xvh(xvh_Vlat) & + ut_DivRad_ut(ut_xvh(xvh_Vrad)) end function ut_Div_xvh_xvh_xvh function xvh_Div_xvh_xvh_xvh(xvh_Vlon,xvh_Vlat,xvh_Vrad) ! ! ベクトル成分である 3 つの格子データに発散を作用させる. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表す. ! ! 極の特異性を回避するためにベクトル場に cosφ/r の重みをかけて ! 計算している. ! ! div V = (r/cosφ)・div (Vcosφ/r) + V_φtanφ/r + V_r/r ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlon !(in) ベクトル場の経度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlat !(in) ベクトル場の緯度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vrad !(in) ベクトル場の動径成分 real(DP), dimension(0:im-1,jc,kc) :: xvh_Div_xvh_xvh_xvh !(out) ベクトル場の発散 xvh_Div_xvh_xvh_xvh & = xvh_Rad/cos(xvh_Lat) & * xvh_ut(ut_Div_xvh_xvh_xvh(xvh_VLon*cos(xvh_Lat)/xvh_Rad, & xvh_VLat*cos(xvh_Lat)/xvh_Rad, & xvh_VRad*cos(xvh_Lat)/xvh_Rad ))& + xvh_VLat*tan(xvh_Lat)/xvh_Rad & + xvh_VRad/xvh_Rad end function xvh_Div_xvh_xvh_xvh function xvh_RotLon_ut_ut(ut_Vrad,ut_Vlat) ! ! ベクトル場の動径成分, 緯度成分である第 1, 2 引数 Vrad, Vlat から ! 回転の経度成分 ! ! 1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r を計算する. ! ! を計算する ! real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vrad !(in) ベクトル場の動径成分 real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vlat !(in) ベクトル場の緯度成分 real(DP), dimension(0:im-1,jc,kc) :: xvh_RotLon_ut_ut !(out) ベクトル場の回転の経度成分 xvh_RotLon_ut_ut = xvh_GradLat_ut(ut_Vrad) & - xvh_ut(ut_RotRad_ut(ut_Vlat)) end function xvh_RotLon_ut_ut function xvh_RotLat_ut_ut(ut_Vlon,ut_Vrad) ! ! ベクトル場の経度成分, 動径成分である第 1, 2 引数 Vlon, Vrad から ! 回転の緯度成分 ! ! 1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ ! ! を計算する. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vlon !(in) ベクトル場の経度成分 real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vrad !(in) ベクトル場の動径成分 real(DP), dimension(0:im-1,jc,kc) :: xvh_RotLat_ut_ut !(out) ベクトル場の回転の緯度成分 xvh_RotLat_ut_ut = xvh_ut(ut_RotRad_ut(ut_Vlon)) & - xvh_GradLon_ut(ut_Vrad) end function xvh_RotLat_ut_ut function ut_RotRad_xvh_xvh(xvh_Vlat,xvh_Vlon) ! ! ベクトルの緯度成分, 経度成分である第 1, 2 引数 Vlat, Vlon に対して ! ベクトル場の回転の動径成分 ! ! 1/rcosφ・∂Vlat/∂λ - 1/rcosφ・∂(Vlon cosφ)/∂φ ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlat !(in) ベクトル場の緯度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlon !(in) ベクトル場の経度成分 real(DP), dimension(nmc,0:lm) :: ut_RotRad_xvh_xvh !(out) ベクトル場の回転の動径成分 ut_RotRad_xvh_xvh = ut_DivLon_xvh(xvh_Vlat) & - ut_DivLat_xvh(xvh_Vlon) end function ut_RotRad_xvh_xvh !--------------- ポロイダル/トロイダルモデル用微分 ----------------- function ut_DRad_ut(ut) ! ! 入力スペクトルデータに動径微分 ∂/∂r を作用する. ! ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を ! 作用させたデータのスペクトル変換のことである. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,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(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:lm) :: ut_DivRad_ut !(out) 発散型動径微分を作用された 2 次元スペクトルデータ ! ut_DivRad_ut = ut_uz(uz_ut(ut_Drad_ut(ut_uz(uz_Rad**2*uz_ut(ut))))/uz_Rad**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(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:lm) :: ut_RotRad_ut !(out) 回転型動径微分を作用された 2 次元スペクトルデータ ! ut_RotRad_ut = ut_uz(uz_ut(ut_Drad_ut(ut_uz(uz_Rad*uz_ut(ut))))/uz_Rad) ut_RotRad_ut = ut_Drad_ut(ut) + ut_uz(1/uz_Rad*uz_ut(ut)) end function ut_RotRad_ut function ua_LaplaH_ua(ua_data) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用する(多層用). ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(DP), intent(in) :: ua_data(:,:) !(in) 入力スペクトルデータ real(DP) :: ua_LaplaH_ua(nmc,size(ua_data,2)) !(out) 入力スペクトルデータのラプラシアン ua_LaplaH_ua = ua_Lapla_ua(ua_data) end function ua_LaplaH_ua function ua_LaplaHInv_ua(ua_data) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用する(多層用). ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(DP), intent(in) :: ua_data(:,:) !(in) 入力スペクトルデータ real(DP) :: ua_LaplaHInv_ua(nmc,size(ua_data,2)) !(out) 入力スペクトルデータのラプラシアン ua_LaplaHInv_ua = ua_LaplaInv_ua(ua_data) end function ua_LaplaHInv_ua 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(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:lm) :: ut_Lapla_ut !(out) ラプラシアンを作用された 2 次元スペクトルデータ ut_Lapla_ut = ut_DivRad_ut(ut_Drad_ut(ut)) & + ut_uz(uz_ut(ua_LaplaH_ua(ut))/uz_Rad**2) end function ut_Lapla_ut !!$ !--------------- ポロイダル/トロイダルモデル用微分 ----------------- !!$ !!$ function wt_KxRGrad_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する. !!$ ! !!$ real(DP), dimension(nc,0:lm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(DP), dimension(nc,0:lm) :: wt_KxRGrad_wt !!$ !(out) 経度微分を作用された 2 次元スペクトルデータ !!$ !!$ wt_KxRGrad_wt = wa_Dlon_wa(wt) !!$ !!$ end function wt_KxRGrad_wt !!$ !!$ function xvz_KGrad_wt(wt) ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r !!$ ! !!$ ! 入力スペクトルデータに対応する格子データに軸方向微分 !!$ ! !!$ ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r !!$ ! !!$ ! を作用させた格子データが返される. !!$ ! ここでベクトル k は球の中心から北極向きの単位ベクトルである. !!$ ! !!$ real(DP), dimension(nc,0:lm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(DP), dimension(0:im-1,jc,0:km) :: xvz_KGrad_wt !!$ !(out) 軸方向微分を作用された 2 次元スペクトルデータ !!$ !!$ xvz_KGrad_wt = cos(xvz_Lat)*xvz_GradLat_wt(wt) & !!$ + sin(xvz_Lat)*xvz_wt(wt_Drad_wt(wt)) !!$ end function xvz_KGrad_wt function ut_L2_ut(ut) ! ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する. ! ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. ! 入力スペクトルデ ータに対応する格子点データに演算子 ! ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用させたデータのスペクトル変換が返される. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:lm) :: ut_L2_ut !(out) L^2 演算子を作用された 2 次元スペクトルデータ ut_L2_ut = -ua_LaplaH_ua(ut) end function ut_L2_ut function ut_L2Inv_ut(ut) ! ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を ! 作用する. ! ! スペクトルデータに L^2 演算子を作用させる関数 ut_L2_ut の逆計算を ! 行う関数である. ! real(DP), dimension(nmc,0:lm), intent(in) :: ut !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(DP), dimension(nmc,0:lm) :: ut_L2Inv_ut !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ ut_L2Inv_ut = -ua_LaplaHInv_ua(ut) end function ut_L2Inv_ut !!$ function wt_QOperatorMPI_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに対応する格子点データに演算子 !!$ ! !!$ ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) !!$ ! !!$ ! を作用させたデータのスペクトル変換が返される. !!$ ! !!$ real(DP), dimension(nc,0:lm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(DP), dimension(nc,0:lm) :: wt_QOperatorMPI_wt !!$ !(out) Q 演算子を作用された 2 次元スペクトルデータ !!$ !!$ wt_QOperatorMPI_wt = & !!$ wt_xvz(xvz_KGrad_wt(wt) - xvz_KGrad_wt(wt_L2_wt(wt))/2) & !!$ - wt_L2_wt(wt_xvz(xvz_KGrad_wt(wt)))/2 !!$ !!$ end function wt_QOperatorMPI_wt function ut_RadRot_xvh_xvh(xvh_VLON,xvh_VLAT) ! r・(▽×v) ! ! ベクトルの渦度と動径ベクトルの内積 r・(▽×v) を計算する. ! ! 第 1, 2 引数(v[λ], v[φ])がそれぞれベクトルの経度成分, 緯度成分を表す. ! ! r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ ! ! のスペクトル データが返される. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLON !(in) ベクトルの経度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLAT !(in) ベクトルの緯度成分 real(DP), dimension(nmc,0:lm) :: ut_RadRot_xvh_xvh !(out) ベクトルの渦度と動径ベクトルの内積 ut_RadRot_xvh_xvh = ut_uz( ua_DivLon_xvb(xvh_VLAT) & - ua_DivLat_xvb(xvh_VLON) ) !!$ ut_RadRot_xvh_xvh = ut_uz(uz_wz(wz_wh(wb_DivLon_xvb(xvh_VLAT) & !!$ - wb_DivLat_xvb(xvh_VLON)))) end function ut_RadRot_xvh_xvh function ut_RadRotRot_xvh_xvh_xvh(xvh_VLON,xvh_VLAT,xvh_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(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLON !(in) ベクトルの経度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLAT !(in) ベクトルの緯度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VRAD !(in) ベクトルの動径成分 real(DP), dimension(nmc,0:lm) :: ut_RadRotRot_xvh_xvh_xvh !(out) ベクトル v の r・(▽×▽×v) ut_RadRotRot_xvh_xvh_xvh = & ut_RotRad_ut(ut_uz( & (ua_DivLon_xvb(xvh_VLON)+ ua_DivLat_xvb(xvh_VLAT)) )) & + ut_L2_ut(ut_xvh(xvh_VRAD/xvh_RAD)) end function ut_RadRotRot_xvh_xvh_xvh subroutine ut_RadRot_RadRotRot( & xvh_VLON,xvh_VLAT,xvh_VRAD, ut_Tor, ut_LaplaPol ) ! ! ベクトル v に対して r・(▽×v) と r・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分, ! 緯度成分, 動径成分を表す. ! ! r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ ! ! r・(▽×▽×v) = 1/r ∂/∂r (r・( 1/cosφ・∂v[λ]/∂λ ! + 1/cosφ・∂(v[φ] cosφ)/∂φ ) ) ! + L^2 v[r]/r ! ! のスペクトルデータが返される. ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回少なくてすむ. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLON !(in) ベクトルの経度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLAT !(in) ベクトルの緯度成分 real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VRAD !(in) ベクトルの動径成分 real(DP), dimension(nmc,0:lm) :: ut_Tor !(out) ベクトル v の r・(▽×v) real(DP), dimension(nmc,0:lm) :: ut_LaplaPol !(out) ベクトル v の r・(▽×▽×v) call ut_Vector2VorDivMPI(xvh_VLon, xvh_VLat, ut_Tor, ut_LaplaPol) ut_LaplaPol = ut_RotRad_ut(ut_LaplaPol) & + ut_L2_ut(ut_xvh(xvh_VRAD/xvh_RAD)) end subroutine ut_RadRot_RadRotRot subroutine ut_Potential2VectorMPI(& xvh_VLON,xvh_VLAT,xvh_VRAD,ut_TORPOT,ut_POLPOT) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψr) + ▽x▽x(Φr) ! ! の各成分を計算する ! real(DP), dimension(0:im-1,jc,kc) :: xvh_VLON !(out) ベクトル場の経度成分 real(DP), dimension(0:im-1,jc,kc) :: xvh_VLAT !(out) ベクトル場の緯度成分 real(DP), dimension(0:im-1,jc,kc) :: xvh_VRAD !(out) ベクトル場の動径成分 real(DP), dimension(nmc,0:lm), intent(in) :: ut_TORPOT !(in) トロイダルポテンシャル real(DP), dimension(nmc,0:lm), intent(in) :: ut_POLPOT !(in) ポロイダルポテンシャル real(DP), dimension(nmc,0:km) :: uz_Tor !(in) トロイダルポテンシャル real(DP), dimension(nmc,0:km) :: uz_DPolDrad !(in) ポロイダルポテンシャル uz_Tor = uz_ut(ut_TORPOT) uz_DPolDRad = uz_ut(ut_RotRad_ut(ut_POLPOT)) call ua_StreamPotential2VectorMPI(-uz_Tor, uz_DPolDRad, xvh_VLON, xvh_VLAT) !!$ xvh_VLON = xvh_VLON/cos(xvh_Lat) !!$ xvh_VLAT = xvh_VLAT/cos(xvh_Lat) xvh_VRAD = xvh_ut(ut_L2_ut(ut_POLPOT))/xvh_RAD end subroutine ut_Potential2VectorMPI subroutine ut_Potential2RotationMPI(& xvh_RotVLON,xvh_RotVLAT,xvh_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(DP), dimension(0:im-1,jc,kc), intent(OUT) :: xvh_RotVLON !(out) 回転の経度成分 real(DP), dimension(0:im-1,jc,kc), intent(OUT) :: xvh_RotVLAT !(out) 回転の緯度成分 real(DP), dimension(0:im-1,jc,kc), intent(OUT) :: xvh_RotVRAD !(out) 回転の動径成分 ! 入力ベクトル場を表すポテンシャル real(DP), dimension(nmc,0:lm), intent(in) :: ut_TORPOT !(in) トロイダルポテンシャル real(DP), dimension(nmc,0:lm), intent(in) :: ut_POLPOT !(in) ポロイダルポテンシャル call ut_Potential2VectorMPI( & xvh_RotVLON,xvh_RotVLAT,xvh_RotVRAD, & -ut_Lapla_ut(ut_POLPOT), ut_TORPOT) end subroutine ut_Potential2RotationMPI subroutine ut_Vector2VorDivMPI(xvh_VLon, xvh_VLat, ut_Vor, ut_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xvh_VLon(0:im-1,jc,kc) !(in) 速度経度成分 real(8), intent(in) :: xvh_VLat(0:im-1,jc,kc) !(in) 速度緯度成分 real(8), intent(out) :: ut_Vor(nmc,0:lm) !(out) 渦度 real(8), intent(out) :: ut_Div(nmc,0:lm) !(out) 発散 real(8) :: wh_Vor(nc,kc) !速度経度成分 real(8) :: wh_Div(nc,kc) !速度緯度成分 integer :: k do k=1,kc call w_Vector2VorDivMPI & (xvh_VLon(:,:,k), xvh_VLat(:,:,k), wh_Vor(:,k), wh_Div(:,k)) enddo ut_Vor = ut_wh(wh_Vor) ut_Div = ut_wh(wh_Div) end subroutine ut_Vector2VorDivMPI !--------------- 積分計算 ----------------- !----(入力データ xvz)--- function vh_IntLon_xvh(xvh) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)積分. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(jc,kc) :: vh_IntLon_xvh !(out) 経度方向(帯状)積分された 2 次元子午面格子点データ integer :: i vh_IntLon_xvh = 0.0d0 do i=0,im-1 vh_IntLon_xvh(:,:) = vh_IntLon_xvh(:,:) & + xvh(i,:,:) * x_Lon_Weight(i) enddo end function vh_IntLon_xvh function xh_IntLat_xvh(xvh) ! ! 3 次元格子点データの緯度方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) cosφ dφ を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(0:im-1,kc) :: xh_IntLat_xvh ! 緯度円格子点データ !(out) 緯度積分された 2 次元緯度動径格子点データ xh_IntLat_xvh = xb_IntLat_xvb(xvh) end function xh_IntLat_xvh function xv_IntRad_xvh(xvh) ! 動径積分 ! ! 3 次元格子点データの動径方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dr を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(0:im-1,jc) :: xv_IntRad_xvh !(out) 動径積分された 2 次元経度緯度(水平, 球面)格子点データ real(DP), dimension(0:im-1,jc) :: xv_IntRadTmp integer :: i,j,k,ierr xv_IntRadTmp = 0.0d0 do k=1,kc do j=1,jc do i=0,im-1 xv_IntRadTmp(i,j) = xv_IntRadTmp(i,j) & + xvh(i,j,k) * h_Rad_Weight(k) enddo enddo end do CALL MPI_ALLREDUCE(xv_IntRadTmp,xv_IntRad_xvh,im*jc,MPI_REAL8, & MPI_SUM,MPI_COMM_V,IERR) end function xv_IntRad_xvh function x_IntLatRad_xvh(xvh) ! ! 3 次元格子点データの緯度動径(子午面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dφdr ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(0:im-1) :: x_IntLatRad_xvh !(out) 緯度動径(子午面)積分された 1 次元経度格子点データ x_IntLatRad_xvh = x_IntRad_xh(xb_IntLat_xvb(xvh)) end function x_IntLatRad_xvh function v_IntLonRad_xvh(xvh) ! ! 3 次元格子点データの経度動径(緯度円)積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dλdr を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(jc) :: v_IntLonRad_xvh !(out) 経度動径(緯度円)積分された 1 次元緯度格子点データ v_IntLonRad_xvh = v_IntLon_xv(xv_IntRad_xvh(xvh)) end function v_IntLonRad_xvh function h_IntLonLat_xvh(xvh) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(kc) :: h_IntLonLat_xvh !(out) 緯度経度(水平, 球面)積分された 1 次元動径格子点データ h_IntLonLat_xvh = b_IntLonLat_xvb(xvh) end function h_IntLonLat_xvh function IntLonLatRad_xvh(xvh) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP) :: IntLonLatRad_xvh !(out) 全球積分値 IntLonLatRad_xvh = IntLonLat_xv(xv_IntRad_xvh(xvh)) end function IntLonLatRad_xvh !----(入力データ vz)--- function h_IntLat_vh(vh) ! 緯度積分 ! ! 2 次元(VZ)格子点データの緯度方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) cosφ dφ を計算する. ! real(DP), dimension(jc,kc), intent(in) :: vh !(in) 2 次元緯度動径(子午面)格子点データ real(DP), dimension(kc) :: h_IntLat_vh !(out) 緯度積分された 1 次元動径格子点データ h_IntLat_vh = b_IntLat_vb(vh) end function h_IntLat_vh function v_IntRad_vh(vh) ! 動径積分 ! ! 2 次元(VH)格子点データの動径方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) r^2dr を計算する. ! real(DP), dimension(jc,kc), intent(in) :: vh !(in) 2 次元緯度動径(子午面)格子点データ real(DP), dimension(jc) :: v_IntRad_vh !(out) 動径積分された 1 次元緯度格子点データ integer :: j do j=1,jc v_IntRad_vh(j) = IntRad_h(vh(j,:)) enddo end function v_IntRad_vh function IntLatRad_vh(vh) ! ! 2 次元(VH)格子点データの緯度動径積分(子午面)および平均 ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2cosφ dφdr を計算する. ! real(DP), dimension(jc,0:km), intent(in) :: vh !(in) 2 次元緯度動径(子午面)格子点データ real(DP) :: IntLatRad_vh !(out) 積分値 !IntLatRad_vh = IntRad_h(b_IntLat_vb(vh)) IntLatRad_vh = IntLat_v(v_IntRad_vh(vh)) end function IntLatRad_vh !----(入力データ xz)--- function h_IntLon_xh(xh) ! ! 2 次元(XH)格子点データの経度方向積分. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ を計算する. ! real(DP), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(DP), 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(DP), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(DP), dimension(0:im-1) :: x_IntRad_xh !(out) 動径積分された 1 次元経度格子点データ integer :: i do i=0,im-1 x_IntRad_xh(i) = IntRad_h(xh(i,:)) enddo end function x_IntRad_xh function IntLonRad_xh(xh) ! 経度動径(緯度円)積分 ! ! 2 次元(XH)格子点データの経度動径積分 ! ! 2 次元データ f(λ,r) に対して∫f(λ,r) r^2dλdr を計算する. ! real(DP), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(DP) :: IntLonRad_xh !(out) 積分値 IntLonRad_xh = IntLon_x(x_IntRad_xh(xh)) end function IntLonRad_xh !----(入力データ z)--- function IntRad_h(h) ! 動径積分 ! ! 1 次元(Z)格子点データの動径方向域積分. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr を計算する. ! real(DP), dimension(kc), intent(in) :: h !(in) 1 次元動径格子点データ real(DP) :: IntRad_h !(out) 積分値 real(DP) :: 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_V,IERR) end function IntRad_h function IntRad_z(z) ! 動径積分 ! ! 1 次元(Z)格子点データの動径方向域積分. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr を計算する. ! real(DP), dimension(0:km), intent(in) :: z !(in) 1 次元動径格子点データ real(DP) :: IntRad_z !(out) 積分値 integer :: k IntRad_z = 0.0d0 do k=0,km IntRad_z = IntRad_z + z(k) * z_Rad_Weight(k) enddo end function IntRad_z !--------------- 平均計算 ----------------- !----(入力データ xvh)--- function vh_AvrLon_xvh(xvh) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ/2π を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(jc,kc) :: vh_AvrLon_xvh !(out) 経度方向(帯状)平均された 2 次元子午面格子点データ vh_AvrLon_xvh = vh_IntLon_xvh(xvh)/sum(x_Lon_Weight) end function vh_AvrLon_xvh function xh_AvrLat_xvh(xvh) ! 緯度積分 ! ! 3 次元格子点データの緯度方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)cosφ dφ/2 を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(im,kc) :: xh_AvrLat_xvh !(out) 緯度平均された 2 次元緯度動径格子点データ xh_AvrLat_xvh = xh_IntLat_xvh(xvh)/sum(y_Lat_Weight) end function xh_AvrLat_xvh function xv_AvrRad_xvh(xvh) ! ! 3 次元格子点データの動径方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dr/((r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(0:im-1,jc) :: xv_AvrRad_xvh ! 水平格子点データ !(out) 動径平均された 2 次元経度緯度(水平, 球面)格子点データ xv_AvrRad_xvh = xv_IntRad_xvh(xvh)/sum(z_Rad_Weight) end function xv_AvrRad_xvh function x_AvrLatRad_xvh(xvh) ! 緯度動径(子午面)積分 ! ! 3 次元格子点データの緯度動径(子午面)平均 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(0:im-1) :: x_AvrLatRad_xvh !(out) 緯度動径(子午面)平均された 1 次元経度格子点データ x_AvrLatRad_xvh = x_IntLatRad_xvh(xvh) & /( sum(y_Lat_Weight)*sum(z_Rad_Weight) ) end function x_AvrLatRad_xvh function v_AvrLonRad_xvh(xvh) ! 経度動径(緯度円)積分 ! ! 3 次元格子点データの経度動径(緯度円)平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(jc) :: v_AvrLonRad_xvh !(out) 経度動径(緯度円)平均された 1 次元緯度格子点データ v_AvrLonRad_xvh = v_IntLonRad_xvh(xvh) & /(sum(x_Lon_Weight)*sum(z_Rad_Weight)) end function v_AvrLonRad_xvh function h_AvrLonLat_xvh(xvh) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ /4π ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP), dimension(kc) :: h_AvrLonLat_xvh !(out) 緯度経度(水平, 球面)平均された 1 次元動径格子点データ h_AvrLonLat_xvh = h_IntLonLat_xvh(xvh) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight)) end function h_AvrLonLat_xvh function AvrLonLatRad_xvh(xvh) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr /(4π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh !(in) 3 次元経度緯度動径格子点データ real(DP) :: AvrLonLatRad_xvh !(out) 全球平均値 AvrLonLatRad_xvh = IntLonLatRad_xvh(xvh) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(z_Rad_Weight)) end function AvrLonLatRad_xvh !----(入力データ vh)--- function h_AvrLat_vh(vh) ! ! 2 次元(VH)格子点データの緯度方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) cosφ dφ/2 を計算する. ! real(DP), dimension(jc,kc), intent(in) :: vh !(in) 2 次元緯度動径(子午面)格子点データ real(DP), dimension(kc) :: h_AvrLat_vh !(out) 緯度平均された 1 次元動径格子点データ h_AvrLat_vh = h_IntLat_vh(vh)/sum(y_Lat_Weight) end function h_AvrLat_vh function v_AvrRad_vh(vh) ! ! 2 次元(VH)格子点データの動径方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2dr /((r[o]^3-r[i]^3)/3) ! を計算する. ! real(DP), dimension(jc,kc), intent(in) :: vh !(in) 2 次元緯度動径(子午面)格子点データ real(DP), dimension(jc) :: v_AvrRad_vh !(out) 動径平均された 1 次元緯度格子点データ v_AvrRad_vh = v_IntRad_vh(vh)/sum(z_Rad_Weight) end function v_AvrRad_vh function AvrLatRad_vh(vh) ! 緯度動径(子午面)積分 ! ! 2 次元(VH)格子点データの緯度動径(子午面)平均 ! ! 2 次元データ f(φ,r) に対して ! ! ∫f(φ,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(DP), dimension(jc,kc), intent(in) :: vh !(in) 2 次元緯度動径(子午面)格子点データ real(DP) :: AvrLatRad_vh !(out) 平均値 AvrLatRad_vh = IntLatRad_vh(vh)/(sum(y_Lat_Weight)*sum(z_Rad_Weight)) end function AvrLatRad_vh !----(入力データ xh)--- function h_AvrLon_xh(xh) ! 経度(帯状)積分 ! ! 2 次元(XH)格子点データの経度方向平均. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ/2π を計算する. ! real(DP), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(DP), 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(DP), dimension(0:im-1,kc), intent(in) :: xh !(in) 2 次元緯度動径格子点データ real(DP), 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(DP), dimension(0:im-1,kc), intent(in) :: xh ! (in) 2 次元格子点データ real(DP) :: AvrLonRad_xh ! 積分値 AvrLonRad_xh = IntLonRad_xh(xh)/(sum(x_Lon_Weight)*sum(z_Rad_Weight)) end function AvrLonRad_xh !----(入力データ z)--- function AvrRad_h(h) ! ! 1 次元(Z)格子点データの動径方向域平均. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr /((r[o]^3-r[i]^3)/3) を ! 計算する. ! real(DP), dimension(kc), intent(in) :: h !(in) 1 次元動径格子点データ real(DP) :: AvrRad_h !(out) 平均値 AvrRad_h = IntRad_h(h)/sum(z_Rad_Weight) end function AvrRad_h function AvrRad_z(z) ! ! 1 次元(Z)格子点データの動径方向域平均. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr /((r[o]^3-r[i]^3)/3) を ! 計算する. ! real(DP), dimension(0:km), intent(in) :: z !(in) 1 次元動径格子点データ real(DP) :: AvrRad_z !(out) 平均値 AvrRad_z = IntRad_z(z)/sum(z_Rad_Weight) end function AvrRad_z !--------------- 補間計算 ----------------- function Interpolate_array000_ut(ut_data,alon,alat,arad) ! ! 緯度 alon, 経度 alat 動径 arad における関数値を ! その球面調和変換係数 wa_data から補間計算する ! real(DP), intent(IN) :: ut_data(nmc,0:lm) ! スペクトルデータ real(DP), intent(IN) :: alon ! 補間する位置(経度) real(DP), intent(IN) :: alat ! 補間する位置(緯度) real(DP), intent(IN) :: arad ! 補間する位置(動径) real(DP) :: Interpolate_array000_ut ! 補間した値 real(DP) :: wa_data(nc,1) ! 補間した値 wa_data = wa_ua(reshape(a_Interpolate_at(ut_data,arad),(/nmc,1/))) Interpolate_array000_ut = Interpolate_w(wa_data(:,1),alon,alat) end function Interpolate_array000_ut function a_Interpolate_array100_ut(ut_data,alon,alat,arad) ! ! 緯度 alon, 経度 alat 動径 arad における関数値を ! その球面調和変換係数 wa_data から補間計算する ! real(DP), intent(IN) :: ut_data(:,:) ! スペクトルデータ real(DP), intent(IN) :: alon(:) ! 補間する位置(経度) real(DP), intent(IN) :: alat ! 補間する位置(緯度) real(DP), intent(IN) :: arad ! 補間する位置(動径) real(DP) :: a_Interpolate_array100_ut(size(alon))! 補間した値 real(DP) :: wa_data(nc,1) ! 補間した値 wa_data = wa_ua(reshape(a_Interpolate_at(ut_data,arad),(/nmc,1/))) a_Interpolate_array100_ut = Interpolate_w(wa_data(:,1),alon,alat) end function a_Interpolate_array100_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(DP), intent(in) :: ut_TORPOT(nmc,0:lm) !(in) 流線関数(スペクトルデータ) real(DP), dimension(0:nm,-nm:nm,0:km) :: nmz_ToroidalEnergySpectrum_ut !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間) real(DP), dimension(0:nm,-nm:nm,0:km) :: nmz_Work ! 作業領域 real(DP), dimension(nc,0:km) ::wz_DATA ! 作業領域 integer :: l, n, m, nm_ary(2) integer :: ierr wz_DATA = wz_wt(wt_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 * wz_Data(l,:)**2 ! + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2 else nmz_Work(n,m,:) = nmz_Work(n,m,:) & + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2 ! + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2 nmz_Work(n,-m,:) = nmz_Work(n,-m,:) & + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2 ! + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_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_H,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(DP), dimension(nmc,0:lm), intent(in) :: ut_TORPOT !(in) トロイダルポテンシャル real(DP), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_ut !(out) エネルギースペクトルトロイダル成分 real(DP), dimension(nc,0:km) :: wz_DATA ! 作業領域 real(DP), dimension(0:nm,0:km) :: nz_Work integer :: l, n, m, nm_ary(2) integer :: ierr real(DP) :: factor nz_Work = 0.0D0 wz_DATA = wz_wt(wt_ut(ut_TORPOT)) do l=1,nc nm_ary = nm_l(l) n = nm_ary(1) ; m = abs(nm_ary(2)) if ( m == 0 ) then factor = 1.0D0 else factor = 2.0D0 endif if ( n .ge. 0 ) then nz_Work(n,:) = nz_Work(n,:) & + factor * 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * wz_Data(l,:)**2 ! + 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * wz_Data(l,:)**2 endif enddo CALL MPI_ALLREDUCE(nz_Work,nz_ToroidalEnergySpectrum_ut,(nm+1)*(km+1), & MPI_REAL8,MPI_SUM,MPI_COMM_H,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(DP), intent(in) :: ut_POLPOT(nmc,0:lm) !(in) 流線関数(スペクトルデータ) real(DP), dimension(0:nm,-nm:nm,0:km) :: nmz_PoloidalEnergySpectrum_ut !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間) real(DP), dimension(0:nm,-nm:nm,0:km) :: nmz_Work ! 作業領域 real(DP), dimension(nc,0:km) ::wz_DATA1 ! 作業領域 real(DP), dimension(nc,0:km) ::wz_DATA2 ! 作業領域 integer :: l, n, m, nm_ary(2) integer :: ierr wz_Data1 = wz_wt(wt_ut(ut_POLPOT)) wz_Data2 = wz_uz(uz_Rad*uz_ut(ut_DRad_ut(ut_POLPOT))) & ! d(rφ)/dr + wz_wt(wt_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) & * ( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 ) ! + 0.5 * n*(n+1) * (4*pi) & ! * ( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 ) else nmz_Work(n,m,:) = nmz_Work(n,m,:) & + 0.5 * n*(n+1) * (4*pi) & * ( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 ) ! + 0.5 * n*(n+1) * (4*pi) & ! * ( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 ) nmz_Work(n,-m,:) = nmz_Work(n,-m,:) & + 0.5 * n*(n+1) * (4*pi) & * ( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 ) ! + 0.5 * n*(n+1) * (4*pi) & ! * ( wz_Data2(l,:)**2 + n*(n+1)*wz_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_H,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(DP), dimension(nmc,0:lm), intent(in) :: ut_POLPOT !(in) ポロイダルポテンシャル real(DP), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_ut !(out) エネルギースペクトルポロイダル成分 real(DP), dimension(nc,0:km) ::wz_DATA1 ! 作業領域 real(DP), dimension(nc,0:km) ::wz_DATA2 ! 作業領域 real(DP), dimension(0:nm,0:km) ::nz_Work integer :: l, n, m, nm_ary(2) integer :: ierr real(DP) :: factor nz_Work = 0.0D0 wz_Data1 = wz_wt(wt_ut(ut_POLPOT)) wz_Data2 = wz_uz(uz_Rad*uz_ut(ut_DRad_ut(ut_POLPOT))) & ! d(rφ)/dr + wz_wt(wt_ut(ut_POLPOT)) ! = rdφ/dr+φ do l=1,nc nm_ary = nm_l(l) n = nm_ary(1) ; m = abs(nm_ary(2)) if ( m == 0 ) then factor = 1.0D0 else factor = 2.0D0 endif if ( n .ge. 0 ) then nz_Work(n,:) = nz_Work(n,:) & + factor * 0.5* n*(n+1)* (4*pi) & *( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 ) endif enddo CALL MPI_ALLREDUCE(nz_Work,nz_PoloidalEnergySpectrum_ut,(nm+1)*(km+1), & MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR) end function nz_PoloidalEnergySpectrum_ut !--------------- 境界値問題 ----------------- subroutine ut_BoundariesTau(ut,values,cond) ! ! スペクトルデータにディリクレ・ノイマン境界条件を適用する ! Chebyshev 空間での境界条件適用(タウ法) ! ! チェビシェフ空間において境界条件を満たすべく高次の係数を ! 定める方法をとっている(タウ法). ! real(DP), dimension(nmc,0:lm),intent(inout) :: ut !(inout) 境界条件を適用するデータ. 修正された値を返す. real(DP), dimension(nmc,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(DP), dimension(nmc,0:lm),intent(inout) :: ut !(inout) 境界条件を適用するデータ. 修正された値を返す. real(DP), dimension(nmc,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(DP), dimension(nmc,0:lm),intent(inout) :: ut_TORPOT !(inout) 境界条件を適用するデータ. 修正された値を返す. real(DP), dimension(nmc,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(DP), dimension(:,:), allocatable :: alu integer, dimension(:), allocatable :: kp real(DP), dimension(0:lm) :: t_data real(DP), dimension(0:km) :: z_data logical :: rigid1=.true. ! 境界条件 logical :: rigid2=.true. ! 境界条件 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(DP), dimension(nmc,0:lm),intent(inout) :: ut_TORPOT !(inout) 境界条件を適用するデータ. 修正された値を返す. real(DP), dimension(nmc,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(DP), dimension(nmc,0:km) :: uz_TORPOT real(DP), dimension(:,:), allocatable :: alu integer, dimension(:), allocatable :: kp real(DP), dimension(0:lm) :: t_data real(DP), dimension(0:km) :: z_data logical :: rigid1=.true. ! 境界条件 logical :: rigid2=.true. ! 境界条件 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Φから計算する. ! ! チェビシェフ格子点空間で境界条件を適用している. ! この関数を用いるためには wt_Initial にて設定する ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく ! しておく必要がある. ! ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は ! ! ▽^2Φ = f ! Φ = const. at boundaries. ! ∂Φ/∂r = 0 at boundaries (粘着条件) ! or ∂^2Φ/∂r^2 = 0 at boundaries (応力なし条件) ! ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される. ! real(DP), dimension(nmc,0:km),intent(in) :: uz !(in) 入力▽^2φ分布 real(DP), dimension(nmc,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(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(nmc,0:km) :: uz_work real(DP), dimension(0:km,0:km) :: gg real(DP), dimension(0:km,0:km) :: gg_work logical :: rigid1=.true. ! 境界条件 logical :: rigid2=.true. ! 境界条件 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_uz','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(nmc,0:km,0:km),kp(nmc,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,nmc 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,nmc 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,nmc 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(DP), dimension(nmc,0:lm),intent(in) :: ut !(in) 入力▽^2φ分布 real(DP), dimension(nmc,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(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(nmc,0:km) :: uz_work real(DP), dimension(nmc,0:lm) :: ut_work logical :: rigid1=.true. ! 境界条件 logical :: rigid2=.true. ! 境界条件 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(nmc,0:lm,0:lm),kp(nmc,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(DP), dimension(nmc,0:lm),intent(in) :: ut !(in) 入力▽^2φ分布 real(DP), dimension(nmc,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(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(nmc,0:km) :: uz_work real(DP), dimension(nmc,0:lm) :: ut_work real(DP), dimension(0:lm,0:lm) :: tt_I real(DP), dimension(0:lm,0:km) :: tz_work logical :: rigid1=.true. ! 境界条件 logical :: rigid2=.true. ! 境界条件 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(nmc,0:km,0:lm),kp(nmc,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,nmc 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,nmc 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,nmc 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(DP), dimension(nmc,0:lm),intent(inout) :: ut_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(:,:), allocatable :: ut_I real(DP), 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(nmc,0:lm,0:lm),kp(nmc,0:lm)) allocate(ut_I(nmc,0:lm),uz_PSI(nmc,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(DP), dimension(nmc,0:lm),intent(inout) :: ut_TOR !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(:,:), allocatable :: ut_I real(DP), dimension(:,:), allocatable :: uz_PSI real(DP), dimension(nmc,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(nmc,0:km,0:lm),kp(nmc,0:lm)) allocate(ut_I(nmc,0:lm),uz_PSI(nmc,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(DP), dimension(nmc,0:lm),intent(inout) :: ut_POL !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(:,:), allocatable :: ut_I real(DP), dimension(:,:), allocatable :: uz_PSI real(DP), dimension(:,:), allocatable :: uz_DPSIDR real(DP), 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(nmc,0:lm,0:lm),kp(nmc,0:lm)) allocate(ut_I(nmc,0:lm),u_n(nmc)) allocate(uz_PSI(nmc,0:km),uz_DPSIDR(nmc,0:km)) do n=1,nmc 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(DP), dimension(nmc,0:lm),intent(inout) :: ut_POL !(inout) 境界条件を適用するデータ. 修正された値を返す. logical, intent(IN), optional :: new !(in) true だと境界条件計算用行列を強制的に新たに作る. ! default は false. real(DP), dimension(:,:,:), allocatable :: alu integer, dimension(:,:), allocatable :: kp real(DP), dimension(:,:), allocatable :: ut_I real(DP), dimension(:,:), allocatable :: uz_PSI real(DP), dimension(:,:), allocatable :: uz_DPSIDR real(DP), dimension(:), allocatable :: u_n real(DP), dimension(nmc,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(nmc,0:lm,0:lm),kp(nmc,0:lm)) allocate(ut_I(nmc,0:lm),u_n(nmc)) allocate(uz_PSI(nmc,0:km),uz_DPSIDR(nmc,0:km)) do n=1,nmc 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 !--------------- 終了処理 ----------------- subroutine ut_mpi_Finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! このサブルーチンを単独で用いるのでなく, ! 上位サブルーチン wa_mpi_Finalize を使用すること. ! deallocate(h_Rad) deallocate(h_Rad_Weight) deallocate(xvh_Lon) deallocate(xvh_Lat) deallocate(xvh_Rad) deallocate(wh_Rad) deallocate(uz_Rad) call ua_mpi_finalize call at_Finalize call MessageNotify('M','ut_mpi_finalize',& 'ut_mpi_module_mint (2021/11/15) is finalized') end subroutine ut_mpi_Finalize end module ut_mpi_module_mint