ae_module.f90

Go to the documentation of this file.

! -*- mode: f90; coding: utf-8 -*-
!-----------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!-----------------------------------------------------------------------
module ae_module
  use dc_message, only: messagenotify
  use dc_types, only: dp
  use spml_utils, only: pi
  implicit none

  private
  public ae_initial
  public ag_ae
  public ae_ag
  public g_e
  public e_g
  public ae_dx_ae
  public e_dx_e
  public a_int_ag
  public int_g
  public a_avr_ag
  public avr_g
  public interpolate_e
  public a_interpolate_ae
  public g_x
  public g_x_weight

  real(DP), allocatable    :: g_x(:)
  real(DP), allocatable    :: g_x_weight(:)

  integer(8)               :: im = 32    ! 格子点の数
  integer(8)               :: km = 10    ! 切断波数
  real(DP)                 :: xl = 1.0d0 ! 領域の大きさ
  integer(8), allocatable  :: it(:)
  real(DP),   allocatable  :: t(:)

  save im, km, it, t, xl, g_x, g_x_weight

contains
  !-----------------------------------------------------------------------
  subroutine ae_initial(i,k,xmin,xmax)

    integer, intent(in) :: i
    integer, intent(in) :: k
    real(DP),intent(in) :: xmin
    real(DP),intent(in) :: xmax
    integer(8)          :: ii

    im = i
    km = k
    xl = xmax - xmin

    if ( im <= 0 .or. km <= 0 ) then
       call messagenotify('E','ae_Initial', &
            'Number of grid points and waves should be positive')
    elseif ( mod(im,2_8) /= 0 ) then
       call messagenotify('E','ae_Initial', &
            'Number of grid points should be even')
    elseif ( km > im/2 ) then
       call messagenotify('E','ae_Initial', &
            'KM shoud be less equal to IM/2')
    endif

    allocate(it(im/2),t(im*3/2))

    call fxrini(im,it,t)

    allocate(g_x(0:im-1))
    do ii=0,im-1
       g_x(ii) = xmin + (xl/im)*ii
    enddo

    allocate(g_x_weight(0:im-1))
    g_x_weight = xl/im

    call messagenotify(&
         'M','ae_initial','ae_module (2019/09/17) is initialized')

  end subroutine ae_initial

  !-----------------------------------------------------------------------
  function ag_ae(ae)

    real(DP), dimension(:,-km:), intent(in)  :: ae
    real(DP), dimension(size(ae,1),0:im-1)   :: ag_ae

    integer(8) :: nm, k
    nm=size(ae,1)

    if ( size(ae,2) < 2*km+1 ) then
       call messagenotify('E','ag_ae', &
            'The Fourier dimension of input data too small.')
    elseif ( size(ae,2) > 2*km+1 ) then
       call messagenotify('W','ag_ae', &
            'The Fourier dimension of input data too large.')
    endif

    ag_ae       = 0.0d0
    ag_ae(:,0)  = ae(:,0)
    do k=1,min(km,im/2-1)
       ag_ae(:,2*k)   = ae(:,k)
       ag_ae(:,2*k+1) = ae(:,-k)
    enddo
    if ( km == im/2 ) then
      ag_ae(:,1) = ae(:,km)
    end if

    call fxrtba(nm,im,ag_ae,it,t)

  end function ag_ae

  !-----------------------------------------------------------------------
  function g_e(e)

    real(DP), dimension(-km:km), intent(in) :: e
    real(DP), dimension(0:im-1)             :: g_e

    real(DP), dimension(1,size(e))  :: ae_work
    real(DP), dimension(1,0:im-1)   :: ag_work

    ae_work(1,:) = e
    ag_work = ag_ae(ae_work)
    g_e = ag_work(1,:)

  end function g_e

  !-----------------------------------------------------------------------
  function ae_ag(ag)

    real(DP), dimension(:,:), intent(in)   :: ag
    real(DP), dimension(size(ag,1),-km:km) :: ae_ag

    real(DP), dimension(size(ag,1),0:im-1) :: ag_work
    integer(8) :: nm, k

    nm = size(ag,1)
    if ( size(ag,2) < im ) then
       call messagenotify('E','ae_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call messagenotify('W','ae_ag', &
            'The Grid points of input data too large.')
    endif
    ag_work = ag

    call fxrtfa(nm,im,ag_work,it,t)

    do k=1,min(km,im/2-1)
       ae_ag(:, k) = ag_work(:,2*k)
       ae_ag(:,-k) = ag_work(:,2*k+1)
    enddo
    ae_ag(:,0) = ag_work(:,0)
    if ( km == im/2 ) then
      ae_ag(:,km) = ag_work(:,1)
    end if

  end function ae_ag

  !-----------------------------------------------------------------------
  function e_g(g)

    real(DP), dimension(0:im-1), intent(in)  :: g
    real(DP), dimension(-km:km)              :: e_g

    real(DP), dimension(1,size(g))        :: ag_work
    real(DP), dimension(1,-km:km)         :: ae_work

    ag_work(1,:) = g
    ae_work = ae_ag(ag_work)
    e_g = ae_work(1,:)

  end function e_g

  !-----------------------------------------------------------------------
  function ae_dx_ae(ae)

    real(DP), dimension(:,-km:), intent(in)  :: ae
    real(DP), dimension(size(ae,1),-km:km)   :: ae_dx_ae

    integer(8) ::  k

    if ( size(ae,2) < 2*km+1 ) then
       call messagenotify('W','ae_Dx_ae', &
            'The Fourier dimension of input data too small.')
    elseif ( size(ae,2) > 2*km+1 ) then
       call messagenotify('W','ae_Dx_ae', &
            'The Fourier dimension of input data too large.')
    endif

    do k=-km,km
       ae_dx_ae(:,k) = -(2*pi*k/xl)*ae(:,-k)
    enddo
  end function ae_dx_ae

  !-----------------------------------------------------------------------
  function e_dx_e(e)

    real(DP), dimension(-km:km), intent(in)     :: e
    real(DP), dimension(-km:km)                 :: e_Dx_e

    real(DP), dimension(1,-km:km)               :: ae_work

    ae_work(1,:) = e
    ae_work = ae_dx_ae(ae_work)
    e_dx_e = ae_work(1,:)

  end function e_dx_e

  !-----------------------------------------------------------------------
  function a_int_ag(ag)

    real(DP), dimension(:,0:), intent(in)     :: ag
    real(DP), dimension(size(ag,1))           :: a_Int_ag

    integer(8) :: i

    if ( size(ag,2) < im ) then
       call messagenotify('E','ae_Int_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call messagenotify('W','ae_Int_ag', &
            'The Grid points of input data too large.')
    endif

    a_int_ag = 0.0d0
    do i=0,im-1
       a_int_ag(:) = a_int_ag(:) + ag(:,i)*g_x_weight(i)
    enddo
  end function a_int_ag

  !-----------------------------------------------------------------------
  function int_g(g)

    real(DP), dimension(0:im-1), intent(in)   :: g
    real(DP)                                  :: Int_g

    int_g = sum(g*g_x_weight)

  end function int_g

  !-----------------------------------------------------------------------
  function a_avr_ag(ag)

    real(DP), dimension(:,0:), intent(in)     :: ag
    real(DP), dimension(size(ag,1))           :: a_Avr_ag

    a_avr_ag = a_int_ag(ag)/sum(g_x_weight)

  end function a_avr_ag

  !-----------------------------------------------------------------------
  function avr_g(g)

    real(DP), dimension(0:im-1), intent(in)   :: g
    real(DP)                                  :: Avr_g

    avr_g = int_g(g)/sum(g_x_weight)

  end function avr_g

  !-----------------------------------------------------------------------
  function interpolate_e(e_data,xval)

    real(DP), dimension(-km:km), intent(IN) :: e_data
    real(DP), intent(in)                    :: xval
    real(DP)                                :: Interpolate_e

    integer(8) :: k

    interpolate_e = e_data(0)
    do k=1,km
       interpolate_e = interpolate_e                             &
            &  + 2.0d0 * (   e_data(k)  * cos(2*pi*k/xl*xval)    &
            &              - e_data(-k) * sin(2*pi*k/xl*xval) )
    end do

  end function interpolate_e

  !-----------------------------------------------------------------------
  function a_interpolate_ae(ae_data,xval)

    real(DP), dimension(:,-km:), intent(IN)   :: ae_data
    real(DP), intent(in)                      :: xval
    real(DP), dimension(size(ae_data,1))      :: a_Interpolate_ae

    integer(8) :: k

    a_interpolate_ae = ae_data(:,0)

    do k=1,km
       a_interpolate_ae = a_interpolate_ae                          &
            &  + 2.0d0 * (   ae_data(:,k)  * cos(2*pi*k/xl*xval)    &
            &              - ae_data(:,-k) * sin(2*pi*k/xl*xval) )
    end do

  end function a_interpolate_ae

end module ae_module