!-- !---------------------------------------------------------------------- ! Copyright (c) 2015 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 w_plane_module ! ! spml/w_module モジュールは 2 次元無限平面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! Fortran90 関数を提供する. ! ! 内部で w_module を呼んでいる. スペクトルデータおよび格子点データの ! 格納方法や変換の詳しい計算法については w_module の下部ルーチンである ! ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! !履歴 2015/12/26 竹広真一 新規作成 ! !++ module w_plane_module2 ! != w_plane_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id:$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/w_plane_module モジュールは 2 次元無限平面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための Fortran90 ! 関数を提供する. ! ! 内部で w_module を呼んでいる. スペクトルデータおよび格子点データの ! 格納方法や変換の詳しい計算法については w_module の下部ルーチンである ! ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (w_, nm_, n_, xy_, x_, y_) は, 返す値の形を示している. ! w_ :: スペクトルデータ ! xy_ :: 2 次元格子点データ ! ! * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv)は, その関数の作用を表している. ! ! * 関数名の最後 (_w_w, _w, _xy) は, 入力変数の形スペクトルデータ ! および格子点データであることを示している. ! _w :: スペクトルデータ ! _xy :: 2 次元格子点データ ! !=== 各データの種類の説明 ! ! * xy : 2 次元格子点データ. ! * 変数の種類と次元は real(8), dimension(0:im-1,1:jm). ! * im, jm はそれぞれ経度, 緯度座標の格子点数であり, サブルーチン ! w_plane_Initial にてあらかじめ設定しておく. ! ! * w : スペクトルデータ. ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1)). ! * nm は球面調和函数の最大全波数であり, サブルーチン w_Initial にて ! あらかじめ設定しておく. ! * スペクトルデータの格納のされ方は関数 l_nm, nm_l によって ! 調べることができる. ! ! * w_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xy_ で始まる関数が返す値は 2 次元格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! w_plane_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 終了処理 ! ! w_plane_Finalize :: モジュールの終了処理(割り付け配列の解放)をおこなう. ! !==== 座標変数 ! ! sp_S, sp_Phi :: 格子点データの 2 次元極座標(格子点データ型 2 次元配列) ! sp_X, sp_Y :: 格子点データの 2 次元 X, Y 座標(格子点データ型 2 次元配列) ! !==== 基本変換 ! ! sp_w :: スペクトルデータから格子データへの変換 ! w_sp :: 格子データからスペクトルデータへの変換 ! !==== 微分 ! ! w_Lapla_w :: スペクトルデータにラプラシアンを作用させる ! w_LaplaInv_w :: スペクトルデータにラプラシアンの逆変換を作用させる ! sp_Dx_w :: スペクトルデータに d/dx を作用させる ! sp_Dy_w :: スペクトルデータに d/dy を作用させる ! w_Dx_sp :: 格子データに d/dx を作用させる ! w_Dy_sp :: 格子データに d/dy を作用させる ! !==== 補間 ! ! Interpolate_w :: スペクトルデータから任意の点での値を求める. ! !==== 積分・平均 ! ! IntSPhi_sp, AvrSPhi_sp :: 2 次元格子点データの全領域積分および平均 ! use dc_message, only : MessageNotify use w_module, only : p_Phi=>x_Lon, s_Lat=>y_Lat, & ps_Non_w=>xy_w, w_Non_ps=>w_xy, & w_Initial, w_Finalize, l_nm, nm_l, & w_Lapla_w, w_LaplaInv_w,w_DPhi_w => w_DLon_w, & xy_GradLambda_w, xy_GradMu_w, w_DivMu_xy, & InterpolateOnLonLat_w => Interpolate_w implicit none integer :: im=64 ! 格子点の設定(東西) integer :: jm=32 ! 格子点の設定(南北) integer :: nm=21 ! 切断波数の設定 real(8) :: Radius ! 球半径 real(8), allocatable :: s_S(:) ! 2 次元極座標動径 real(8), allocatable :: ps_S(:,:), ps_Phi(:,:) real(8), allocatable :: s_Mu(:) ! 球面緯度(sin φ) real(8), allocatable :: ps_Mu(:,:) ! 球面緯度(sin φ) real(8), allocatable :: ps_X(:,:), ps_Y(:,:) ! 2 次元直交座標 private public w_plane_Initial ! 初期化 public w_plane_Finalize ! 終了処理 public s_S, p_Phi ! 格子座標(im,jm) public ps_S, ps_Phi ! 格子座標(im,jm) public ps_X, ps_Y ! 格子座標(im,jm) public ps_w, w_ps ! 変換関数 public l_nm, nm_l ! スペクトル成分位置情報 public ps_Lapla_w, w_LaplaInv_ps ! ラプラシアンと逆演算 public w_DPhi_w ! 方位角微分 public ps_DS_w, ps_DPhi_w ! 微分(極座標) public w_DPhi_ps ! 微分(極座標) public ps_DX_w, ps_DY_w ! 微分(直交座標) public w_Lapla_w, w_LaplaInv_w ! 球面上でのラプラシアンと逆演算 !!$ public w_DX_ps, w_DY_ps ! 微分(直交座標) public Interpolate_w ! ps 座標での補間 public InterpolateOnXY_w ! xy 座標での補間 public xy_w ! xy 座標での補間 save im, jm, nm, Radius, s_S, ps_S, ps_Phi, s_Mu, ps_Mu, ps_X, ps_Y contains !--------------- 初期化 ----------------- subroutine w_plane_initial(n_in,i_in,j_in,r_in,np_in) ! ! スペクトル変換の格子点数, 波数および OPENMP 使用時の ! 最大スレッド数を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! ! np_in に 1 より大きな値を指定すれば ISPACK の球面調和函数変換 ! OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには, ! 実行時に環境変数 OMP_NUM_THREADS を np_in 以下の数字に設定する等の ! システムに応じた準備が必要となる. ! ! np_in に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない. ! integer,intent(in) :: i_in !(in) 格子点数(方位角) integer,intent(in) :: j_in !(in) 格子点数(動径) integer,intent(in) :: n_in !(in) 切断波数の設定 real(8),intent(in) :: r_in !(in) 変換パラメター(球の半径) integer,intent(in), optional :: np_in !(in) OPENMP での最大スレッド数 im = i_in ; jm = j_in ; nm = n_in ; Radius = r_in if ( present (np_in) )then call w_initial(n_in,i_in,j_in,np_in) else call w_initial(n_in,i_in,j_in) endif allocate(s_S(jm)) allocate(s_Mu(jm),ps_Mu(0:im-1,jm)) allocate(ps_S(0:im-1,jm),ps_Phi(0:im-1,jm)) allocate(ps_X(0:im-1,jm),ps_Y(0:im-1,jm)) s_Mu = sin(s_Lat) ; ps_Mu = spread(s_Mu,1,im) s_S = 2*Radius*sqrt((1+s_Mu)/(1-s_Mu)) ps_S = spread(s_S,1,im) ps_Phi = spread(p_Phi,2,jm) ps_X = ps_S * cos(ps_Phi) ; ps_Y = ps_S * sin(ps_Phi) call MessageNotify('M','w_plane_initial',& 'w_plane_module2 (2016/01/15) is initialized') end subroutine w_plane_initial !--------------- 終了処理 ----------------- subroutine w_plane_Finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! 解像度を変更する際にはこのサブルーチンを呼んで終了処理を ! おこなったのちに再度 w_Initial で初期設定しなければ ! ならない. ! call w_Finalize call MessageNotify('M','w_plane_Finalize',& 'w_plane_module2 (2015/12/26) is finalized') end subroutine w_plane_Finalize !--------------- 基本変換 ----------------- function w_ps(ps_Data) real(8), intent(IN) :: ps_Data(0:im-1,jm) real(8) :: w_ps((nm+1)**2) w_ps = w_Non_ps(ps_Data/(1-ps_Mu)) end function w_ps function ps_w(w_Data) real(8), intent(IN) :: w_Data((nm+1)**2) real(8) :: ps_w(0:im-1,jm) ps_w = (1-ps_Mu)*ps_Non_w(w_Data) end function ps_w !--------------- 微分計算 ----------------- function ps_Lapla_w(w_data) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = (1/s)(d/ds)[s(d/ds)] + 1/s^2 (d^2/d\lambda^2) ! ! を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8) :: ps_Lapla_w(0:im-1,jm) !(out) 入力スペクトルデータのラプラシアン real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ ps_Lapla_w = (1-ps_Mu)**2/(4*Radius**2) & * ( ps_Non_w(w_Lapla_w(w_Non_ps(ps_w(w_data))))) !!$ ps_Lapla_w = (1-ps_Mu)**3/(4*Radius**2) * ps_Non_w(w_Lapla_w(w_data)) & !!$ - (1-ps_Mu)**2/(4*Radius**2) & !!$ * ps_Non_w(w_DivMu_xy((1-ps_Mu**2)*ps_Non_w(w_data))) end function ps_Lapla_w function w_LaplaInv_ps(ps_data) ! ! 入力スペクトルデータに逆ラプラシアン ! ! ▽^{-2} ! =[(1/s)(d/ds)[s(d/ds)] + 1/s^2 (d^2/d\lambda^2)]^{-1} ! ! を作用する. ! ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8) :: w_LaplaInv_ps((nm+1)*(nm+1)) !(out) スペクトルデータの逆ラプラシアン real(8), intent(in) :: ps_data(0:im-1,jm) !(in) 入力格子データ w_LaplaInv_ps = w_ps(ps_Non_w(w_LaplaInv_w(w_Non_ps(ps_data*4*Radius**2/(1-ps_Mu)**2)))) end function w_LaplaInv_ps !--------------- 微分計算 (S,Phi座標系用) ----------------- function ps_DS_w(w_data) ! ! スペクトルデータに動径微分 d/ds ! を作用させて格子点データに変換して返す. ! real(8) :: ps_Ds_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型緯度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ ps_Ds_w = xy_GradMu_w(w_data)/ps_S * (1-ps_Mu) & - ps_Non_w(w_data) * (1-ps_Mu**2)/ps_S end function ps_DS_w function ps_DPhi_w(w_data) ! ! スペクトルデータに動径微分 d/ds ! を作用させて格子点データに変換して返す. ! real(8) :: ps_DPhi_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型緯度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ ps_DPhi_w = xy_GradLambda_w(w_data) * (1-ps_Mu) end function ps_DPhi_w function w_DPhi_ps(ps_data) ! ! 格子点データに方位角微分 d/d\lambda ! を作用させてスペクトルデータに変換して返す(1 層用). ! real(8) :: w_DPhi_ps((nm+1)*(nm+1)) !(out) 格子点データを発散型経度微分したスペクトルデータ real(8), intent(in) :: ps_data(0:im-1,1:jm) !(in) 入力格子点データ w_DPhi_ps = w_Non_ps(ps_data/(1-ps_Mu),ipow=0,iflag=-1) end function w_DPhi_ps !!$ function w_sDS_ps(ps_data) !!$ ! !!$ ! 格子点データに動径微分 s(d/ds) を作用させて !!$ ! スペクトルデータに変換して返す. !!$ ! !!$ real(8) :: w_sDS_ps((nm+1)*(nm+1)) !!$ !(out) 格子点データを発散型緯度微分したスペクトルデータ !!$ !!$ real(8), intent(in) :: ps_data(0:im-1,1:jm) !!$ !(in) 入力格子点データ !!$ !!$ w_sDS_ps = w_ps(ps_data,ipow=2,iflag=1) !!$ !!$ end function w_sDS_ps !--------------- 微分計算 (xy 座標系用) ----------------- function ps_DX_w(w_data) ! ! スペクトルデータに微分 d/dx ! を作用させて格子点データに変換して返す. ! real(8) :: ps_DX_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型緯度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ ps_DX_w = cos(ps_Phi)*ps_DS_w(w_data) & - sin(ps_Phi)/ps_S * ps_DPhi_w(w_data) end function ps_DX_w function ps_DY_w(w_data) ! ! スペクトルデータに微分 d/dy ! を作用させて格子点データに変換して返す. ! real(8) :: ps_DY_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型緯度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ ps_DY_w = sin(ps_Phi)*ps_DS_w(w_data) & + cos(ps_Phi)/ps_S * ps_DPhi_w(w_data) end function ps_DY_w !--------------- 補間計算 ----------------- function Interpolate_w(w_data,pval,sval) real(8), intent(IN) :: w_data((nm+1)*(nm+1)) ! real(8), intent(IN) :: pval ! 補間する位置(方位角) real(8), intent(IN) :: sval ! 補間する位置(動径) real(8) :: Interpolate_w ! 補間値 real(8) :: lat, t t = (sval/(2*Radius))**2 lat = asin((t-1)/(t+1)) Interpolate_w = InterpolateOnLonLat_w(w_data,pval,lat)*(1-sin(lat)) end function Interpolate_w function InterpolateOnXY_w(w_data,xval,yval) real(8), intent(IN) :: w_data((nm+1)*(nm+1)) ! real(8), intent(IN) :: xval ! 補間する位置(X) real(8), intent(IN) :: yval ! 補間する位置(Y) real(8) :: InterpolateOnXY_w ! 補間値 real(8) :: lon, lat real(8) :: sval, t real(8), parameter :: PI=3.141592653589793 sval = sqrt(xval**2+yval**2) t = (sval/(2*Radius))**2 lat = asin((t-1)/(t+1)) if ( xval .eq. 0.0D0 .AND. yval .eq. 0.0D0 ) then lon = 0.0D0 else if ( xval .eq. 0.0D0 .AND. yval .gt. 0.0D0 ) then lon = PI/2.0D0 else if ( xval .eq. 0.0D0 .AND. yval .lt. 0.0D0 ) then lon = 3.0D0*PI/2.0D0 else if ( xval .gt. 0.0D0 .AND. yval .ne. 0.0D0 ) then lon = atan(yval/xval) else if ( xval .lt. 0.0D0 .AND. yval .ne. 0.0D0 ) then lon = atan(yval/xval) + PI endif if ( -1.0D0/2.0D0*PI .lt. lon .AND. lon .lt. 0.0 ) then lon = lon + PI endif InterpolateOnXY_w = InterpolateOnLonLat_w(w_data,lon,lat)*(1-sin(lat)) end function InterpolateOnXY_w subroutine xy_w( w_data, imx, jmy, xmin, xmax, ymin, ymax, & xy_Data, x_X, y_Y) real(8), intent(IN) :: w_data((nm+1)*(nm+1)) ! 入力スペクトルデータ integer, intent(IN) :: imx ! X 格子点数 integer, intent(IN) :: jmy ! Y 格子点数 real(8), intent(IN) :: xmin, xmax ! X 座標範囲 real(8), intent(IN) :: ymin, ymax ! Y 座標 real(8), intent(OUT) :: xy_Data(imx,jmy) ! XY データ real(8), intent(OUT) :: x_X(imx) ! X 格子点座標値 real(8), intent(OUT) :: y_Y(jmy) ! Y 格子点座標値 integer :: i,j do i=1,imx x_X(i) = xmin + (xmax-xmin)/(imx-1)*(i-1) end do do j=1,jmy y_Y(j) = ymin + (ymax-ymin)/(jmy-1)*(j-1) end do do j=1,jmy do i=1,imx xy_Data(i,j) = InterpolateOnXY_w(w_data,x_X(i),y_Y(j)) enddo enddo end subroutine xy_w end module w_plane_module2