!-- !---------------------------------------------------------------------- ! Copyright(c) 2017 SPMDODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 scee_module_fftj ! ! spml/scee_module_fftj モジュールは平行平板間での 3 次元流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向にフーリエ数変換および上下の境界壁を扱うための ! sin/cos 変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ee_module_fftj, asc_module を用いている. 最下部ではフーリエ ! およびチェビシェフ変換のエンジンとして ISPACK の Fortran77 ! サブルーチンを用いている. ! ! !履歴 2017/05/26 竹広真一 作成 ! !凡例 ! データ種類と index ! x : 水平(X) y : 水平(Y) z : 鉛直 ! e : フーリエ変換スペクトル ! l : フーリエ関数スペクトル(X 方向波数) ! m : フーリエ関数スペクトル(Y 方向波数) ! s : sin 関数スペクトル ! c : cos 関数スペクトル ! a : 任意の次元 ! ! zyx : 3 次元格子点データ ! yx : 水平 2 次元格子点データ ! zy : 鉛直 2 次元格子点データ ! zx : 鉛直 2 次元格子点データ ! ! zee : 水平スペクトル鉛直格子点データ ! see : スペクトルデータ鉛直 sin 変換格子 ! cee : スペクトルデータ鉛直 cos 変換格子 ! !++ module scee_module_fftj ! != scee_module_fftj ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/scee_module_fftj モジュールは平行平板間での 3 次元流体運動を ! スペクトル法によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向にフーリエ数変換および上下の境界壁を扱うための ! sin/cos 変換を用いる場合のスペクトル計算のためのさまざまな ! 関数を提供する. ! ! 内部で ee_module_fftj, asc_module を用いている. 最下部ではフーリエ ! およびチェビシェフ変換のエンジンとして ISPACK の Fortran77 ! サブルーチンを用いている. ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (see_, cee_, zyx_, zee_, ee_, yx_, x_, y_, z_, a_) は, ! 返す値の形を示している. ! see_ :: スペクトルデータ(2 重フーリエ・sin変換) ! cee_ :: スペクトルデータ(2 重フーリエ・cos変換) ! zyx_ :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・) ! zee_ :: 水平スペクトル, 鉛直格子点データ ! e2a_ :: 1 次元化した水平スペクトル, 任意座標データ ! aee_ :: 任意座標データ, 水平スペクトルデータ ! ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 (see_, cee_, zyx_, zee_, ee_, yx_, x_, y_, z_, a_) は, ! 入力変数の形がスペクトルデータおよび格子点データであることを示している. ! _see :: スペクトルデータ(2 重フーリエ・sin変換) ! _cee :: スペクトルデータ(2 重フーリエ・cos変換) ! _zyx :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・) ! _zee :: 水平スペクトル, 鉛直格子点データ ! _e2a :: 1 次元化した水平スペクトル, 任意座標データ ! _aee :: 任意座標データ, 水平スペクトルデータ ! !=== 各データの種類の説明 ! ! * zyx : 3 次元格子点データ(鉛直, 水平 2 次元) ! * 変数の種類と次元は real(8), dimension(0:km,0:jm-1,0:im-1). ! * im, jm, km はそれぞれ水平 X, Y, 鉛直 Z 座標の格子点数であり, ! サブルーチン tee_Initial にてあらかじめ設定しておく. ! ! * see : スペクトルデータ ! * 変数の種類と次元は real(8), dimension(nm,-mm:mm,-lm:lm). ! * lm, mm は X,Y 方向最大波数, nm は sin 級数の最大次数 ! であり, サブルーチン tee_Initial にてあらかじめ設定しておく. ! ! * cee : スペクトルデータ ! * 変数の種類と次元は real(8), dimension(0:nm,-mm:mm,-lm:lm). ! * lm, mm は X,Y 方向最大波数, nm は cos 級数の最大次数 ! であり, サブルーチン tee_Initial にてあらかじめ設定しておく. ! ! * zee : 水平スペクトル, 鉛直格子点データ. ! * 変数の種類と次元は real(8), dimension(0:km,-mm:mm,-lm:lm). ! ! * e2a : 1 次元化した水平スペクトル, 任意座標データ ! * 変数の種類と次元は real(8), dimension((2*mm+1)*(2*lm+1),:) ! ! * aee : 任意座標データ, 水平スペクトルデータ ! * 変数の種類と次元は real(8), dimension(:,-mm:mm,-lm:lm). ! ! * see_, cee_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * zyx_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * zee_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! scee_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_X, y_Y, z_Z :: 格子点座標(水平 X,Y, 鉛直 Z 座標)を ! 格納した1 次元配列 ! x_X_Weight, y_X_Weight, z_Z_Weight :: 重み座標を格納した 1 次元配列 ! zyx_X, zyx_Y, zyx_Z :: 格子点データの水平鉛直座標(X,Y,Z) ! (格子点データ型 3 次元配列) ! yx_X, yx_Y :: 格子点データの水平座標(X,Y) ! zy_Z, zy_Y :: 格子点データの鉛直水平座標(Z,Y) ! zx_Z, zx_X :: 格子点データの鉛直水平座標(Z,X) ! (格子点データ型 2 次元配列) ! !==== 基本変換 ! ! zyx_see, see_zyx :: スペクトルデータと 3 次元格子データの間の変換 ! (2 重フーリエ, sin 変換) ! zyx_cee, cee_zyx :: スペクトルデータと 3 次元格子データの間の変換 ! (2 重フーリエ, cos 変換) ! zyx_zee, zee_zyx :: 3 次元格子データと水平スペクトル・鉛直格子データとの間 ! の変換 (2 重フーリエ変換) ! zee_see, see_zee :: スペクトルデータと水平スペクトル・鉛直格子データとの間 ! の変換 (sin変換) ! zee_cee, cee_zee :: スペクトルデータと水平スペクトル・鉛直格子データとの間 ! の変換 (cos変換) ! ee_yx, yx_ee :: スペクトルデータと 2 次元水平格子データの間の変換 ! (2 重フーリエ変換) ! az_as, as_az :: 格子データと sin 変換データの間の変換を同時に複数個行う ! az_ac, ac_az :: 格子データと cos 変換データの間の変換を同時に複数個行う ! e2a_aee, aee_e2a :: 水平スペクトル軸を 1 次元化し転置, 転置 ２次元化する. ! !==== 微分 ! ! see_Dx_see :: スペクトルデータに動径微分∂/∂x を作用させる ! cee_Dx_cee :: スペクトルデータに動径微分∂/∂x を作用させる ! see_Dy_see :: スペクトルデータに動径微分∂/∂y を作用させる ! cee_Dy_cee :: スペクトルデータに動径微分∂/∂y を作用させる ! see_Dz_cee :: スペクトルデータに動径微分∂/∂z を作用させる ! cee_Dz_see :: スペクトルデータに動径微分∂/∂z を作用させる ! ! see_Lapla_see :: スペクトルデータにラプラシアンを作用させる ! cee_Lapla_cee :: スペクトルデータにラプラシアンを作用させる ! see_LaplaInv_see :: スペクトルデータに逆ラプラシアンを作用させる ! cee_LaplaInv_cee :: スペクトルデータに逆ラプラシアンを作用させる ! see_LaplaH_see :: スペクトルデータに水平ラプラシアンを作用させる ! cee_LaplaH_cee :: スペクトルデータに水平ラプラシアンを作用させる ! see_LaplaHInv_see :: スペクトルデータに逆水平ラプラシアンを作用させる ! cee_LaplaHInv_cee :: スペクトルデータに逆水平ラプラシアンを作用させる ! ! see_Div_zyx_zyx_zyx :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! cee_Div_zyx_zyx_zyx :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! !==== トロイダルポロイダル計算用微分 ! ! see_ZRot_zyx_zyx :: ベクトル v の渦度と動径ベクトル r の内積 ! z・(▽×v) を計算する ! cee_ZRot_zyx_zyx :: ベクトル v の渦度と動径ベクトル r の内積 ! z・(▽×v) を計算する ! see_ZRotRot_zyx_zyx_zyx :: ベクトルの v の z・(▽×▽×v) を計算する ! cee_ZRotRot_zyx_zyx_zyx :: ベクトルの v の z・(▽×▽×v) を計算する ! scee_Potential2Vector_see_cee :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! scee_Potential2Vector_cee_see :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! scee_Potential2Rotation_see_cee :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! scee_Potential2Rotation_cee_see :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! zee_ToroidalEnergySpectrum_cee, zee_ToroidalEnergySpectrum_see :: ! トロイダルポテンシャルからエネルギーのフーリエ各成分を計算する ! zee_PoloidalEnergySpectrum_tee, zee_PoloidalEnergySpectrum_cee :: ! ポロイダルポテンシャルからエネルギーのフーリエ各成分を計算する ! !==== 積分・平均(3 次元データ) ! ! IntZYX_zyx, AvrZYX_zyx :: 3 次元格子点データの全領域積分および平均 ! z_IntYX_zyx, z_AvrYX_zyx :: 3 次元格子点データの水平積分および平均 ! y_IntZX_zyx, y_AvrZX_zyx :: 3 次元格子点データのZX積分および平均 ! z_IntZY_zyx, z_AvrZY_zyx :: 3 次元格子点データのZY積分および平均 ! zy_IntX_zyx, zy_AvrX_zyx :: 3 次元格子点データの水平X方向積分および平均 ! zx_IntY_zyx, zx_AvrY_zyx :: 3 次元格子点データの水平Y方向積分および平均 ! zx_IntZ_zyx, zx_AvrZ_zyx :: 3 次元格子点データの鉛直方向積分および平均 ! !==== 積分・平均(2 次元データ) ! ! IntYX_yx, AvrYX_yx :: 2 次元格子点データの水平積分および平均 ! IntZX_zx, AvrZX_zx :: 2 次元(ZX)格子点データのZX積分および平均 ! IntZY_zy, AvrZY_zy :: 2 次元(ZY)格子点データのZY積分および平均 ! y_IntX_yx, y_AvrX_yx :: 水平 2 次元格子点データのX方向積分および平均 ! x_IntY_yx, x_AvrY_yx :: 水平2 次元格子点データのY方向積分および平均 ! z_IntX_zx, z_AvrX_zx :: 2 次元(ZX)格子点データのX方向積分および平均 ! x_IntZ_zx, x_AvrZ_zx :: 2 次元(ZX)格子点データのZ方向積分および平均 ! z_IntY_zy, z_AvrY_zy :: 2 次元(ZY)格子点データのY方向積分および平均 ! y_IntZ_zy, y_AvrZ_zy :: 2 次元(ZY)格子点データのZ方向積分および平均 ! !==== 積分・平均(1 次元データ) ! ! IntX_x, AvrX_x :: 1 次元(X)格子点データのX方向積分および平均 ! IntY_y, AvrY_y :: 1 次元(Y)格子点データのY方向積分および平均 ! IntZ_z, AvrZ_z :: 1 次元(Z)格子点データのZ方向積分および平均 ! !==== 補間計算 ! ! Interpolate_cee :: スペクトルデータから任意の点の値を補間する. ! Interpolate_see :: スペクトルデータから任意の点の値を補間する. ! use dc_message use lumatrix use ee_module_fftj use asc_module, z_Z => g_X, z_Z_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, & c_Dz_s => c_Dx_s, ac_Dz_as => ac_Dx_as, & s_Dz_c => s_Dx_c, as_Dz_ac => as_Dx_ac implicit none private public scee_Initial public x_X, x_X_Weight public y_Y, y_Y_Weight public z_Z, z_Z_Weight public yx_X, yx_Y, zy_Z, zy_Y, zx_X, zx_Z public zyx_X, zyx_Y, zyx_Z public zee_Z public scee_VMiss public ee_yx, yx_ee public ac_Dz_as, as_Dz_ac, c_Dz_s, s_Dz_c public az_ac, az_as, ac_az, as_az, z_c, z_s, c_z, s_z public zyx_see, see_zyx, zee_see, see_zee public zyx_cee, cee_zyx, zee_cee, cee_zee public zyx_zee, zee_zyx public e2z_e2s, e2s_e2z public e2z_e2c, e2c_e2z public aee_e2a, e2a_aee, ee_e2, e2_ee public cee_Dx_cee, see_Dx_see, cee_Dy_cee, see_Dy_see, see_Dz_cee, cee_Dz_see public cee_Lapla_cee, see_Lapla_see public cee_LaplaInv_cee, see_LaplaInv_see public cee_LaplaH_cee, see_LaplaH_see, cee_LaplaHInv_cee, see_LaplaHInv_see public cee_Div_zyx_zyx_zyx, see_Div_zyx_zyx_zyx public zy_IntX_zyx, zx_IntY_zyx, yx_IntZ_zyx public x_IntZY_zyx, y_IntZX_zyx, z_IntYX_zyx public IntZYX_zyx public x_IntY_yx, y_IntX_yx, IntYX_yx public z_IntY_zy, y_IntZ_zy, IntZY_zy public z_IntX_zx, x_IntZ_zx, IntZX_zx public IntX_x, IntY_y, IntZ_z public zy_AvrX_zyx, zx_AvrY_zyx, yx_AvrZ_zyx public x_AvrZY_zyx, y_AvrZX_zyx, z_AvrYX_zyx public AvrZYX_zyx public x_AvrY_yx, y_AvrX_yx, AvrYX_yx public z_AvrY_zy, y_AvrZ_zy, AvrZY_zy public z_AvrX_zx, x_AvrZ_zx, AvrZX_zx public AvrX_x, AvrY_y, AvrZ_z public cee_ZRot_zyx_zyx, see_ZRot_zyx_zyx public cee_ZRotRot_zyx_zyx_zyx, see_ZRotRot_zyx_zyx_zyx public scee_Potential2vector_cee_see, scee_Potential2vector_see_cee public scee_Potential2Rotation_cee_see, scee_Potential2Rotation_see_cee public Interpolate_cee, Interpolate_see public zee_ToroidalEnergySpectrum_cee, zee_ToroidalEnergySpectrum_see ! nz_ToroidalEnergySpectrum_wt public zee_PoloidalEnergySpectrum_cee, zee_PoloidalEnergySpectrum_see ! nz_PoloidalEnergySpectrum_wt integer :: im=32, jm=32, km=16 ! 格子点の設定(X,Y,Z) integer :: lm=10, mm=10, nm=16 ! 切断波数の設定(水平X,Y, 鉛直Z) real(8) :: xl, yl, zl ! 領域の大きさ(水平X,Y, 鉛直Z) real(8), parameter :: pi=3.1415926535897932385D0 real(8), dimension(:,:,:), allocatable :: zyx_X, zyx_Y, zyx_Z ! 座標 real(8), dimension(:,:), allocatable :: zy_Z, zy_Y ! 座標 real(8), dimension(:,:), allocatable :: zx_Z, zx_X ! 座標 real(8), dimension(:,:,:), allocatable :: zee_Z ! 座標 real(8) :: scee_VMiss = -999.0 ! 欠損値 save im, jm, km, lm, mm, nm, xl, yl, zl contains !--------------- 初期化 ----------------- subroutine scee_Initial(i,j,k,l,m,n,xmin,xmax,ymin,ymax,zmin,zmax) ! ! スペクトル変換の格子点数, 波数, 各座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! integer,intent(in) :: i ! 格子点数(水平X) integer,intent(in) :: j ! 格子点数(水平Y) integer,intent(in) :: k ! 格子点数(鉛直Z) integer,intent(in) :: l ! 切断波数(水平X波数) integer,intent(in) :: m ! 切断波数(水平Y波数) integer,intent(in) :: n ! 切断波数(鉛直Z波数) real(8),intent(in) :: xmin, xmax ! X 方向領域 real(8),intent(in) :: ymin, ymax ! Y 方向領域 real(8),intent(in) :: zmin, zmax ! Z 方向領域 integer :: id im = i ; jm = j ; km = k lm = l ; mm = m ; nm = n xl = xmax - xmin yl = ymax - ymin zl = zmax - zmin call ee_Initial(im,jm,lm,mm,xmin,xmax,ymin,ymax,id) call asc_Initial(km,nm,zmin,zmax) allocate(zyx_X(0:km,0:jm-1,0:im-1)) allocate(zyx_Y(0:km,0:jm-1,0:im-1)) allocate(zyx_Z(0:km,0:jm-1,0:im-1)) allocate(zy_Z(0:km,0:jm-1)) allocate(zy_Y(0:km,0:jm-1)) allocate(zx_Z(0:km,0:im-1)) allocate(zx_X(0:km,0:im-1)) allocate(zee_Z(0:km,-mm:mm,-lm:lm)) zyx_X = spread(yx_X,1,km+1) zyx_Y = spread(yx_Y,1,km+1) zyx_Z = spread(spread(z_Z,2,jm),3,im) zy_Z = spread(z_Z,2,jm) zy_Y = spread(y_Y,1,km+1) zx_Z = spread(z_Z,2,im) zx_X = spread(x_X,1,km+1) zee_Z = spread(spread(z_Z,2,2*mm+1),3,2*lm+1) call MessageNotify('M','tee_initial', & 'scee_module_fftj (2017/05/28) is initialized') end subroutine scee_Initial !--------------- 基本変換 ----------------- function zyx_see(see) ! ! スペクトルデータから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 2 重フーリエ sin スペクトルデータ real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_see !(out) 3 次元格子点データ zyx_see = zyx_zee(zee_see(see)) end function zyx_see function see_zyx(zyx) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_zyx !(out) 2 重フーリエ sin スペクトルデータ see_zyx = see_zee(zee_zyx(zyx)) end function see_zyx function zyx_cee(cee) ! ! スペクトルデータから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 2 重フーリエ cos スペクトルデータ real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_cee !(out) 3 次元格子点データ zyx_cee = zyx_zee(zee_cee(cee)) end function zyx_cee function cee_zyx(zyx) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_zyx !(out) 2 重フーリエ cos スペクトルデータ cee_zyx = cee_zee(zee_zyx(zyx)) end function cee_zyx function zyx_zee(zee) ! ! 水平スペクトル・鉛直格子点データから 3 次元格子点データへ(逆)変換する. ! real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee !(in) 2 次元水平スペクトル・鉛直格子点データ real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_zee !(out) 3 次元格子点データ real(8), dimension(-mm:mm,-lm:lm) :: ee real(8), dimension(0:jm-1,0:im-1) :: yx integer :: k do k = 0, km !zyx_zee(k,:,:) = yx_ee(zee(k,:,:)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! ee = zee(k,-mm:mm,-lm:lm) yx = yx_ee(ee) zyx_zee(k,:,:) = yx enddo end function zyx_zee function zee_zyx(zyx) ! ! 3 次元格子データから水平スペクトル・鉛直格子点データへ(正)変換する. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_zyx !(out) 2 次元スペクトル・鉛直格子点データ real(8), dimension(-mm:mm,-lm:lm) :: ee real(8), dimension(0:jm-1,0:im-1) :: yx integer :: k do k = 0, km !zee_zyx(k,:,:) = ee_yx(zyx(k,:,:)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! yx = zyx(k,0:jm-1,0:im-1) ee = ee_yx(yx) zee_zyx(k,:,:) = ee enddo end function zee_zyx function zee_see(see) ! ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 2 次元水平スペクトル鉛直 sin スペクトルデータ real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_see !(out) 2 次元水平スペクトル・鉛直格子点データ zee_see = aee_e2a(e2z_e2s(e2a_aee(see))) end function zee_see function see_zee(zee) ! ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee !(in) 2 次元水平スペクトル・鉛直格子点データ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_zee !(out) 2 次元水平鉛直 sin スペクトルデータ see_zee = aee_e2a(e2s_e2z(e2a_aee(zee))) end function see_zee function zee_cee(cee) ! ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 2 次元水平スペクトル鉛直 cos スペクトルデータ real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_cee !(out) 2 次元水平スペクトル・鉛直格子点データ zee_cee = aee_e2a(e2z_e2c(e2a_aee(cee))) end function zee_cee function cee_zee(zee) ! ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee !(in) 2 次元水平スペクトル・鉛直格子点データ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_zee !(out) 2 次元水平鉛直 cos スペクトルデータ cee_zee = aee_e2a(e2c_e2z(e2a_aee(zee))) end function cee_zee function e2z_e2s(e2s) ! ! 鉛直スペクトルを格子点に変換する ! real(8), dimension((2*lm+1)*(2*mm+1),nm) :: e2s !(in) 1 次元化された水平スペクトル・鉛直スペクトル real(8), dimension((2*lm+1)*(2*mm+1),0:km) :: e2z_e2s !(out) 1 次元化された水平スペクトル・鉛直格子座標 e2z_e2s = az_as(e2s) end function e2z_e2s function e2s_e2z(e2z) ! ! 鉛直格子点をスペクトルに変換する ! real(8), dimension((2*lm+1)*(2*mm+1),0:km) :: e2z !(in) 1 次元化された水平スペクトル・鉛直格子座標 real(8), dimension((2*lm+1)*(2*mm+1),nm) :: e2s_e2z !(out) 1 次元化された水平スペクトル・鉛直スペクトル e2s_e2z = as_az(e2z) end function e2s_e2z function e2z_e2c(e2c) ! ! 鉛直スペクトルを格子点に変換する ! real(8), dimension((2*lm+1)*(2*mm+1),0:nm) :: e2c !(in) 1 次元化された水平スペクトル・鉛直スペクトル real(8), dimension((2*lm+1)*(2*mm+1),0:km) :: e2z_e2c !(out) 1 次元化された水平スペクトル・鉛直格子座標 e2z_e2c = az_ac(e2c) end function e2z_e2c function e2c_e2z(e2z) ! ! 鉛直格子点をスペクトルに変換する ! real(8), dimension((2*lm+1)*(2*mm+1),0:km) :: e2z !(in) 1 次元化された水平スペクトル・鉛直格子座標 real(8), dimension((2*lm+1)*(2*mm+1),0:nm) :: e2c_e2z !(out) 1 次元化された水平スペクトル・鉛直スペクトル e2c_e2z = ac_az(e2z) end function e2c_e2z function e2a_aee(aee) ! ! 水平スペクトルを 1 次元化転置する. ! real(8), dimension(:,-mm:,-lm:), intent(in) :: aee !(in) 任意座標・水平 2 次元スペクトルデータ dimension(:,-mm:mm,-lm:lm) real(8), dimension((2*lm+1)*(2*mm+1),size(aee,1)) :: e2a_aee !(out) 1 次元化された水平スペクトル・任意座標データ integer :: j,k,l,m ! !e2a_aee = transpose(reshape(aee,(/size(aee,1),(2*lm+1)*(2*mm+1)/))) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do k=1,size(aee,1) do l=-lm,lm do m=-mm,mm j = (m+mm+1)+ (2*mm+1)*(l+lm) e2a_aee(j,k) = aee(k,m,l) enddo enddo enddo end function e2a_aee function aee_e2a(e2a) ! ! 水平スペクトルを転置展開する. ! real(8), dimension(:,:),intent(in) :: e2a !(in) 1 次元化された水平スペクトル・任意座標データ ! dimmension((2*mm*1)*(2*lm*1),:) real(8), dimension(size(e2a,2),-mm:mm,-lm:lm) :: aee_e2a !(out) 任意座標・水平 2 次元スペクトルデータ integer :: j,k,l,m ! !aee_e2a = reshape(transpose(e2a),(/size(e2a,2),2*mm+1,2*lm+1/)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do k=1,size(e2a,2) do l=-lm,lm do m=-mm,mm j = (m+mm+1) + (2*mm+1)*(l+lm) aee_e2a(k,m,l) = e2a(j,k) enddo enddo enddo end function aee_e2a function e2_ee(ee) ! ! 水平スペクトルを 1 次元化転置する. ! real(8), dimension(-mm:,-lm:), intent(in) :: ee !(in) 任意座標・水平 2 次元スペクトルデータ dimension(-mm:mm,-lm:lm) real(8), dimension((2*lm+1)*(2*mm+1)) :: e2_ee !(out) 1 次元化された水平スペクトル・任意座標データ integer :: j,l,m !e2_ee = reshape(e2a_aee(reshape(ee,(/1,2*mm+1,2*lm+1/))), & ! (/(2*lm+1)*(2*mm+1)/)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do l=-lm,lm do m=-mm,mm j = (m+mm+1)+ (2*mm+1)*(l+lm) e2_ee(j) = ee(m,l) enddo enddo end function e2_ee function ee_e2(e2) ! ! 水平スペクトルを転置展開する. ! real(8), dimension(:),intent(in) :: e2 !(in) 1 次元化された水平スペクトル・任意座標データ ! dimmension((2*mm*1)*(2*lm*1),:) real(8), dimension(-mm:mm,-lm:lm) :: ee_e2 !(out) 任意座標・水平 2 次元スペクトルデータ integer :: j,l,m !ee_e2 = reshape(aee_e2a(reshape(e2,(/(2*mm+1)*(2*lm+1),1/))), & ! (/2*mm+1,2*lm+1/)) ! ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す ! do l=-lm,lm do m=-mm,mm j = (m+mm+1) + (2*mm+1)*(l+lm) ee_e2(m,l) = e2(j) enddo enddo end function ee_e2 !--------------- 微分計算 ----------------- function cee_Dx_cee(cee) ! ! 入力スペクトルデータに水平微分 ∂/∂x を作用する. ! ! スペクトルデータのX微分とは, 対応する格子点データにX微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_Dx_cee !(in) X微分されたスペクトルデータ integer :: n do n=0,nm cee_Dx_cee(n,:,:) = ee_Dx_ee(cee(n,:,:)) enddo end function cee_Dx_cee function see_Dx_see(see) ! ! 入力スペクトルデータに水平微分 ∂/∂x を作用する. ! ! スペクトルデータのX微分とは, 対応する格子点データにX微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 入力スペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_Dx_see !(in) X微分されたスペクトルデータ integer :: n do n=1,nm see_Dx_see(n,:,:) = ee_Dx_ee(see(n,:,:)) enddo end function see_Dx_see function cee_Dy_cee(cee) ! ! 入力スペクトルデータに水平微分 ∂/∂y を作用する. ! ! スペクトルデータのY微分とは, 対応する格子点データにY微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_Dy_cee !(in) X微分されたスペクトルデータ integer :: n do n=0,nm cee_Dy_cee(n,:,:) = ee_Dy_ee(cee(n,:,:)) enddo end function cee_Dy_cee function see_Dy_see(see) ! ! 入力スペクトルデータに水平微分 ∂/∂y を作用する. ! ! スペクトルデータのY微分とは, 対応する格子点データにY微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 入力スペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_Dy_see !(in) X微分されたスペクトルデータ integer :: n do n=1,nm see_Dy_see(n,:,:) = ee_Dy_ee(see(n,:,:)) enddo end function see_Dy_see function see_Dz_cee(cee) ! ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する. ! ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 入力スペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_Dz_cee !(in) 鉛直微分されたスペクトルデータ see_Dz_cee = aee_e2a(e2s_Dz_e2c(e2a_aee(cee))) end function see_Dz_cee function cee_Dz_see(see) ! ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する. ! ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_Dz_see !(in) 鉛直微分されたスペクトルデータ cee_Dz_see = aee_e2a(e2c_Dz_e2s(e2a_aee(see))) end function cee_Dz_see function e2s_Dz_e2c(e2c) ! ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する. ! ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension((2*lm+1)*(2*mm+1),0:nm), intent(IN) :: e2c !(in) 入力スペクトルデータ real(8), dimension((2*lm+1)*(2*mm+1),nm) :: e2s_Dz_e2c !(in) 鉛直微分されたスペクトルデータ e2s_Dz_e2c = as_Dz_ac(e2c) end function e2s_Dz_e2c function e2c_Dz_e2s(e2s) ! ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する. ! ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を ! 作用させたデータのスペクトル変換のことである. ! real(8), dimension((2*lm+1)*(2*mm+1),nm), intent(IN) :: e2s !(in) 入力スペクトルデータ real(8), dimension((2*lm+1)*(2*mm+1),0:nm) :: e2c_Dz_e2s !(in) 鉛直微分されたスペクトルデータ e2c_Dz_e2s = ac_Dz_as(e2s) end function e2c_Dz_e2s function cee_Lapla_cee(cee) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2 ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_Lapla_cee !(out) ラプラシアンを作用された 2 次元スペクトルデータ cee_Lapla_cee = cee_LaplaH_cee(cee) + cee_Dz_see(see_Dz_cee(cee)) end function cee_Lapla_cee function see_Lapla_see(see) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2 ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 2 次元球面調和函数 sin スペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_Lapla_see !(out) ラプラシアンを作用された 2 次元スペクトルデータ see_Lapla_see = see_LaplaH_see(see) + see_Dz_cee(cee_Dz_see(see)) end function see_Lapla_see function cee_LaplaInv_cee(cee) ! ! スペクトルデータに逆ラプラシアンを作用させる. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_LaplaInv_cee !(out) ラプラシアンを作用された 2 次元スペクトルデータ real(8), allocatable, dimension(:,:,:) :: cee_LaplaInvFactor logical :: first=.true. save first, cee_LaplaInvFactor if ( first ) then first = .false. allocate(cee_LaplaInvFactor(0:nm,-mm:mm,-lm:lm)) cee_LaplaInvFactor = 1.0D0 cee_LaplaInvFactor = 1.0D0/cee_Lapla_cee(cee_LaplaInvFactor) cee_LaplaInvFactor(0,0,0) = 0.0D0 end if cee_LaplaInv_cee = cee_LaplaInvFactor * cee end function cee_LaplaInv_cee function see_LaplaInv_see(see) ! ! スペクトルデータに逆ラプラシアンを作用させる. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 2 次元球面調和函数 cos スペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_LaplaInv_see !(out) ラプラシアンを作用された 2 次元スペクトルデータ real(8), allocatable, dimension(:,:,:) :: see_LaplaInvFactor logical :: first=.true. save first, see_LaplaInvFactor if ( first ) then first = .false. allocate(see_LaplaInvFactor(nm,-mm:mm,-lm:lm)) see_LaplaInvFactor = 1.0D0 see_LaplaInvFactor = 1.0D0/see_Lapla_see(see_LaplaInvFactor) end if see_LaplaInv_see = see_LaplaInvFactor * see end function see_LaplaInv_see function cee_LaplaH_cee(cee) ! ! 入力スペクトルデータに水平ラプラシアン ! ! ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2 ! ! を作用する. ! ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 入力スペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_LaplaH_cee !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ integer :: n do n=0,nm cee_LaplaH_cee(n,:,:) = ee_Lapla_ee(cee(n,:,:)) enddo end function cee_LaplaH_cee function see_LaplaH_see(see) ! ! 入力スペクトルデータに水平ラプラシアン ! ! ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2 ! ! を作用する. ! ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 入力スペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_LaplaH_see !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ integer :: n do n=1,nm see_LaplaH_see(n,:,:) = ee_Lapla_ee(see(n,:,:)) enddo end function see_LaplaH_see function cee_LaplaHInv_cee(cee) ! ! 入力スペクトルデータに逆水平ラプラシアン ! ! ▽^-2_H = (∂^2/∂X^2 + ∂^2/∂Y^2)^-1 ! ! を作用する. ! ! スペクトルデータの逆水平ラプラシアンとは, 対応する格子点データに ! 逆水平ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_LaplaHInv_cee !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ integer :: n do n=0,nm cee_LaplaHInv_cee(n,:,:) = ee_LaplaInv_ee(cee(n,:,:)) enddo end function cee_LaplaHInv_cee function see_LaplaHInv_see(see) ! ! 入力スペクトルデータに逆水平ラプラシアン ! ! ▽^-2_H = (∂^2/∂X^2 + ∂^2/∂Y^2)^-1 ! ! を作用する. ! ! スペクトルデータの逆水平ラプラシアンとは, 対応する格子点データに ! 逆水平ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(nm,-mm:mm,-lm:lm) :: see_LaplaHInv_see !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ integer :: n do n=1,nm see_LaplaHInv_see(n,:,:) = ee_LaplaInv_ee(see(n,:,:)) enddo end function see_LaplaHInv_see function see_Div_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ) ! ! べクトル成分である 3 つの格子データに発散を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX !(in) ベクトル場のX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY !(in) ベクトル場のY成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ !(in) ベクトル場のZ成分 real(8), dimension(nm,-mm:mm,-lm:lm) :: see_Div_zyx_zyx_zyx !(out) ベクトル場の発散 see_Div_zyx_zyx_zyx = see_Dx_see(see_zyx(zyx_VX)) & + see_Dy_see(see_zyx(zyx_VY)) & + see_Dz_cee(cee_zyx(zyx_VZ)) end function see_Div_zyx_zyx_zyx function cee_Div_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ) ! ! べクトル成分である 3 つの格子データに発散を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX !(in) ベクトル場のX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY !(in) ベクトル場のY成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ !(in) ベクトル場のZ成分 real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_Div_zyx_zyx_zyx !(out) ベクトル場の発散 cee_Div_zyx_zyx_zyx = cee_Dx_cee(cee_zyx(zyx_VX)) & + cee_Dy_cee(cee_zyx(zyx_VY)) & + cee_Dz_see(see_zyx(zyx_VZ)) end function cee_Div_zyx_zyx_zyx !--------------- 積分計算 ----------------- !----(入力データ zyx)--- function zy_IntX_zyx(zyx) ! X積分 ! ! 3 次元格子点データのX方向積分. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:km,0:jm-1) :: zy_IntX_zyx !(out) X方向積分された 2 次元ZY格子点データ integer :: i zy_IntX_zyx = 0.0d0 do i=0,im-1 zy_IntX_zyx(:,:) = zy_IntX_zyx(:,:) & + zyx(:,:,i) * x_X_Weight(i) enddo end function zy_IntX_zyx function zx_IntY_zyx(zyx) ! ! 3 次元格子点データのY方向域積分. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:km,0:im-1) :: zx_IntY_zyx !(out) Y積分された 2 次元ZY格子点データ. integer :: j zx_IntY_zyx = 0.0d0 do j=0,jm-1 zx_IntY_zyx(:,:) = zx_IntY_zyx(:,:) & + zyx(:,j,:) * y_Y_Weight(j) enddo end function zx_IntY_zyx function yx_IntZ_zyx(zyx) ! Z積分 ! ! 3 次元格子点データのZ方向積分. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:jm-1,0:im-1) :: yx_IntZ_zyx !(out) Z積分された 2 次元YX(水平)格子点データ integer :: k yx_IntZ_zyx = 0.0d0 do k=0,km yx_IntZ_zyx(:,:) = yx_IntZ_zyx(:,:) & + zyx(k,:,:) * z_Z_Weight(k) enddo end function yx_IntZ_zyx function x_IntZY_zyx(zyx) ! ! 3 次元格子点データのZY積分 ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元ZYX格子点データ real(8), dimension(0:im-1) :: x_IntZY_zyx !(out) ZY(子午面)積分された 1 次元X格子点データ integer :: j, k x_IntZY_zyx = 0.0D0 do j=0,jm-1 do k=0,km x_IntZY_zyx = x_IntZY_zyx & + zyx(k,j,:) * y_Y_Weight(j) * z_Z_Weight(k) enddo enddo end function x_IntZY_zyx function y_IntZX_zyx(zyx) ! ! 3 次元格子点データのZX積分. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元ZYX格子点データ real(8), dimension(0:jm-1) :: y_IntZX_zyx !(out) ZX積分された 1 次元Y格子点データ integer :: i, k y_IntZX_zyx = 0 do i=0,im-1 do k=0,km y_IntZX_zyx = y_IntZX_zyx & + zyx(k,:,i) * x_X_Weight(i) * z_Z_Weight(k) enddo enddo end function y_IntZX_zyx function z_IntYX_zyx(zyx) ! YX(水平)積分 ! ! 3 次元格子点データのYX(水平)積分 ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元ZYX格子点データ real(8), dimension(0:km) :: z_IntYX_zyx !(out) YX(水平)積分された 1 次元Z格子点データ integer :: i, j z_IntYX_zyx = 0 do j=0,jm-1 do i=0,im-1 z_IntYX_zyx = z_IntYX_zyx & + zyx(:,j,i) * x_X_Weight(i) * y_Y_Weight(j) enddo enddo end function z_IntYX_zyx function IntZYX_zyx(zyx) ! ZYX(全領域)積分 ! ! 3 次元格子点データのZYX(全領域)積分 ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8) :: IntZYX_zyx !(out) 全領域積分値 integer :: i, j, k IntZYX_zyx = 0 do i=0,im-1 do j=0,jm-1 do k=0,km IntZYX_zyx = IntZYX_zyx & + zyx(k,j,i) * x_X_Weight(i) & * y_Y_Weight(j) * z_Z_Weight(k) enddo enddo enddo end function IntZYX_zyx !----(入力データ zy)--- function z_IntY_zy(zy) ! Y積分 ! ! 2 次元(ZY)格子点データのY方向積分. ! real(8), dimension(0:km,0:jm-1), intent(in) :: zy !(in) 2 次元ZY(子午面)格子点データ real(8), dimension(0:km) :: z_IntY_zy !(out) Y積分された 1 次元Z格子点データ integer :: j z_IntY_zy = 0.0d0 do j=0,jm-1 z_IntY_zy(:) = z_IntY_zy(:) + zy(:,j) * y_Y_Weight(j) enddo end function z_IntY_zy function y_IntZ_zy(zy) ! Z積分 ! ! 2 次元(ZY)格子点データのZ方向域積分. ! real(8), dimension(0:km,0:jm-1), intent(in) :: zy !(in) 2 次元ZY格子点データ real(8), dimension(0:jm-1) :: y_IntZ_zy !(out) Z積分された 1 次元Y格子点データ integer :: k y_IntZ_zy = 0.0d0 do k=0,km y_IntZ_zy(:) = y_IntZ_zy(:) & + zy(k,:) * z_Z_Weight(k) enddo end function y_IntZ_zy function IntZY_zy(zy) ! ! 2 次元(ZY)格子点データのZY積分 ! real(8), dimension(0:km,0:jm-1), intent(in) :: zy !(in) 2 次元ZY(子午面)格子点データ real(8) :: IntZY_zy !(out) 積分値 integer :: j, k IntZY_zy = 0 do j=0,jm-1 do k=0,km IntZY_zy = IntZY_zy & + zy(k,j) * y_Y_Weight(j) * z_Z_Weight(k) enddo enddo end function IntZY_zy !----(入力データ zx)--- function z_IntX_zx(zx) ! ! 2 次元(ZX)格子点データのX方向積分. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZY格子点データ real(8), dimension(0:km) :: z_IntX_zx !(out) X積分された 1 次元Z格子点データ integer :: i z_IntX_zx = 0.0d0 do i=0,im-1 z_IntX_zx(:) = z_IntX_zx(:) + zx(:,i) * x_X_Weight(i) enddo end function z_IntX_zx function x_IntZ_zx(zx) ! ! 2 次元(ZX)格子点データのZ方向積分. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZY格子点データ real(8), dimension(0:im-1) :: x_IntZ_zx !(out) Z積分された 1 次元X格子点データ integer :: k x_IntZ_zx = 0.0d0 do k=0,km x_IntZ_zx(:) = x_IntZ_zx(:) & + zx(k,:) * z_Z_Weight(k) enddo end function x_IntZ_zx function IntZX_zx(zx) ! ZX積分 ! ! 2 次元(ZX)格子点データのZX積分 ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZY格子点データ real(8) :: IntZX_zx !(out) 積分値 integer :: i, k IntZX_zx = 0 do i=0,im-1 do k=0,km IntZX_zx = IntZX_zx & + zx(k,i) * x_X_Weight(i) * z_Z_Weight(k) enddo enddo end function IntZX_zx !----(入力データ z)--- function IntZ_z(z) ! Z積分 ! ! 1 次元(Z)格子点データのZ方向積分. ! real(8), dimension(0:km), intent(in) :: z !(in) 1 次元Z格子点データ real(8) :: IntZ_z !(out) 積分値 integer :: k IntZ_z = 0.0d0 do k=0,km IntZ_z = IntZ_z + z(k) * z_Z_Weight(k) enddo end function IntZ_z !--------------- 平均計算 ----------------- !----(入力データ zyx)--- function zy_AvrX_zyx(zyx) ! X積分 ! ! 3 次元格子点データのX方向平均. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元ZYX格子点データ real(8), dimension(0:km,0:jm-1) :: zy_AvrX_zyx !(out) X方向平均された 2 次元子午面格子点データ zy_AvrX_zyx = zy_IntX_zyx(zyx)/sum(x_X_Weight) end function zy_AvrX_zyx function zx_AvrY_zyx(zyx) ! Y平均 ! ! 3 次元格子点データのY方向域平均. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:km,0:im-1) :: zx_AvrY_zyx !(out) Y平均された 2 次元ZX格子点データ zx_AvrY_zyx = zx_IntY_zyx(zyx)/sum(y_Y_Weight) end function zx_AvrY_zyx function yx_AvrZ_zyx(zyx) ! ! 3 次元格子点データのZ方向平均. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:jm-1,0:im-1) :: yx_AvrZ_zyx !(out) Z平均された 2 次元YX(水平)格子点データ yx_AvrZ_zyx = yx_IntZ_zyx(zyx)/sum(z_Z_Weight) end function yx_AvrZ_zyx function x_AvrZY_zyx(zyx) ! ZY積分 ! ! 3 次元格子点データのZY平均 ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元ZYX格子点データ real(8), dimension(0:im-1) :: x_AvrZY_zyx !(out) ZY平均された 1 次元X格子点データ x_AvrZY_zyx = x_IntZY_zyx(zyx) & /( sum(y_Y_Weight)*sum(z_Z_Weight) ) end function x_AvrZY_zyx function y_AvrZX_zyx(zyx) ! ZX積分 ! ! 3 次元格子点データのZX平均. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:jm-1) :: y_AvrZX_zyx !(out) ZX平均された 1 次元Y格子点データ y_AvrZX_zyx = y_IntZX_zyx(zyx) & /(sum(x_X_Weight)*sum(z_Z_Weight)) end function y_AvrZX_zyx function z_AvrYX_zyx(zyx) ! YX(水平)積分 ! ! 3 次元格子点データのYX(水平)積分 ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元格子点データ real(8), dimension(0:km) :: z_AvrYX_zyx !(out) YX(水平)平均された 1 次元Z格子点データ z_AvrYX_zyx = z_IntYX_zyx(zyx) & /(sum(x_X_Weight)*sum(y_Y_Weight)) end function z_AvrYX_zyx function AvrZYX_zyx(zyx) ! ZYX(全領域)積分 ! ! 3 次元格子点データのZYX(全領域)積分 ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx !(in) 3 次元ZYX格子点データ real(8) :: AvrZYX_zyx !(out) 全領域平均値 AvrZYX_zyx = IntZYX_zyx(zyx) & /(sum(x_X_Weight)*sum(y_Y_Weight) * sum(z_Z_Weight)) end function AvrZYX_zyx !----(入力データ zy)--- function z_AvrY_zy(zy) ! ! 2 次元(ZY)格子点データのY方向平均. ! real(8), dimension(0:km,0:jm-1), intent(in) :: zy !(in) 2 次元ZY格子点データ real(8), dimension(0:km) :: z_AvrY_zy !(out) Y平均された 1 次元Z格子点データ z_AvrY_zy = z_IntY_zy(zy)/sum(y_Y_Weight) end function z_AvrY_zy function y_AvrZ_zy(zy) ! ! 2 次元(ZY)格子点データのZ方向平均. ! real(8), dimension(0:km,0:jm-1), intent(in) :: zy !(in) 2 次元ZY格子点データ real(8), dimension(0:jm-1) :: y_AvrZ_zy !(out) Z平均された 1 次元Y格子点データ y_AvrZ_zy = y_IntZ_zy(zy)/sum(z_Z_Weight) end function y_AvrZ_zy function AvrZY_zy(zy) ! ZY平均 ! ! 2 次元(ZY)格子点データのZY平均 ! real(8), dimension(0:km,0:jm-1), intent(in) :: zy !(in) 2 次元ZY(子午面)格子点データ real(8) :: AvrZY_zy !(out) 平均値 AvrZY_zy = IntZY_zy(zy)/(sum(y_Y_Weight)*sum(z_Z_Weight)) end function AvrZY_zy !----(入力データ zx)--- function z_AvrX_zx(zx) ! X積分 ! ! 2 次元(ZX)格子点データのX方向平均. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZY格子点データ real(8), dimension(0:km) :: z_AvrX_zx !(out) X平均された 1 次元Z格子点データ z_AvrX_zx = z_IntX_zx(zx)/sum(x_X_Weight) end function z_AvrX_zx function x_AvrZ_zx(zx) ! Z平均 ! ! 2 次元(ZX)格子点データのZ方向平均. ! real(8), dimension(0:km,0:im-1), intent(in) :: zx !(in) 2 次元ZY格子点データ real(8), dimension(0:im-1) :: x_AvrZ_zx !(out) Z平均された 1 次元X格子点データ x_AvrZ_zx = x_IntZ_zx(zx)/sum(z_Z_Weight) end function x_AvrZ_zx function AvrZX_zx(zx) ! ZX積分 ! ! 2 次元(ZX)格子点データのZX平均 ! real(8), dimension(0:km,0:im-1), intent(in) :: zx ! (in) 2 次元格子点データ real(8) :: AvrZX_zx ! 平均値 AvrZX_zx = IntZX_zx(zx)/(sum(x_X_Weight)*sum(z_Z_Weight)) end function AvrZX_zx !----(入力データ z)--- function AvrZ_z(z) ! ! 1 次元(Z)格子点データのZ方向平均. ! real(8), dimension(0:km), intent(in) :: z !(in) 1 次元Z格子点データ real(8) :: AvrZ_z !(out) 平均値 AvrZ_z = IntZ_z(z)/sum(z_Z_Weight) end function AvrZ_z !--------------- ポロイダル/トロイダルモデル用微分 ----------------- !!$ function wt_KxRGrad_wt(wt) !!$ ! !!$ ! 入力スペクトルデータにX微分 k×r・▽ = ∂/∂λを作用する. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_KxRGrad_wt !!$ !(out) X微分を作用された 2 次元スペクトルデータ !!$ !!$ wt_KxRGrad_wt = wa_Dlon_wa(wt) !!$ !!$ end function wt_KxRGrad_wt !!$ !!$ function zyx_KGrad_wt(wt) ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r !!$ ! !!$ ! 入力スペクトルデータに対応する格子データに軸方向微分 !!$ ! !!$ ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r !!$ ! !!$ ! を作用させた格子データが返される. !!$ ! ここでベクトル k は球の中心から北極向きの単位ベクトルである. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(0:km,0:jm-1,0:im-1) :: yzx_KGrad_wt !!$ !(out) 軸方向微分を作用された 2 次元スペクトルデータ !!$ !!$ xzy_KGrad_wt = cos(xyz_Y)*xyz_GradY_wt(wt) & !!$ + sin(xyz_Y)*xyz_wt(wt_Drad_wt(wt)) !!$ !!$ end function zyx_KGrad_wt !!$ !!$ function wt_L2_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する. !!$ ! !!$ ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる. !!$ ! 入力スペクトルデ ータに対応する格子点データに演算子 !!$ ! !!$ ! L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ) !!$ ! !!$ ! を作用させたデータのスペクトル変換が返される. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_L2_wt !!$ !(out) L^2 演算子を作用された 2 次元スペクトルデータ !!$ !!$ wt_L2_wt = -wa_Lapla_wa(wt) !!$ !!$ end function wt_L2_wt !!$ !!$ function wt_L2Inv_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を !!$ ! 作用する. !!$ ! !!$ ! スペクトルデータに L^2 演算子を作用させる関数 wt_L2_wt の逆計算を !!$ ! 行う関数である. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_L2Inv_wt !!$ !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ !!$ !!$ wt_L2Inv_wt = -wa_LaplaInv_wa(wt) !!$ !!$ end function wt_L2Inv_wt !!$ !!$ function wt_QOperator_wt(wt) !!$ ! !!$ ! 入力スペクトルデータに対応する格子点データに演算子 !!$ ! !!$ ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) !!$ ! !!$ ! を作用させたデータのスペクトル変換が返される. !!$ ! !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt !!$ !(in) 2 次元球面調和函数チェビシェフスペクトルデータ !!$ !!$ real(8), dimension(-lm:lm,-mm,mm,0:nm) :: wt_QOperator_wt !!$ !(out) Q 演算子を作用された 2 次元スペクトルデータ !!$ !!$ wt_QOperator_wt = & !!$ wt_xyz(xyz_KGrad_wt(wt) - xyz_KGrad_wt(wt_L2_wt(wt))/2) & !!$ - wt_L2_wt(wt_xyz(xyz_KGrad_wt(wt)))/2 !!$ !!$ end function wt_QOperator_wt function cee_ZRot_zyx_zyx(zyx_VX,zyx_VY) ! z・(▽×v) ! ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する. ! ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す. ! ! z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y ! ! のスペクトル データが返される. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX !(in) ベクトルのX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY !(in) ベクトルのY成分 real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_ZRot_zyx_zyx !(out) ベクトルの渦度とZベクトルの内積 cee_ZRot_zyx_zyx = cee_DX_cee(cee_zyx(zyx_VY)) & - cee_DY_cee(cee_zyx(zyx_VX)) end function cee_ZRot_zyx_zyx function see_ZRot_zyx_zyx(zyx_VX,zyx_VY) ! z・(▽×v) ! ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する. ! ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す. ! ! z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y ! ! のスペクトル データが返される. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX !(in) ベクトルのX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY !(in) ベクトルのY成分 real(8), dimension(nm,-mm:mm,-lm:lm) :: see_ZRot_zyx_zyx !(out) ベクトルの渦度とZベクトルの内積 see_ZRot_zyx_zyx = see_DX_see(see_zyx(zyx_VY)) & - see_DY_see(see_zyx(zyx_VX)) end function see_ZRot_zyx_zyx function cee_ZRotRot_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ) ! ! ベクトル v に対して z・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分, ! Y成分, Z成分を表す. ! ! z・(▽×▽×v) = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) ) ! - ▽_H^2 v[z] ! ! のスペクトルデータが返される. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX !(in) ベクトルのX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY !(in) ベクトルのY成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ !(in) ベクトルのZ成分 real(8), dimension(0:nm,-mm:mm,-lm:lm) :: cee_ZRotRot_zyx_zyx_zyx !(out) ベクトル v の z・(▽×▽×v) cee_ZRotRot_zyx_zyx_zyx = & cee_DZ_see( see_DX_see(see_zyx(zyx_VX)) & + see_DY_see(see_zyx(zyx_VY)) ) & - cee_LaplaH_cee(cee_zyx(zyx_VZ)) end function cee_ZRotRot_zyx_zyx_zyx function see_ZRotRot_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_VZ) ! ! ベクトル v に対して z・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分, ! Y成分, Z成分を表す. ! ! z・(▽×▽×v) = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) ) ! - ▽_H^2 v[z] ! ! のスペクトルデータが返される. ! real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VX !(in) ベクトルのX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VY !(in) ベクトルのY成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx_VZ !(in) ベクトルのZ成分 real(8), dimension(nm,-mm:mm,-lm:lm) :: see_ZRotRot_zyx_zyx_zyx !(out) ベクトル v の z・(▽×▽×v) see_ZRotRot_zyx_zyx_zyx = & see_DZ_cee( cee_DX_cee(cee_zyx(zyx_VX)) & + cee_DY_cee(cee_zyx(zyx_VY)) ) & - see_LaplaH_see(see_zyx(zyx_VZ)) end function see_ZRotRot_zyx_zyx_zyx subroutine scee_Potential2Vector_cee_see(& zyx_VX,zyx_VY,zyx_VZ,cee_Torvel,see_Polvel) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψz) + ▽x▽x(Φz) ! ! の各成分を計算する ! real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_VX !(out) ベクトル場のX成分 real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_VY !(out) ベクトル場のY成分 real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_VZ !(out) ベクトル場のZ成分 real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Torvel !(in) トロイダルポテンシャル real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_Polvel !(in) ポロイダルポテンシャル zyx_VX = zyx_cee( cee_DY_cee(cee_Torvel) & + cee_DX_cee(cee_DZ_see(see_Polvel)) ) zyx_VY = zyx_cee( - cee_DX_cee(cee_Torvel) & + cee_DY_cee(cee_DZ_see(see_Polvel)) ) zyx_VZ = -zyx_see(see_LaplaH_see(see_Polvel)) end subroutine scee_Potential2Vector_cee_see subroutine scee_Potential2Vector_see_cee(& zyx_VX,zyx_VY,zyx_VZ,see_Torvel,cee_Polvel) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψz) + ▽x▽x(Φz) ! ! の各成分を計算する ! real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_VX !(out) ベクトル場のX成分 real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_VY !(out) ベクトル場のY成分 real(8), dimension(0:km,0:jm-1,0:im-1) :: zyx_VZ !(out) ベクトル場のZ成分 real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_Torvel !(in) トロイダルポテンシャル real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Polvel !(in) ポロイダルポテンシャル zyx_VX = zyx_see( see_DY_see(see_Torvel) & + see_DX_see(see_DZ_cee(cee_Polvel)) ) zyx_VY = zyx_see( - see_DX_see(see_Torvel) & + see_DY_see(see_DZ_cee(cee_Polvel)) ) zyx_VZ = -zyx_cee(cee_LaplaH_cee(cee_Polvel)) end subroutine scee_Potential2Vector_see_cee subroutine scee_Potential2Rotation_cee_see(& zyx_RotVX,zyx_RotVY,zyx_RotVZ,cee_Torvel,see_Polvel) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψz) + ▽x▽x(Φz) ! ! に対して, その回転 ! ! ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z) ! ! を計算する. ! ベクトル場の回転 real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVX !(out) 回転のX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVY !(out) 回転のY成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVZ !(out) 回転のZ成分 ! 入力ベクトル場を表すポテンシャル real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Torvel !(in) トロイダルポテンシャル real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_Polvel !(in) ポロイダルポテンシャル call scee_Potential2Vector_see_cee(& zyx_RotVX,zyx_RotVY,zyx_RotVZ,-see_Lapla_see(see_Polvel),cee_Torvel) end subroutine scee_Potential2Rotation_cee_see subroutine scee_Potential2Rotation_see_cee(& zyx_RotVX,zyx_RotVY,zyx_RotVZ,see_Torvel,cee_Polvel) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψz) + ▽x▽x(Φz) ! ! に対して, その回転 ! ! ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z) ! ! を計算する. ! ベクトル場の回転 real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVX !(out) 回転のX成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVY !(out) 回転のY成分 real(8), dimension(0:km,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVZ !(out) 回転のZ成分 ! 入力ベクトル場を表すポテンシャル real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_Torvel !(in) トロイダルポテンシャル real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_Polvel !(in) ポロイダルポテンシャル call scee_Potential2Vector_cee_see(& zyx_RotVX,zyx_RotVY,zyx_RotVZ,-cee_Lapla_cee(cee_Polvel),see_Torvel) end subroutine scee_Potential2Rotation_see_cee !--------------- 補間計算 ----------------- function Interpolate_cee(cee_data,x,y,z) ! ! (x,y,z) における関数値を ! そのスペクトル係数 tee_data から補間計算する ! real(8), intent(IN) :: cee_data(0:nm,-mm:mm,-lm:lm) ! スペクトルデータ real(8), intent(IN) :: x ! 補間する位置(X) real(8), intent(IN) :: y ! 補間する位置(Y) real(8), intent(IN) :: z ! 補間する位置(Z) real(8) :: Interpolate_cee ! 補間した値 Interpolate_cee = & Interpolate_ee(ee_e2(a_Interpolate_ac(e2a_aee(cee_data),z)),x,y) end function Interpolate_cee function Interpolate_see(see_data,x,y,z) ! ! (x,y,z) における関数値を ! そのスペクトル係数 tee_data から補間計算する ! real(8), intent(IN) :: see_data(nm,-mm:mm,-lm:lm) ! スペクトルデータ real(8), intent(IN) :: x ! 補間する位置(X) real(8), intent(IN) :: y ! 補間する位置(Y) real(8), intent(IN) :: z ! 補間する位置(Z) real(8) :: Interpolate_see ! 補間した値 Interpolate_see = & Interpolate_ee(ee_e2(a_Interpolate_as(e2a_aee(see_data),z)),x,y) end function Interpolate_see !--------------- ポロイダル/トロイダルモデル用スペクトル解析 ---------------- function zee_ToroidalEnergySpectrum_cee(cee_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 水平 X 波数 l, Y 波数 m の各成分を計算する ! ! * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー ! スペクトルは (1/2)(l**2+m**2)|ψ(l,m,z)|^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での ! 水平平均エネルギーに等しい. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_TORPOT !(in) トロイダルポテンシャル real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_ToroidalEnergySpectrum_cee !(out) エネルギースペクトルトロイダル成分 real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Work integer :: l, m zee_Work = zee_cee(cee_Torpot) do l=-lm,lm do m=-mm,mm zee_ToroidalEnergySpectrum_cee(:,m,l) & = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) & * ( zee_Work(:,m,l)**2 + zee_Work(:,-m,-l)**2 ) enddo enddo end function zee_ToroidalEnergySpectrum_cee function zee_ToroidalEnergySpectrum_see(see_TORPOT) ! ! トロイダルポテンシャルから, トロイダルエネルギーの ! 水平 X 波数 l, Y 波数 m の各成分を計算する ! ! * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分 ! ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー ! スペクトルは (1/2)(l**2+m**2)|ψ(l,m,z)|^2 と計算される. ! ! * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での ! 水平平均エネルギーに等しい. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_TORPOT !(in) トロイダルポテンシャル real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_ToroidalEnergySpectrum_see !(out) エネルギースペクトルトロイダル成分 real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Work integer :: l, m zee_Work = zee_see(see_Torpot) do l=-lm,lm do m=-mm,mm zee_ToroidalEnergySpectrum_see(:,m,l) & = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) & * ( zee_Work(:,m,l)**2 + zee_Work(:,-m,-l)**2 ) enddo enddo end function zee_ToroidalEnergySpectrum_see !!$ function nz_ToroidalEnergySpectrum_wt(wt_TORPOT) !!$ ! !!$ ! トロイダルポテンシャルから, トロイダルエネルギーの !!$ ! 球面調和函数全波数の各成分を計算する. !!$ ! !!$ ! * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分 !!$ ! ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは !!$ ! Σ[m=-n]^n(1/2)n(n+1)4πr^2ψ(n,m,r)^2 と計算される. !!$ ! !!$ ! * 全てのエネルギースペクトル成分の和をZ積分したもの(r^2の重み無し) !!$ ! が球殻内での全エネルギーに等しい. !!$ ! !!$ real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_TORPOT !!$ !(in) トロイダルポテンシャル !!$ !!$ real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wt !!$ !(out) エネルギースペクトルトロイダル成分 !!$ !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA ! 作業領域 !!$ integer :: n, m !!$ !!$ wz_DATA = wz_wt(wt_TORPOT) !!$ do n=0,nm !!$ nz_ToroidalEnergySpectrum_wt(n,:) & !!$ = 0.5 * n*(n+1)* (4*pi) * z_Z**2 & !!$ * sum(wz_Data(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) !!$ enddo !!$ !!$ end function nz_ToroidalEnergySpectrum_wt function zee_PoloidalEnergySpectrum_cee(cee_POLPOT) ! ! ポロイダルポテンシャルから, ポロイダルエネルギーの ! 水平 X 波数 l, Y 波数 m の各成分を計算する. ! ! * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分 ! φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー ! スペクトルは ! ! (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2} ! ! と計算される. ! ! * 全てのエネルギースペクトル成分の和をZ積分したもの ! が水平平均エネルギーに等しい. ! real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: cee_POLPOT !(in) ポロイダルポテンシャル real(8), dimension(0:km,-nm:nm,-lm:lm) :: zee_PoloidalEnergySpectrum_cee !(out) エネルギースペクトルポロイダル成分 real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Data ! 作業領域 real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_DData ! 作業領域 integer :: l, m zee_Data = zee_cee(cee_POLPOT) zee_DData = zee_see(see_DZ_cee(cee_POLPOT)) do l=-lm,lm do m=-mm,mm zee_PoloidalEnergySpectrum_cee(:,m,l) = & + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2) & *( zee_DData(:,m,l)**2 + zee_DData(:,-m,-l)**2 & + ((2*pi*l/xl)**2+(2*pi*m/yl)**2) & *( zee_Data(:,m,l)**2 + zee_Data(:,-m,-l)**2)) enddo enddo end function zee_PoloidalEnergySpectrum_cee function zee_PoloidalEnergySpectrum_see(see_POLPOT) ! ! ポロイダルポテンシャルから, ポロイダルエネルギーの ! 水平 X 波数 l, Y 波数 m の各成分を計算する. ! ! * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分 ! φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー ! スペクトルは ! ! (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2} ! ! と計算される. ! ! * 全てのエネルギースペクトル成分の和をZ積分したもの ! が水平平均エネルギーに等しい. ! real(8), dimension(nm,-mm:mm,-lm:lm), intent(in) :: see_POLPOT !(in) ポロイダルポテンシャル real(8), dimension(0:km,-nm:nm,-lm:lm) :: zee_PoloidalEnergySpectrum_see !(out) エネルギースペクトルポロイダル成分 real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Data ! 作業領域 real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_DData ! 作業領域 integer :: l, m zee_Data = zee_see(see_POLPOT) zee_DData = zee_cee(cee_DZ_see(see_POLPOT)) do l=-lm,lm do m=-mm,mm zee_PoloidalEnergySpectrum_see(:,m,l) = & + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2) & *( zee_DData(:,m,l)**2 + zee_DData(:,-m,-l)**2 & + ((2*pi*l/xl)**2+(2*pi*m/yl)**2) & *( zee_Data(:,m,l)**2 + zee_Data(:,-m,-l)**2)) enddo enddo end function zee_PoloidalEnergySpectrum_see !!$ function nz_PoloidalEnergySpectrum_wt(wt_POLPOT) !!$ ! !!$ ! ポロイダルポテンシャルから, ポロイダルエネルギーの !!$ ! 球面調和函数全波数の各成分を計算する !!$ ! !!$ ! * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分 !!$ ! φ(n,m,r)から全波数 n 成分のポロイダルエネルギースペクトルは !!$ ! !!$ ! Σ[m=-n]^n ((1/2)n(n+1)4πr^2{[d(rφ(n,m,r))/dr]^2 !!$ ! + n(n+1)φ(n,m,r)^2} !!$ ! !!$ ! と計算される. !!$ ! !!$ ! * 全ての全波数に対してのエネルギースペクトル成分の和をZ積分したもの !!$ ! (r^2の重み無し)が球殻内での全エネルギーに等しい. !!$ ! !!$ real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_POLPOT !!$ !(in) ポロイダルポテンシャル !!$ !!$ real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_wt !!$ !(out) エネルギースペクトルポロイダル成分 !!$ !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1 ! 作業領域 !!$ real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2 ! 作業領域 !!$ integer :: n, m !!$ !!$ wz_Data1 = wz_wt(wt_POLPOT) !!$ wz_Data2 = wz_Z*wz_wt(wt_DZ_wt(wt_POLPOT)) & ! d(rφ)/dr !!$ + wz_wt(wt_POLPOT) ! = rdφ/dr+φ !!$ !!$ do n=0,nm !!$ nz_PoloidalEnergySpectrum_wt(n,:) = & !!$ + 0.5* n*(n+1)* (4*pi) & !!$ *( sum(wz_Data2(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) & !!$ + n*(n+1)*sum(wz_Data1(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) ) !!$ enddo !!$ !!$ end function nz_PoloidalEnergySpectrum_wt end module scee_module_fftj