! -*- mode: f90; coding: utf-8 -*- !------------------------------------------------------------------------- ! Copyright (c) 2019-2020 SPMODEL Development Group. All rights reserved. ! ! 履歴 2019/03/25 竹広真一 ! 2019/10/05 佐々木洋平 ! 2020/08/12 竹広真一 ! 2020/11/07 竹広真一 セクター計算オプション導入 ! !------------------------------------------------------------------------- !> @copyright Copyright (C) SPMODEL Development Group, 2019. !> License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail !> @author Shin-ichi Takehiro, Youhei SASAKI !> @brief !> @ref w_mpi_module の下部モジュール: 基本変換 !> !> @warning !> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい． !> ユーザは w_mpi_module を呼ぶこと． !> module w_mpi_module_base_mint use dc_types, only: DP use dc_message, only: MessageNotify use spml_utils, only: pi use iso_c_binding, only: c_ptr, c_f_pointer use mpi implicit none private public im, jm, nn, mm, nm, mi ! 格子点数, 切断波数 public it, t, p, r, jc ! 変換用作業配列 public j1, j2 ! 分散格子データの大きさ public jl ! 分散スペクトルデータの大きさ public nc public np, ip public c 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 !> 格子点の設定(東西) integer(8) :: im=64 !> 格子点の設定(南北) integer(8) :: jm=32 !> 切断波数(東西波数)の設定 integer(8) :: mm !> 切断波数(全波数)の設定 integer(8) :: nn !> 計算する最大の全波数の設定 integer(8) :: nm !> 経度方向対称性 integer(8) :: mi !> 経度座標格子点 real(DP), allocatable :: x_Lon(:) !> 経度座標格子点重み real(DP), allocatable :: x_Lon_Weight(:) !> 緯度座標分散格子点 real(DP), allocatable :: v_Lat(:) !> 緯度座標分散格子点重み real(DP), allocatable :: v_Lat_Weight(:) !> 緯度座標格子点 real(DP), allocatable :: y_Lat(:) !> 緯度座標格子点重み real(DP), allocatable :: y_Lat_Weight(:) !> 分散座標経度格子点, 2次元配列 real(DP), allocatable :: xv_Lon(:,:) !> 分散座標緯度格子点, 2次元配列 real(DP), allocatable :: xv_Lat(:,:) !> 座標経度格子点, 2次元配列 real(DP), allocatable :: xy_Lon(:,:) !> 座標緯度格子点, 2次元配列 real(DP), allocatable :: xy_Lat(:,:) !> コミュニケータ integer(8) :: icom !> 変換用配列 integer(8), allocatable :: it(:) !> 変換用配列 integer(8), allocatable :: jc(:) !> 変換用配列 real(DP), allocatable :: t(:) !> 変換用配列 real(DP), allocatable :: r(:) !> 変換用配列 real(DP), pointer :: p(:,:) !> 分散配置用変数 integer(8) :: j1 !> 分散配置用変数 integer(8) :: j2 !> 分散格子データの大きさ integer :: jl !> 分散スペクトルデータの大きさ integer :: nc !> MPI プロセス数 integer :: np !> MPI プロセス番号 integer :: ip !> グリッド種類 integer(8) :: ig=1 !> 変換用配列 real(DP), pointer :: g(:,:) !> 変換用配列 real(DP), pointer :: w(:) !> 変換用配列 real(DP), pointer :: g1(:,:) !> 変換用配列 real(DP), pointer :: g2(:,:) !> 変換用配列 real(DP), pointer :: ww(:) !> 変換用配列 type(c_ptr) :: pp !> 変換用配列 type(c_ptr) :: pg !> 変換用配列 type(c_ptr) :: pw !> 変換用配列 type(c_ptr) :: pww !> 変換用配列 type(c_ptr) :: pg1 !> 変換用配列 type(c_ptr) :: pg2 !> 作業配列 real(DP), allocatable :: c(:) !> 分散配置用変数 integer, allocatable :: j1p(:) !> 分散配置用変数 integer, allocatable :: j2p(:) !> 分散格子データの大きさ integer, allocatable :: jlp(:) !> 分散スペクトルデータの大きさ integer, allocatable :: ncp(:) !> 分散/全データ変換用配列 integer, allocatable :: snd_rank(:) !> 分散/全データ変換用配列 integer, allocatable :: rcv_rank(:) !> ベクトル長 integer(8) :: jv !> 初期化フラグ logical :: w_base_mpi_initialize=.false. !--------------------------------------------------------------------- !> 波数 (n,m) よりスペクトルデータの格納位置 l を返す !> interface l_nm module procedure l_nm_array00 module procedure l_nm_darray00 module procedure l_nm_array01 module procedure l_nm_darray01 module procedure l_nm_array10 module procedure l_nm_darray10 module procedure l_nm_array11 module procedure l_nm_darray11 end interface l_nm !--------------------------------------------------------------------- !> スペクトルデータの位置(l) より，対応する波数 (n,m) を返す !> interface nm_l module procedure nm_l_int module procedure nm_l_dint module procedure nm_l_array module procedure nm_l_darray end interface nm_l save im, jm, nm, mm, nn ! 格子点数, 切断波数を記憶 save it, t, p, r, jc ! 変換用配列を記憶 save icom save mi 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 !--------------------------------------------------------------------- !> スペクトル変換の格子点数, 波数および OPENMP 使用時の !> 最大スレッド数を設定する. !> !> 実際の使用には上位サブルーチン w_mpi_initial() を用いること. !> subroutine w_base_mpi_Initial(n_in,i_in,j_in,mpi_comm,mint) !> 格子点数(東西), 2の巾乗(<=2048) integer,intent(in) :: i_in !> 格子点数(南北), 4 の倍数 integer,intent(in) :: j_in !> 切断全波数 integer,intent(in) :: n_in !> コミュニケーター integer,intent(in), optional :: mpi_comm !> 経度方向対称性 integer,intent(in), optional :: mint integer :: i, j integer :: ierr ! 三角波数切断 im = i_in jm = j_in nn = n_in nm = n_in+1 mm = n_in if ( present(mpi_comm) ) then icom = mpi_comm else icom = MPI_COMM_WORLD endif if ( present(mint) )then mi = mint else mi = 1_8 end if 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*mi)+1)*(2*nm-mm/(np*mi)*np*mi)/4+mm/(np*mi)+1)) allocate(jc((mm/(np*mi)+1)*(2*nm-mm/(np*mi)*np*mi)/16+mm/(np*mi)+1)) call mxallc(pp,jm/2*(5+2*(mm/(np*mi)+1))) call c_f_pointer(pp,p,[jm/2,5+2*(mm/(np*mi)+1)]) ! 変換用作業配列 allocate(c((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi))) call syinm1(mm,nm,im,it,t,r,icom,mi) call syinm2(mm,nm,jm,ig,p,r,jc,icom,mi) call syinmc(mm,nn,c,icom,mi) ! 変換用変数割付 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*mi)+1)*np*mi*2) call c_f_pointer(pw, w, [2*jv*((jm/jv-1)/np+1)*(mm/(np*mi)+1)*np*mi*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*mi)+1)*np*mi*2*2) call c_f_pointer(pww, ww, [2*jv*((jm/jv-1)/np+1)*(mm/(np*mi)+1)*np*mi*2*2]) ! 各プロセスの担当緯度, 波数情報 call syqrnj(jm,jv,j1,j2,icom) ! 分散格子データの大きさ jl = int(j2-j1+1) ! 分散スペクトルデータの大きさ nc = int((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) 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( int(j1), 1, MPI_INTEGER, & j1p(0), 1, MPI_INTEGER, int(icom), ierr) call MPI_ALLGATHER( int(j2), 1, MPI_INTEGER, & j2p(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/mi ! 経度座標 x_Lon_Weight(i) = 2*(pi/im)/mi ! 経度座標重み end do 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) end if end do 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_mpi_module_base_mint (2020/11/16) is initialized') end subroutine w_base_mpi_Initial !--------------------------------------------------------------------- !> スペクトルデータから格子データへ変換する(1 層用). !> function xv_w(w_data,ipow,iflag) !> スペクトルデータ real(DP), intent(in) :: w_data(nc) !> 作用させる \f$1/\cosφ\f$ の次数. 省略時は 0. integer, intent(in), optional :: ipow !>変換の種類: 省略時は 0. !> * 0 : 通常の逆変換 !> * -1 : 経度微分を作用させた逆変換 !> * 1 : 緯度微分 \f$\cosφ・∂/∂φ\f$ を作用させた逆変換 !> * 2 : \f$\sinφ\f$を作用させた逆変換 integer, intent(in), optional :: iflag !> 格子点データ real(DP) :: xv_w(0:im-1,1:jl) integer, parameter :: ipow_default = 0 integer, parameter :: iflag_default = 0 integer(8) :: ipval, ifval ! 作業用スペクトルデータ(SYCSMX 出力用) real(DP) :: w_Xdata((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SYCSMY 出力用) real(DP) :: w_Ydata((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','xv_w',& 'w_base_mpi_module 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 sytsmg(mm,nm,nn,im,jm,jv,w_data,g,it,t,p,r,jc,w,ipval,icom,mi) xv_w = g(:,1:jl) else if( ifval == -1 ) then call sycsmx(mm,nn,w_data,w_Xdata,icom,mi) call sytsmg(mm,nm,nn,im,jm,jv,w_Xdata,g,it,t,p,r,jc,w,ipval,icom,mi) xv_w = g(:,1:jl) else if( ifval == 1 ) then call sycsmy(mm,nn,w_data,w_Ydata,c,icom,mi) call sytsmg(mm,nm,nm,im,jm,jv,w_Ydata,g,it,t,p,r,jc,w,ipval,icom,mi) xv_w = g(:,1:jl) else if( ifval == 2 ) then call sytsmg(mm,nm,nn,im,jm,jv,w_data,g,it,t,p,r,jc,w,ipval,icom,mi) xv_w = g(:,1:jl)*sin(xv_Lat) else call MessageNotify('E','xv_w','invalid value of iflag') endif end function xv_w !--------------------------------------------------------------------- !> 格子データからスペクトルデータへ(正)変換する(1 層用). !> function w_xv(xv_data,ipow,iflag) !> 格子点データ real(DP), intent(in) :: xv_data(0:im-1,1:jl) !> 変換時に同時に作用させる \f$1/\cosφ\f$ の次数. 省略時は 0. integer, intent(in), optional :: ipow !> 変換の種類. 省略時は 0. !> * 0 : 通常の正変換 !> * -1 : 経度微分を作用させた正変換 !> * 1 : 緯度微分 \f$1/\cosφ・∂(f \cos^2φ)/∂φ\f$ !> を作用させた正変換 !> * 2 : \f$\sinφ\f$を作用させた正変換 integer, intent(in), optional :: iflag !> スペクトルデータ real(DP) :: w_xv(nc) integer, parameter :: ipow_default = 0 ! スイッチデフォルト値 integer, parameter :: iflag_default = 0 ! スイッチデフォルト値 integer(8) :: ipval, ifval ! 作業用スペクトルデータ(SYCSMX 出力用) real(DP) :: w_Xdata((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SYCSMY 出力用) real(DP) :: w_Ydata((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','xv_w',& 'w_base_mpi_module 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 sytgms(mm,nm,nn,im,jm,jv,w_xv,g,it,t,p,r,jc,w,ipval,icom,mi) else if ( ifval == -1 ) then g = 0.0D0 ; g(:,1:jl) = xv_data call sytgms(mm,nm,nn,im,jm,jv,w_Xdata,g,it,t,p,r,jc,w,ipval,icom,mi) call sycsmx(mm,nn,w_Xdata,w_xv,icom,mi) else if ( ifval == 1 ) then g = 0.0D0 ; g(:,1:jl) = xv_data call sytgms(mm,nm,nm,im,jm,jv,w_Ydata,g,it,t,p,r,jc,w,ipval,icom,mi) call sycyms(mm,nn,w_Ydata,w_xv,c,icom,mi) else if ( ifval == 2 ) then g = 0.0D0 ; g(:,1:jl) = xv_data*sin(xv_Lat) call sytgms(mm,nm,nn,im,jm,jv,w_xv,g,it,t,p,r,jc,w,ipval,icom,mi) end if end function w_xv !--------------------------------------------------------------------- !> 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に !> (逆)変換する(1 層用, MPI) !> !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. !> \f{align*}{ !> u \cosφ &= ∂χ/∂λ - \cosφ∂ψ/∂φ, \\ !> v \cosφ &= \cosφ∂χ/∂φ + ∂ψ/∂λ !> \f} !> subroutine w_StreamPotential2VectorMPI(w_Psi, w_Chi, xv_U, xv_V) !> 流線関数 real(DP), intent(in) :: w_Psi(nc) !> 速度ポテンシャル real(DP), intent(in) :: w_Chi(nc) !> 速度経度成分 real(DP), intent(out) :: xv_U(0:im-1,1:jl) !> 速度緯度成分 real(DP), intent(out) :: xv_V(0:im-1,1:jl) real(DP) :: w_Rdata1((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) real(DP) :: w_Rdata2((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SYTS2G 出力用) real(DP) :: w_Xdata((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SYCSMX 出力用) real(DP) :: w_Ydata((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SYCSMY 出力用) if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','w_StreamPotential2VectorMPI',& 'w_base_module not initialize yet.') endif ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ ! call sycsmx(mm,nn,w_Chi,w_Xdata,icom,mi) call sycsmy(mm,nn,w_Psi,w_Ydata,c,icom,mi) call sycmpk(mm,nn,nm,w_Xdata,w_Rdata1,icom,mi) w_Rdata1 = w_Rdata1 - w_Ydata ! ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ ! call sycsmy(mm,nn,w_Chi,w_Ydata,c,icom,mi) call sycsmx(mm,nn,w_Psi,w_Xdata,icom,mi) call sycmpk(mm,nn,nm,w_Xdata,w_Rdata2,icom,mi) w_Rdata2= w_Rdata2 + w_Ydata ! call sytsmv(mm,nm,nm,im,jm,jv, & & w_Rdata1,w_Rdata2,g1,g2,it,t,p,r,jc,ww,1_8,icom,mi) ! u xv_U = g1(:,1:jl) ! v xv_V = g2(:,1:jl) end subroutine w_StreamPotential2VectorMPI !--------------------------------------------------------------------- !> 速度場(格子データ)から渦度・発散(スペクトルデータ)に !> (正)変換する(1 層用, MPI) !> !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. !> \f{align*}{ !> ζ & = 1/\cosφ∂v/∂λ - 1/\cosφ ∂(u \cosφ)/∂φ, \\ !> D & = 1/\cosφ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ !> \f} !> subroutine w_Vector2VorDivMPI(xv_U, xv_V, w_Vor, w_Div) !> 速度経度成分 real(DP), intent(in) :: xv_U(0:im-1,1:jl) !> 速度緯度成分 real(DP), intent(in) :: xv_V(0:im-1,1:jl) !> 流線関数 real(DP), intent(out) :: w_Vor(nc) !> 速度ポテンシャル real(DP), intent(out) :: w_Div(nc) real(DP) :: w_Xdata((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SUCS2X 出力用) real(DP) :: w_Ydata1((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) real(DP) :: w_Ydata2((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SUCS2Y 出力用) real(DP) :: w_Data1((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) real(DP) :: w_Data2((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) 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 sytvms(mm,nm,nn+1,im,jm,jv,& & w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,1_8,icom,mi) call sycmpk(mm,nm,nn,w_Ydata1,w_Xdata,icom,mi) call sycsmx(mm,nn,w_Xdata,w_Div,icom,mi) call sycyms(mm,nn,w_Ydata2,w_Data2,c,icom,mi) ! ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ ! call sycmpk(mm,nm,nn,w_Ydata2,w_Xdata,icom,mi) call sycsmx(mm,nn,w_Xdata,w_Vor,icom,mi) call sycyms(mm,nn,w_Ydata1,w_Data1,c,icom,mi) ! ! 渦度・発散の計算 ! ! ζ = 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 !--------------------------------------------------------------------- !> 速度場(格子データ)から渦度・発散(スペクトルデータ)に !> (正)変換する(1 層用, MPI) !> !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. !> \f{align*}{ !> ζ & = 1/\cosφ∂v/∂λ - 1/\cosφ ∂(u \cosφ)/∂φ, \\ !> D & = 1/\cosφ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ !> \f} !> subroutine w_VectorCosLat2VorDivMPI(xv_UCosLat, xv_VCosLat, w_Vor, w_Div) !> 速度経度成分 * cos(lat) real(DP), intent(in) :: xv_UCosLat(0:im-1,1:jl) !> 速度緯度成分 * cos(lat) real(DP), intent(in) :: xv_VCosLat(0:im-1,1:jl) !> 流線関数 real(DP), intent(out) :: w_Vor(nc) !> 速度ポテンシャル real(DP), intent(out) :: w_Div(nc) real(DP) :: w_Xdata((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SUCS2X 出力用) real(DP) :: w_Ydata1((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) real(DP) :: w_Ydata2((mm/(np*mi)+1)*(2*(nm+1)-mm/(np*mi)*np*mi)) ! 作業用スペクトルデータ(SUCS2Y 出力用) real(DP) :: w_Data1((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) real(DP) :: w_Data2((mm/(np*mi)+1)*(2*(nn+1)-mm/(np*mi)*np*mi)) 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 sytvms(mm,nm,nn+1,im,jm,jv,w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,2_8,icom,mi) call sycmpk(mm,nm,nn,w_Ydata1,w_Xdata,icom,mi) call sycsmx(mm,nn,w_Xdata,w_Div,icom,mi) call sycyms(mm,nn,w_Ydata2,w_Data2,c,icom,mi) ! ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算 ! call sycmpk(mm,nm,nn,w_Ydata2,w_Xdata,icom,mi) call sycsmx(mm,nn,w_Xdata,w_Vor,icom,mi) call sycyms(mm,nn,w_Ydata1,w_Data1,c,icom,mi) ! ! 渦度・発散の計算 ! ! ζ = 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 !--------------------------------------------------------------------- !> 分割緯度データを集める. !> function ay_av(av) !> 緯度格子点データ(分割) (:,jl) real(DP), dimension(:,:), intent(in) :: av !> 緯度格子点データ real(DP), dimension(size(av,1),jm) :: ay_av 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) !> 緯度格子点データ(分割) (:,jc) real(DP), dimension(:,:), intent(in) :: ay !> 緯度格子点データ real(DP), dimension(size(ay,1),jl) :: av_ay av_ay = ay(:,j1:j2) end function av_ay !--------------------------------------------------------------------- !> モジュールの終了処理(割り付け配列の解放)をおこなう. !> !> 実際の使用には上位サブルーチン w_mpi_finalize() を用いること. !> subroutine w_base_mpi_Finalize if ( .not. w_base_mpi_initialize ) then call MessageNotify('W','w_base_mpi_Finalize',& 'w_mpi_module_base 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_mpi_module_base_mint (2020/11/16) is finalized') end subroutine w_base_mpi_Finalize !---------------------------------------------------------------------! !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 引数 n,m がともに整数値の場合, 整数値を返す. !> プロセスに n,m が存在しない場合 0 を返す !> function l_nm_array00(n_in,m_in) !> 全波数 integer, intent(in) :: n_in !> 帯状波数 integer, intent(in) :: m_in !> スペクトルデータの格納位置 integer :: l_nm_array00 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 not initialize yet. Use synmml routine in ISPACK3 directly.') end if n = n_in m = m_in call synmml(nn,n,m,ip_w,l,icom,mi) if ( ip_w /= ip ) then l_nm_array00 = 0 else l_nm_array00 = int(l) endif end function l_nm_array00 !---------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 引数 n,m がともに整数値の場合, 整数値を返す. !> プロセスに n,m が存在しない場合 0 を返す !> function l_nm_darray00(n_in,m_in) !> 全波数 integer(8), intent(in) :: n_in !> 帯状波数 integer(8), intent(in) :: m_in !> スペクトルデータの格納位置 integer :: l_nm_darray00 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 not initialize yet. Use synmml routine in ISPACK3 directly.') end if n = n_in m = m_in call synmml(nn,n,m,ip_w,l,icom,mi) if ( ip_w /= ip ) then l_nm_darray00 = 0 else l_nm_darray00 = int(l) endif end function l_nm_darray00 !--------------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, !> marray と同じ大きさの 1 次元整数配列を返す. !> function l_nm_array01(n,marray) !> 全波数 integer, intent(in) :: n !> 帯状波数 integer, intent(in) :: marray(:) !> スペクトルデータ位置 integer :: l_nm_array01(size(marray)) integer :: i do i=1, size(marray) l_nm_array01(i) = l_nm_array00(n,marray(i)) enddo end function l_nm_array01 !--------------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, !> marray と同じ大きさの 1 次元整数配列を返す. !> function l_nm_darray01(n,marray) !> 全波数 integer(8), intent(in) :: n !> 帯状波数 integer(8), intent(in) :: marray(:) !> スペクトルデータ位置 integer :: l_nm_darray01(size(marray)) integer :: i do i=1, size(marray) l_nm_darray01(i) = l_nm_darray00(n,marray(i)) enddo end function l_nm_darray01 !--------------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 第 1 引数 narray が整数 1 次元配列, 第 2 引数 m が整数の場合, !> narray と同じ大きさの 1 次元整数配列を返す. !> function l_nm_array10(narray,m) !> 全波数 integer, intent(in) :: narray(:) !> 帯状波数 integer, intent(in) :: m !> スペクトルデータ位置 integer :: l_nm_array10(size(narray)) integer :: i do i=1, size(narray) l_nm_array10(i) = l_nm_array00(narray(i),m) enddo end function l_nm_array10 !--------------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 第 1 引数 narray が整数 1 次元配列, 第 2 引数 m が整数の場合, !> narray と同じ大きさの 1 次元整数配列を返す. !> function l_nm_darray10(narray,m) !> 全波数 integer(8), intent(in) :: narray(:) !> 帯状波数 integer(8), intent(in) :: m !> スペクトルデータ位置 integer :: l_nm_darray10(size(narray)) integer :: i do i=1, size(narray) l_nm_darray10(i) = l_nm_darray00(narray(i),m) enddo end function l_nm_darray10 !--------------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, !> narray, marray と同じ大きさの 1 次元整数配列を返す. !> narray, marray は同じ大きさでなければならない. !> function l_nm_array11(narray,marray) !> 全波数 integer, intent(in) :: narray(:) !> 帯状波数 integer, intent(in) :: marray(:) !> スペクトルデータ位置 integer :: l_nm_array11(size(narray)) 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 !--------------------------------------------------------------------- !> 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す. !> !> 第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, !> narray, marray と同じ大きさの 1 次元整数配列を返す. !> narray, marray は同じ大きさでなければならない. !> function l_nm_darray11(narray,marray) !> 全波数 integer(8), intent(in) :: narray(:) !> 帯状波数 integer(8), intent(in) :: marray(:) !> スペクトルデータ位置 integer :: l_nm_darray11(size(narray)) 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_darray11(i) = l_nm_darray00(narray(i),marray(i)) enddo end function l_nm_darray11 !--------------------------------------------------------------------- !> スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す. !> !> 引数 l が整数値の場合, 対応する全波数と帯状波数を !> 長さ 2 の 1 次元整数値を返す. !> nm_l(1) が全波数, nm_l(2) が帯状波数である. !> function nm_l_int(l) !> スペクトルデータの格納位置 integer, intent(in) :: l !> 全波数, 帯状波数 integer :: nm_l_int(2) integer(8) :: n,m ! 全波数, 帯状波数 if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','nm_l_int', & 'w_base_mpi_module not initialize yet. Use sylmnm routine in ISPACK3 directly.') endif call sylmnm(mm,nn,int(l,kind=8),n,m,icom,mi) nm_l_int = (/int(n),int(m)/) end function nm_l_int !--------------------------------------------------------------------- !> スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す. !> !> 引数 l が整数値の場合, 対応する全波数と帯状波数を !> 長さ 2 の 1 次元整数値を返す. !> nm_l(1) が全波数, nm_l(2) が帯状波数である. !> function nm_l_dint(l) !> スペクトルデータの格納位置 integer(8), intent(in) :: l !> 全波数, 帯状波数 integer :: nm_l_dint(2) integer(8) :: n,m ! 全波数, 帯状波数 if ( .not. w_base_mpi_initialize ) then call MessageNotify('E','nm_l_int', & 'w_base_mpi_module not initialize yet. Use sylmnm routine in ISPACK3 directly.') endif call sylmnm(mm,nn,l,n,m,icom,mi) nm_l_dint = (/int(n),int(m)/) end function nm_l_dint !--------------------------------------------------------------------- !> スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す. !> !> 引数 larray が整数 1 次元配列の場合, !> larray に対応する n, m を格納した 2 次元整数配列を返す. !> nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である. !> function nm_l_array(larray) !> 全波数, 帯状波数 integer, intent(in) :: larray(:) !> スペクトルデータの格納位置 integer :: nm_l_array(size(larray),2) integer :: i do i=1, size(larray) nm_l_array(i,:) = nm_l_int(larray(i)) enddo end function nm_l_array !--------------------------------------------------------------------- !> スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す. !> !> 引数 larray が整数 1 次元配列の場合, !> larray に対応する n, m を格納した 2 次元整数配列を返す. !> nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である. !> function nm_l_darray(larray) !> 全波数, 帯状波数 integer(8), intent(in) :: larray(:) !> スペクトルデータの格納位置 integer :: nm_l_darray(size(larray),2) integer :: i do i=1, size(larray) nm_l_darray(i,:) = nm_l_dint(larray(i)) enddo end function nm_l_darray end module w_mpi_module_base_mint