!-- !---------------------------------------------------------------------- ! Copyright (c) 2016 SPMDOEL Development Group !---------------------------------------------------------------------- !表題 ua_base_mpi_module_svpack ! ! spml/ua_base_mpi_module_svpack モジュールは球面上での流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール ua_mpi_module_svpack の下部モジュールであり, ! スペクトル計算の基本的な Fortran90 関数を提供する. ! ! 球面上の 1 層モデル用 w_base_module_svpack モジュールを多層モデル用に ! 並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに ! 対する変換が行える. ! ! 内部で ISPACK2 の SVPACK/SUPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK2 の SVPACK/SUPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_mpi_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! ! !履歴 2016/08/15 竹広真一 wa_base_module_svpack より多層 MPI 対応改造 ! !++ module ua_base_mpi_module_svpack ! != ua_base_mpi_module_svpack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/ua_base_mpi_module_svpack モジュールは球面上での流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール ua_mpi_module_svpack の下部モジュールであり, ! スペクトル計算の基本的な Fortran90 関数を提供する. ! ! 球面上の 1 層モデル用 w_base_module_svpack モジュールを多層モデル用に ! 並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに ! 対する変換が行える. ! ! 内部で ISPACK2 の SVPACK/SUPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK2 の SVPACK/SUPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_mpi_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! use dc_message use w_base_module_svpack, only : im, jm, nn, mm, xy_w, w_xy, & w_StreamPotential2Vector, & w_Vector2VorDiv, w_VectorCosLat2VorDiv use w_base_mpi_module_supack, only : w_base_mpi_initial, nc, l_nm, nm_l use mpi implicit none integer :: km=16 ! 同時に処理する最大データ数(層の数). integer :: kc integer :: k1, k2 integer :: np integer :: ip private public km ! 層数 public ua_base_mpi_Initial ! 初期化サブルーチン public ua_base_mpi_Finalize ! 終了処理 public nc public kc public k1, k2 public l_nm, nm_l public kp_kl public xyb_ua, ua_xyb ! 変換関数 public wb_ua, ua_wb ! 変換関数 public ua_StreamPotential2VectorMPI ! 流線ポテンシャルから速度場計算 public ua_Vector2VorDivMPI ! 速度場から渦度発散を計算 public ua_VectorCosLat2VorDivMPI ! 速度場から渦度発散を計算 save km ! 最大データ数(層数)を記憶 save np, ip save kc, k1, k2 contains !--------------- 初期化 ----------------- subroutine ua_base_mpi_initial(n_in,i_in,j_in,k_in) ! ! ua_base_mpi_module 初期化 ! integer,intent(in) :: n_in !(in) 切断波数 integer,intent(in) :: i_in !(in) 格子点数(東西) integer,intent(in) :: j_in !(in) 格子点数(南北) integer,intent(in) :: k_in !(in) 最大データ数(層数) integer :: ierr km = k_in CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR) call set_layerdiv(ip,k1,k2) kc = k2-k1+1 ! 分散多層データの大きさ call w_base_mpi_initial(n_in,i_in,j_in) call MessageNotify('M','ua_base_mpi_initial', & 'ua_base_mpi_module_svpack (2016/08/15) is initialized') end subroutine ua_base_mpi_initial subroutine set_layerdiv(ip,kk1,kk2) integer, intent(IN) :: ip ! プロセス番号 integer, intent(OUT) :: kk1 ! 担当する層開始番号 integer, intent(OUT) :: kk2 ! 担当する層終了番号 integer :: ki ki=(km-1)/np+1 if(ki*ip.lt.km) then kk1=1+ki*ip kk2=min(km,ki*(ip+1)) else ! 当該プロセスでは層データを担当しない場合 kk1=0 kk2=-1 end if end subroutine set_layerdiv subroutine layer2procinfo(k,ipk,kp) integer, intent(IN) :: k ! 層番号 integer, intent(OUT) :: ipk ! 担当するプロセス番号 integer, intent(OUT) :: kp ! プロセス内での分割後層番号 integer :: iip, kk1, kk2 do iip=0,np-1 call set_layerdiv(iip,kk1,kk2) if ( kk1 .le. k .AND. k .le. kk2 ) then ipk = iip kp = k - kk1 + 1 goto 999 end if end do call MessageNotify('E','layer2porcinfo','Process not found.') 999 continue end subroutine layer2procinfo !--------------- 基本関数 ----------------- function kp_kl(kl) integer, intent(IN) :: kl ! 層番号(全体) integer :: kp_kl ! 層番号(分割プロセス中) integer :: ipk call layer2procinfo(kl,ipk,kp_kl) if ( ip /= ipk ) kp_kl = 0 end function kp_kl !--------------- データ入れ換え ----------------- function wb_ua(ua_data) real(8), intent(IN) :: ua_data(nc,km) real(8) :: wb_ua((mm+1)**2,kc) real(8) :: databuf((MM/NP+1)*(2*(NN+1)-MM/NP*NP)*NP*KM) integer :: nb, kk1, kk2, kp, m, ipk, k, kk, ns, ipdest, l, la, ierr NB=(MM/NP+1)*(2*(NN+1)-MM/NP*NP) do ipk=0,np-1 call set_layerdiv(ipk,kk1,kk2) kp = kk2-kk1+1 CALL MPI_GATHER(ua_data(:,kk1),NB*kp,MPI_REAL8, & databuf,NB*kp,MPI_REAL8,& IPK,MPI_COMM_WORLD,IERR) enddo !$omp parallel do default(shared) private(m,kk,ns,ipdest,l,la) do k=1,kc CALL SUCOPY(nn+1,databuf(nb*(k-1)+1),wb_ua(1,k)) DO M=1,MM KK=M/NP NS=KK*(2*(NN+1)-(KK-1)*NP)+1 CALL SUNM2L(NN,M,M,IPDEST,L) CALL SVNM2L(NN,M,M,LA) CALL SUCOPY(2*(NN+1-M),databuf(ns+nb*kc*ipdest+nb*(k-1)),wb_ua(la,k)) enddo END DO !$omp end parallel do end function wb_ua function ua_wb(wb_data) real(8), intent(IN) :: wb_data((mm+1)**2,kc) real(8) :: ua_wb(nc,km) real(8) :: databuf( (mm/np+1)*(2*(nn+1)-mm/np*np)*np*km ) integer :: ipk, kk1, kk2, kp, k, m, kk, NS, NB, IPDEST, L, LA, IERR NB=(MM/NP+1)*(2*(NN+1)-MM/NP*NP) !$omp parallel do default(shared) private(m,kk,ns,ipdest,l,la) do k=1,kc CALL SUCOPY(nn+1,wb_data(1,k),databuf(nb*(k-1)+1)) DO M=1,MM KK=M/NP NS=KK*(2*(NN+1)-(KK-1)*NP)+1 CALL SUNM2L(NN,M,M,IPDEST,L) CALL SVNM2L(NN,M,M,LA) CALL SUCOPY(2*(NN+1-M),wb_data(la,k),databuf(ns+nb*kc*ipdest+nb*(k-1))) enddo END DO !$omp end parallel do do ipk=0,np-1 call set_layerdiv(ipk,kk1,kk2) kp = kk2-kk1+1 if ( kp == 0 ) then CALL MPI_SCATTER(databuf,NB*KP,MPI_REAL8,& ua_wb(1,1),NB*KP,MPI_REAL8, & ipk,MPI_COMM_WORLD,IERR) else CALL MPI_SCATTER(databuf,NB*KP,MPI_REAL8,& ua_wb(1,kk1),NB*KP,MPI_REAL8, & ipk,MPI_COMM_WORLD,IERR) endif enddo end function ua_wb !--------------- 基本変換 ----------------- function xyb_ua(ua_data,ipow,iflag) ! 球面調和関数スペクトル -> 格子点 ! ! スペクトルデータから格子データへ変換する(多層用). ! real(8), intent(in) :: ua_data(nc,km) !(in) スペクトルデータ ! real(8) :: xyb_ua(0:im-1,1:jm,kc) !(out) 格子点データ ! 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) :: wb_data((mm+1)**2,kc) integer :: k if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif wb_data = wb_ua(ua_data) do k=1,kc xyb_ua(:,:,k) = xy_w(wb_data(:,k),iflag=ifval,ipow=ipval) enddo end function xyb_ua function ua_xyb(xyb_data,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル ! ! 格子データからスペクトルデータへ(正)変換する(多層用). ! real(8), intent(in) :: xyb_data(0:im-1,jm,kc) !(in) 格子点データ real(8) :: ua_xyb(nc,km) !(out) スペクトルデータ 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) :: wb_data((mm+1)**2,kc) integer k if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif do k=1,kc wb_data(:,k) = w_xy(xyb_data(:,:,k),iflag=ifval,ipow=ipval) enddo ua_xyb = ua_wb(wb_data) end function ua_xyb !----------- 速度, 渦度・発散, 流線・ポテンシャル計算 ------------- subroutine ua_StreamPotential2VectorMPI(ua_Psi, ua_Chi, xyb_U, xyb_V) ! ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に ! (逆)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ, ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ ! real(8), intent(in) :: ua_Psi(nc,km) !(in) 流線関数((nm+1)*(nm+1),:) real(8), intent(in) :: ua_Chi(nc,km) !(in) 速度ポテンシャル((nm+1)*(nm+1),:) real(8), intent(out) :: xyb_U(0:im-1,jm,kc) !(out) 速度経度成分(0:im-1,1:jm,:) real(8), intent(out) :: xyb_V(0:im-1,jm,kc) !(out) 速度緯度成分(0:im-1,1:jm,:) real(8) :: wb_Psi((mm+1)**2,kc) !(in) 流線関数((nm+1)*(nm+1),:) real(8) :: wb_Chi((mm+1)**2,kc) !(in) 速度ポテンシャル((nm+1)*(nm+1),:) integer :: k wb_Psi = wb_ua(ua_Psi) wb_Chi = wb_ua(ua_Chi) do k=1,kc call w_StreamPotential2Vector & ( wb_Psi(:,k), wb_Chi(:,k), xyb_U(:,:,k), xyb_V(:,:,k) ) enddo end subroutine ua_StreamPotential2VectorMPI subroutine ua_Vector2VorDivMPI(xyb_U, xyb_V, ua_Vor, ua_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xyb_U(0:im-1,jm,kc) !(in) 速度経度成分(0:im-1,1:jm,:) real(8), intent(in) :: xyb_V(0:im-1,jm,kc) !(in) 速度緯度成分(0:im-1,1:jm,:) real(8), intent(out) :: ua_Vor(nc,km) !(out) 渦度 real(8), intent(out) :: ua_Div(nc,km) !(out) 発散 real(8) :: wb_Vor((mm+1)**2,kc) ! 渦度 real(8) :: wb_Div((mm+1)**2,kc) ! 発散 integer :: k do k=1,kc call w_Vector2VorDiv & (xyb_U(:,:,k), xyb_V(:,:,k), wb_Vor(:,k), wb_Div(:,k)) enddo ua_Vor = ua_wb(wb_Vor) ua_Div = ua_wb(wb_Div) end subroutine ua_Vector2VorDivMPI subroutine ua_VectorCosLat2VorDivMPI(xyb_UCosLat, xyb_VCosLat, ua_Vor, ua_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xyb_UCosLat(0:im-1,jm,kc) !(in) 速度経度成分 * cos(lat) (0:im-1,1:jm,:) real(8), intent(in) :: xyb_VCosLat(0:im-1,jm,kc) !(in) 速度緯度成分 * cos(lat) (0:im-1,1:jm,:) real(8), intent(out) :: ua_Vor(nc,km) !(out) 流線関数 real(8), intent(out) :: ua_Div(nc,km) !(out) 速度ポテンシャル real(8) :: wb_Vor((mm+1)**2,kc) ! 渦度 real(8) :: wb_Div((mm+1)**2,kc) ! 発散 integer :: k do k=1,kc call w_VectorCosLat2VorDiv & (xyb_UCosLat(:,:,k), xyb_VCosLat(:,:,k), wb_Vor(:,k), wb_Div(:,k)) enddo ua_Vor = ua_wb(wb_Vor) ua_Div = ua_wb(wb_Div) end subroutine ua_VectorCosLat2VorDivMPI !--------------- 終了処理 ----------------- subroutine ua_base_mpi_Finalize ! ! ! このサブルーチンを単独で用いるのでなく, ! 上位サブルーチン wa_Finalize を使用すること. ! call MessageNotify('M','ua_base_mpi_finalize',& 'ua_base_mpi_module_svpack is finalized (dummy, 2016/08/15) ') end subroutine ua_base_mpi_Finalize end module ua_base_mpi_module_svpack