!-- !---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- !表題 w_wave_base_module_svpack ! ! spml/w_base_module_svpack モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するためのモジュール ! w_module_svpack の下部モジュールであり, スペクトル計算の基本的な ! Fortran90 関数を提供する. ! ! 内部で ISPACK の LJPACK(SVPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算 ! 法については ISPACK/SVPACK のマニュアルを参照されたい. ! !== 履歴 ! ! 2017/03/09 竹広真一 w_wave_base_module_sjpack より改造 ! ! 制限 ! ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち ! ・波数切断の仕方は三角波数切断に決めうち. ! !++ module w_wave_base_module_svpack ! != w_wave_base_module_svpack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id: w_base_module_svpack.f90 586 2013-05-23 17:42:12Z uwabami $ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要. ! ! spml/w_base_module_svpack モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するためのモジュール ! w_module_svpack の下部モジュールであり, スペクトル計算の基本的な ! Fortran90 関数を提供する. ! ! 内部で ISPACK の LJPACK(SVPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算 ! 法については ISPACK/SVPACK のマニュアルを参照されたい. ! use dc_message, only : MessageNotify use iso_c_binding implicit none integer :: im=2 ! 東西波数成分 integer :: jm=32 ! 格子点の設定(南北) integer :: nm=21 ! 計算する最大の全波数の設定 integer :: nn=22 ! 切断波数(全波数)の設定 integer :: mm=21 ! 切断波数(東西波数)の設定 integer :: m = 1 ! 東西波数の設定 real(8), allocatable :: p(:,:), r(:) ! 変換用配列 integer, allocatable :: jc(:) ! 変換用配列 real(8), allocatable :: c(:) ! 作業配列 real(8), allocatable :: y_Lat(:) ! 緯度経度 real(8), allocatable :: y_Lat_Weight(:) ! 座標重み real(8), allocatable :: xy_Lat(:,:) real(8), pointer :: xy_pdata(:,:),q(:) type(C_PTR) :: pg, pq logical :: w_wave_base_initialize=.false. ! 初期化フラッグ real(8), parameter :: pi=3.1415926535897932385D0 private public im, jm, nn, mm, nm, m ! 格子点数, 切断波数 public p, r ! 変換用作業配列 public w_base_Initial ! 初期化サブルーチン public w_base_Finalize ! 終了サブルーチン public y_Lat ! 格子座標 public y_Lat_Weight ! 格子座標重み public xy_Lat ! 格子座標(im,jm) public l_nm, nm_l ! 波数格納位置 public xy_w, w_xy ! 変換関数 public w_DLon_w ! 経度微分 interface l_nm module procedure l_nm_array00 module procedure l_nm_array01 module procedure l_nm_array10 module procedure l_nm_array11 end interface interface nm_l module procedure nm_l_int module procedure nm_l_array end interface save im, jm, nm, mm, nn, m ! 格子点数, 切断波数を記憶 save p, r ! 変換用配列を記憶 save jc ! 変換用配列の大きさ save c ! 変換用配列の大きさ save xy_pdata, q save w_wave_base_initialize ! 初期化フラグ contains !--------------- 初期化 ----------------- subroutine w_base_Initial(n_in,m_in,i_in,j_in) ! ! スペクトル変換の格子点数, 波数を設定する. ! ! 実際の使用には上位サブルーチン w_wave_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) :: m_in !(in) 東西波数 integer :: i, j real(8), allocatable :: wp(:) ! 作業領域 if ( m <= 0 ) then call MessageNotify('E','w_wave_base_module_initial', & 'longitudinal wavenumber m must be positive') endif im = i_in ; jm = j_in nn = n_in ; nm = n_in+1 ; mm = n_in ! default は三角波数切断 m = m_in ! 東西波数は 1 成分 allocate(p(jm/2,5+2*mm)) ! 変換用配列 allocate(r(((mm+1)*(2*nm-mm-1)+1)/4*3+(2*nm-mm)*(mm+1)/2)) ! 変換用配列 allocate(jc(mm*(2*nm-mm-1)/8+mm)) ! 微分用配列 allocate(c((nn+1)*(nn+1))) ! 変換用配列 allocate(wp(jm/2*mm)) ! 作業領域 allocate(y_Lat(1:jm)) allocate(y_Lat_Weight(1:jm)) ! 格子点座標格納配列 allocate(xy_Lat(0:im-1,1:jm)) ! 格子点座標格納配列 call mlallc(pg,jm*im) call mlallc(pq,jm/2*7) call c_f_pointer(pg, xy_pdata, [IM,JM]) call c_f_pointer(pq, q, [jm/2*7]) call lvinit(mm,nm,jm,p,r,jc,wp) call svinic(nn,c) 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 i=0,im-1 xy_Lat(i,:) = y_Lat enddo w_wave_base_initialize = .true. call MessageNotify('M','w_wave_base_initial',& 'w_wave_base_module_svpack (2017/03/09) for AVX is initialized') end subroutine w_base_Initial !--------------- 基本関数 ----------------- function l_nm_array00(n_in,m_in) ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 引数 n,m がともに整数値の場合, 整数値を返す. ! integer :: l_nm_array00 !(out) スペクトルデータの格納位置 integer, intent(in) :: n_in !(in) 全波数 integer, intent(in) :: m_in !(in) 帯状波数 if ( m_in >= 0 ) then l_nm_array00 = 2*(n_in-m_in+1)-1 else l_nm_array00 = 2*(n_in+m_in+1) endif end function l_nm_array00 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(narray,m_in) ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 第 1 引数 narray が整数 1 次元配列, 第 2 引数 m が整数の場合, ! narray と同じ大きさの 1 次元整数配列を返す. ! integer, intent(in) :: narray(:) !(in) 全波数 integer, intent(in) :: m_in !(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_in) enddo end function l_nm_array10 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(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 ( mod(l,2) == 0 ) then nm_l_int(1) = l/2 + m-1 nm_l_int(2) = -m else nm_l_int(1) = (l+1)/2 + m-1 nm_l_int(2) = m endif end function nm_l_int 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 !--------------- 基本変換 ----------------- function w_DLon_w(w_data) ! ! 入力スペクトルデータに経度微分∂/∂λを作用する. ! ! スペクトルデータの経度微分とは, 対応する格子点データに ! 経度微分を作用させたデータのスペクトル変換のことである. ! real(8) :: w_DLon_w(2*(nn-m+1)) !(out) 入力スペクトルデータのラプラシアン real(8), intent(in) :: w_data(2*(nn-m+1)) !(in) 入力スペクトルデータ integer :: n, n1, n2 do n=m,nn n1 = 2*(n-m+1)-1 n2 = 2*(n-m+1) w_DLon_w(n2) = m * w_data(n1) w_DLon_w(n1) = - m * w_data(n2) end do end function w_DLon_w function xy_w(w_data,ipow,iflag) ! ! スペクトルデータから格子データへ変換する(1 層用). ! real(8) :: xy_w(0:im-1,1:jm) !(out) 格子点データ real(8), intent(in) :: w_data(2*(nn-m+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_Ydata(2*(nn-m+2)) ! 作業用スペクトルデータ(lJCS2Y 出力用) real(8) :: ws(2*(nn-m+1)) ! 変換用作業配列 real(8) :: ws2(2*(nm-m+1)) ! 変換用作業配列 if ( .not. w_wave_base_initialize ) then call MessageNotify('E','xy_w',& 'w_wave_base_module_svpack 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 lvtswg(nm,nn,jm,m,w_Data,xy_pdata,p,q,r,jc,ws,ipval) xy_w = xy_pdata else if( ifval==-1 ) then call lvtswg(nm,nn,jm,m,w_DLon_w(w_data),xy_pdata,p,q,r,jc,ws,ipval) xy_w = xy_pdata else if( ifval==1 ) then call svcswy(nn,m,w_data,w_Ydata,c) call lvtswg(nm,nm,jm,m,w_Ydata,xy_pdata,p,q,r,jc,ws2,ipval) xy_w = xy_pdata else if( ifval==2 ) then call lvtswg(nm,nn,jm,m,w_data,xy_pdata,p,q,r,jc,ws,ipval) xy_w = xy_pdata * sin(xy_Lat) else call MessageNotify('E','xy_w','invalid value of iflag') endif end function xy_w function w_xy(xy_data,ipow,iflag) ! ! 格子データからスペクトルデータへ(正)変換する(1 層用). ! real(8) :: w_xy(2*(nn-m+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 real(8) :: w_Ydata(2*(nn-m+2)) ! 作業用スペクトルデータ(LJCY2S 出力用) real(8) :: ws(2*(nn-m+1)) ! 変換用作業配列 real(8) :: ws2(2*(nm-m+1)) ! 変換用作業配列 if ( .not. w_wave_base_initialize ) then call MessageNotify('E','xy_w',& 'w_wave_base_module_svpack 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 xy_pdata = xy_data call lvtgws(nm,nn,jm,m,w_xy,xy_pdata,p,q,r,jc,ws,ipval) else if ( ifval == -1 ) then xy_pdata = xy_data call lvtgws(nm,nn,jm,m,w_xy,xy_pdata,p,q,r,jc,ws,ipval) w_xy = w_Dlon_w(w_xy) else if ( ifval == 1 ) then xy_pdata = xy_data call lvtgws(nm,nm,jm,m,w_Ydata,xy_pdata,p,q,r,jc,ws2,ipval) call svcyws(nn,m,w_Ydata,w_xy,c) else if ( ifval == 2 ) then xy_pdata = xy_data*sin(xy_Lat) call lvtgws(nm,nn,jm,m,w_xy,xy_pdata,p,q,r,jc,ws,ipval) end if end function w_xy !--------------- 終了処理 ----------------- subroutine w_base_Finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! 実際の使用には上位サブルーチン w_Finalize を用いること. ! if ( .not. w_wave_base_initialize ) then call MessageNotify('W','w_base_Finalize',& 'w_wave_base_module not initialized yet') return endif deallocate(p) ! 変換用配列 deallocate(r) ! 変換用配列 deallocate(jc) ! 変換用作業配列 deallocate(c) ! 変換用作業配列 deallocate(y_Lat) deallocate(y_Lat_Weight) ! 格子点座標格納配列 deallocate(xy_Lat) ! 格子点座標格納配列 nullify(xy_pdata) ! 変換用作業配列 nullify(q) ! 変換用作業配列 call mlfree(pg) ! 変換用作業配列 call mlfree(pq) ! 変換用作業配列 w_wave_base_initialize = .false. call MessageNotify('M','w_base_Finalize',& 'w_wave_base_module_svpack (2017/03/09) for AVX is finalized') end subroutine w_base_Finalize end module w_wave_base_module_svpack