Class | w_base_module_sjpack_cuda |
In: |
libsrc/w_module_sjpack_cuda/w_base_module_sjpack_cuda.f90
|
Authors: | Shin-ichi Takehiro, Youhei SASAKI |
Version: | $Id: w_base_module_sjpack_cuda.f90 590 2013-08-19 08:48:21Z uwabami $ |
Copyright&License: | See COPYRIGHT |
spml/w_base_module_cuda モジュールは球面上での 2 次元流体運動を 球面調和函数を用いたスペクトル法によって数値計算するための モジュール w_module_sjpack_cuda の下部モジュールであり, スペクトル法の 基本的な Fortran90 関数を提供する.
内部で ISPACK の sjpack-cuda Fortran77 サブルーチンを呼んでいる. スペクトルデータおよび格子点データの格納方法や変換の 詳しい計算法については ISPACK/sjpack-cuda のマニュアルを 参照されたい.
Function : | |||
l_nm_array00 : | integer
| ||
n : | integer, intent(in)
| ||
m : | integer, intent(in)
|
全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
引数 n,m がともに整数値の場合, 整数値を返す.
function l_nm_array00(n,m) ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 引数 n,m がともに整数値の場合, 整数値を返す. ! integer :: l_nm_array00 !(out) スペクトルデータの格納位置 integer, intent(in) :: n !(in) 全波数 integer, intent(in) :: m !(in) 帯状波数 if ( .not. w_base_initialize ) then call MessageNotify('E','l_nm_array00', 'w_base_module not initialize yet. Use sjnm2l routine in ISPACK directly.') else call sjnm2l(nn,n,m,l_nm_array00) endif end function l_nm_array00
Function : | |||
l_nm_array01(size(marray)) : | integer
| ||
n : | integer, intent(in)
| ||
marray(:) : | integer, intent(in)
|
スペクトルデータの格納位置
全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, marray と同じ大きさの 1 次元整数配列を返す.
function l_nm_array01(n,marray) ! スペクトルデータの格納位置 ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, ! marray と同じ大きさの 1 次元整数配列を返す. ! integer, intent(in) :: n !(in) 全波数 integer, intent(in) :: marray(:) !(in) 帯状波数 integer :: l_nm_array01(size(marray)) !(out) スペクトルデータ位置 integer :: i do i=1, size(marray) l_nm_array01(i) = l_nm_array00(n,marray(i)) enddo end function l_nm_array01
Function : | |||
l_nm_array10(size(narray)) : | integer
| ||
narray(:) : | integer, intent(in)
| ||
m : | integer, intent(in)
|
全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
第 1 引数 narray が整数 1 次元配列, 第 2 引数 m が整数の場合, narray と同じ大きさの 1 次元整数配列を返す.
function l_nm_array10(narray,m) ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 第 1 引数 narray が整数 1 次元配列, 第 2 引数 m が整数の場合, ! narray と同じ大きさの 1 次元整数配列を返す. ! integer, intent(in) :: narray(:) !(in) 全波数 integer, intent(in) :: m !(in) 帯状波数 integer :: l_nm_array10(size(narray)) !(out) スペクトルデータ位置 integer :: i do i=1, size(narray) l_nm_array10(i) = l_nm_array00(narray(i),m) enddo end function l_nm_array10
Function : | |||
l_nm_array11(size(narray)) : | integer
| ||
narray(:) : | integer, intent(in)
| ||
marray(:) : | integer, intent(in)
|
全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, narray, marray と同じ大きさの 1 次元整数配列を返す. narray, marray は同じ大きさでなければならない.
function l_nm_array11(narray,marray) ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, ! narray, marray と同じ大きさの 1 次元整数配列を返す. ! narray, marray は同じ大きさでなければならない. ! integer, intent(in) :: narray(:) !(in) 全波数 integer, intent(in) :: marray(:) !(in) 帯状波数 integer :: l_nm_array11(size(narray)) !(out) スペクトルデータ位置 integer :: i if ( size(narray) .ne. size(marray) ) then call MessageNotify('E','l_nm_array11', 'dimensions of input arrays n and m are different.') endif do i=1, size(narray) l_nm_array11(i) = l_nm_array00(narray(i),marray(i)) enddo end function l_nm_array11
Function : | |||
nm_l_int(2) : | integer
| ||
l : | integer, intent(in)
|
スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
引数 l が整数値の場合, 対応する全波数と帯状波数を 長さ 2 の 1 次元整数値を返す. nm_l(1) が全波数, nm_l(2) が帯状波数である.
function nm_l_int(l) ! ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す. ! ! 引数 l が整数値の場合, 対応する全波数と帯状波数を ! 長さ 2 の 1 次元整数値を返す. ! nm_l(1) が全波数, nm_l(2) が帯状波数である. ! integer :: nm_l_int(2) !(out) 全波数, 帯状波数 integer, intent(in) :: l !(in) スペクトルデータの格納位置 if ( .not. w_base_initialize ) then call MessageNotify('E','nm_l_int', 'w_base_module not initialize yet. Use sjl2nm routine in ISPACK directly.') else call sjl2nm(nn,l,nm_l_int(1),nm_l_int(2)) endif end function nm_l_int
Function : | |||
nm_l_array(size(larray),2) : | integer
| ||
larray(:) : | integer, intent(in)
|
スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
引数 larray が整数 1 次元配列の場合, larray に対応する n, m を格納した 2 次元整数配列を返す. nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である.
function nm_l_array(larray) ! ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す. ! ! 引数 larray が整数 1 次元配列の場合, ! larray に対応する n, m を格納した 2 次元整数配列を返す. ! nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である. ! integer, intent(in) :: larray(:) !(out) 全波数, 帯状波数 integer :: nm_l_array(size(larray),2) !(in) スペクトルデータの格納位置 integer :: i do i=1, size(larray) nm_l_array(i,:) = nm_l_int(larray(i)) enddo end function nm_l_array
Subroutine : | |||
w_Psi((mm+1)*(mm+1)) : | real(8), intent(in)
| ||
w_Chi((mm+1)*(mm+1)) : | real(8), intent(in)
| ||
xy_U(0:im-1,1:jm) : | real(8), intent(out)
| ||
xy_V(0:im-1,1:jm) : | real(8), intent(out)
|
流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に (逆)変換する(1 層用)
スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ, v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ
subroutine w_StreamPotential2Vector(w_Psi, w_Chi, xy_U, xy_V) ! ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に ! (逆)変換する(1 層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ, ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ ! real(8), intent(in) :: w_Psi((mm+1)*(mm+1)) !(in) 流線関数 real(8), intent(in) :: w_Chi((mm+1)*(mm+1)) !(in) 速度ポテンシャル real(8), intent(out) :: xy_U(0:im-1,1:jm) !(out) 速度経度成分 real(8), intent(out) :: xy_V(0:im-1,1:jm) !(out) 速度緯度成分 real(8) :: w_Rdata((mm+4)*mm+2) ! 作業用スペクトルデータ(SJTS2G 出力用) real(8) :: w_Xdata((mm+1)*(mm+1)) ! 作業用スペクトルデータ(SJCS2X 出力用) real(8) :: w_Ydata((mm+4)*mm+2) ! 作業用スペクトルデータ(SJCS2Y 出力用) real(8) :: ws2((2*nm+1-mm)*mm+nm+1) real(8) :: wg(im*jm) ! 変換用作業配列 real(8) :: w(jm*im) ! 変換用作業配列 if ( .not. w_base_initialize ) then call MessageNotify('E','w_StreamPotential2Vector', 'w_base_module_sjpack_cuda not initialize yet.') endif ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ の計算 ! call sjcs2x(mm,w_Chi,w_Xdata) call sjcs2y(mm,w_Psi,w_Ydata,c) call sjcrup(mm,nm,w_Xdata,w_Rdata) w_Rdata = w_Rdata - w_Ydata ! ! u の計算 ! call sjvs2g(mm,nm,nm,im,jm,w_Rdata,xy_U, it,t,p,r,ws2,wg,w,1,ip) ! ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ の計算 ! call sjcs2y(mm,w_Chi,w_Ydata,c) call sjcs2x(mm,w_Psi,w_Xdata) call sjcrup(mm,nm,w_Xdata,w_Rdata) w_Rdata= w_Rdata + w_Ydata ! ! v の計算 ! call sjvs2g(mm,nm,nm,im,jm,w_Rdata,xy_V, it,t,p,r,ws2,wg,w,1,ip) end subroutine w_StreamPotential2Vector
Subroutine : | |||
xy_U(0:im-1,1:jm) : | real(8), intent(in)
| ||
xy_V(0:im-1,1:jm) : | real(8), intent(in)
| ||
w_Vor((mm+1)*(mm+1)) : | real(8), intent(out)
| ||
w_Div((mm+1)*(mm+1)) : | real(8), intent(out)
|
速度場(格子データ)から渦度・発散(スペクトルデータ)に (正)変換する(1 層用)
スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
subroutine w_Vector2VorDiv(xy_U, xy_V, w_Vor, w_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(1 層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xy_U(0:im-1,1:jm) !(in) 速度経度成分 real(8), intent(in) :: xy_V(0:im-1,1:jm) !(in) 速度緯度成分 real(8), intent(out) :: w_Vor((mm+1)*(mm+1)) !(out) 流線関数 real(8), intent(out) :: w_Div((mm+1)*(mm+1)) !(out) 速度ポテンシャル real(8) :: w_Xdata((mm+1)*(mm+1)) ! 作業用スペクトルデータ(SJCS2X 出力用) real(8) :: w_Ydata((mm+4)*nm+2) ! 作業用スペクトルデータ(SJCY2S 出力用) real(8) :: w_Data1((mm+1)*(mm+1)) real(8) :: w_Data2((mm+1)*(mm+1)) real(8) :: ws2((2*nm+1-mm)*mm+nm+1) real(8) :: wg(im*jm) ! 変換用作業配列 real(8) :: w(jm*im) ! 変換用作業配列 if ( .not. w_base_initialize ) then call MessageNotify('E','w_Vector2VorDiv', 'w_base_module_sjpack_cuda not initialize yet.') endif ! ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ の計算 ! call sjvg2s(mm,nm,nm,im,jm,w_Ydata,xy_U, it,t,p,r,ws2,wg,w,1,ip) call sjcrdn(mm,nm,w_Ydata,w_Xdata) call sjcs2x(mm,w_Xdata,w_Div) call sjcy2s(mm,w_Ydata,w_Data1,c) ! ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算 ! call sjvg2s(mm,nm,nm,im,jm,w_Ydata,xy_V, it,t,p,r,ws2,wg,w,1,ip) call sjcrdn(mm,nm,w_Ydata,w_Xdata) call sjcs2x(mm,w_Xdata,w_Vor) call sjcy2s(mm,w_Ydata,w_Data2,c) ! ! 渦度・発散の計算 ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! w_Vor = w_Vor - w_Data1 w_Div = w_Div + w_Data2 end subroutine w_Vector2VorDiv
Subroutine : | |||
n_in : | integer,intent(in)
| ||
i_in : | integer,intent(in)
| ||
j_in : | integer,intent(in)
| ||
np_in : | integer,intent(in), optional
|
スペクトル変換の格子点数, 波数および OPENMP 使用時の 最大スレッド数を設定する.
実際の使用には上位サブルーチン w_Initial を用いること.
subroutine w_base_Initial(n_in,i_in,j_in,np_in) ! ! スペクトル変換の格子点数, 波数および OPENMP 使用時の ! 最大スレッド数を設定する. ! ! 実際の使用には上位サブルーチン w_Initial を用いること. ! integer,intent(in) :: i_in !(in) 格子点数(東西), 2の巾乗(<=2048) integer,intent(in) :: j_in !(in) 格子点数(南北), 4 の倍数 integer,intent(in) :: n_in !(in) 切断全波数 integer,intent(in), optional :: np_in !(in) OPENMP での最大スレッド数(ダミー) integer :: i, j im = i_in jm = j_in nn = n_in nm = n_in+1 mm = n_in ! default は三角波数切断 if ( present(np_in) )then if ( np_in .gt. 1 ) then call MessageNotify('W','w_base_Initial', 'OpenMP computation not implemented in w_module_sjpack_cuda.') endif endif if ( w_base_initialize == .true. ) then call MessageNotify('W','w_base_Initial', 'w_base_Initial called again.') call sjvcls(ip) deallocate(p,r,t,c) deallocate(x_Lon,x_Lon_Weight) deallocate(y_Lat,y_Lat_Weight) deallocate(xy_Lon,xy_Lat) endif allocate(p(jm/2,mm+4)) ! 変換用配列 allocate(r((mm+1)*(2*nm-mm-1)+1)) ! 変換用配列 allocate(t(im*6)) ! 変換用配列 allocate(c((mm+1)*(mm+1))) ! 変換用作業配列 allocate(x_Lon(0:im-1)) ! 格子点座標格納配列(経度) allocate(x_Lon_Weight(0:im-1)) allocate(xy_Lon(0:im-1,1:jm)) allocate(y_Lat(1:jm)) allocate(y_Lat_Weight(1:jm)) ! 格子点座標格納配列 allocate(xy_Lat(0:im-1,1:jm)) ! 格子点座標格納配列 call sjinit(mm,nm,jm,im,p,r,it,t) call sjvopn(mm,nm,jm,im,p,r,ip) call sjinic(mm,c) do i=0,im-1 x_Lon(i) = 2*pi/im*i ! 経度座標 x_Lon_Weight(i) = 2*pi/im ! 経度座標重み enddo do j=1,jm/2 y_Lat(jm/2+j) = asin(p(j,1)) ! 緯度座標 y_Lat(jm/2-j+1) = -asin(p(j,1)) ! 緯度座標 y_Lat_Weight(jm/2+j) = 2*p(j,2) ! 緯度重み(Gauss grid) y_Lat_Weight(jm/2-j+1) = 2*p(j,2) ! 緯度重み(Gauss grid) enddo do j=1,jm xy_Lon(:,j) = x_Lon enddo do i=0,im-1 xy_Lat(i,:) = y_Lat enddo w_base_initialize = .true. call MessageNotify('M','w_base_initial', 'w_base_module_sjpack_cuda (2013/02/12) is initialized') end subroutine w_base_Initial
Function : | |||
w_xy((mm+1)*(mm+1)) : | real(8)
| ||
xy_data(0:im-1,1:jm) : | real(8), intent(in)
| ||
ipow : | integer, intent(in), optional
| ||
iflag : | integer, intent(in), optional
|
格子データからスペクトルデータへ(正)変換する(1 層用).
function w_xy(xy_data,ipow,iflag) ! ! 格子データからスペクトルデータへ(正)変換する(1 層用). ! real(8) :: w_xy((mm+1)*(mm+1)) !(out) スペクトルデータ real(8), intent(in) :: xy_data(0:im-1,1:jm) !(in) 格子点データ integer, intent(in), optional :: ipow !(in) 変換時に同時に作用させる 1/cosφ の次数. 省略時は 0. integer, intent(in), optional :: iflag ! 変換の種類 ! 0 : 通常の正変換 ! -1 : 経度微分を作用させた正変換 ! 1 : 緯度微分 1/cosφ・∂(f cos^2φ)/∂φ を作用させた正変換 ! 2 : sinφを作用させた正変換 ! 省略時は 0. integer, parameter :: ipow_default = 0 ! スイッチデフォルト値 integer, parameter :: iflag_default = 0 ! スイッチデフォルト値 integer ipval, ifval !, i, j real(8) :: w_Rdata((2*nn+1-mm)*mm+nn+1) ! 作業用スペクトルデータ(SJTS2G 出力用) real(8) :: w_Xdata((mm+1)*(mm+1)) ! 作業用スペクトルデータ(SJCS2X 出力用) real(8) :: w_Ydata((mm+4)*nm+2) ! 作業用スペクトルデータ(SJCY2S 出力用) !! real(8) :: ws(2*(nn+1)) ! 変換用作業配列 real(8) :: ws((2*nn+1-mm)*mm+nn+1) real(8) :: ws2((2*nm+1-mm)*mm+nm+1) real(8) :: wg(im*jm) ! 変換用作業配列 real(8) :: w(jm*im) ! 変換用作業配列 if ( .not. w_base_initialize ) then call MessageNotify('E','xy_w', 'w_base_module_sjpack_cuda not initialize yet.') endif if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif if ( ifval == 0 ) then call sjvg2s(mm,nm,nn,im,jm,w_Rdata,xy_data, it,t,p,r,ws,wg,w,ipval,ip) call sjcrdn(mm,nn,w_Rdata,w_xy) else if ( ifval == -1 ) then call sjvg2s(mm,nm,nn,im,jm,w_Rdata,xy_data, it,t,p,r,ws,wg,w,ipval,ip) call sjcrdn(mm,nn,w_Rdata,w_Xdata) call sjcs2x(mm,w_Xdata,w_xy) else if ( ifval == 1 ) then call sjvg2s(mm,nm,nm,im,jm,w_Ydata,xy_data, it,t,p,r,ws2,wg,w,ipval,ip) call sjcy2s(mm,w_Ydata,w_xy,c) else if ( ifval == 2 ) then call sjvg2s(mm,nm,nn,im,jm,w_Rdata,xy_data*sin(xy_Lat), it,t,p,r,ws,wg,w,ipval,ip) call sjcrdn(mm,nn,w_Rdata,w_xy) end if end function w_xy
Function : | |||
xy_w(0:im-1,1:jm) : | real(8)
| ||
w_data((mm+1)*(mm+1)) : | real(8), intent(in)
| ||
ipow : | integer, intent(in), optional
| ||
iflag : | integer, intent(in), optional
|
スペクトルデータから格子データへ変換する(1 層用).
function xy_w(w_data,ipow,iflag) ! ! スペクトルデータから格子データへ変換する(1 層用). ! real(8) :: xy_w(0:im-1,1:jm) !(out) 格子点データ real(8), intent(in) :: w_data((mm+1)*(mm+1)) !(in) スペクトルデータ integer, intent(in), optional :: ipow !(in) 作用させる 1/cosφ の次数. 省略時は 0. integer, intent(in), optional :: iflag !(in) 変換の種類 ! 0 : 通常の正変換 ! -1 : 経度微分を作用させた逆変換 ! 1 : 緯度微分 cosφ・∂/∂φ を作用させた逆変換 ! 2 : sinφを作用させた逆変換 ! 省略時は 0. ! integer, parameter :: ipow_default = 0 integer, parameter :: iflag_default = 0 integer ipval, ifval real(8) :: w_Rdata((2*nn+1-mm)*mm+nn+1) ! 作業用スペクトルデータ(SJTS2G 出力用) real(8) :: w_Xdata((mm+1)*(mm+1)) ! 作業用スペクトルデータ(SJCS2X 出力用) real(8) :: w_Ydata((mm+4)*mm+2) ! 作業用スペクトルデータ(SJCS2Y 出力用) !!$ real(8) :: ws(2*(nn+1)) ! 変換用作業配列 real(8) :: ws((2*nn+1-mm)*mm+nn+1) real(8) :: ws2((2*nm+1-mm)*mm+nm+1) real(8) :: wg(im*jm) ! 変換用作業配列 real(8) :: w(jm*im) ! 変換用作業配列 if ( .not. w_base_initialize ) then call MessageNotify('E','xy_w', 'w_base_module_sjpack_cuda not initialize yet.') endif if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif if ( ifval==0 ) then call sjcrup(mm,nn,w_data,w_Rdata) call sjvs2g(mm,nm,nn,im,jm,w_Rdata,xy_w, it,t,p,r,ws,wg,w,ipval,ip) else if( ifval==-1 ) then call sjcs2x(mm,w_data,w_Xdata) call sjcrup(mm,nn,w_Xdata,w_Rdata) call sjvs2g(mm,nm,nn,im,jm,w_Rdata,xy_w, it,t,p,r,ws,wg,w,ipval,ip) else if( ifval==1 ) then call sjcs2y(mm,w_data,w_Ydata,c) call sjvs2g(mm,nm,nm,im,jm,w_Ydata,xy_w, it,t,p,r,ws2,wg,w,ipval,ip) else if( ifval==2 ) then call sjcrup(mm,nn,w_data,w_Rdata) call sjvs2g(mm,nm,nn,im,jm,w_Rdata,xy_w, it,t,p,r,ws,wg,w,ipval,ip) xy_w = xy_w * sin(xy_Lat) else call MessageNotify('E','xy_w','invalid value of iflag') endif end function xy_w