!-- !---------------------------------------------------------------------- ! Copyright (c) 2019 SPMDOEL Development Group !---------------------------------------------------------------------- !表題 ua_base_mpi_module_sypack ! ! spml/ua_base_mpi_module_sypack モジュールは球面上での流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール ua_mpi_module_sypack の下部モジュールであり, ! スペクトル計算の基本的な Fortran90 関数を提供する. ! ! 球面上の 1 層モデル用 u_base_module_sypack モジュールを多層モデル用に ! 並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに ! 対する変換が行える. ! ! 内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK3 の SYPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって u_mpi_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! ! !履歴 2016/08/15 竹広真一 wa_base_module_sypack より多層 MPI 対応改造 ! !++ module ua_base_mpi_module_sypack ! != ua_base_mpi_module_sypack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/ua_base_mpi_module_sypack モジュールは球面上での流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール ua_mpi_module_sypack の下部モジュールであり, ! スペクトル計算の基本的な Fortran90 関数を提供する. ! ! 球面上の 1 層モデル用 w_base_module_sypack モジュールを多層モデル用に ! 並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに ! 対する変換が行える. ! ! 内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK3 の SYPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_mpi_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! use dc_message use w_base_mpi_module_sypack, only : im, jm, jc=>jl, nn, mm, & xv_Lon, xv_Lat, x_Lon, v_Lat, x_Lon_Weight, v_Lat_Weight, & xv_w, w_xv, & w_StreamPotential2VectorMPI, w_Vector2VorDivMPI, & w_VectorCosLat2VorDivMPI, & w_base_mpi_initial, w_base_mpi_finalize, nc, l_nm, nm_l use mpi implicit none integer :: km=16 ! 同時に処理する最大データ数(層の数). integer :: nproc, nproc_h, nproc_v integer :: iproc, iproc_h, iproc_v integer :: ks, ke, kc integer, allocatable :: ksp(:), kep(:), kcp(:) integer :: nms, nme, nmc integer, allocatable :: nmsp(:), nmep(:), nmcp(:) integer :: MPI_COMM_W integer :: MPI_COMM_H integer :: MPI_COMM_V logical :: first_wb_ua=.true. logical :: first_ua_wb=.true. logical :: first_wa_ua=.true. logical :: first_aa_ab=.true. private !public km ! 層数 public ua_base_mpi_Initial ! 初期化サブルーチン public ua_base_mpi_Finalize ! 終了処理 public MPI_COMM_W, MPI_COMM_H, MPI_COMM_V public nproc, nproc_h, nproc_v public iproc, iproc_h, iproc_v public jc, nc public km public ks, ke, kc public nms, nme, nmc public l_nm, nm_l public lu_nm, nm_lu public kp_kl public xv_Lon, xv_Lat public x_Lon, v_Lat public x_Lon_Weight, v_Lat_Weight public xvb_ua, ua_xvb ! 変換関数 public wb_ua, ua_wb ! 変換関数 public wa_ua, ua_wa ! 変換関数 public aa_ab, ab_aa ! 変換関数 public ua_StreamPotential2VectorMPI ! 流線ポテンシャルから速度場計算 public ua_Vector2VorDivMPI ! 速度場から渦度発散を計算 public ua_VectorCosLat2VorDivMPI ! 速度場から渦度発散を計算 save km ! 最大データ数(層数)を記憶 save nproc, iproc, nproc_h, iproc_h, nproc_v, iproc_v save kc, ks, ke save MPI_COMM_W, MPI_COMM_H, MPI_COMM_V save first_wb_ua, first_ua_wb, first_wa_ua, first_aa_ab contains !--------------- 初期化 ----------------- subroutine ua_base_mpi_initial(n_in,i_in,j_in,k_in,npv,mpi_comm) ! ! 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,intent(in),optional :: npv !(in) 層方向分割数 integer,intent(in),optional :: mpi_comm !(in) 全体コミュニケーター integer :: ierr, iv km = k_in if ( present(npv) ) then nproc_v = npv else nproc_v = 1 endif if ( present(mpi_comm) ) then MPI_COMM_W = mpi_comm else MPI_COMM_W = MPI_COMM_WORLD endif CALL MPI_COMM_SIZE(MPI_COMM_W,NPROC,IERR) CALL MPI_COMM_RANK(MPI_COMM_W,IPROC,IERR) NPROC_H = NPROC/NPROC_V IPROC_V = IPROC/NPROC_H IPROC_H = MOD(IPROC,NPROC_H) CALL MPI_COMM_SPLIT(MPI_COMM_W, iproc_v, iproc_h, MPI_COMM_H, IERR ) CALL MPI_COMM_SPLIT(MPI_COMM_W, iproc_h, iproc_v, MPI_COMM_V, IERR ) !!$ CALL MPI_COMM_RANK( MPIA_COMM_H, MYID_H_TMP , NERR ) !!$ CALL MPI_COMM_SIZE( MPIA_COMM_H, NPROC_H_TMP, NERR ) !!$ !!$ CALL MPI_COMM_RANK( MPIA_COMM_V, MYID_V_TMP , NERR ) !!$ CALL MPI_COMM_SIZE( MPIA_COMM_V, NPROC_V_TMP, NERR ) !!$ call set_layerdiv(npv,iproc_v,ks,ke) !!$ kc = ke-ks+1 ! 分散多層データの大きさ call w_base_mpi_initial(n_in,i_in,j_in,mpi_comm_h) allocate(ksp(0:NPROC_V-1)) allocate(kep(0:NPROC_V-1)) allocate(kcp(0:NPROC_V-1)) call mpidevide(km,nproc_v,ksp,kep,kcp) ksp = ksp ; kep = kep ks = ksp(iproc_v) ; ke = kep(iproc_v) ; kc = kcp(iproc_v) IF(iproc.EQ.0)THEN WRITE(6,*)' RANK KSTART KEND KC ' WRITE(6,*)'------------------------------' DO IV=0,NPROC_V-1 WRITE(6,'(4i7)') IV, KSP(IV), KEP(IV), KCP(IV) ENDDO WRITE(6,*)'' END IF allocate(nmsp(0:NPROC_V-1)) allocate(nmep(0:NPROC_V-1)) allocate(nmcp(0:NPROC_V-1)) call mpidevide(nc,nproc_v,nmsp,nmep,nmcp) nms = nmsp(iproc_v) ; nme = nmep(iproc_v) ; nmc = nmcp(iproc_v) IF(iproc.EQ.0)THEN WRITE(6,*)' RANK NSTART NEND NC ' WRITE(6,*)'------------------------------' DO IV=0,NPROC_V-1 WRITE(6,'(4i7)') IV, NMSP(IV), NMEP(IV), NMCP(IV) ENDDO WRITE(6,*)'' ENDIF call MessageNotify('M','ua_base_mpi_initial', & 'ua_base_mpi_module_sypack (2019/04/26) 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)/npv+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,npv-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 subroutine mpidevide(mm,np,msp,mep,mcp) ! ! 鉛直分割情報 ! integer, intent(IN) :: mm ! 全格子点数 integer, intent(IN) :: np ! 分割プロセス数 integer, intent(OUT) :: msp(0:np-1) ! 分割開始インデックス integer, intent(OUT) :: mep(0:np-1) ! 分割終了インデックス integer, intent(OUT) :: mcp(0:np-1) ! 分割インデックス数 integer :: i, mint, mmod mint = mm/np mmod = mod(mm,np) mcp = mint do i = 0,mmod-1 mcp(i) = mcp(i)+1 enddo msp(0)=1 ; mep(0) = msp(0) + mcp(0) -1 do i = 1,np-1 msp(i) = msp(i-1) + mcp(i-1) mep(i) = msp(i) + mcp(i) - 1 end do end subroutine mpidevide !--------------- 基本関数 ----------------- function kp_kl(kl) integer, intent(IN) :: kl ! 層番号(全体) integer :: kp_kl ! 層番号(分割プロセス中) if ( ks .le. kl .AND. kl .le. ke ) then kp_kl = kl - ks + 1 else kp_kl = 0 endif end function kp_kl function lu_nm(n,m) integer, intent(IN) :: n ! 全波数 integer, intent(IN) :: m ! 帯状波数 integer :: lu_nm ! 位置情報(分割プロセス中) integer :: ll ! 位置情報(全体) ll = l_nm(n,m) if ( nms .le. ll .AND. ll .le. nme ) then lu_nm = ll-nms+1 else lu_nm = 0 endif end function lu_nm function nm_lu(l) integer, intent(IN) :: l ! 位置情報(分割プロセス中) integer :: nm_lu(2) ! 全波数, 帯状波数 nm_lu = nm_l(l+nms-1) end function nm_lu !--------------- データ入れ換え ----------------- function wb_ua(ua) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(nmc,km), intent(in) :: ua !(in) 2 次元球面調和函数(水平再分割) (nc,:) real(8), dimension(nc,kc) :: wb_ua !(out) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), allocatable :: ua_recv(:,:) integer :: IRQSTS integer :: IST( MPI_STATUS_SIZE ) integer :: ierr, iv integer, parameter :: itag=110 integer, allocatable ::snd_rank_v(:) integer, allocatable ::rcv_rank_v(:) integer, allocatable ::rcv_idx(:) save snd_rank_v, rcv_rank_v, rcv_idx if ( first_wb_ua ) then first_wb_ua = .false. if ( allocated(snd_rank_v) ) deallocate(snd_rank_v) if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v) if ( allocated(rcv_idx) ) deallocate(rcv_idx) allocate(snd_rank_v(1:nproc_v-1)) allocate(rcv_rank_v(1:nproc_v-1)) allocate(rcv_idx(0:nproc_v-1)) do iv=1,nproc_v-1 rcv_rank_v(iv) = mod(iproc_v+iv,nproc_v) snd_rank_v(iv) = mod(iproc_v-iv+nproc_v,nproc_v) enddo rcv_idx(0) = 1 do iv=1,nproc_v-1 rcv_idx(iv) = rcv_idx(iv-1)+nmcp(iv-1) enddo endif call MPI_Barrier(MPI_COMM_V,IERR) do iv=1,nproc_v-1 call MPI_ISEND(ua(1,ksp(rcv_rank_v(iv))),nmc*kcp(rcv_rank_v(iv)), & MPI_REAL8, rcv_rank_v(iv),& itag,MPI_COMM_V,IRQSTS,IERR) allocate(ua_recv(nmcp(snd_rank_v(iv)),kc)) call MPI_RECV(ua_recv,kc*nmcp(snd_rank_v(iv)), MPI_REAL8,& snd_rank_v(iv), itag, MPI_COMM_V, IST, IERR) wb_ua(nmsp(snd_rank_v(iv)):nmep(snd_rank_v(iv)),:) = ua_recv deallocate(ua_recv) CALL MPI_WAIT(IRQSTS,IST,IERR) enddo wb_ua(nms:nme,:) = ua(:,ks:ke) end function wb_ua function ua_wb(wb) ! ! 水平スペクトル・層分割データを集める. ! real(8), dimension(nc,kc), intent(in) :: wb !(in) 動径格子点データ(分割) real(8), dimension(nmc,km) :: ua_wb !(out) 動径格子点データ real(8), allocatable :: wb_send(:,:) integer :: IRQSTS integer :: IST( MPI_STATUS_SIZE ) integer :: ierr, iv integer, parameter :: itag=111 integer, allocatable ::snd_rank_v(:) integer, allocatable ::rcv_rank_v(:) integer, allocatable ::rcv_idx(:) save snd_rank_v, rcv_rank_v, rcv_idx if ( first_ua_wb ) then first_ua_wb = .false. if ( allocated(snd_rank_v) ) deallocate(snd_rank_v) if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v) if ( allocated(rcv_idx) ) deallocate(rcv_idx) allocate(snd_rank_v(1:nproc_v-1)) allocate(rcv_rank_v(1:nproc_v-1)) allocate(rcv_idx(0:nproc_v-1)) do iv=1,nproc_v-1 rcv_rank_v(iv) = mod(iproc_v+iv,nproc_v) snd_rank_v(iv) = mod(iproc_v-iv+nproc_v,nproc_v) enddo rcv_idx(0) = 1 do iv=1,nproc_v-1 rcv_idx(iv) = rcv_idx(iv-1)+kcp(iv-1) enddo endif call MPI_Barrier(MPI_COMM_V,IERR) do iv=1,nproc_v-1 allocate(wb_send(nmcp(rcv_rank_v(iv)),kc)) wb_send = wb(nmsp(rcv_rank_v(iv)):nmep(rcv_rank_v(iv)),:) call MPI_ISEND(wb_send, nmcp(rcv_rank_v(iv))*kc,& MPI_REAL8, rcv_rank_v(iv),& itag,MPI_COMM_V,IRQSTS,IERR) deallocate(wb_send) call MPI_RECV(ua_wb(1,rcv_idx(snd_rank_v(iv))),& nmc*kcp(snd_rank_v(iv)), MPI_REAL8, & snd_rank_v(iv), itag, MPI_COMM_V, IST, IERR) CALL MPI_WAIT(IRQSTS,IST,IERR) enddo ua_wb(:,ks:ke) = wb(nms:nme,:) end function ua_wb function wa_ua(ua) ! ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension(:,:), intent(in) :: ua !(in) 2 次元球面調和函数(水平再分割) (nc,:) real(8), dimension(nc,size(ua,2)) :: wa_ua !(out) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(size(ua,2),nmc) :: au real(8), dimension(size(ua,2),nc) :: aw integer :: a_size integer :: IRQSTS integer :: IST( MPI_STATUS_SIZE ) integer :: ierr, ip_v integer, parameter :: itag=110 ! integer, allocatable ::snd_rank(:) integer, allocatable ::snd_rank_v(:) ! integer, allocatable ::rcv_rank(:) integer, allocatable ::rcv_rank_v(:) integer, allocatable ::rcv_idx(:) save snd_rank_v, rcv_rank_v, rcv_idx !, snd_rank, rcv_rank if ( first_wa_ua ) then first_wa_ua = .false. if ( allocated(snd_rank_v) ) deallocate(snd_rank_v) if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v) if ( allocated(rcv_idx) ) deallocate(rcv_idx) allocate(snd_rank_v(1:nproc_v-1)) allocate(rcv_rank_v(1:nproc_v-1)) allocate(rcv_idx(0:nproc_v-1)) do ip_v=1,nproc_v-1 rcv_rank_v(ip_v) = mod(iproc_v+ip_v,nproc_v) snd_rank_v(ip_v) = mod(iproc_v-ip_v+nproc_v,nproc_v) enddo rcv_idx(0) = 1 do ip_v=1,nproc_v-1 rcv_idx(ip_v) = rcv_idx(ip_v-1)+nmcp(ip_v-1) enddo endif au = transpose(ua) a_size = size(ua,2) call MPI_Barrier(MPI_COMM_V,IERR) do ip_v=1,nproc_v-1 call MPI_ISEND(au, a_size*nmc, MPI_REAL8, rcv_rank_v(ip_v),& itag,MPI_COMM_V,IRQSTS,IERR) call MPI_RECV(aw(1,rcv_idx(snd_rank_v(ip_v))), & a_size*nmcp(snd_rank_v(ip_v)), MPI_REAL8, & snd_rank_v(ip_v), itag, MPI_COMM_V, IST, IERR) CALL MPI_WAIT(IRQSTS,IST,IERR) enddo aw(:,rcv_idx(iproc_v):rcv_idx(iproc_v)+nmcp(iproc_v)-1) = au wa_ua = transpose(aw) end function wa_ua function ua_wa(wa) ! ! 水平スペクトルデータを再分割する ! real(8), dimension(:,:), intent(in) :: wa !(in) 2 次元球面調和函数データ (nc,:) real(8), dimension(nmc,size(wa,2)) :: ua_wa !(out) 2 次元球面調和函数データ(水平再分割) ua_wa = wa(nms:nme,:) end function ua_wa !!$ function wb_ua(ua_data) !!$ real(8), intent(IN) :: ua_data(nmc,km) !!$ real(8) :: wb_ua(nc,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 aa_ab(ab) ! ! 水平スペクトル・動径格子点データを集める. ! real(8), dimension(:,:), intent(in) :: ab !(in) 動径格子点データ(分割) real(8), dimension(size(ab,1),km) :: aa_ab !(out) 動径格子点データ integer :: a_size integer :: IRQSTS integer :: IST( MPI_STATUS_SIZE ) integer :: ierr, ip_v integer, parameter :: itag=111 integer, allocatable ::snd_rank_v(:) integer, allocatable ::rcv_rank_v(:) integer, allocatable ::rcv_idx(:) save snd_rank_v, rcv_rank_v, rcv_idx if ( first_aa_ab ) then first_aa_ab = .false. if ( allocated(snd_rank_v) ) deallocate(snd_rank_v) if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v) if ( allocated(rcv_idx) ) deallocate(rcv_idx) allocate(snd_rank_v(1:nproc_v-1)) allocate(rcv_rank_v(1:nproc_v-1)) allocate(rcv_idx(0:nproc_v-1)) do ip_v=1,nproc_v-1 rcv_rank_v(ip_v) = mod(iproc_v+ip_v,nproc_v) snd_rank_v(ip_v) = mod(iproc_v-ip_v+nproc_v,nproc_v) enddo rcv_idx(0) = 0 do ip_v=1,nproc_v-1 rcv_idx(ip_v) = rcv_idx(ip_v-1)+kcp(ip_v-1) enddo endif a_size = size(ab,1) call MPI_Barrier(MPI_COMM_V,IERR) do ip_v=1,nproc_v-1 call MPI_ISEND(ab, a_size*kc, MPI_REAL8, rcv_rank_v(ip_v),& itag,MPI_COMM_V,IRQSTS,IERR) call MPI_RECV(aa_ab(1,rcv_idx(snd_rank_v(ip_v))),& a_size*kcp(snd_rank_v(ip_v)), MPI_REAL8, & snd_rank_v(ip_v), itag, MPI_COMM_V, IST, IERR) CALL MPI_WAIT(IRQSTS,IST,IERR) enddo aa_ab(:,rcv_idx(iproc_v):rcv_idx(iproc_v)+kcp(iproc_v)-1) = ab end function aa_ab function ab_aa(aa) ! ! 水平スペクトル・動径格子点データを分割する. ! real(8), dimension(:,:), intent(in) :: aa !(in) 動径格子点データ(:,0:km) real(8), dimension(size(aa,1),kc) :: ab_aa !(out) 動径格子点データ(分割) ab_aa = aa(:,ks:ke) end function ab_aa !--------------- 基本変換 ----------------- function xvb_ua(ua_data,ipow,iflag) ! 球面調和関数スペクトル -> 格子点 ! ! スペクトルデータから格子データへ変換する(多層用). ! real(8), intent(in) :: ua_data(nmc,km) !(in) スペクトルデータ ! real(8) :: xvb_ua(0:im-1,1:jc,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(nc,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 xvb_ua(:,:,k) = xv_w(wb_data(:,k),iflag=ifval,ipow=ipval) enddo end function xvb_ua function ua_xvb(xvb_data,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル ! ! 格子データからスペクトルデータへ(正)変換する(多層用). ! real(8), intent(in) :: xvb_data(0:im-1,jc,kc) !(in) 格子点データ real(8) :: ua_xvb(nmc,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(nc,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_xv(xvb_data(:,:,k),iflag=ifval,ipow=ipval) enddo ua_xvb = ua_wb(wb_data) end function ua_xvb !----------- 速度, 渦度・発散, 流線・ポテンシャル計算 ------------- subroutine ua_StreamPotential2VectorMPI(ua_Psi, ua_Chi, xvb_U, xvb_V) ! ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に ! (逆)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ, ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ ! real(8), intent(in) :: ua_Psi(nmc,km) !(in) 流線関数((nm+1)*(nm+1),:) real(8), intent(in) :: ua_Chi(nmc,km) !(in) 速度ポテンシャル((nm+1)*(nm+1),:) real(8), intent(out) :: xvb_U(0:im-1,jc,kc) !(out) 速度経度成分(0:im-1,1:jm,:) real(8), intent(out) :: xvb_V(0:im-1,jc,kc) !(out) 速度緯度成分(0:im-1,1:jm,:) real(8) :: wb_Psi(nc,kc) !(in) 流線関数((nm+1)*(nm+1),:) real(8) :: wb_Chi(nc,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_StreamPotential2VectorMPI & ( wb_Psi(:,k), wb_Chi(:,k), xvb_U(:,:,k), xvb_V(:,:,k) ) enddo end subroutine ua_StreamPotential2VectorMPI subroutine ua_Vector2VorDivMPI(xvb_U, xvb_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) :: xvb_U(0:im-1,jc,kc) !(in) 速度経度成分(0:im-1,1:jm,:) real(8), intent(in) :: xvb_V(0:im-1,jc,kc) !(in) 速度緯度成分(0:im-1,1:jm,:) real(8), intent(out) :: ua_Vor(nmc,km) !(out) 渦度 real(8), intent(out) :: ua_Div(nmc,km) !(out) 発散 real(8) :: wb_Vor(nc,kc) ! 渦度 real(8) :: wb_Div(nc,kc) ! 発散 integer :: k do k=1,kc call w_Vector2VorDivMPI & (xvb_U(:,:,k), xvb_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(xvb_UCosLat, xvb_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) :: xvb_UCosLat(0:im-1,jc,kc) !(in) 速度経度成分 * cos(lat) (0:im-1,1:jm,:) real(8), intent(in) :: xvb_VCosLat(0:im-1,jc,kc) !(in) 速度緯度成分 * cos(lat) (0:im-1,1:jm,:) real(8), intent(out) :: ua_Vor(nmc,km) !(out) 流線関数 real(8), intent(out) :: ua_Div(nmc,km) !(out) 速度ポテンシャル real(8) :: wb_Vor(nc,kc) ! 渦度 real(8) :: wb_Div(nc,kc) ! 発散 integer :: k do k=1,kc call w_VectorCosLat2VorDivMPI & (xvb_UCosLat(:,:,k), xvb_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 を使用すること. ! deallocate(ksp) deallocate(kep) deallocate(kcp) deallocate(nmsp) deallocate(nmep) deallocate(nmcp) call w_base_mpi_Finalize first_wb_ua=.true. first_ua_wb=.true. first_wa_ua=.true. first_aa_ab=.true. call MessageNotify('M','ua_base_mpi_finalize',& 'ua_base_mpi_module_sypack is finalized (2019/05/29) ') end subroutine ua_base_mpi_Finalize end module ua_base_mpi_module_sypack