!-- !---------------------------------------------------------------------- ! Copyright (c) 2016 SPMODEL Development Group. !---------------------------------------------------------------------- !表題 wa_interpolate_mpi_module_supack ! ! spml/wa_interpolate_mpi_module_supack モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール wa_mpi_module_supack の下部モジュールであり, スペクトル法による ! 補間計算のための Fortran90 関数を提供する. ! ! 補間計算の方法については doc/wa_module.tex を参照のこと. ! 内部で ISPACK2 の SUPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK2/SUPACK のマニュアルを参照されたい. ! !履歴 2016/03/16 竹広真一 w_interpolate_mpi_module_supack より改造 ! !++ module wa_interpolate_mpi_module_supack ! != wa_interpolate_mpi_module_supack ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id: wa_interpolate_mpi_module_supack.f90 590 2013-08-19 08:48:21Z uwabami $ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/wa_interplata_mpi_module_supack モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール wa_module の下部モジュールであり, スペクトル法による ! 補間計算のための Fortran90 関数を提供する. ! ! 内部で ISPACK2 の SUPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK2/SUPACK のマニュアルを参照されたい. ! use dc_message, only : MessageNotify use w_base_mpi_module_supack, only : nm=>nn, nc, l_nm use mpi implicit none private public a_Interpolate_wa ! 補間関数 interface a_Interpolate_wa ! ! 緯度 alon, 経度 alat における関数値を ! その球面調和変換係数 wa_data から補間計算する ! ! 入力する緯度経度座標は, ! １点, 経度1点緯度複数点, 経度複数点緯度1点, 経度緯度複数点 ! の 4 種類 ! module procedure a_Interpolate_array00_wa end interface interface alpha module procedure alpha_array0 end interface interface Pmm module procedure Pmm_array0 end interface contains !--------------- 補間計算 ----------------- function a_Interpolate_array00_wa(wa_data,alon,alat) ! ! 緯度 alat, 経度 alon における関数値を ! その球面調和変換係数 wa_data から補間計算する ! real(8), intent(IN) :: wa_data(:,:) ! スペクトルデータ real(8), intent(IN) :: alon ! 補間する位置(経度) real(8), intent(IN) :: alat ! 補間する位置(緯度) real(8) :: a_Interpolate_array00_wa(size(wa_data,2)) ! 補間した値 real(8) :: mu real(8), dimension(size(wa_data,2)) :: y0, y1, y2, AnmPnm real(8), dimension(size(wa_data,2)) :: gdata logical :: LFlagAnd, LFlagOr integer :: k,m integer :: ierr logical,parameter :: wdata_check = .false. mu = sin(alat) gdata = 0.0D0 if ( wdata_check ) then !---- 波数チェック LFlagAnd = .true. ; LFlagOr = .false. do k=nm,0,-1 LFlagAnd = (l_nm(k,0) /= 0) .and. LFlagAnd LFlagOr = (l_nm(k,0) /= 0) .or. LFlagOr enddo if ( (.not. LFlagAnd) .and. LFlagOr ) then call MessageNotify & ('E','wa_Interpolate_array00','Invalid spectral data for m=0') endif do m=1,nm LFlagAnd = .true. ; LFlagOr = .false. do k=nm,m,-1 LFlagAnd = (l_nm(k,m) /= 0) .and. LFlagAnd LFlagOr = (l_nm(k,m) /= 0) .or. LFlagOr enddo if ( (.not. LFlagAnd) .and. LFlagOr ) then call MessageNotify & ('E','wa_Interpolate_array00','Invalid spectral data for m nonzedo') endif end do do m=1,nm LFlagAnd = .true. ; LFlagOr = .false. do k=nm,m,-1 LFlagAnd = (l_nm(k,-m) /= 0) .and. LFlagAnd LFlagOr = (l_nm(k,-m) /= 0) .or. LFlagOr enddo if ( (.not. LFlagAnd) .and. LFlagOr ) then call MessageNotify & ('E','wa_Interpolate_array00','Invalid spectral data for m nonzedo') endif end do endif !---- Σa_n^0 P_n^0 の計算 if ( l_nm(0,0) /= 0 ) then y2 = 0 ; y1 = 0 do k=nm,1,-1 y0 = alpha(k,0,mu) * y1 + beta(k+1,0)*y2 + wa_data(l_nm(k,0),:) y2 = y1 ; y1 = y0 enddo gdata = gdata + ( beta(1,0) * y2 + mu*sqrt(3.0D0) * y1 & + wa_data(l_nm(0,0),:) ) * Pmm(0,mu) endif !---- Σ Re[s_n^m] P_n^m exp(imλ)の計算 do m=1,nm if ( l_nm(m,m) /= 0 ) then y2 = 0 ; y1 = 0 do k=nm,m+1,-1 y0 = alpha(k,m,mu) * y1 + beta(k+1,m) * y2 + wa_data(l_nm(k,m),:) y2 = y1 ; y1 = y0 enddo AnmPnm =(wa_data(l_nm(m,m),:) + beta(m+1,m)*y2 & + mu*sqrt(2.0D0*m+3)*y1 ) * Pmm(m,mu) gdata = gdata + AnmPnm * 2 * cos(m*alon) endif end do !---- Σ Im[s_n^m] P_n^m exp(imλ)の計算 do m=1,nm if ( l_nm(m,-m) /= 0 ) then y2 = 0 ; y1 = 0 do k=nm,m+1,-1 y0 = alpha(k,m,mu) * y1 + beta(k+1,m)*y2 + wa_data(l_nm(k,-m),:) y2 = y1 ; y1 = y0 enddo AnmPnm =(wa_data(l_nm(m,-m),:) + beta(m+1,m)*y2 & + mu*sqrt(2.0D0*m+3)*y1 ) * Pmm(m,mu) gdata = gdata - AnmPnm * 2 * sin(m*alon) endif end do CALL MPI_ALLREDUCE(gdata,a_Interpolate_array00_wa,size(wa_data,2),& MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR) end function a_Interpolate_array00_wa !--------------- 下部ルーチン ----------------- function alpha_array0(n,m,x) ! ! 漸化式の P_n^m の係数 ! integer, intent(IN) :: n,m real(8), intent(IN) :: x real(8) :: alpha_array0 alpha_array0 = sqrt( (2.0D0*n+3)*(2.0D0*n+1)/((n-m+1)*(n+m+1)) ) * x end function alpha_array0 function beta(n,m) ! ! 漸化式の P_{n-1}^m の係数 ! integer, intent(IN) :: n,m real(8) :: beta beta = - sqrt( (2.0D0*n+3)*(n+m)*(n-m)/((2*n-1)*(n+m+1)*(n-m+1)) ) end function beta function Pmm_array0(m,x) ! ! ルジャンドル陪函数 P_m^m(x) の計算 ! sqrt( factrl(2*m+1) )/(2.0D0**m * factrl(m)) & ! * (1-x**2)**(m/2.0D0) ! ! 係数が発散しないように対数で計算する ! integer, intent(IN) :: m ! ルジャンドル陪函数の次数 real(8), intent(IN) :: x ! 引き数 real(8) :: Pmm_array0 ! ルジャンドル陪函数の値 real(8) :: gammaln ! Γ関数の対数 if ( m < 0 ) call MessageNotify('E','Pmm in wa_Intepolate_module',& 'order m should be larger equal to zero') Pmm_array0 = 1.0 if ( m > 0 )then Pmm_array0 = & exp(gammaln(2.0D0*m+2.0)/2.0D0 - m*log(2.0D0) - gammaln(m+1.0D0)) & * (1-x**2)**(m/2.0D0) endif end function Pmm_array0 end module wa_interpolate_mpi_module_supack