! 内部重力波の鉛直速度分布を求める
! 準圧縮系の場合
! 途中で分散関係を解く
! 地形: cos 型

!-- 実行時の確認事項---

!-  * 計算領域幅, グリッド間隔,
!-  * 変数 Temp, N

!----------------------

program analytical_solution

  use gtool_history                             ! モジュール指定

  use fft8_module

  implicit none

  !-- 
  real(8), parameter     :: N    = 2.0d-2       ! 浮力振動数
  real(8), parameter     :: U    = 1.0d+1       ! 平均流の大きさ
!!$  real(8), parameter     :: Cp   = 
  real(8), parameter     :: Grav = 9.8d0
  real(8), parameter     :: GasC   = 8.314d0 / 28.0d-3
  real(8), parameter     :: Temp0 = 231.0d0    ! Grav/(Cp*N**2)
  real(8), parameter     :: Press0 = 965.0d2
  real(8), parameter     :: Dens00 = Press0 / ( GasC * Temp0 )
  real(8), parameter     :: ScaleHeight = GasC * Temp0 / Grav
  !-- 
  !--
  integer                :: nx, nz
!!$  real(8), parameter     :: zrange = 10.0d+3
  real(8), parameter     :: zrange = 10.0d+2
  !
!!$  real(8)            :: mountheight  =  10.0d0
  real(8)            :: mountheight  = 100.0d0
  ! A1
!!$  real(8)            :: mounthalfwid = 5000.0d0
!!$  real(8), parameter     :: xrange = 200.0d+3
!!$  real(8), parameter     :: dx = 100            ! グリッド間隔
!!$  real(8), parameter     :: dz = 100            ! グリッド間隔
  ! A2
!!$  real(8)            :: mounthalfwid = 500.0d0
!!$  real(8), parameter     :: xrange = 200.0d+3
!!$  real(8), parameter     :: dx = 100            ! グリッド間隔
!!$  real(8), parameter     :: dz = 100            ! グリッド間隔
  ! A3
  real(8)            :: mounthalfwid = 100.0d0
  real(8), parameter     :: xrange = 40.0d+3
  real(8), parameter     :: dx = 20            ! グリッド間隔
  real(8), parameter     :: dz = 20            ! グリッド間隔
  ! A4
!!$  real(8)            :: mounthalfwid = 50.0d0
!!$  real(8), parameter     :: xrange = 10.0d+3
!!$  real(8), parameter     :: dx = 5            ! グリッド間隔
!!$  real(8), parameter     :: dz = 5            ! グリッド間隔
  ! D1
!!$  real(8)            :: mounthalfwid = 500.0d0
  !--
  ! D1
!!$  real(8), parameter     :: xrange = 100.0d+3
!!$  real(8), parameter     :: dx = 50            ! グリッド間隔
!!$  real(8), parameter     :: dz = 50            ! グリッド間隔

  real(8), allocatable   :: x(:), z(:)
  real(8)                :: val                 ! 計算のための変数
  integer                :: i, j, l             ! ループのための添え字
  integer                :: lneg
  integer                :: lmin
  integer                :: lmax

  complex(8)             :: m                     ! z 方向の波数
  real(8), parameter     :: pi = acos(-1.0d0)

  !-- 計算する物理量
  real(8), allocatable   :: wave_w        (:, :)             ! 鉛直速度
  real(8), allocatable   :: wave_w_sub    (:, :)
  real(8), allocatable   :: wave_exner    (:, :)
  real(8), allocatable   :: wave_exner_sub(:, :)
  real(8), allocatable   :: wave_u        (:, :)
  real(8), allocatable   :: wave_u_sub    (:, :)
  real(8), allocatable   :: wave_pt       (:, :)
  real(8), allocatable   :: wave_pt_sub   (:, :)
  real(8), allocatable   :: wave_uwbar    (:)


  real(8), allocatable   :: x_hs(:)             ! 鉛直速度

  real(8)            :: mountcenterX
  real(8)            :: dwn
  real(8)            :: mountamp
  real(8)            :: mountwn

  complex(8) :: AmpU
  complex(8) :: AmpW
  real(8)              :: anal_UWbar
  real(8), allocatable :: w_wn        (:)
  real(8), allocatable :: w_anal_UWbar(:)

  real(8), allocatable :: z_Dens0(:)

  real(8), allocatable   :: xz_ur(:,:)
  real(8), allocatable   :: xz_ui(:,:)
  real(8), allocatable   :: xz_wr(:,:)
  real(8), allocatable   :: xz_wi(:,:)
  integer                :: ita( 5, 2 )
  real(8), allocatable   :: dta(:,:)
  real(8), allocatable   :: w_fft_UWbar(:)
  real(8), allocatable   :: fft_UWbar


  !-- グリッド数を計算

  nx = xrange/dx
  nz = zrange/dz

  !-- allocate 文

  allocate( x(nx) )
  allocate( z(nz) )
  allocate( wave_w        (nx,nz) )
  allocate( wave_w_sub    (nx,nz) )
  allocate( wave_exner    (nx,nz) )
  allocate( wave_exner_sub(nx,nz) )
  allocate( wave_u        (nx,nz) )
  allocate( wave_u_sub    (nx,nz) )
  allocate( wave_pt       (nx,nz) )
  allocate( wave_pt_sub   (nx,nz) )
  allocate( wave_uwbar       (nz) )
  allocate( z_Dens0          (nz) )

  allocate( x_hs(nx) )

  !--------------------------------------


  !-- 絵を描くための格子を作る
  !--

  do i = 1, nx
!!$     x(i) = dx*(i-1) + dx/2 !!数値解の出力と格子点をそろえるため
     x(i) = dx*(i-1)
  end do

  do j = 1, nz
!!$     z(j) = dz*(j-1) + dz/2 !!数値解の出力と格子点をそろえるため
     z(j) = dz*(j-1)
  end do

  !--------------------------------------

  mountcenterX = dx * dble( nx/2+1 - 1 )

  do i = 1, nx
     x_hs(i) = mountheight * mounthalfwid**2 / ( mounthalfwid**2 + ( X(i) - mountcenterX )**2 )
  end do

!!$  do i = 1, nx
!!$     write( 60, * ) X(i), x_hs(i)
!!$  end do



  z_Dens0 = Dens00 * exp( - z(:) / ScaleHeight )


  ! !-- 準圧縮系での鉛直速度の解析解
  ! !--

  wave_w     = 0.0d0
  wave_exner = 0.0d0
  wave_u     = 0.0d0
  wave_pt    = 0.0d0


  x_hs = 0.0d0


!!$  dwn = 5.0d-7
!!$  ! parameters for dwn = 5e-7
!!$  ! propagating : k = 0 - 1.999e-3
!!$  lmin = 1 ; lmax = 3999

!!$  dwn = 2.0d-7
!!$  ! parameters for dwn = 2e-7
!!$  ! propagating : k = 0 - 1.999e-3
!!$  lmin = 1 ; lmax = 4999


  dwn = 1.0d-6

  ! N = 2e-2, U = 10.0 === > K = 2e-3
  !
  ! parameters for dwn = 1e-6
  ! full : k = 0 - 1e-2
!!$  lmin = 1 ; lmax = 10000
  ! propagating : k = 0 - 1.999e-3
!  lmin = 1 ; lmax = 1999
!  lmin = 1 ; lmax = 2200
!  lmin = 1 ; lmax = 5000
  lmin = 1 ; lmax = 10000
  ! evanescent : k = 2e-3 - 1e-2
!!$  lmin = 2000 ; lmax = 10000
  !
  ! propagating : k = 1.733e-3 - 1e-2 : theta =  - deg.
!!$  lmin = 1466 ; lmax = 1732
  ! propagating : k = 1.733e-3 - 1e-2 : theta =  - deg.
!!$  lmin = 1733 ; lmax = 1999
  ! propagating : k = 0 - 1.732e-3 : theta = 30 - 90 deg. ???
!!$  lmin = 1 ; lmax = 1732
  ! propagating : k = 0 - 0.999e-3 : theta = 60 - 90 deg. ???
!!$  lmin = 1 ; lmax =  999
  !
  anal_UWbar = 0.0d0
  allocate( w_wn        (lmin:lmax) )
  allocate( w_anal_UWbar(lmin:lmax) )
  !
  do l = lmin, lmax
    write( 6, * ) lmin, l, lmax
    mountwn     = dwn * dble(l-1)
    mountamp    = mountheight*mounthalfwid**2 * pi/mounthalfwid / pi * exp( - mounthalfwid * mountwn )
    w_wn(l) = mountwn
    do i = 1, nx
      x_hs(i) = x_hs(i) &
        & + mountamp * cos( mountwn * ( X(i) - mountcenterX ) ) &
        &   * dwn
    end do
!!$    if ( l > 1 ) then
!!$      write( 62, * ) mountwn, 2.0d0*pi/mountwn, mountamp
!!$    end if
    if ( l > 1 ) then
     call gw( &
        & N, U, &
        & nx, nz, &
        & x, z, &
        & mountamp, 2.0d0*pi/mountwn, mountcenterx, &
        & wave_w_sub, wave_exner_sub, wave_u_sub, wave_pt_sub, &
        & AmpU, AmpW &
        & )
     wave_w     = wave_w     + wave_w_sub     * dwn
     wave_exner = wave_exner + wave_exner_sub * dwn
     wave_u     = wave_u     + wave_u_sub     * dwn
     wave_pt    = wave_pt    + wave_pt_sub    * dwn

     w_anal_UWbar(l) = ( AmpU * conjg(AmpW) + conjg(AmpU) * AmpW ) / 4.0d0
    else
     w_anal_UWbar(l) = 0.0d0
    end if
  end do
  w_anal_UWbar = w_anal_UWbar * dwn*dwn * 2.0d0*pi/dwn / xrange
  w_anal_UWbar = w_anal_UWbar * Dens00

  anal_UWbar = 0.0d0
  do l = lmin, lmax
    anal_UWbar = anal_UWbar + w_anal_UWbar(l)
  end do
  write( 6, * ) '1 : UWbar = ', anal_UWbar
  open( 71, file = 'momfluxspe_analytic.txt', status = 'unknown' )
  write( 71, * ) '# mountheight = ', mountheight
  write( 71, * ) '# <index>   <wn (m-1)>   <uw spe. dens. (kg s-1)>'
  do l = lmin, lmax
    mountwn     = dwn * dble(l-1)
    write( 71, * ) l, mountwn, w_anal_UWbar(l) / dwn
!!$    write( 71, * ) l, mountwn, w_anal_UWbar(l)
  end do
  close( 71 )


  !------
  wave_uwbar = 0.0d0
  do j = 1, nz
     do i = 1, nx
        wave_uwbar(j) = wave_uwbar(j) + wave_u(i,j) * wave_w(i,j) / dble(nx)
     end do
  end do
  wave_uwbar(:) = wave_uwbar(:) * z_Dens0
  !-----

  wave_u = wave_u + U

  !-----


  allocate( xz_ur(nx,nz) )
  allocate( xz_ui(nx,nz) )
  allocate( xz_wr(nx,nz) )
  allocate( xz_wi(nx,nz) )
  xz_ur = wave_u
  xz_ui = 0.0d0
  xz_wr = wave_w
  xz_wi = 0.0d0
  allocate( dta( nx, 2 ) )
  call fft8_fft_init( nx, ita, dta )
  call fft8_fft( nx, nz, ita, dta, xz_ur, xz_ui, 1 )
  call fft8_fft( nx, nz, ita, dta, xz_wr, xz_wi, 1 )

  allocate( w_fft_UWbar(nx) )
  j = 1
  l = 1
  w_fft_UWbar(l) = 0.0d0
  do l = 1+1, nx
    lneg = nx-(l-2)
    w_fft_UWbar(l) = xz_ur(l,j) * xz_wr(lneg,j) - xz_ui(l,j) * xz_wi(lneg,j)
  end do
  w_fft_UWbar = w_fft_UWbar * z_Dens0(j)
!!$  w_fft_UWbar = w_fft_UWbar &
!!$    & / ( (2.0d0*pi/xrange)*(2.0d0*pi/xrange) * 2.0d0*pi/(2.0d0*pi/xrange) / xrange )
  fft_UWbar = 0.0d0
  do l = 1, nx
!!$    fft_UWbar = fft_UWbar &
!!$      & + w_fft_UWbar(l) &
!!$      &   * ( (2.0d0*pi/xrange)*(2.0d0*pi/xrange) * 2.0d0*pi/(2.0d0*pi/xrange) / xrange )
    fft_UWbar = fft_UWbar + w_fft_UWbar(l)
  end do
  write( 6, * ) '2 : UWbar = ', fft_UWbar
  do l = 1, nx
    mountwn     = (2.0d0*pi/xrange)*(l-1)
    write( 72, * ) l, mountwn, w_fft_UWbar(l) / (2.0d0*pi/xrange)
  end do


!!$  anal_UWbar = 0.0d0
!!$  fft_UWbar  = 0.0d0
!!$  do l = 1, 1999
!!$!    write( 6, * ) l, w_anal_UWbar(l), w_fft_UWbar(l)
!!$    anal_UWbar = anal_UWbar + w_anal_UWbar(l)
!!$    fft_UWbar  = fft_UWbar  + w_fft_UWbar(l)
!!$    write( 6, * ) l, anal_UWbar, fft_UWbar
!!$  end do
!!$  stop

!!$  do i = 1, nx
!!$     write( 61, * ) x(i), x_hs(i)
!!$  end do


  !-- 絵を描くための準備
  !--
  call HistoryCreate( &                             ! ヒストリー作成
    & file='analytical_solution.nc', title='Diffusion equation', &
    & source='Sample program of gtool_history/gtool5',   &
    & institution='GFD_Dennou Club davis project',       &
    & dims=(/'x ','z ', 'wn'/), dimsizes=(/nx,nz,(lmax-lmin+1)/), &
    & longnames=(/'X-coordinate','Z-coordinate', 'wavenumber  '/), &
    & units=(/'m  ','m  ', 'm-1'/) )

  call HistoryPut('x' ,x)                           ! 次元変数出力
  call HistoryPut('z' ,z)                           ! 次元変数出力
  call HistoryPut('wn',w_wn)                        ! 次元変数出力
 
  call HistoryAddVariable( &                       ! 変数定義
    & varname='VelZ', dims=(/'x','z'/), &
    & longname='vertical velocity', units='m/s', xtype='double')
  call HistoryAddVariable( &                       ! 変数定義
    & varname='Exner', dims=(/'x','z'/), &
    & longname='Exner function', units='1', xtype='double')
  call HistoryAddVariable( &                       ! 変数定義
    & varname='VelX', dims=(/'x','z'/), &
    & longname='x velocity', units='m/s', xtype='double')
  call HistoryAddVariable( &                       ! 変数定義
    & varname='PTemp', dims=(/'x','z'/), &
    & longname='potential temperature', units='K', xtype='double')
  call HistoryAddVariable( &                       ! 変数定義
    & varname='MfluxSpe', dims=(/'wn'/), &
    & longname='momentum flux spectrum', units='kg s-1', xtype='double')
  call HistoryAddVariable( &                       ! 変数定義
    & varname='Mflux', dims=(/''/), &
    & longname='momentum flux', units='kg m-1 s-1', xtype='double')

  call HistoryPut('VelZ' ,wave_w    )
  call HistoryPut('Exner',wave_exner)
  call HistoryPut('VelX' ,wave_u    )
  call HistoryPut('PTemp',wave_pt   )
  call HistoryPut('MfluxSpe',w_anal_UWbar/dwn )
  call HistoryPut('Mflux',anal_UWbar )

  call HistoryClose

  write( 6, * ) 'UWbar = ', anal_UWbar
!!$  do j = 1, nz
  do j = 1, 5
     write( 6, * ) j, wave_uwbar(j)
  end do

  stop
end program analytical_solution


subroutine gw( &
     & N, U, &
     & nx, nz, &
     & x, z, &
     & h_0, lambda, mountcenterx, &
     & wave_w, wave_exner, wave_u, wave_pt, &
     & AmpU, AmpW &
     & )

  implicit none

  real(8), intent(in ) :: N
  real(8), intent(in ) :: U
  integer, intent(in ) :: nx, nz
  real(8), intent(in ) :: x(nx), z(nz)
  real(8), intent(in ) :: h_0
  real(8), intent(in ) :: lambda
  real(8), intent(in ) :: mountcenterx
  real(8), intent(out) :: wave_w    (nx, nz)
  real(8), intent(out) :: wave_exner(nx, nz)
  real(8), intent(out) :: wave_u    (nx, nz)
  real(8), intent(out) :: wave_pt   (nx, nz)
  complex(8), intent(out) :: AmpU
  complex(8), intent(out) :: AmpW


  real(8)                :: val                 ! 計算のための変数
  integer                :: i, j, l             ! ループのための添え字
  complex(8)             :: c = (0.0,1.0)       ! 虚数 i=c=(実部, 虚部)

  real(8)                :: k                     ! x 方向の波数
  complex(8)             :: m                     ! z 方向の波数
  real(8), parameter     :: pi = acos(-1.0d0)

  !-- 物理定数
  real(8), parameter     :: Grav     = 9.8d0    ! 重力加速度
  real(8), parameter     :: GasRdry  = 296.7d0  ! 乾燥大気の気体定数
  real(8), parameter     :: GasCpdry = 1039.6d0 ! 乾燥大気の定圧比熱
  real(8), parameter     :: GasCvdry = 742.9d0  ! 乾燥大気の定積比熱
  real(8)                :: H                   ! 高さのスケールハイト
  real(8)                :: Hstar               ! 高さのスケールハイト
  real(8)                :: c_s                 ! 音速

  ! 実験設定によって変える量
  !real(8), parameter     :: Temp = 300.0d0      ! 大気の温度(等温大気を仮定)
  real(8)                :: Temp

  complex(8)                :: AmpExner
  complex(8)                :: AmpPTemp
  complex(8)                :: Wave

  !-- 計算する物理量

  ! 温度を定義
  Temp = Grav**2/(GasCpdry*N**2)

  ! スケールハイトを定義
  H     = Temp*GasRdry /Grav
  Hstar = Temp*GasCpdry/Grav

  ! 音速を定義
  c_s = sqrt( GasCpdry*GasRdry*Temp/GasCvdry )

  ! 水平波数を定義
  k = 2*pi/lambda


  !--------------------------------------

  !-- 分散関係式を解く
  !--

  val=(N/U)**2-(k)**2-1/(2*H)**2 + (U*k)**2/c_s**2
  if (val >= 0) then
     val = sqrt(val)
     m = cmplx(val,0.0)
  else
     val = sqrt(-val)
     m = cmplx(0.0,val)
  end if

  print *, k, m

  ! !--------------------------------------

  ! !-- 準圧縮系での鉛直速度の解析解
  ! !--

  AmpW = U*h_0*c*k
  AmpPTemp = -c/(-u*k)*Temp/Hstar &
       & * AmpW
  AmpExner = &
       & c/(GasCpDry*Temp*(c*m+(-1.0d0/Hstar+1.0d0/(2.0d0*H)))) &
       & * ( -U*k-Grav/Hstar/(-U*k) ) &
       & * AmpW
!!$  AmpU = &
!!$       & k*GasCpDry*Temp/(-U*k) &
!!$       & * AmpExner

  AmpU = &
       & - k &
       & * 1.0d0 / ( c*m - ( 1.0d0/Hstar - 1.0d0/(2.0d0*H) ) ) &
       & * ( U*k + Grav/Hstar/(-U*k) ) &
       & * h_0
  AmpU = &
       & - c_s**2 &
       & * ( c*m + 1.0d0/Hstar - 1.0d0/(2.0d0*H) ) &
       & / ( - U + c_s**2/U ) &
       & * h_0

  do j = 1, nz
     do i = 1, nx
        Wave = exp( c*( k*(x(i) - mountcenterx )+m*z(j) ) ) * exp(z(j)/(2*H))
        wave_w    (i, j) = AmpW     * Wave
        wave_exner(i, j) = AmpExner * Wave
        wave_u    (i, j) = AmpU     * Wave
        wave_pt   (i, j) = AmpPTemp * Wave
     end do
  end do


end subroutine gw
