!-- !---------------------------------------------------------------------- ! Copyright(c) 2017 SPMDODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 wsc_module_svpack ! ! spml/wsc_module_svpack モジュールは球面上および球殻内での流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向に球面調和函数変換および鉛直方向に sin/cos 関数変換を ! 用いる場合のスペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wa_module_svpack, asc_module を用いている. 最下部では ! 球面調和変換およびチェビシェフ変換のエンジンとして ISPACK の ! Fortran77 サブルーチンを用いている. ! ! 関数, サブルーチンの名前と機能は wt_module のものと同じである. ! したがって use 文を wsc_module から wsc_module_svpack に ! 変更するだけで SVPACK の機能が使えるようになる. ! ! ただし l_nm, nm_l の使い方には注意されたい. wt_module の l_nm, nm_l は ! wt_Initial で初期化しなくとも用いることができる(結果が切断波数に依らない)が, ! wsc_module_svpack のものは初期化したのちにしか使うことができない. ! ! !履歴 2017/05/21 竹広真一 ! !凡例 ! データ種類と index ! x : 経度 y : 緯度 z : 動径 ! w : 球面調和関数スペクトル ! n : 球面調和関数スペクトル(水平全波数) ! m : 球面調和関数スペクトル(帯状波数) ! s : sin 関数スペクトル ! c : cos 関数スペクトル ! a : 任意の次元 ! ! xyz : 3 次元格子点データ ! xy : 水平 2 次元格子点データ ! yz : 子午面 2 次元格子点データ ! xz : 緯度面 2 次元格子点データ ! ! wz : 水平スペクトル動径格子点データ ! ws : スペクトルデータ(鉛直 sin) ! wc : スペクトルデータ(鉛直 cos) ! !++ module wsc_module_svpack ! != wsc_module_svpack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/wsc_module_svpack モジュールは球面上および球殻内での流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向に球面調和函数変換および上下の境界壁を扱うための ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で wa_module_sjpack, asc_module を用いている. 最下部では ! 球面調和変換およびチェビシェフ変換のエンジンとして ISPACK の ! Fortran77 サブルーチンを用いている. ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (ws_, wc_, nmz_,nz_,xyz_,wz_,w_,xy_,x_,y_,s_,c_,a_) は, ! 返す値の形を示している. ! ws_,wc_ :: スペクトルデータ(球面調和函数・三角関数変換) ! nmz_ :: 水平スペクトルデータ(全波数 n, 帯状波数各成分, 動径) ! nz_ :: 水平スペクトルデータ(全波数 n, 動径) ! xyz_ :: 3 次元格子点データ(経度・緯度・動径) ! wz_ :: 水平スペクトル, 動径格子点データ ! ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 (ws_,wc_,xyz_,wz_,w_,xy_,x_,y_,z_,a_) は, 入力変数の ! 形がスペクトルデータおよび格子点データであることを示している. ! _ws,_wc :: スペクトルデータ ! _xyz :: 3 次元格子点データ ! _xyz_xyz :: 2 つの3 次元格子点データ, ... ! !=== 各データの種類の説明 ! ! * xyz : 3 次元格子点データ(経度・緯度・動径) ! * 変数の種類と次元は real(8), dimension(0:im-1,1:jm,0:km). ! * im, jm, km はそれぞれ経度, 緯度, 動径座標の格子点数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * ws : スペクトルデータ(鉛直 sin) ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),lm). ! * nm は球面調和函数の最大全波数, lm は三角関数展開の最大波数 ! であり, サブルーチン wsc_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! ! * nmz : 水平スペクトルデータの並んだ 3 次元配列. ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,0:km). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wt_Initial にて ! あらかじめ設定しておく. ! ! * nz : スペクトルデータの並んだ 2 次元配列. ! * 変数の種類と次元は real(8), dimension(0:nm,0:km). ! * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり, ! サブルーチン wt_Initial にてあらかじめ設定しておく. ! ! * wz : 水平スペクトル, 動径格子点データ. ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),0:km). ! ! * ws_, wc_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyz_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * wz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! wsc_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_Lon, y_Lat, z_Rad :: 格子点座標(緯度, 経度, 動径座標)を ! 格納した1 次元配列 ! x_Lon_Weight, y_Lat_Weight, z_Rad_Weight :: 重み座標を格納した 1 次元配列 ! xyz_Lon, xyz_Lat, xyz_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) ! (格子点データ型 3 次元配列) ! !==== 基本変換 ! ! xyz_ws, ws_xyz :: スペクトルデータと 3 次元格子データの間の変換 ! (球面調和函数, sin 変換) ! xyz_wc, wc_xyz :: スペクトルデータと 3 次元格子データの間の変換 ! (球面調和函数, cos 変換) ! xyz_wz, wz_xyz :: 3 次元格子データと水平スペクトル・動径格子データとの間 ! の変換 (球面調和函数) ! wz_ws, ws_wz :: スペクトルデータと水平スペクトル・動径格子データとの間 ! の変換 (sin 変換) ! wz_wc, wc_wz :: スペクトルデータと水平スペクトル・動径格子データとの間 ! の変換 (cos 変換) ! w_xy, xy_w :: スペクトルデータと 2 次元水平格子データの間の変換 ! (球面調和函数変換) ! az_as, as_az :: 格子データと sin 変換スペクトルデータの間の変換を ! 同時に複数個行う ! az_as, as_az :: 格子データと cos 変換スペクトル間の変換を ! 同時に複数個行う ! l_nm, nm_l :: スペクトルデータの格納位置と全波数・帯状波数の変換 ! !==== 微分 ! ! ws_DRad_wc, wc_DRad_ws :: スペクトルデータに動径微分∂/∂r を作用させる ! ws_DivRad_wc, wc_DivRad_ws :: スペクトルデータに発散型動径微分 ! ∂/∂r ! ws_RotRad_wc, wc_RotRad_ws :: スペクトルデータに回転型動径微分 ! ∂/∂r を作用させる ! ws_Lapla_ws, wc_Lapla_wc :: スペクトルデータにラプラシアンを作用させる ! ws_LaplaInv_ws, wc_LaplaInv_wc :: スペクトルデータに逆ラプラシアンを作用させる ! xyz_GradLon_ws, xyz_GradLon_wc :: スペクトルデータに勾配型経度微分 ! 1/(a cosφ)・∂/∂λを作用させる ! xyz_GradLat_wt, xyz_GradLat_wc :: スペクトルデータに勾配型緯度微分 ! 1/a・∂/∂φを作用させる ! ws_DivLon_xyz, wc_DivLon_xyz :: 格子データに発散型経度微分 ! 1/(a cosφ)・∂/∂λを作用させる ! ws_DivLat_xyz, wc_DivLat_xyz :: 格子データに発散型緯度微分 ! (1/a cosφ)・∂(g cosφ)/∂φを作用させる ! wc_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! ws_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyz_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! wc_RotRad_xyz_xyz :: ベクトル場の回転の動径成分を計算する ! !==== トロイダルポロイダル計算用微分 ! ! ws_KxRGrad_ws, wc_KxRGrad_wc :: スペクトルデータに ! 経度微分 k×r・▽ = ∂/∂λを作用させる ! ws_L2_ws, wc_L2_wc :: スペクトルデータに ! L2 演算子 = -水平ラプラシアンを作用させる ! ws_L2Inv_ws, ws_L2Inv_ws :: スペクトルデータに ! L2 演算子の逆 = -逆水平ラプラシアンを作用させる ! wc_RadRot_xyz_xyz :: ベクトル v の渦度と動径ベクトル r の内積 ! r・(▽×v) を計算する ! ws_RadRotRot_xyz_xyz_xyz :: ベクトルの v の r・(▽×▽×v) を計算する ! wsc_Potential2Vector :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! wsc_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! nmz_ToroidalEnergySpectrum_wc, nz_ToroidalEnergySpectrum_wc :: ! トロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! nmz_PoloidalEnergySpectrum_ws, nz_PoloidalEnergySpectrum_ws :: ! ポロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! !==== 積分・平均(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 次元格子点データの ! 動径方向積分および平均 ! !==== 積分・平均(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)格子点データの動径方向積分および ! 平均 ! !==== 積分・平均(1 次元データ) ! ! IntLon_x, AvrLon_x :: 1 次元(X)格子点データの経度方向積分および平均 ! IntLat_y, AvrLat_y :: 1 次元(Y)格子点データの緯度方向積分および平均 ! IntRad_z, AvrRad_z :: 1 次元(Z)格子点データの動径方向積分および平均 ! !==== 補間計算 ! ! Interpolate_ws, Interpolate_wc :: スペクトルデータから任意の点の値を補間する. ! use dc_message use lumatrix use wa_module_svpack use asc_module, z_RAD => g_X, z_RAD_WEIGHT => g_X_WEIGHT, & as_az => as_ag, az_as => ag_as, & ac_az => ac_ag, az_ac => ag_ac, & s_z => s_g, z_s => g_s, c_z => c_g, z_c => g_c, & s_Dr_c => s_Dx_c, as_Dr_ac => as_Dx_ac, & c_Dr_s => s_Dx_c, ac_Dr_as => ac_Dx_as implicit none private public wsc_Initial public x_Lon, x_Lon_Weight public y_Lat, y_Lat_Weight public z_Rad, z_Rad_Weight public l_nm, nm_l public xy_Lon, xy_Lat public xyz_Lon, xyz_Lat, xyz_Rad public wz_Rad public wsc_VMiss public w_xy, xy_w public ac_Dr_as, c_Dr_s, az_as, as_az public as_Dr_ac, s_Dr_c, az_ac, ac_az public xyz_ws, ws_xyz, wz_ws, ws_wz public xyz_wc, wc_xyz, wz_wc, wc_wz public xyz_wz, wz_xyz public ws_DRad_wc, wc_DRad_ws, wc_DivRad_ws, ws_DivRad_wc public wc_RotRad_ws, ws_RotRad_wc public ws_Lapla_ws, wc_Lapla_wc, ws_LaplaInv_ws, wc_LaplaInv_wc public xyz_GradLon_ws, xyz_GradLon_wc, xyz_gradlat_ws, xyz_gradlat_wc public wc_DivLon_xyz, ws_DivLon_xyz, wc_DivLat_xyz, ws_DivLat_xyz public wc_Div_xyz_xyz_xyz, ws_Div_xyz_xyz_xyz, xyz_Div_xyz_xyz_xyz public wc_RotRad_xyz_xyz ! xyz_RotLon_ws_wc, xyz_RotLat_wc_ws public yz_IntLon_xyz, xz_IntLat_xyz, xy_IntRad_xyz public x_IntLatRad_xyz, y_IntLonRad_xyz, z_IntLonLat_xyz public IntLonLatRad_xyz public x_IntLat_xy, y_IntLon_xy, IntLonLat_xy public z_IntLat_yz, y_IntRad_yz, IntLatRad_yz public z_IntLon_xz, x_IntRad_xz, IntLonRad_xz public IntLon_x, IntLat_y, IntRad_z public yz_AvrLon_xyz, xz_AvrLat_xyz, xy_AvrRad_xyz public x_AvrLatRad_xyz, y_AvrLonRad_xyz, z_AvrLonLat_xyz public AvrLonLatRad_xyz public x_AvrLat_xy, y_AvrLon_xy, AvrLonLat_xy public z_AvrLat_yz, y_AvrRad_yz, AvrLatRad_yz public z_AvrLon_xz, x_AvrRad_xz, AvrLonRad_xz public AvrLon_x, AvrLat_y, AvrRad_z public wc_KxRGrad_wc, ws_KxRGrad_ws public wc_L2_wc, ws_L2_ws, wc_L2Inv_wc, ws_L2Inv_ws public wc_RadRot_xyz_xyz, ws_RadRotRot_xyz_xyz_xyz public wsc_Potential2vector, wsc_Potential2Rotation public Interpolate_wc, Interpolate_ws public nmz_ToroidalEnergySpectrum_wc, nz_ToroidalEnergySpectrum_wc public nmz_PoloidalEnergySpectrum_ws, nz_PoloidalEnergySpectrum_ws 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 :: xyz_LON, xyz_LAT, xyz_RAD ! 座標 real(8), dimension(:,:), allocatable :: wz_RAD ! 座標 real(8) :: wsc_VMiss = -999.0 ! 欠損値 save im, jm, km, nm, lm, ri, ro contains !--------------- 初期化 ----------------- subroutine wsc_Initial(i,j,k,n,l,r_in,r_out,np,wa_init) ! ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! ! np に 1 より大きな値を指定すれば ISPACK の球面調和函数変換 ! OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには, ! 実行時に環境変数 OMP_NUM_THREADS を np 以下の数字に設定する等の ! システムに応じた準備が必要となる. ! ! np に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない. ! ! 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 :: np ! OPENMP での最大スレッド数 logical,intent(in), optional :: wa_init ! wa_initial スイッチ logical :: wa_initialize=.true. ! wa_initial スイッチ im = i ; jm = j ; km = k nm = n ; lm = l ri = r_in ; ro = r_out if ( present(wa_init) ) then wa_initialize = wa_init else wa_initialize = .true. endif if ( wa_initialize ) then if ( present(np) ) then call wa_Initial(nm,im,jm,km+1,np) else call wa_Initial(nm,im,jm,km+1) endif endif call asc_Initial(km,lm,r_in,r_out) allocate(xyz_Lon(0:im-1,1:jm,0:km)) allocate(xyz_Lat(0:im-1,1:jm,0:km)) allocate(xyz_Rad(0:im-1,1:jm,0:km)) allocate(wz_Rad((nm+1)*(nm+1),0:km)) xyz_Lon = spread(xy_Lon,3,km+1) xyz_Lat = spread(xy_Lat,3,km+1) xyz_Rad = spread(spread(z_Rad,1,jm),1,im) wz_Rad = spread(z_Rad,1,(nm+1)*(nm+1)) z_Rad_Weight = z_Rad_Weight * Ro**2 ! r^2 dr の積分重み call MessageNotify('M','wsc_initial', & 'wsc_module_svpack (2017/05/21) is initialized') end subroutine wsc_Initial !--------------- 基本変換 ----------------- function xyz_ws(ws) ! ! スペクトルデータから 3 次元格子点データへ(逆)変換する. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_ws !(out) 3 次元経度緯度動径格子点データ xyz_ws = xya_wa(wz_ws(ws)) end function xyz_ws function ws_xyz(xyz) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),lm) :: ws_xyz !(out) 2 次元球面調和函数 sin スペクトルデータ ws_xyz = ws_wz(wa_xya(xyz)) end function ws_xyz function xyz_wc(wc) ! ! スペクトルデータから 3 次元格子点データへ(逆)変換する. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_wc !(out) 3 次元経度緯度動径格子点データ xyz_wc = xya_wa(wz_wc(wc)) end function xyz_wc function wc_xyz(xyz) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_xyz !(out) 2 次元球面調和函数 cos スペクトルデータ wc_xyz = wc_wz(wa_xya(xyz)) end function wc_xyz function xyz_wz(wz) ! ! 水平スペクトル・動径格子点データから 3 次元格子点データへ(逆)変換する. ! real(8), dimension((nm+1)*(nm+1),0:km), intent(in) :: wz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_wz !(out) 3 次元経度緯度動径格子点データ xyz_wz = xya_wa(wz) end function xyz_wz function wz_xyz(xyz) ! ! 3 次元格子データから水平スペクトル・動径格子点データへ(正)変換する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),0:km) :: wz_xyz !(out) 2 次元球面調和函数スペクトル・動径格子点データ wz_xyz = wa_xya(xyz) end function wz_xyz function wz_ws(ws) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:km) :: wz_ws !(out) 2 次元球面調和函数スペクトル・動径格子点データ wz_ws = az_as(ws) end function wz_ws function ws_wz(wz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension((nm+1)*(nm+1),0:km), intent(in) :: wz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(8), dimension((nm+1)*(nm+1),lm) :: ws_wz !(out) 2 次元球面調和函数 sin スペクトルデータ ws_wz = as_az(wz) end function ws_wz function wz_wc(wc) ! ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:km) :: wz_wc !(out) 2 次元球面調和函数スペクトル・動径格子点データ wz_wc = az_ac(wc) end function wz_wc function wc_wz(wz) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension((nm+1)*(nm+1),0:km), intent(in) :: wz !(in) 2 次元球面調和函数スペクトル・動径格子点データ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_wz !(out) 2 次元球面調和函数 cos スペクトルデータ wc_wz = ac_az(wz) end function wc_wz !--------------- 微分計算 ----------------- function wc_DRad_ws(ws) ! ! 入力スペクトルデータに動径微分 ∂/∂r を作用する. ! ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_DRad_ws !(in) 動径微分された2 次元球面調和函数 cos スペクトルデータ wc_DRad_ws = ac_Dr_as(ws) end function wc_DRad_ws function ws_DRad_wc(wc) ! ! 入力スペクトルデータに動径微分 ∂/∂r を作用する. ! ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_DRad_wc !(in) 動径微分された2 次元球面調和函数 sin スペクトルデータ ws_DRad_wc = as_Dr_ac(wc) end function ws_DRad_wc function wc_DivRad_ws(ws) ! ! 入力スペクトルデータに近似的発散型動径微分 ∂/∂r を作用する. ! ! スペクトルデータの発散型動径微分とは, 対応する格子点データに ! 発散型動径微分を作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_DivRad_ws !(out) 発散型動径微分を作用された 2 次元スペクトルデータ wc_DivRad_ws = wc_Drad_ws(ws) end function wc_DivRad_ws function ws_DivRad_wc(wc) ! ! 入力スペクトルデータに近似的発散型動径微分 ∂/∂r を作用する. ! ! スペクトルデータの発散型動径微分とは, 対応する格子点データに ! 発散型動径微分を作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_DivRad_wc !(out) 発散型動径微分を作用された 2 次元スペクトルデータ ws_DivRad_wc = ws_Drad_wc(wc) end function ws_DivRad_wc function wc_RotRad_ws(ws) ! ! 入力スペクトルデータに近似的回転型動径微分 ∂/∂r を作用する. ! ! スペクトルデータの回転型動径微分とは, 対応する格子点データに ! 回転型動径微分を作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_RotRad_ws !(out) 回転型動径微分を作用された 2 次元スペクトルデータ wc_RotRad_ws = wc_Drad_ws(ws) end function wc_RotRad_ws function ws_RotRad_wc(wc) ! ! 入力スペクトルデータに近似的回転型動径微分 ∂/∂r を作用する. ! ! スペクトルデータの回転型動径微分とは, 対応する格子点データに ! 回転型動径微分を作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_RotRad_wc !(out) 回転型動径微分を作用された 2 次元スペクトルデータ ws_RotRad_wc = ws_Drad_wc(wc) end function ws_RotRad_wc function ws_Lapla_ws(ws) ! 入力スペクトルデータに近似的ラプラシアン ! ! ▽^2 = 1/a^2 cos^2φ・∂^2/∂λ^2 ! + 1/a^2 cosφ・∂/∂φ(cosφ∂/∂φ) ! + ∂^2/∂r^2 ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_Lapla_ws !(out) ラプラシアンを作用された 2 次元スペクトルデータ ws_Lapla_ws = ws_Drad_wc(wc_Drad_ws(ws)) + wa_Lapla_wa(ws)/Ro**2 end function ws_Lapla_ws function wc_Lapla_wc(wc) ! 入力スペクトルデータに近似的ラプラシアン ! ! ▽^2 = 1/a^2 cos^2φ・∂^2/∂λ^2 ! + 1/a^2 cosφ・∂/∂φ(cosφ∂/∂φ) ! + ∂^2/∂r^2 ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_Lapla_wc !(out) ラプラシアンを作用された 2 次元スペクトルデータ wc_Lapla_wc = wc_DRad_ws(ws_DRad_wc(wc)) + wa_Lapla_wa(wc)/Ro**2 end function wc_Lapla_wc function wc_LaplaInv_wc(wc) ! ! スペクトルデータに逆ラプラシアンを作用させる. ! real(8), dimension((nm+1)*(nm+1),0:lm),intent(in) :: wc real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_LaplaInv_wc real(8), allocatable, dimension(:,:) :: wc_LaplaInvFactor logical :: first=.true. integer :: n, m, l save first, wc_LaplaInvFactor if ( first ) then first = .false. allocate(wc_LaplaInvFactor((nm+1)*(nm+1),0:lm)) do n=0,nm do m=-n,n do l=0,lm if ( n==0 .and. l==0 ) then wc_LaplaInvFactor(l_nm(n,0),l) = 0.0D0 else wc_LaplaInvFactor(l_nm(n,m),l) & = -1.0D0/( n*(n+1)/Ro**2 + (PI/(Ro-Ri)*l)**2 ) end if end do end do end do end if wc_LaplaInv_wc = wc_LaplaInvFactor * wc end function wc_LaplaInv_wc function ws_LaplaInv_ws(ws) ! ! スペクトルデータに逆ラプラシアンを作用させる. ! real(8), dimension((nm+1)*(nm+1),lm),intent(in) :: ws real(8), dimension((nm+1)*(nm+1),lm) :: ws_LaplaInv_ws real(8), allocatable, dimension(:,:) :: ws_LaplaInvFactor logical :: first=.true. integer :: n, m, l save first, ws_LaplaInvFactor if ( first ) then first = .false. allocate(ws_LaplaInvFactor((nm+1)*(nm+1),lm)) do n=0,nm do m=-n,n do l=1,lm ws_LaplaInvFactor(l_nm(n,m),l) & = -1.0D0/( n*(n+1)/Ro**2 + (PI/(Ro-Ri)*l)**2 ) end do end do end do end if ws_LaplaInv_ws = ws_LaplaInvFactor * ws end function ws_LaplaInv_ws function xyz_GradLon_wc(wc) ! ! スペクトルデータに近似的勾配型経度微分 1/(a cosφ)・∂/∂λ ! を作用させる. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_GradLon_wc !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ xyz_GradLon_wc = xya_GradLon_wa(wz_wc(wc))/Ro end function xyz_GradLon_wc function xyz_GradLon_ws(ws) ! ! スペクトルデータに近似的勾配型経度微分 1/(a cosφ)・∂/∂λ ! を作用させる. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_GradLon_ws !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ xyz_GradLon_ws = xya_GradLon_wa(wz_ws(ws))/Ro end function xyz_GradLon_ws function xyz_GradLat_wc(wc) ! ! スペクトルデータに近似的勾配型経度微分 1/a ∂/∂φ を作用させる. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_GradLat_wc !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ xyz_GradLat_wc = xya_GradLat_wa(wz_wc(wc))/Ro end function xyz_GradLat_wc function xyz_GradLat_ws(ws) ! ! スペクトルデータに近似的勾配型経度微分 1/a ∂/∂φ を作用させる. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_GradLat_ws !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ xyz_GradLat_ws = xya_GradLat_wa(wz_ws(ws))/Ro end function xyz_GradLat_ws function wc_DivLon_xyz(xyz) ! ! 格子点データに近似的発散型経度微分 1/(a cosφ)・∂/∂λ を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_DivLon_xyz !(out) 発散型経度微分を作用された cos スペクトルデータ wc_DivLon_xyz = wc_wz(wa_DivLon_xya(xyz/Ro)) end function wc_DivLon_xyz function ws_DivLon_xyz(xyz) ! ! 格子点データに近似的発散型経度微分 1/(a cosφ)・∂/∂λ を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),lm) :: ws_DivLon_xyz !(out) 発散型経度微分を作用された sin スペクトルデータ ws_DivLon_xyz = ws_wz(wa_DivLon_xya(xyz/Ro)) end function ws_DivLon_xyz function wc_DivLat_xyz(xyz) ! ! 格子データに近似的発散型緯度微分 1/(a cosφ)・∂(f cosφ)/∂φ を ! 作用させたスペクトルデータを返す. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_DivLat_xyz !(out) 発散型緯度微分を作用された cos スペクトルデータ wc_DivLat_xyz = wc_wz(wa_divlat_xya(xyz/Ro)) end function wc_DivLat_xyz function ws_DivLat_xyz(xyz) ! ! 格子データに近似的発散型緯度微分 1/(a cosφ)・∂(f cosφ)/∂φ を ! 作用させたスペクトルデータを返す. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),lm) :: ws_DivLat_xyz !(out) 発散型緯度微分を作用された cos スペクトルデータ ws_DivLat_xyz = ws_wz(wa_divlat_xya(xyz/Ro)) end function ws_DivLat_xyz function wc_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad) ! ! べクトル成分である 3 つの格子データに近似的発散を作用させた ! スペクトルデータを返す. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表し, 近似的発散は ! ! 1/(a cosφ)・∂u/∂λ + 1/(a cosφ)・∂(v cosφ)/∂φ + ∂w/∂r ! ! と計算される. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vrad !(in) ベクトル場の動径成分 real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_Div_xyz_xyz_xyz !(out) ベクトル場の発散 wc_Div_xyz_xyz_xyz = wc_DivLon_xyz(xyz_Vlon) & + wc_DivLat_xyz(xyz_Vlat) & + wc_DivRad_ws(ws_xyz(xyz_Vrad)) end function wc_Div_xyz_xyz_xyz function ws_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad) ! ! べクトル成分である 3 つの格子データに近似的発散を作用させた ! スペクトルデータを返す. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表し, 近似的発散は ! ! 1/(a cosφ)・∂u/∂λ + 1/(a cosφ)・∂(v cosφ)/∂φ + ∂w/∂r ! ! と計算される. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vrad !(in) ベクトル場の動径成分 real(8), dimension((nm+1)*(nm+1),lm) :: ws_Div_xyz_xyz_xyz !(out) ベクトル場の発散 ws_Div_xyz_xyz_xyz = ws_DivLon_xyz(xyz_Vlon) & + ws_DivLat_xyz(xyz_Vlat) & + ws_DivRad_wc(wc_xyz(xyz_Vrad)) end function ws_Div_xyz_xyz_xyz function xyz_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad) ! ! ベクトル成分である 3 つの格子データに近似的発散を作用させる. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表す. ! ! 極の特異性を回避するためにベクトル場に cosφの重みをかけて ! 計算している. ! ! div V = (1/cosφ)・div (Vcosφ) + V_φtanφ/a ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vrad !(in) ベクトル場の動径成分 real(8), dimension(0:im-1,1:jm,0:km) :: xyz_Div_xyz_xyz_xyz !(out) ベクトル場の発散 xyz_Div_xyz_xyz_xyz & = 1.0D0/cos(xyz_Lat) & * xyz_wc(wc_Div_xyz_xyz_xyz(xyz_VLon*cos(xyz_Lat), & xyz_VLat*cos(xyz_Lat), & xyz_VRad*cos(xyz_Lat))) & + xyz_VLat*tan(xyz_Lat)/Ro end function xyz_Div_xyz_xyz_xyz !!$ function xyz_RotLon_ws_wc(ws_Vrad,wc_Vlat) !!$ ! !!$ ! ベクトル場の動径成分, 緯度成分である第 1, 2 引数 Vrad, Vlat から !!$ ! 回転の近似的経度成分 !!$ ! !!$ ! 1/a ∂Vrad/∂φ-∂Vlat/∂r を計算する. !!$ ! !!$ ! を計算する !!$ ! !!$ real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws_Vrad !!$ !(in) ベクトル場の動径成分 !!$ !!$ real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc_Vlat !!$ !(in) ベクトル場の緯度成分 !!$ !!$ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_RotLon_ws_wc !!$ !(out) ベクトル場の回転の経度成分 !!$ !!$ xyz_RotLon_ws_wc = xyz_GradLat_ws(ws_Vrad) & !!$ - xyz_ws(ws_RotRad_wc(wc_Vlat)) !!$ !!$ end function xyz_RotLon_ws_wc !!$ !!$ function xyz_RotLat_wc_ws(wc_Vlon,ws_Vrad) !!$ ! !!$ ! ベクトル場の経度成分, 動径成分である第 1, 2 引数 Vlon, Vrad から !!$ ! 回転の緯度成分 !!$ ! !!$ ! 1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ !!$ ! !!$ ! を計算する. !!$ ! !!$ real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc_Vlon !!$ !(in) ベクトル場の経度成分 !!$ !!$ real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws_Vrad !!$ !(in) ベクトル場の動径成分 !!$ !!$ real(8), dimension(0:im-1,1:jm,0:km) :: xyz_RotLat_wc_ws !!$ !(out) ベクトル場の回転の緯度成分 !!$ !!$ xyz_RotLat_wc_ws = xyz_ws(ws_RotRad_wc(wc_Vlon)) & !!$ - xyz_GradLon_ws(ws_Vrad) !!$ !!$ end function xyz_RotLat_wc_ws function wc_RotRad_xyz_xyz(xyz_Vlat,xyz_Vlon) ! ! ベクトルの緯度成分, 経度成分である第 1, 2 引数 Vlat, Vlon に対して ! ベクトル場の近似的回転の動径成分 ! ! 1/(a cosφ)・∂Vlat/∂λ - 1/(a cosφ)・∂(Vlon cosφ)/∂φ ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_Vlon !(in) ベクトル場の経度成分 real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_RotRad_xyz_xyz !(out) ベクトル場の回転の動径成分 wc_RotRad_xyz_xyz = wc_DivLon_xyz(xyz_Vlat) - wc_DivLat_xyz(xyz_Vlon) end function wc_RotRad_xyz_xyz !--------------- 積分計算 ----------------- !----(入力データ xyz)--- function yz_IntLon_xyz(xyz) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)積分. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm,0:km) :: yz_IntLon_xyz !(out) 経度方向(帯状)積分された 2 次元子午面格子点データ integer :: i yz_IntLon_xyz = 0.0d0 do i=0,im-1 yz_IntLon_xyz(:,:) = yz_IntLon_xyz(:,:) & + xyz(i,:,:) * x_Lon_Weight(i) enddo end function yz_IntLon_xyz function xz_IntLat_xyz(xyz) ! ! 3 次元格子点データの緯度方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) cosφ dφ を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,0:km) :: xz_IntLat_xyz !(out) 緯度積分された 2 次元緯度動径格子点データ. ! 緯度円格子点データ integer :: j xz_IntLat_xyz = 0.0d0 do j=1,jm xz_IntLat_xyz(:,:) = xz_IntLat_xyz(:,:) & + xyz(:,j,:) * y_Lat_Weight(j) enddo end function xz_IntLat_xyz function xy_IntRad_xyz(xyz) ! 動径積分 ! ! 3 次元格子点データの動径方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dr を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,1:jm) :: xy_IntRad_xyz !(out) 動径積分された 2 次元経度緯度(水平, 球面)格子点データ integer :: k xy_IntRad_xyz = 0.0d0 do k=0,km xy_IntRad_xyz(:,:) = xy_IntRad_xyz(:,:) & + xyz(:,:,k) * z_Rad_Weight(k) enddo end function xy_IntRad_xyz function x_IntLatRad_xyz(xyz) ! ! 3 次元格子点データの緯度動径(子午面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dφdr ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1) :: x_IntLatRad_xyz !(out) 緯度動径(子午面)積分された 1 次元経度格子点データ integer :: j, k x_IntLatRad_xyz = 0.0D0 do k=0,km do j=1,jm x_IntLatRad_xyz = x_IntLatRad_xyz & + xyz(:,j,k) * y_Lat_Weight(j) * z_Rad_Weight(k) enddo enddo end function x_IntLatRad_xyz function y_IntLonRad_xyz(xyz) ! ! 3 次元格子点データの経度動径(緯度円)積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dλdr を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm) :: y_IntLonRad_xyz !(out) 経度動径(緯度円)積分された 1 次元緯度格子点データ integer :: i, k y_IntLonRad_xyz = 0 do k=0,km do i=0,im-1 y_IntLonRad_xyz = y_IntLonRad_xyz & + xyz(i,:,k) * x_Lon_Weight(i) * z_Rad_Weight(k) enddo enddo end function y_IntLonRad_xyz function z_IntLonLat_xyz(xyz) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:km) :: z_IntLonLat_xyz !(out) 緯度経度(水平, 球面)積分された 1 次元動径格子点データ integer :: i, j z_IntLonLat_xyz = 0 do j=1,jm do i=0,im-1 z_IntLonLat_xyz = z_IntLonLat_xyz & + xyz(i,j,:) * x_Lon_Weight(i) * y_Lat_Weight(j) enddo enddo end function z_IntLonLat_xyz function IntLonLatRad_xyz(xyz) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8) :: IntLonLatRad_xyz !(out) 全球積分値 integer :: i, j, k IntLonLatRad_xyz = 0 do k=0,km do j=1,jm do i=0,im-1 IntLonLatRad_xyz = IntLonLatRad_xyz & + xyz(i,j,k) * x_Lon_Weight(i) & * y_Lat_Weight(j) * z_Rad_Weight(k) enddo enddo enddo end function IntLonLatRad_xyz !----(入力データ yz)--- function z_IntLat_yz(yz) ! 緯度積分 ! ! 2 次元(YZ)格子点データの緯度方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) cosφ dφ を計算する. ! real(8), dimension(jm,0:km), intent(in) :: yz !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(0:km) :: z_IntLat_yz !(out) 緯度積分された 1 次元動径格子点データ integer :: j z_IntLat_yz = 0.0d0 do j=1,jm z_IntLat_yz(:) = z_IntLat_yz(:) + yz(j,:) * y_Lat_Weight(j) enddo end function z_IntLat_yz function y_IntRad_yz(yz) ! 動径積分 ! ! 2 次元(YZ)格子点データの動径方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) r^2dr を計算する. ! real(8), dimension(1:jm,0:km), intent(in) :: yz !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(1:jm) :: y_IntRad_yz !(out) 動径積分された 1 次元緯度格子点データ integer :: k y_IntRad_yz = 0.0d0 do k=0,km y_IntRad_yz(:) = y_IntRad_yz(:) & + yz(:,k) * z_Rad_Weight(k) enddo end function y_IntRad_yz function IntLatRad_yz(yz) ! ! 2 次元(YZ)格子点データの緯度動径積分(子午面)および平均 ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2cosφ dφdr を計算する. ! real(8), dimension(1:jm,0:km), intent(in) :: yz !(in) 2 次元緯度動径(子午面)格子点データ real(8) :: IntLatRad_yz !(out) 積分値 integer :: j, k IntLatRad_yz = 0 do k=0,km do j=1,jm IntLatRad_yz = IntLatRad_yz & + yz(j,k) * y_Lat_Weight(j) * z_Rad_Weight(k) enddo enddo end function IntLatRad_yz !----(入力データ xz)--- function z_IntLon_xz(xz) ! ! 2 次元(XZ)格子点データの経度方向積分. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ を計算する. ! real(8), dimension(0:im-1,0:km), intent(in) :: xz !(in) 2 次元緯度動径格子点データ real(8), dimension(0:km) :: z_IntLon_xz !(out) 経度積分された 1 次元動径格子点データ integer :: i z_IntLon_xz = 0.0d0 do i=0,im-1 z_IntLon_xz(:) = z_IntLon_xz(:) + xz(i,:) * x_Lon_Weight(i) enddo end function z_IntLon_xz function x_IntRad_xz(xz) ! ! 2 次元(XZ)格子点データの動径方向域積分. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r) r^2dr を計算する. ! real(8), dimension(0:im-1,0:km), intent(in) :: xz !(in) 2 次元緯度動径格子点データ real(8), dimension(0:im-1) :: x_IntRad_xz !(out) 動径積分された 1 次元経度格子点データ integer :: k x_IntRad_xz = 0.0d0 do k=0,km x_IntRad_xz(:) = x_IntRad_xz(:) & + xz(:,k) * z_Rad_Weight(k) enddo end function x_IntRad_xz function IntLonRad_xz(xz) ! 経度動径(緯度円)積分 ! ! 2 次元(XZ)格子点データの経度動径積分 ! ! 2 次元データ f(λ,r) に対して∫f(λ,r) r^2dλdr を計算する. ! real(8), dimension(0:im-1,0:km), intent(in) :: xz !(in) 2 次元緯度動径格子点データ real(8) :: IntLonRad_xz !(out) 積分値 integer :: i, k IntLonRad_xz = 0 do k=0,km do i=0,im-1 IntLonRad_xz = IntLonRad_xz & + xz(i,k) * x_Lon_Weight(i) * z_Rad_Weight(k) enddo enddo end function IntLonRad_xz !----(入力データ z)--- function IntRad_z(z) ! 動径積分 ! ! 1 次元(Z)格子点データの動径方向域積分. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr を計算する. ! real(8), dimension(0:km), intent(in) :: z !(in) 1 次元動径格子点データ real(8) :: 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 !--------------- 平均計算 ----------------- !----(入力データ xyz)--- function yz_AvrLon_xyz(xyz) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ/2π を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm,0:km) :: yz_AvrLon_xyz !(out) 経度方向(帯状)平均された 2 次元子午面格子点データ yz_AvrLon_xyz = yz_IntLon_xyz(xyz)/sum(x_Lon_Weight) end function yz_AvrLon_xyz function xz_AvrLat_xyz(xyz) ! 緯度積分 ! ! 3 次元格子点データの緯度方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)cosφ dφ/2 を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,0:km) :: xz_AvrLat_xyz !(out) 緯度平均された 2 次元緯度動径格子点データ xz_AvrLat_xyz = xz_IntLat_xyz(xyz)/sum(y_Lat_Weight) end function xz_AvrLat_xyz function xy_AvrRad_xyz(xyz) ! ! 3 次元格子点データの動径方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dr/((r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,1:jm) :: xy_AvrRad_xyz !(out) 動径平均された 2 次元経度緯度(水平, 球面)格子点データ ! 水平格子点データ xy_AvrRad_xyz = xy_IntRad_xyz(xyz)/sum(z_Rad_Weight) end function xy_AvrRad_xyz function x_AvrLatRad_xyz(xyz) ! 緯度動径(子午面)積分 ! ! 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,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1) :: x_AvrLatRad_xyz !(out) 緯度動径(子午面)平均された 1 次元経度格子点データ x_AvrLatRad_xyz = x_IntLatRad_xyz(xyz) & /( sum(y_Lat_Weight)*sum(z_Rad_Weight) ) end function x_AvrLatRad_xyz function y_AvrLonRad_xyz(xyz) ! 経度動径(緯度円)積分 ! ! 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,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jm) :: y_AvrLonRad_xyz !(out) 経度動径(緯度円)平均された 1 次元緯度格子点データ y_AvrLonRad_xyz = y_IntLonRad_xyz(xyz) & /(sum(x_Lon_Weight)*sum(z_Rad_Weight)) end function y_AvrLonRad_xyz function z_AvrLonLat_xyz(xyz) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ /4π ! ! を計算する. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:km) :: z_AvrLonLat_xyz !(out) 緯度経度(水平, 球面)平均された 1 次元動径格子点データ z_AvrLonLat_xyz = z_IntLonLat_xyz(xyz) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight)) end function z_AvrLonLat_xyz function AvrLonLatRad_xyz(xyz) ! 緯度経度動径(全球)積分 ! ! 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,0:km), intent(in) :: xyz !(in) 3 次元経度緯度動径格子点データ real(8) :: AvrLonLatRad_xyz !(out) 全球平均値 AvrLonLatRad_xyz = IntLonLatRad_xyz(xyz) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(z_Rad_Weight)) end function AvrLonLatRad_xyz !----(入力データ yz)--- function z_AvrLat_yz(yz) ! ! 2 次元(YZ)格子点データの緯度方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) cosφ dφ/2 を計算する. ! real(8), dimension(1:jm,0:km), intent(in) :: yz !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(0:km) :: z_AvrLat_yz !(out) 緯度平均された 1 次元動径格子点データ z_AvrLat_yz = z_IntLat_yz(yz)/sum(y_Lat_Weight) end function z_AvrLat_yz function y_AvrRad_yz(yz) ! ! 2 次元(YZ)格子点データの動径方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2dr /((r[o]^3-r[i]^3)/3) ! を計算する. ! real(8), dimension(1:jm,0:km), intent(in) :: yz !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(1:jm) :: y_AvrRad_yz !(out) 動径平均された 1 次元緯度格子点データ y_AvrRad_yz = y_IntRad_yz(yz)/sum(z_Rad_Weight) end function y_AvrRad_yz function AvrLatRad_yz(yz) ! 緯度動径(子午面)積分 ! ! 2 次元(YZ)格子点データの緯度動径(子午面)平均 ! ! 2 次元データ f(φ,r) に対して ! ! ∫f(φ,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(1:jm,0:km), intent(in) :: yz !(in) 2 次元緯度動径(子午面)格子点データ real(8) :: AvrLatRad_yz !(out) 平均値 AvrLatRad_yz = IntLatRad_yz(yz)/(sum(y_Lat_Weight)*sum(z_Rad_Weight)) end function AvrLatRad_yz !----(入力データ xz)--- function z_AvrLon_xz(xz) ! 経度(帯状)積分 ! ! 2 次元(XZ)格子点データの経度方向平均. ! ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ/2π を計算する. ! real(8), dimension(0:im-1,0:km), intent(in) :: xz !(in) 2 次元緯度動径格子点データ real(8), dimension(0:km) :: z_AvrLon_xz !(out) 経度平均された 1 次元動径格子点データ z_AvrLon_xz = z_IntLon_xz(xz)/sum(x_Lon_Weight) end function z_AvrLon_xz function x_AvrRad_xz(xz) ! 動径積分 ! ! 2 次元(XZ)格子点データの動径方向域平均. ! ! 2 次元データ f(λ,r) に対して ! ! ∫f(λ,r) r^2dr /((r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,0:km), intent(in) :: xz !(in) 2 次元緯度動径格子点データ real(8), dimension(0:im-1) :: x_AvrRad_xz !(out) 動径平均された 1 次元経度格子点データ x_AvrRad_xz = x_IntRad_xz(xz)/sum(z_Rad_Weight) end function x_AvrRad_xz function AvrLonRad_xz(xz) ! 経度動径(緯度円)積分 ! ! 2 次元(XZ)格子点データの経度動径平均 ! ! 2 次元データ f(λ,r) に対して ! ! ∫f(λ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,0:km), intent(in) :: xz ! (in) 2 次元格子点データ real(8) :: AvrLonRad_xz ! 積分値 AvrLonRad_xz = IntLonRad_xz(xz)/(sum(x_Lon_Weight)*sum(z_Rad_Weight)) end function AvrLonRad_xz !----(入力データ z)--- function AvrRad_z(z) ! ! 1 次元(Z)格子点データの動径方向域平均. ! ! 1 次元データ f(r) に対して ∫f(r) r^2dr /((r[o]^3-r[i]^3)/3) を ! 計算する. ! real(8), dimension(0:km), intent(in) :: z !(in) 1 次元動径格子点データ real(8) :: AvrRad_z !(out) 平均値 AvrRad_z = IntRad_z(z)/sum(z_Rad_Weight) end function AvrRad_z !--------------- ポロイダル/トロイダルモデル用微分 ----------------- function wc_KxRGrad_wc(wc) ! ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_KxRGrad_wc !(out) 経度微分を作用された 2 次元スペクトルデータ wc_KxRGrad_wc = wa_Dlon_wa(wc) end function wc_KxRGrad_wc function ws_KxRGrad_ws(ws) ! ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_KxRGrad_ws !(out) 経度微分を作用された 2 次元スペクトルデータ ws_KxRGrad_ws = wa_Dlon_wa(ws) end function ws_KxRGrad_ws function wc_L2_wc(wc) ! ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する. ! ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. ! 入力スペクトルデ ータに対応する格子点データに演算子 ! ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用させたデータのスペクトル変換が返される. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_L2_wc !(out) L^2 演算子を作用された 2 次元スペクトルデータ wc_L2_wc = -wa_Lapla_wa(wc) end function wc_L2_wc function ws_L2_ws(ws) ! ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する. ! ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. ! 入力スペクトルデ ータに対応する格子点データに演算子 ! ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用させたデータのスペクトル変換が返される. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_L2_ws !(out) L^2 演算子を作用された 2 次元スペクトルデータ ws_L2_ws = -wa_Lapla_wa(ws) end function ws_L2_ws function wc_L2Inv_wc(wc) ! ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を ! 作用する. ! ! スペクトルデータに L^2 演算子を作用させる関数 wc_L2_wc の逆計算を ! 行う関数である. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_L2Inv_wc !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ wc_L2Inv_wc = -wa_LaplaInv_wa(wc) end function wc_L2Inv_wc function ws_L2Inv_ws(ws) ! ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を ! 作用する. ! ! スペクトルデータに L^2 演算子を作用させる関数 ws_L2_ws の逆計算を ! 行う関数である. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension((nm+1)*(nm+1),lm) :: ws_L2Inv_ws !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ ws_L2Inv_ws = -wa_LaplaInv_wa(ws) end function ws_L2Inv_ws function wc_RadRot_xyz_xyz(xyz_VLON,xyz_VLAT) ! r・(▽×v) ! ! ベクトルの渦度と動径ベクトルの内積 r・(▽×v) を計算する. ! ! 第 1, 2 引数(v[λ], v[φ])がそれぞれベクトルの経度成分, 緯度成分を表す. ! ! r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ ! ! のスペクトル データが返される. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_VLAT !(in) ベクトルの緯度成分 real(8), dimension((nm+1)*(nm+1),0:lm) :: wc_RadRot_xyz_xyz !(out) ベクトルの渦度と動径ベクトルの内積 wc_RadRot_xyz_xyz = wc_wz(wa_DivLon_xya(xyz_VLAT) & - wa_DivLat_xya(xyz_VLON)) end function wc_RadRot_xyz_xyz function ws_RadRotRot_xyz_xyz_xyz(xyz_VLON,xyz_VLAT,xyz_VRAD) ! ! ベクトル v に対して近似的な r・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分, ! 緯度成分, 動径成分を表す. ! ! r・(▽×▽×v) = ∂/∂r (( 1/cosφ・∂v[λ]/∂λ ! + 1/cosφ・∂(v[φ] cosφ)/∂φ ) ) ! + L^2 v[r]/a ! ! のスペクトルデータが返される. ! real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_VLAT !(in) ベクトルの緯度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(in) :: xyz_VRAD !(in) ベクトルの動径成分 real(8), dimension((nm+1)*(nm+1),lm) :: ws_RadRotRot_xyz_xyz_xyz !(out) ベクトル v の r・(▽×▽×v) ws_RadRotRot_xyz_xyz_xyz = & ws_RotRad_wc(wc_wz( & (wa_DivLon_xya(xyz_VLON)+ wa_DivLat_xya(xyz_VLAT)))) & + ws_L2_ws(ws_xyz(xyz_VRAD/Ro)) end function ws_RadRotRot_xyz_xyz_xyz subroutine wsc_Potential2Vector(& xyz_VLON,xyz_VLAT,xyz_VRAD,wc_TORPOT,ws_POLPOT) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψr) + ▽x▽x(Φr) ! ! の各成分を計算する ! real(8), dimension(0:im-1,1:jm,0:km) :: xyz_VLON !(out) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jm,0:km) :: xyz_VLAT !(out) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jm,0:km) :: xyz_VRAD !(out) ベクトル場の動径成分 real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc_TORPOT !(in) トロイダルポテンシャル real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws_POLPOT !(in) ポロイダルポテンシャル xyz_VLON = Ro * xyz_GradLat_wc(wc_TORPOT) & + xya_GradLon_wa(wz_wc(wc_RotRad_ws(ws_POLPOT))) xyz_VLAT = - Ro * xyz_GradLon_wc(wc_TORPOT) & + xya_GradLat_wa(wz_wc(wc_RotRad_ws(ws_POLPOT))) xyz_VRAD = xyz_ws(ws_L2_ws(ws_POLPOT))/Ro end subroutine wsc_Potential2Vector subroutine wsc_Potential2Rotation(& xyz_RotVLON,xyz_RotVLAT,xyz_RotVRAD,wc_TORPOT,ws_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,1:jm,0:km), intent(OUT) :: xyz_RotVLON !(out) 回転の経度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(OUT) :: xyz_RotVLAT !(out) 回転の緯度成分 real(8), dimension(0:im-1,1:jm,0:km), intent(OUT) :: xyz_RotVRAD !(out) 回転の動径成分 ! 入力ベクトル場を表すポテンシャル real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc_TORPOT !(in) トロイダルポテンシャル real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws_POLPOT !(in) ポロイダルポテンシャル xyz_RotVLON = Ro * xyz_GradLat_ws(-ws_Lapla_ws(ws_POLPOT)) & + xya_GradLon_wa(wz_ws(ws_RotRad_wc(wc_TORPOT))) xyz_RotVLAT = - Ro * xyz_GradLon_ws(-ws_Lapla_ws(ws_POLPOT)) & + xya_GradLat_wa(wz_ws(ws_RotRad_wc(wc_TORPOT))) xyz_RotVRAD = xyz_wc(wc_L2_wc(wc_TORPOT))/Ro end subroutine wsc_Potential2Rotation !--------------- 補間計算 ----------------- function Interpolate_wc(wc_data,alon,alat,arad) ! ! 緯度 alon, 経度 alat 動径 arad における関数値を ! その球面調和変換係数 wa_data から補間計算する ! real(8), intent(IN) :: wc_data((nm+1)**2,0:lm) ! スペクトルデータ real(8), intent(IN) :: alon ! 補間する位置(経度) real(8), intent(IN) :: alat ! 補間する位置(緯度) real(8), intent(IN) :: arad ! 補間する位置(動径) real(8) :: Interpolate_wc ! 補間した値 Interpolate_wc = & Interpolate_w(a_Interpolate_ac(wc_data,arad),alon,alat) end function Interpolate_wc function Interpolate_ws(ws_data,alon,alat,arad) ! ! 緯度 alon, 経度 alat 動径 arad における関数値を ! その球面調和変換係数 wa_data から補間計算する ! real(8), intent(IN) :: ws_data((nm+1)**2,lm) ! スペクトルデータ real(8), intent(IN) :: alon ! 補間する位置(経度) real(8), intent(IN) :: alat ! 補間する位置(緯度) real(8), intent(IN) :: arad ! 補間する位置(動径) real(8) :: Interpolate_ws ! 補間した値 Interpolate_ws = & Interpolate_w(a_Interpolate_as(ws_data,arad),alon,alat) end function Interpolate_ws !--------------- ポロイダル/トロイダルモデル用スペクトル解析 ---------------- function nmz_ToroidalEnergySpectrum_wc(wc_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する ! ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(n,m,r)から全波数 n, 帯状波数 m 成分のトロイダルエネルギー ! スペクトルは (1/2)n(n+1)4πa^2|ψ(n,m,r)|^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和を動径積分したもの(a^2の重み無し) ! が球殻内での全エネルギーに等しい. ! ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される. ! wsc_VMiss によって設定できる (初期値は -999.0) ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc_TORPOT !(in) トロイダルポテンシャル real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_ToroidalEnergySpectrum_wc !(out) エネルギースペクトルトロイダル成分 real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA ! 作業領域 integer :: n, m nmz_ToroidalEnergySpectrum_wc = wsc_VMiss wz_DATA = wz_wc(wc_TORPOT) do n=0,nm nmz_ToroidalEnergySpectrum_wc(n,0,:) & = 0.5 * n*(n+1)* (4*pi) * Ro**2 & * wz_DATA(l_nm(n,0),:)**2 do m=1,n nmz_ToroidalEnergySpectrum_wc(n,m,:) & = 0.5 * n*(n+1)* (4*pi) * Ro**2 & * (wz_DATA(l_nm(n,m),:)**2+wz_DATA(l_nm(n,-m),:)**2) nmz_ToroidalEnergySpectrum_wc(n,-m,:) & = nmz_ToroidalEnergySpectrum_wc(n,m,:) enddo enddo end function nmz_ToroidalEnergySpectrum_wc function nz_ToroidalEnergySpectrum_wc(wc_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 球面調和函数全波数の各成分を計算する. ! ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは ! Σ[m=-n]^n(1/2)n(n+1)4πa^2|ψ(n,m,r)|^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和を動径積分したもの(a^2の重み無し) ! が球殻内での全エネルギーに等しい. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wc_TORPOT !(in) トロイダルポテンシャル real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wc !(out) エネルギースペクトルトロイダル成分 real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA ! 作業領域 integer :: n, m wz_DATA = wz_wc(wc_TORPOT) do n=0,nm nz_ToroidalEnergySpectrum_wc(n,:) & = 0.5 * n*(n+1)* (4*pi) * Ro**2 * wz_Data(l_nm(n,0),:)**2 do m=1,n nz_ToroidalEnergySpectrum_wc(n,:) & = nz_ToroidalEnergySpectrum_wc(n,:) & + 0.5 * n*(n+1)* (4*pi) * Ro**2 & * 2* (wz_Data(l_nm(n,m),:)**2 + wz_Data(l_nm(n,-m),:)**2) enddo enddo end function nz_ToroidalEnergySpectrum_wc function nmz_PoloidalEnergySpectrum_ws(ws_POLPOT) ! ! ポロイダルポテンシャルから, ポロイダルエネルギーの ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する. ! ! * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分 ! φ(n,m,r)から全波数 n, 帯状波数 m 成分のポロイダルエネルギー ! スペクトルは ! ! (1/2)n(n+1)4π{|d(aφ(n,m,r))/dr|^2 + n(n+1)|φ(n,m,r)|^2} ! ! と計算される. ! ! * 全てのエネルギースペクトル成分の和を動径積分したもの(a^2の重み無し) ! が球殻内での全エネルギーに等しい. ! ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される. ! 欠損値の値はモジュール変数 wsc_VMiss によって設定できる ! (初期値は -999.0) ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws_POLPOT !(in) ポロイダルポテンシャル real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_PoloidalEnergySpectrum_ws !(out) エネルギースペクトルポロイダル成分 real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1 ! 作業領域 real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2 ! 作業領域 integer :: n, m nmz_PoloidalEnergySpectrum_ws = wsc_VMiss wz_Data1 = wz_ws(ws_POLPOT) wz_Data2 = Ro*wz_wc(wc_DRad_ws(ws_POLPOT)) ! d(rφ)/dr ~ a dφ/dr do n=0,nm nmz_PoloidalEnergySpectrum_ws(n,0,:) = & + 0.5* n*(n+1)* (4*pi) & *( wz_Data2(l_nm(n,0),:)**2 & + n*(n+1)*wz_Data1(l_nm(n,0),:)**2 ) do m=1,n nmz_PoloidalEnergySpectrum_ws(n,m,:) = & + 0.5* n*(n+1)* (4*pi) & *( wz_Data2(l_nm(n,m),:)**2 + wz_Data2(l_nm(n,-m),:)**2 & + n*(n+1)* ( wz_Data1(l_nm(n,m),:)**2 + wz_Data1(l_nm(n,-m),:)**2)) nmz_PoloidalEnergySpectrum_ws(n,-m,:) = & nmz_PoloidalEnergySpectrum_ws(n,m,:) enddo enddo end function nmz_PoloidalEnergySpectrum_ws function nz_PoloidalEnergySpectrum_ws(ws_POLPOT) ! ! ポロイダルポテンシャルから, ポロイダルエネルギーの ! 球面調和函数全波数の各成分を計算する ! ! * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分 ! φ(n,m,r)から全波数 n 成分のポロイダルエネルギースペクトルは ! ! Σ[m=-n]^n ((1/2)n(n+1)4π{|d(rφ(n,m,r))/dr|^2 ! + n(n+1)|φ(n,m,r)|^2} ! ! と計算される. ! ! * 全ての全波数に対してのエネルギースペクトル成分の和を動径積分したもの ! (r^2の重み無し)が球殻内での全エネルギーに等しい. ! real(8), dimension((nm+1)*(nm+1),lm), intent(in) :: ws_POLPOT !(in) ポロイダルポテンシャル real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_ws !(out) エネルギースペクトルポロイダル成分 real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1 ! 作業領域 real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2 ! 作業領域 integer :: n, m wz_Data1 = wz_ws(ws_POLPOT) wz_Data2 = Ro * wz_wc(wc_DRad_ws(ws_POLPOT)) ! d(rφ)/dr ~ a dφ/dr do n=0,nm nz_PoloidalEnergySpectrum_ws(n,:) & = 0.5* n*(n+1)* (4*pi) & *( wz_Data2(l_nm(n,0),:)**2 + n*(n+1)*wz_Data1(l_nm(n,0),:)**2 ) do m=1,n nz_PoloidalEnergySpectrum_ws(n,:) & = nz_PoloidalEnergySpectrum_ws(n,:) & + 2 * 0.5* n*(n+1)* (4*pi) & *( wz_Data2(l_nm(n,m),:)**2 + wz_Data2(l_nm(n,-m),:)**2 & + n*(n+1)*(wz_Data1(l_nm(n,m),:)**2 +wz_Data1(l_nm(n,-m),:)**2)) enddo enddo end function nz_PoloidalEnergySpectrum_ws end module wsc_module_svpack