!---------------------------------------------------------------------- ! Copyright (c) 2019 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 w_base_mpi_module_sypack ! ! spml3/w_base_mpi_module_sypack モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法と MPI によって数値計算するための ! モジュール w_mpi_module_sypack の下部モジュールであり, スペクトル計算の ! 基本的な Fortran90 関数を提供する. ! ! 内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK3/SYPACK のマニュアルを参照されたい. ! ! !履歴 2019/03/25 竹広真一 ! module w_base_mpi_module_sypack ! ! w_base_mpi_module_sypack ! ! spml3/w_base_mpi_module_sypack モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法と MPI によって数値計算するための ! モジュール w_mpi_module_sypack の下部モジュールであり, スペクトル法の ! 基本的なな Fortran90 関数を提供する. ! ! 内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK3/SYPACK のマニュアルを参照されたい. ! use dc_message use iso_c_binding use mpi implicit none integer(8) :: im=64 ! 格子点の設定(東西) integer(8) :: jm=32 ! 格子点の設定(南北) integer(8) :: mm ! 切断波数(東西波数)の設定 integer(8) :: nn ! 切断波数(全波数)の設定 integer(8) :: nm ! 計算する最大の全波数の設定 integer(8) :: ig=1 ! グリッド種類 integer(8) :: icom ! コミュニケータ integer(8), allocatable :: it(:) ! 変換用配列 integer(8), allocatable :: jc(:) ! 変換用配列 real(8), allocatable :: t(:) ! 変換用配列 real(8), allocatable :: r(:) ! 変換用配列 real(8), pointer :: p(:,:) ! 変換用配列 real(8), pointer :: g(:,:) ! 変換用配列 real(8), pointer :: w(:) ! 変換用配列 real(8), pointer :: g1(:,:),g2(:,:) ! 変換用配列 real(8), pointer :: ww(:) ! 変換用配列 type(c_ptr) :: pp ! 変換用配列 type(c_ptr) :: pg ! 変換用配列 type(c_ptr) :: pw ! 変換用配列 type(c_ptr) :: pww,pg1,pg2 ! 変換用配列 real(8), allocatable :: c(:) ! 作業配列 integer(8) :: j1, j2 ! 分散配置用変数 integer :: jl ! 分散格子データの大きさ integer :: nc ! 分散スペクトルデータの大きさ integer, allocatable :: j1p(:), j2p(:) ! 分散配置用変数 integer, allocatable :: jlp(:) ! 分散格子データの大きさ integer, allocatable :: ncp(:) ! 分散スペクトルデータの大きさ integer, allocatable :: snd_rank(:) ! 分散/全データ変換用配列 integer, allocatable :: rcv_rank(:) ! 分散/全データ変換用配列 integer(8) :: jv ! ベクトル長 integer :: np ! MPI プロセス数 integer :: ip ! MPI プロセス番号 real(8), allocatable :: x_Lon(:), x_Lon_Weight(:) ! 経度 real(8), allocatable :: v_Lat(:), v_Lat_Weight(:) ! 緯度 real(8), allocatable :: y_Lat(:), y_Lat_Weight(:) ! 緯度 real(8), allocatable :: xv_Lon(:,:), xv_Lat(:,:) real(8), allocatable :: xy_Lon(:,:), xy_Lat(:,:) logical :: w_base_mpi_initialize=.false. ! 初期化フラッグ real(8), parameter :: pi=3.1415926535897932385D0 private public im, jm, nn, mm, nm ! 格子点数, 切断波数 public it, t, p, r, jc ! 変換用作業配列 public j1, j2 ! 分散格子データの大きさ public jl ! 分散スペクトルデータの大きさ public nc public np, ip public icom ! コミュニケーター public w_base_mpi_Initial ! 初期化サブルーチン public w_base_mpi_Finalize ! 終了サブルーチン public x_Lon ! 格子座標 public x_Lon_Weight ! 格子座標重み public v_Lat ! 格子座標 public v_Lat_Weight ! 格子座標重み public y_Lat ! 格子座標 public y_Lat_Weight ! 格子座標重み public xv_Lon, xv_Lat ! 格子座標(im,jl) public xy_Lon, xy_Lat ! 格子座標(im,jm) public l_nm, nm_l ! 波数格納位置 public xv_w, w_xv ! 変換関数 public w_StreamPotential2VectorMPI ! 流線ポテンシャルから速度場計算 public w_Vector2VorDivMPI ! 速度場から渦度発散を計算 public w_VectorCosLat2VorDivMPI ! 速度場から渦度発散を計算 public ay_av, av_ay 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 ! 格子点数, 切断波数を記憶 save it, t, p, r, jc ! 変換用配列を記憶 save icom save jv save j1, j2 save jl, nc save j1p, j2p save jlp, ncp save np, ip save snd_rank, rcv_rank save c ! 変換用配列 save w, g, pw, pg save ww, g1, g2, pww, pg1, pg2 save w_base_mpi_initialize ! 初期化フラグ contains !--------------- 初期化 ----------------- subroutine w_base_mpi_Initial(n_in,i_in,j_in,mpi_comm) ! ! スペクトル変換の格子点数, 波数および 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 :: mpi_comm !(in) コミュニケーター integer :: i, j integer :: ierr ! 三角波数切断 im = i_in ; jm = j_in nn = n_in ; nm = n_in+1 ; mm = n_in ! default は三角波数切断 if ( present(mpi_comm) ) then icom = mpi_comm else icom = MPI_COMM_WORLD endif CALL MPI_COMM_SIZE(int(ICOM),NP,IERR) CALL MPI_COMM_RANK(int(ICOM),IP,IERR) allocate(t(im*3/2)) allocate(it(im/2)) allocate(r(5*(mm/np+1)*(2*nm-mm/np*np)/4+mm/np+1)) allocate(jc((mm/np+1)*(2*nm-mm/np*np)/16+mm/np+1)) call mxallc(pp,jm/2*(5+2*(mm/np+1))) call c_f_pointer(pp,p,[jm/2,5+2*(mm/np+1)]) allocate(c((mm/np+1)*(2*(nn+1)-mm/np*np))) ! 変換用作業配列 call syini1(mm,nm,im,it,t,r,icom) call syini2(mm,nm,jm,ig,p,r,jc,icom) call syinic(mm,nn,c,icom) ! 変換用変数割付 call syqrjv(jm,jv) call mxallc(pg,im*((jm/jv-1)/np+1)*jv) call c_f_pointer(pg, g, [im,((jm/jv-1)/np+1)*jv]) call mxallc(pw,2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2) call c_f_pointer(pw, w, [2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2]) call mxallc(pg1,im*((jm/jv-1)/np+1)*jv) call c_f_pointer(pg1, g1, [im,((jm/jv-1)/np+1)*jv]) call mxallc(pg2,im*((jm/jv-1)/np+1)*jv) call c_f_pointer(pg2, g2, [im,((jm/jv-1)/np+1)*jv]) call mxallc(pww,2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2*2) call c_f_pointer(pww, ww, [2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2*2]) ! 各プロセスの担当緯度, 波数情報 call syqrnj(jm,jv,j1,j2,icom) jl = int(j2-j1+1) ! 分散格子データの大きさ nc = int((mm/np+1)*(2*(nn+1)-mm/np*np)) ! 分散スペクトルデータの大きさ allocate(j1p(0:np-1),j2p(0:np-1)) allocate(jlp(0:np-1)) allocate(ncp(0:np-1)) call MPI_ALLGATHER( jl, 1, MPI_INTEGER, & jlp(0), 1, MPI_INTEGER, int(icom), ierr) call MPI_ALLGATHER( j1, 1, MPI_INTEGER, & j1p(0), 1, MPI_INTEGER, int(icom), ierr) call MPI_ALLGATHER( j2, 1, MPI_INTEGER, & j2p(0), 1, MPI_INTEGER, int(icom), ierr) call MPI_ALLGATHER( nc, 1, MPI_INTEGER, & ncp(0), 1, MPI_INTEGER, int(icom), ierr) allocate(snd_rank(1:np-1)) allocate(rcv_rank(1:np-1)) do i=1,np-1 rcv_rank(i) = mod(ip+i,np) snd_rank(i) = mod(ip+np-i,np) enddo ! 座標変数設定 allocate(x_Lon(0:im-1)) ! 格子点座標格納配列(経度) allocate(x_Lon_Weight(0:im-1)) allocate(v_Lat(jl),v_Lat_Weight(jl)) ! 格子点座標格納配列 allocate(y_Lat(jm),y_Lat_Weight(jm)) ! 格子点座標格納配列 allocate(xv_Lon(0:im-1,jl),xv_Lat(0:im-1,jl)) ! 格子点座標格納配列 allocate(xy_Lon(0:im-1,jm),xy_Lat(0:im-1,jm)) ! 格子点座標格納配列 do i=0,int(im)-1 x_Lon(i) = 2*pi/im*i ! 経度座標 x_Lon_Weight(i) = 2*pi/im ! 経度座標重み enddo do j=int(j1),int(j2) if ( j .gt. jm/2 ) then v_Lat(j-j1+1) = asin(p(j-jm/2,1)) ! 緯度座標 v_Lat_Weight(j-j1+1) = 2*p(j-jm/2,2) ! 緯度重み(Gauss grid) else v_Lat(j-j1+1) = -asin(p(jm/2-j+1,1)) ! 緯度座標 v_Lat_Weight(j-j1+1) = 2*p(jm/2-j+1,2) ! 緯度重み(Gauss grid) endif enddo do i=0,int(im)-1 xv_Lat(i,:) = v_Lat enddo xy_Lat = ay_av(xv_Lat) ; y_Lat = xy_Lat(0,:) y_Lat_Weight = reshape(ay_av(reshape(v_Lat_Weight, (/1,jl/))),(/jm/)) do j=1,int(jl) xv_Lon(:,j) = x_Lon enddo do j=1,int(jm) xy_Lon(:,j) = x_Lon enddo w_base_mpi_initialize = .true. call MessageNotify('M','w_base_mpi_initial',& 'w_base_mpi_module_sypack (2019/03/25) is initialized') end subroutine w_base_mpi_Initial !--------------- 基本関数 ----------------- function l_nm_array00(n_in,m_in) ! ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. ! ! 引数 n,m がともに整数値の場合, 整数値を返す. ! プロセスに n,m が存在しない場合 0 を返す ! integer :: l_nm_array00 !(out) スペクトルデータの格納位置 integer, intent(in) :: n_in !(in) 全波数 integer, intent(in) :: m_in !(in) 帯状波数 integer(8) :: ip_w integer(8) :: l, n, m if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','l_nm_array00',& 'w_base_mpi_module_sypack not initialize yet. Use synm2l routine in ISPACK3 directly.') else n = n_in ; m = m_in call synm2l(nn,n,m,ip_w,l,icom) if ( ip_w /= ip ) then l_nm_array00 = 0 else l_nm_array00 = int(l) endif 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) ! ! 全波数(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(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) スペクトルデータの格納位置 integer(8) :: n,m ! 全波数, 帯状波数 if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','nm_l_int',& 'w_base_mpi_module_sypack not initialize yet. Use syl2nm routine in ISPACK3 directly.') else call syl2nm(mm,nn,int(l,kind=8),n,m,icom) nm_l_int = (/int(n),int(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 xv_w(w_data,ipow,iflag) ! ! スペクトルデータから格子データへ変換する(1 層用). ! real(8) :: xv_w(0:im-1,1:jl) !(out) 格子点データ real(8), intent(in) :: w_data(nc) !(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(8) :: ipval, ifval real(8) :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np)) ! 作業用スペクトルデータ(SYCS2X 出力用) real(8) :: w_Ydata((mm/np+1)*(2*(nm+1)-mm/np*np)) ! 作業用スペクトルデータ(SYCS2Y 出力用) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','xv_w',& 'w_base_mpi_module_sypack 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 syts2g(mm,nm,nn,im,jm,jv,w_data,g,it,t,p,r,jc,w,ipval,icom) xv_w = g(:,1:jl) else if( ifval==-1 ) then call sycs2x(mm,nn,w_data,w_Xdata,icom) call syts2g(mm,nm,nn,im,jm,jv,w_Xdata,g,it,t,p,r,jc,w,ipval,icom) xv_w = g(:,1:jl) else if( ifval==1 ) then call sycs2y(mm,nn,w_data,w_Ydata,c,icom) call syts2g(mm,nm,nm,im,jm,jv,w_Ydata,g,it,t,p,r,jc,w,ipval,icom) xv_w = g(:,1:jl) else if( ifval==2 ) then call syts2g(mm,nm,nn,im,jm,jv,w_data,g,it,t,p,r,jc,w,ipval,icom) xv_w = g(:,1:jl)*sin(xv_Lat) else call MessageNotify('E','xv_w','invalid value of iflag') endif end function xv_w function w_xv(xv_data,ipow,iflag) ! ! 格子データからスペクトルデータへ(正)変換する(1 層用). ! real(8) :: w_xv(nc) !(out) スペクトルデータ real(8), intent(in) :: xv_data(0:im-1,1:jl) !(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(8) :: ipval, ifval real(8) :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np)) ! 作業用スペクトルデータ(SYCS2X 出力用) real(8) :: w_Ydata((mm/np+1)*(2*(nm+1)-mm/np*np)) ! 作業用スペクトルデータ(SYCS2Y 出力用) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','xv_w',& 'w_base_mpi_module_sypack 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 g = 0.0D0 ; g(:,1:jl) = xv_data call sytg2s(mm,nm,nn,im,jm,jv,w_xv,g,it,t,p,r,jc,w,ipval,icom) else if ( ifval == -1 ) then g = 0.0D0 ; g(:,1:jl) = xv_data call sytg2s(mm,nm,nn,im,jm,jv,w_Xdata,g,it,t,p,r,jc,w,ipval,icom) call sycs2x(mm,nn,w_Xdata,w_xv,icom) else if ( ifval == 1 ) then g = 0.0D0 ; g(:,1:jl) = xv_data call sytg2s(mm,nm,nm,im,jm,jv,w_Ydata,g,it,t,p,r,jc,w,ipval,icom) call sycy2s(mm,nn,w_Ydata,w_xv,c,icom) else if ( ifval == 2 ) then g = 0.0D0 ; g(:,1:jl) = xv_data*sin(xv_Lat) call sytg2s(mm,nm,nn,im,jm,jv,w_xv,g,it,t,p,r,jc,w,ipval,icom) end if 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(nc) !(in) 流線関数 real(8), intent(in) :: w_Chi(nc) !(in) 速度ポテンシャル real(8), intent(out) :: xv_U(0:im-1,1:jl) !(out) 速度経度成分 real(8), intent(out) :: xv_V(0:im-1,1:jl) !(out) 速度緯度成分 real(8) :: w_Rdata1((mm/np+1)*(2*(nm+1)-mm/np*np)) real(8) :: w_Rdata2((mm/np+1)*(2*(nm+1)-mm/np*np)) ! 作業用スペクトルデータ(SYTS2G 出力用) real(8) :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np)) ! 作業用スペクトルデータ(SYCS2X 出力用) real(8) :: w_Ydata((mm/np+1)*(2*(nm+1)-mm/np*np)) ! 作業用スペクトルデータ(SYCS2Y 出力用) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','w_StreamPotential2VectorMPI',& 'w_base_module not initialize yet.') endif ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ の計算 ! call sycs2x(mm,nn,w_Chi,w_Xdata,icom) call sycs2y(mm,nn,w_Psi,w_Ydata,c,icom) call sycrpk(mm,nn,nm,w_Xdata,w_Rdata1,icom) w_Rdata1 = w_Rdata1 - w_Ydata ! ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ の計算 ! call sycs2y(mm,nn,w_Chi,w_Ydata,c,icom) call sycs2x(mm,nn,w_Psi,w_Xdata,icom) call sycrpk(mm,nn,nm,w_Xdata,w_Rdata2,icom) w_Rdata2= w_Rdata2 + w_Ydata ! call syts2v(mm,nm,nm,im,jm,jv,w_Rdata1,w_Rdata2,g1,g2,it,t,p,r,jc,ww,1_8,icom) ! ! u の計算 ! xv_U = g1(:,1:jl) ! ! v の計算 ! xv_V = g2(:,1:jl) 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:jl) !(in) 速度経度成分 real(8), intent(in) :: xv_V(0:im-1,1:jl) !(in) 速度緯度成分 real(8), intent(out) :: w_Vor(nc) !(out) 流線関数 real(8), intent(out) :: w_Div(nc) !(out) 速度ポテンシャル real(8) :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np)) ! 作業用スペクトルデータ(SUCS2X 出力用) real(8) :: w_Ydata1((mm/np+1)*(2*(nm+1)-mm/np*np)) real(8) :: w_Ydata2((mm/np+1)*(2*(nm+1)-mm/np*np)) ! 作業用スペクトルデータ(SUCS2Y 出力用) real(8) :: w_Data1((mm/np+1)*(2*(nn+1)-mm/np*np)) real(8) :: w_Data2((mm/np+1)*(2*(nn+1)-mm/np*np)) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','w_Vector2VorDivMPI',& 'w_base_module not initialize yet.') endif ! ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ の計算 ! g1 = 0.0D0 ; g1(:,1:jl) = xv_U g2 = 0.0D0 ; g2(:,1:jl) = xv_V call sytv2s(mm,nm,nn+1,im,jm,jv,w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,1_8,icom) call sycrpk(mm,nm,nn,w_Ydata1,w_Xdata,icom) call sycs2x(mm,nn,w_Xdata,w_Div,icom) call sycy2s(mm,nn,w_Ydata2,w_Data2,c,icom) ! ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算 ! call sycrpk(mm,nm,nn,w_Ydata2,w_Xdata,icom) call sycs2x(mm,nn,w_Xdata,w_Vor,icom) call sycy2s(mm,nn,w_Ydata1,w_Data1,c,icom) ! ! 渦度・発散の計算 ! ! ζ = 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_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:jl) !(in) 速度経度成分 * cos(lat) real(8), intent(in) :: xv_VCosLat(0:im-1,1:jl) !(in) 速度緯度成分 * cos(lat) real(8), intent(out) :: w_Vor(nc) !(out) 流線関数 real(8), intent(out) :: w_Div(nc) !(out) 速度ポテンシャル real(8) :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np)) ! 作業用スペクトルデータ(SUCS2X 出力用) real(8) :: w_Ydata1((mm/np+1)*(2*(nm+1)-mm/np*np)) real(8) :: w_Ydata2((mm/np+1)*(2*(nm+1)-mm/np*np)) ! 作業用スペクトルデータ(SUCS2Y 出力用) real(8) :: w_Data1((mm/np+1)*(2*(nn+1)-mm/np*np)) real(8) :: w_Data2((mm/np+1)*(2*(nn+1)-mm/np*np)) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','w_VectorCosLat2VorDivMPI',& 'w_base_module not initialize yet.') endif ! ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ の計算 ! g1 = 0.0D0 ; g1(:,1:jl) = xv_UCosLat g2 = 0.0D0 ; g2(:,1:jl) = xv_VCosLat call sytv2s(mm,nm,nn+1,im,jm,jv,w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,2_8,icom) call sycrpk(mm,nm,nn,w_Ydata1,w_Xdata,icom) call sycs2x(mm,nn,w_Xdata,w_Div,icom) call sycy2s(mm,nn,w_Ydata2,w_Data2,c,icom) ! ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算 ! call sycrpk(mm,nm,nn,w_Ydata2,w_Xdata,icom) call sycs2x(mm,nn,w_Xdata,w_Vor,icom) call sycy2s(mm,nn,w_Ydata1,w_Data1,c,icom) ! ! 渦度・発散の計算 ! ! ζ = 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_VectorCosLat2VorDivMPI !----------- y 方向結合・分割 ------------- function ay_av(av) ! ! 分割緯度データを集める. ! real(8), dimension(:,:), intent(in) :: av !(in) 緯度格子点データ(分割) (:,jl) real(8), dimension(size(av,1),jm) :: ay_av !(out) 緯度格子点データ integer :: a_size integer :: IRQSTS integer :: IST( MPI_STATUS_SIZE ) integer :: ierr, i integer, parameter :: itag=103 a_size = size(av,1) call MPI_Barrier(int(icom),IERR) do i=1,np-1 call MPI_ISEND(av, a_size*jl, MPI_REAL8, rcv_rank(i),& itag,int(icom),IRQSTS,IERR) if ( j1p(snd_rank(i)) .ge. 1 ) then call MPI_RECV(ay_av(1,j1p(snd_rank(i))),& a_size*jlp(snd_rank(i)), MPI_REAL8, & snd_rank(i), itag, int(icom), IST, IERR) endif CALL MPI_WAIT(IRQSTS,IST,IERR) enddo ay_av(:,j1:j2) = av end function ay_av function av_ay(ay) ! ! 緯度データを分割する ! real(8), dimension(:,:), intent(in) :: ay !(in) 緯度格子点データ(分割) (:,jc) real(8), dimension(size(ay,1),jl) :: av_ay !(out) 緯度格子点データ av_ay = ay(:,j1:j2) end function av_ay !--------------- 終了処理 ----------------- subroutine w_base_mpi_Finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! 実際の使用には上位サブルーチン w_Finalize を用いること. ! if ( .not. w_base_mpi_initialize ) then call MessageNotify('W','w_base_mpi_Finalize',& 'w_base_mpi_module_sypack not initialized yet') return endif deallocate(it) ! 変換用配列 deallocate(t) ! 変換用配列 deallocate(r) ! 変換用配列 deallocate(jc) ! 変換用配列 deallocate(c) ! 変換用作業配列 call mxfree(pp) ! 変換用作業配列 call mxfree(pg) ! 変換用作業配列 call mxfree(pg1) ! 変換用作業配列 call mxfree(pg2) ! 変換用作業配列 call mxfree(pw) ! 変換用作業配列 call mxfree(pww) ! 変換用作業配列 deallocate(j1p,j2p) deallocate(jlp) deallocate(ncp) deallocate(snd_rank) deallocate(rcv_rank) deallocate(x_Lon) ! 格子点座標格納配列(経度) deallocate(x_Lon_weight) deallocate(xv_Lon) deallocate(v_Lat) deallocate(v_Lat_Weight) ! 格子点座標格納配列 deallocate(xv_Lat) ! 格子点座標格納配列 deallocate(y_Lat,y_Lat_Weight) ! 格子点座標格納配列 deallocate(xy_Lon,xy_Lat) ! 格子点座標格納配列 w_base_mpi_initialize = .false. call MessageNotify('M','w_base_mpi_Finalize',& 'w_base_mpi_module_sypack (2019/03/25) is finalized') end subroutine w_base_mpi_Finalize end module w_base_mpi_module_sypack