Class | at_module |
In: |
src/at_module.f90
|
a_Int_ag : | real(8), dimension(size(ag,1)) |
ag : | real(8), dimension(:,0:), intent(in) |
————— 積分計算 ——————
function a_Int_ag(ag) real(8), dimension(:,0:), intent(in) :: ag real(8), dimension(size(ag,1)) :: a_Int_ag integer :: i if ( size(ag,2) < im+1 ) then call MessageNotify('E','ae_ag', 'The Grid points of input data too small.') elseif ( size(ag,2) > im+1 ) then call MessageNotify('W','ae_ag', 'The Grid points of input data too large.') endif a_Int_ag = 0.0d0 do i=0,im a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_X_Weight(i) enddo end function a_Int_ag
ag_at : | double precision, dimension(size(at_data,1),0:im) |
at_data : | double precision, dimension(:,:), intent(in) |
—- 逆変換 —-
function ag_at(at_data) ! スペクトル -> 台形格子(2 次元) double precision, dimension(:,:), intent(in) :: at_data double precision, dimension(size(at_data,1),0:im) :: ag_at real(8), dimension(size(at_data,1)*im) :: y integer :: m m = size(at_data,1) if ( size(at_data,2)-1 < km ) then call MessageNotify('E','ag_at', 'The Chebyshev dimension of input data too small.') elseif ( size(at_data,2)-1 > km ) then call MessageNotify('W','ag_at', 'The Chebyshev dimension of input data too large.') endif ag_at = 0.0D0 ag_at(:,0:km)=at_data call fttctb(m,im,ag_at,y,it,t) end function ag_at
at_Dx_at : | real(8), dimension(size(at_data,1),0:size(at_data,2)-1) |
at_data : | real(8), dimension(:,0:), intent(in) |
—- 微分計算 —-
function at_Dx_at(at_data) real(8), dimension(:,0:), intent(in) :: at_data real(8), dimension(size(at_data,1),0:size(at_data,2)-1) :: at_Dx_at integer :: m, k integer :: nm, kmax nm=size(at_data,1) kmax=size(at_data,2)-1 if ( kmax < km ) then call MessageNotify('W','at_Dx_at', 'The Chebyshev dimension of input data too small.') elseif ( kmax > km ) then call MessageNotify('E','at_Dx_at', 'The Chebyshev dimension of input data too large.') endif if ( kmax == im ) then do m=1,nm at_Dx_at(m,kmax) = 0. at_Dx_at(m,kmax-1) = 2 * km * at_data(m,kmax) /2 enddo else do m=1,nm at_Dx_at(m,kmax) = 0. at_Dx_at(m,kmax-1) = 2 * km * at_data(m,kmax) ! スタートはグリッド対応最大波数未満. Factor 1/2 不要 enddo endif do k=kmax-2,0,-1 do m=1,nm at_Dx_at(m,k) = at_Dx_at(m,k+2) + 2*(k+1)*at_data(m,k+1) enddo enddo do k=0,kmax do m=1,nm at_Dx_at(m,k) = 2/xl * at_Dx_at(m,k) enddo enddo end function at_Dx_at
i : | integer, intent(in) |
k : | integer, intent(in) |
xmin : | real(8), intent(in) |
xmax : | real(8), intent(in) |
—- 初期化 —-
subroutine at_Initial(i,k,xmin,xmax) integer, intent(in) :: i, k real(8), intent(in) :: xmin,xmax integer :: ii,kk im=i km=k xl = xmax-xmin if ( im <= 0 .or. km <= 0 ) then call MessageNotify('E','at_initial', 'Number of grid points and waves should be positive') elseif ( mod(im,2) /= 0 ) then call MessageNotify('E','at_initial','Number of grid points should be even') elseif ( km > im ) then call MessageNotify('E','at_initial','KM shoud be less equal IM') endif allocate(t(3*im)) call fttcti(im,it,t) allocate(g_X(0:im)) do ii=0,im g_X(ii) = (xmax+xmin)/2 + xl/2 * cos(pi*ii/im) enddo allocate(g_X_Weight(0:im)) do ii=0,im g_X_Weight(ii) = 1.0 do kk=2,km,2 g_X_Weight(ii) = g_X_Weight(ii) + 2/(1D0-kk**2) * cos(kk*ii*pi/im) enddo if ( (km == im) .and. (mod(im,2)==0) ) then ! 最後の和は factor 1/2. g_X_Weight(ii) = g_X_Weight(ii) - 1/(1D0-km**2)* cos(km*ii*pi/im) endif g_X_Weight(ii) = 2D0/im * g_X_Weight(ii) * xl/2 enddo g_X_Weight(0) = g_X_Weight(0) / 2 g_X_Weight(im) = g_X_Weight(im) / 2 end subroutine at_Initial