!---------------------------------------------------------------------- ! Copyright (c) 2008-2013 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 w_base_mpi_module ! ! spml/w_base_mpi_module モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法と MPI によって数値計算するための ! モジュール w_mpi_module の下部モジュールであり, スペクトル計算の ! 基本的な Fortran90 関数を提供する. ! ! 内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! ! !履歴 2008/05/26 竹広真一 w_base_module を MPI 化 ! 2010/01/07 佐々木洋平 RDoc 用のドキュメント修正, ! 2012/03/30 竹広真一 コメント修正 ! 2013/02/12 竹広真一 w_StreamPotential2VectorMPI, ! w_Vector2VorDivMPI 導入 ! 2013/02/15 竹広真一 w_VectorCosLat2VorDivMPI 導入 ! 2013/02/23 竹広真一 w_base_mpi_Finalize 導入 ! module w_base_mpi_module ! ! w_base_mpi_module ! ! spml/w_base_mpi_module モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法と MPI によって数値計算するための ! モジュール w_mpi_module の下部モジュールであり, スペクトル法の ! 基本的なな Fortran90 関数を提供する. ! ! 内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! use dc_message use w_base_module, only : im, jm, nm, x_Lon implicit none integer :: it(6) ! 変換用配列(分散格子点用) real(8), allocatable :: t(:) ! 変換用配列(分散格子点用) integer, allocatable :: ip(:) ! 変換用配列(分散格子点用) real(8), allocatable :: p(:), r(:) ! 変換用配列(分散格子点用) integer, allocatable :: ia(:) ! 変換用配列(分散格子点用) real(8), allocatable :: a(:) ! 変換用配列(分散格子点用) real(8), allocatable :: y(:) ! 変換用配列(分散格子点用) integer :: jc ! 分散配置用変数 real(8), allocatable :: yy(:,:) ! 変換用配列 real(8), allocatable :: q(:) ! 作業配列 real(8), allocatable :: ww(:), ws(:) ! 作業配列 real(8), allocatable :: w(:) ! 作業配列 real(8), allocatable :: v_Lat(:),v_Lat_Weight(:) ! 緯度経度 real(8), allocatable :: xv_Lon(:,:), xv_Lat(:,:) real(8), allocatable :: xv_work(:,:) ! w_xv,xv_w 変換用配列 integer :: id=65, jd=33 ! xv_work の大きさ real(8), parameter :: pi=3.1415926535897932385D0 logical :: w_base_initialize=.false. ! 初期化フラッグ private private im, jm, nm ! Intel Fortran 対策 public it, t, y, ip, p, r, ia, a ! 変換用作業配列 public id, jd ! 作業用配列の大きさ public jc ! 分散配置情報 public w_base_mpi_Initial ! 初期化サブルーチン public w_base_mpi_Finalize ! 終了サブルーチン public v_Lat, v_Lat_Weight ! 経度分散格子座標・重み public xv_Lon, xv_Lat ! 分散格子座標(im,jc) public xv_w, w_xv ! 変換関数 public w_StreamPotential2VectorMPI ! 流線ポテンシャルから速度場計算 public w_Vector2VorDivMPI ! 速度場から渦度発散を計算 public w_VectorCosLat2VorDivMPI ! 速度場から渦度発散を計算 save it, t, y, ip, p, r, ia, a ! 変換用配列を記憶 save jc ! 分散格子点配列の大きさ save id, jd ! 変換用配列の大きさ save w_base_initialize ! 初期化フラグ contains !--------------- 初期化 ----------------- subroutine w_base_mpi_Initial ! ! スペクトル変換の格子点数, 波数を設定する. ! ! 実際の使用には上位サブルーチン w_mpi_Initial を用いること. ! integer :: iw, i, j allocate(t(im*2)) ! 変換用配列(分散配置) allocate(ip(((nm+1)/2+nm+1)*2)) ! 変換用配列(分散配置) allocate(p(((nm+1)/2+nm+1)*jm)) ! 変換用配列(分散配置) allocate(r(((nm+1)/2*2+3)*(nm/2+1))) ! 変換用配列(分散配置) allocate(ia((nm+1)*(nm+1)*4)) ! 変換用配列(分散配置) allocate(a((nm+1)*(nm+1)*6)) ! 変換用配列(分散配置) allocate(y(jm*2)) ! 変換用配列(分散配置) ! 注意 : 別ルーチンによって w_base_Initial が呼んであることを仮定 call snmini(nm,im,jm,jc,it,t,y,ip,p,r,ia,a) if ( im/2*2 .eq. im ) then id = im+1 else id = im endif if ( jc/2*2 .eq. jc ) then jd = jc+1 else jd = jc endif allocate(xv_work(id,jd)) ! 変換用配列 allocate(q(((nm+1)/2+nm+1)*jm)) ! 作業配列 iw=max((nm+4)*(nm+3),jd*3*(nm+1),jd*im) allocate(ws(iw),ww(iw), w((nm+1)*(nm+1))) ! 作業用配列 allocate(yy(jc/2,4)) ! 変換用配列 allocate(v_Lat(jc),v_Lat_Weight(jc)) ! 格子点座標格納配列 allocate(xv_Lon(0:im-1,jc),xv_Lat(0:im-1,jc)) ! 格子点座標格納配列 yy = reshape(y(1:2*jc),(/jc/2,4/)) do j=1,jc/2 v_Lat(jc/2+j) = asin(yy(j,1)) ! 緯度座標 v_Lat(jc/2-j+1) = -asin(yy(j,1)) ! 緯度座標 v_Lat_Weight(jc/2+j) = 2*yy(j,2) ! 緯度重み(Gauss grid) v_Lat_Weight(jc/2-j+1) = 2*yy(j,2) ! 緯度重み(Gauss grid) enddo do j=1,jc xv_Lon(:,j) = x_Lon enddo do i=0,im-1 xv_Lat(i,:) = v_Lat enddo w_base_initialize = .true. call MessageNotify('M','w_base_mpi_initial',& 'w_base_mpi_module (2013/02/23) is initialized') end subroutine w_base_mpi_Initial !--------------- 基本変換(格子点分散配置) ----------------- function xv_w(w_data,ipow,iflag) ! ! スペクトルデータから分散格子データへ変換する(1 層用). ! real(8) :: xv_w(0:im-1,jc) !(out) 格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) スペクトルデータ integer, intent(in), optional :: ipow !(in) 作用させる 1/cosφ の次数. 省略時は 0. integer, intent(in), optional :: iflag !(in) 変換の種類 ! 0 : 通常の正変換 ! 1 : 経度微分を作用させた正変換 ! -1 : 緯度微分を作用させた正変換 ! 2 : sinφを作用させた正変換 ! 省略時は 0. ! integer, parameter :: ipow_default = 0 integer, parameter :: iflag_default = 0 integer ipval, ifval if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif call snts2g(nm,im,id,jc,jd,1,w_data,xv_work,& it,t,y,ip,p,r,ia,a,q,ws,ww,ipval,ifval) xv_w=xv_work(1:im,1:jc) end function xv_w function w_xv(xv_data,ipow,iflag) ! ! 格子データからスペクトルデータへ(正)変換する(1 層用). ! real(8) :: w_xv((nm+1)*(nm+1)) !(out) スペクトルデータ real(8), intent(in) :: xv_data(0:im-1,jc) !(in) 格子点データ integer, intent(in), optional :: ipow !(in) 変換時に同時に作用させる 1/cosφ の次数. 省略時は 0. integer, intent(in), optional :: iflag ! 変換の種類 ! 0 : 通常の正変換 ! 1 : 経度微分を作用させた正変換 ! -1 : 緯度微分を作用させた正変換 ! 2 : sinφを作用させた正変換 ! 省略時は 0. integer, parameter :: ipow_default = 0 ! スイッチデフォルト値 integer, parameter :: iflag_default = 0 ! スイッチデフォルト値 integer ipval, ifval if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif xv_work(1:im,1:jc)=xv_data call sntgms(nm,im,id,jc,jd,1,xv_work,w_xv,& it,t,y,ip,p,r,ia,a,q,ws,ww,ipval,ifval,w) end function w_xv !----------- 速度, 渦度・発散, 流線・ポテンシャル計算 ------------- subroutine w_StreamPotential2VectorMPI(w_Psi, w_Chi, xv_U, xv_V) ! ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に ! (逆)変換する(1 層用, MPI) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ, ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ ! real(8), intent(in) :: w_Psi((nm+1)*(nm+1)) !(in) 流線関数 real(8), intent(in) :: w_Chi((nm+1)*(nm+1)) !(in) 速度ポテンシャル real(8), intent(out) :: xv_U(0:im-1,1:jc) !(out) 速度経度成分 real(8), intent(out) :: xv_V(0:im-1,1:jc) !(out) 速度緯度成分 if ( .not. w_base_initialize ) then call MessageNotify('E','w_StreamPotential2VectorMPI',& 'w_base_module not initialize yet.') endif ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ の計算 ! xv_U = xv_w(w_Chi,ipow=1,iflag=-1) - xv_w(w_Psi,ipow=1,iflag=1) ! ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ の計算 ! xv_V = xv_w(w_Chi,ipow=1,iflag=1) + xv_w(w_Psi,ipow=1,iflag=-1) end subroutine w_StreamPotential2VectorMPI subroutine w_Vector2VorDivMPI(xv_U, xv_V, w_Vor, w_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(1 層用, MPI) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xv_U(0:im-1,1:jm) !(in) 速度経度成分 real(8), intent(in) :: xv_V(0:im-1,1:jm) !(in) 速度緯度成分 real(8), intent(out) :: w_Vor((nm+1)*(nm+1)) !(out) 流線関数 real(8), intent(out) :: w_Div((nm+1)*(nm+1)) !(out) 速度ポテンシャル if ( .not. w_base_initialize ) then call MessageNotify('E','w_Vector2VorDivMPI',& 'w_base_module not initialize yet.') endif ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! w_Vor = w_xv(xv_V,ipow=1,iflag=-1) - w_xv(xv_U,ipow=1,iflag=1) ! ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! w_Div = w_xv(xv_U,ipow=1,iflag=-1) + w_xv(xv_V,ipow=1,iflag=1) end subroutine w_Vector2VorDivMPI subroutine w_VectorCosLat2VorDivMPI(xv_UCosLat, xv_VCosLat, w_Vor, w_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(1 層用, MPI) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xv_UCosLat(0:im-1,1:jm) !(in) 速度経度成分 * cos(lat) real(8), intent(in) :: xv_VCosLat(0:im-1,1:jm) !(in) 速度緯度成分 * cos(lat) real(8), intent(out) :: w_Vor((nm+1)*(nm+1)) !(out) 流線関数 real(8), intent(out) :: w_Div((nm+1)*(nm+1)) !(out) 速度ポテンシャル if ( .not. w_base_initialize ) then call MessageNotify('E','w_VectorCosLat2VorDivMPI',& 'w_base_module not initialize yet.') endif ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! w_Vor = w_xv(xv_VCosLat,ipow=2,iflag=-1) & - w_xv(xv_UCosLat,ipow=2,iflag=1) ! ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! w_Div = w_xv(xv_UCosLat,ipow=2,iflag=-1) & + w_xv(xv_VCosLat,ipow=2,iflag=1) end subroutine w_VectorCosLat2VorDivMPI !--------------- 終了処理 ----------------- subroutine w_base_mpi_Finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! 実際の使用には上位サブルーチン w_Finalize を用いること. ! if ( .not. w_base_initialize ) then call MessageNotify('W','w_base_mpi_Finalize',& 'w_base_mpi_module not initialized yet') return endif deallocate(t) ! 変換用配列(分散配置) deallocate(ip) ! 変換用配列(分散配置) deallocate(p) ! 変換用配列(分散配置) deallocate(r) ! 変換用配列(分散配置) deallocate(ia) ! 変換用配列(分散配置) deallocate(a) ! 変換用配列(分散配置) deallocate(y) ! 変換用配列(分散配置) deallocate(xv_work) ! 変換用配列 deallocate(q) ! 作業配列 deallocate(ws,ww,w) ! 作業用配列 deallocate(yy) ! 変換用配列 deallocate(v_Lat,v_Lat_Weight) ! 格子点座標格納配列 deallocate(xv_Lon,xv_Lat) ! 格子点座標格納配列 w_base_initialize = .false. call MessageNotify('M','w_base_mpi_Finalize',& 'w_base_mpi_module (2013/02/23) is finalized') end subroutine w_base_mpi_Finalize end module w_base_mpi_module