! 内部重力波の中層大気伝搬の数値シミュレーション(f90)
! 境界条件:   上端: すべてゼロ(粘性によるスポンジ層を設定)
!            下端: u, Tは鉛直微分ゼロ, w=0
!         水平方向:周期境界条件 
! 
! 空間差分: 中央差分
! 時間差分: Runge-Kutta-Gill 法
! 2/26: 連続の式,x方向の運動方程式 の空間差分の式作成
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! プログラム主要部分
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
program main

  use gtool_history
! 物理定数: g, c_v, R, p_s
  use phys_const
  ! 時間空間解像度, 空間サイズ: grd_x, grd_z, L_x, L_z, dx,dz, dt
  use CalcParameter
  ! 境界条件
  use boundary_condition
  ! 初期値出力
  use Goya_Miyahara_basic_state
  ! 空間微分
  use gravity_wave_eq
  ! 時間発展
  use leap_frog
  implicit none

  !書き出し用の次元配列
  real(8), dimension(1:grd_x,1:grd_z) :: x
  real(8), dimension(1:grd_x,1:grd_z) :: z
  real(8), dimension(1:grd_x,1:grd_z) :: dw            !とりあえず鉛直速度
  real(8), dimension(1:grd_x,1:grd_z) :: du            !とりあえず鉛直速度
  real(8), dimension(1:grd_x,1:grd_z) :: temp          !とりあえず温度？ 
  real(8), dimension(1:grd_x,1:grd_z) :: dens
  real(8), dimension(1:grd_x,1:grd_z) :: pres
  real(8), dimension(1:grd_x,1:grd_z) :: theta
  real(8), dimension(1:grd_x,1:grd_z) :: theta_B
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! 予報変数
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  real(8) :: &
       & rho(-1:grd_x+2,0:grd_z+1),    &    !方程式の解(密度)
       & drho_u(-1:grd_x+2,0:grd_z+1), &    !方程式の解(密度*x速度)
       & rho_w(-1:grd_x+2,0:grd_z+1),  &    !方程式の解(密度*z速度)
       & drho_T(-1:grd_x+2,0:grd_z+1), &    !方程式の解(密度*温度)
       & p(-1:grd_x+2,0:grd_z+1),      &    !圧力
       & drho(-1:grd_x+2,0:grd_z+1),   &   
       & dp(-1:grd_x+2,0:grd_z+1)          

  real(8) :: &
       & rho_dt(1:grd_x, 1:grd_z),      &
       & drho_u_dt(1:grd_x, 1:grd_z),   &
       & drho_w_dt(1:grd_x, 1:grd_z), &
       & drho_T_dt(1:grd_x, 1:grd_z)

  real(8) :: &  ! leap flog のための退避場所
       & rho_b1(-1:grd_x+2,0:grd_z+1),    &    !方程式の解(密度)
       & drho_u_b1(-1:grd_x+2,0:grd_z+1), &    !方程式の解(密度*x速度)
       & rho_w_b1(-1:grd_x+2,0:grd_z+1),  &    !方程式の解(密度*z速度)
       & drho_T_b1(-1:grd_x+2,0:grd_z+1), &       !方程式の解(密度*温度)
!
       & rho_b2(-1:grd_x+2,0:grd_z+1),    &    !方程式の解(密度)
       & drho_u_b2(-1:grd_x+2,0:grd_z+1), &    !方程式の解(密度*x速度)
       & rho_w_b2(-1:grd_x+2,0:grd_z+1),  &    !方程式の解(密度*z速度)
       & drho_T_b2(-1:grd_x+2,0:grd_z+1)       !方程式の解(密度*温度)

  
  real(8) :: tim              !時間
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!初期値
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  real(8),dimension(0:grd_z+1) :: &
       & u_B, T_B, rho_B, p_B, rho_u_B, rho_T_B
  real(8) :: &
       & K_h(1:grd_z),  &
       & K_vu(1:grd_z), &
       & K_vw(1:grd_z), &
       & Rf(1:grd_z)
  ! real(8) :: & 
  !      & u_B(-1:grd_x+2,0:grd_z+1),    &
  !      & T_B(-1:grd_x+2,0:grd_z+1),    &
  !      & rho_B(-1:grd_x+2,0:grd_z+1),  &
  !      & p_B(-1:grd_x+2,0:grd_z+1),    &
  !      & rho_u_B(-1:grd_x+2,0:grd_z+1),&
  !      & rho_T_B(-1:grd_x+2,0:grd_z+1)
       
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! その他の作業領域
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  integer :: i, j,  ii, ifn

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! 基本場と初期値
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  ! 基本場
  tim = dble(0)
  call basic_state( &
       & u_B, T_B, p_B, rho_B, K_h, K_vu, K_vw, Rf)

  
  call basic_boundary(u_B)
  call basic_boundary(T_B)
  call basic_boundary(rho_B)
 
  do j = 1,grd_z
     rho_u_B(j) = u_B(j)*(rho_B(j)+rho_B(j+1))*0.5d+0
     rho_T_B(j) = T_B(j)*rho_B(j)
     p_B(j)     = rho_T_B(j)*R
  end do
 
  call basic_boundary(rho_u_B)
  call basic_boundary(rho_T_B)
  call basic_boundary(p_B)

  ! 初期値設定
  drho(-1:grd_x+2,0:grd_z+1)   = dble(0)
  drho_u(-1:grd_x+2,0:grd_z+1) = dble(0)
  rho_w(-1:grd_x+2,0:grd_z+1)  = dble(0) 
  drho_T(-1:grd_x+2,0:grd_z+1) = dble(0)
  dp(-1:grd_x+2,0:grd_z+1)     = dble(0)

!  基本場の生成
  do j = 1, grd_z
     do i = 1, grd_x
        rho(i,j)     = rho_B(j)
        p(i,j)       = p_B(j)
     end do
  end do
  
  call p_boundary(rho,grd_z+1)
  call p_boundary(p,grd_z+1)


  

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !RKGの値の初期化
 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
! 出力用の値を設定
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  do i=1,grd_x
     do j=1,grd_z
! !        dw(i,j)   = rho_w(i,j-1)
!         dw(i,j)   = dble(0)
!         du(i,j)   = drho_u(i,j)
!         temp(i,j) = T_B(j)
!         pres(i,j) = p(i,j)
!         dens(i,j) = rho(i,j)

        !        dw(i,j)   = rho_w(i,j-1)
        dw(i,j)   = rho_w(i,j-1)
        du(i,j)   = drho_u(i,j)
        temp(i,j) = T_B(j)
        pres(i,j) = p(i,j)
        dens(i,j) = rho(i,j)
        theta(i,j)= T_B(j)*((p_s/p(i,j))**(R/(c_v+R)))
        theta_B(i,j) = theta(i,j)
     end do
  end do

  ! temp = 0.0d+0
  ! pres = 0.0d+0
  ! dens = 0.0d+0
!  theta = 0.0d+0
  
  call dw_boundary(dw,grd_z)  
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! gtool5 ヒストリ作成
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ! 書き出し用の次元変数出力
  do i = 1, grd_x
     x(i,1:grd_z) = (i-1)*dx
  end do
  do i = 1, grd_z
     z(1:grd_x,i) = (i-1)*dz
  end do
  call HistoryCreate( &                        ! ヒストリー作成
     & file='data_2d_gravity_wave.nc', title='gravity_wave',             &
     & source='2d_gravity_wave.f90 gtool5 test',                         &
     & institution='GFD_Dennou Club davis project',                      &
     & dims=(/'x','z','t'/), dimsizes=(/grd_x,grd_z,0/),               & 
     & longnames=(/'X-coordinate','Z-coordinate','time        '/),       &
     & units=(/'m','m','s'/),                                            &
     & origin=real(tim), interval=real(Intv*dt) )

  call HistoryPut('x',x(1:grd_x,1))                       ! 次元変数出力
  call HistoryPut('z',z(1,1:grd_z))                       ! 次元変数出力
  
  call HistoryAddVariable( &                   ! 変数定義
       & varname='dw', dims=(/'x','z','t'/), &
       & longname='vertical velocity', units='m/s', xtype='double')

  call HistoryAddVariable( &                   ! 変数定義
       & varname='du', dims=(/'x','z','t'/), &
       & longname='horizontal velocity', units='m/s', xtype='double')
  
  call HistoryAddVariable( &                   ! 変数定義
       & varname='temp', dims=(/'x','z','t'/), &
       & longname='Temperautre', units='K', xtype='double')
  
  call HistoryAddVariable( &                   ! 変数定義
       & varname='dens', dims=(/'x','z','t'/), &
       & longname='density', units='kg/m^3', xtype='double')

  call HistoryAddVariable( &                   ! 変数定義
       & varname='pres', dims=(/'x','z','t'/), &
       & longname='pressure', units='Pa', xtype='double')

  call HistoryAddVariable( &                   ! 変数定義
       & varname='theta', dims=(/'x','z','t'/), &
       & longname='potential temperature', units='K', xtype='double')

  
  call HistoryPut('du',du)
  call HistoryPut('dw',dw)
  call HistoryPut('temp',temp)                 !変数(の初期値)出力
  call HistoryPut('dens',dens)
  call HistoryPut('pres',pres)
  call HistoryPut('theta',theta)

  do i=-1,grd_x+2
     if (mod(i,grd_x) == grd_x/2) then
        rho_w(i,1)  = PrimeDist
        ! rho_w(i,1)  = 1.0d-5
     end if
  end do
  
  do IFN=1, idmax
     
! RKG 1step %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     call gravity_wave( &
          & rho, drho_u, rho_w, drho_T,          &
          & rho_B, p_B, rho_u_B, rho_T_B, p,     &
          & K_h, K_vu, K_vw, Rf,                 &
          & rho_dt, drho_u_dt, drho_w_dt, drho_T_dt )
     
     call leap_frog_2D_filter(rho, rho_b1, rho_b2, rho_dt, grd_x, grd_z, ifn)
     call p_boundary(rho,grd_z+1)

     call leap_frog_2D_filter(drho_u, drho_u_b1,drho_u_b2, drho_u_dt, grd_x, grd_z, ifn)
     call uT_boundary(drho_u,grd_z+1)   

     call leap_frog_2D_filter(rho_w, rho_w_b1,rho_w_b2, drho_w_dt, grd_x, grd_z, ifn)
     call dw_boundary(rho_w,grd_z+1)

     call leap_frog_2D_filter(drho_T, drho_T_b1,drho_T_b2, drho_T_dt, grd_x, grd_z, ifn)
     call uT_boundary(drho_T,grd_z+1)
        
     do j = 1, grd_z
        do i = 1, grd_x
           p(i,j) = R*(drho_T(i,j)+rho_T_B(j))
        end do
     end do
     call p_boundary(p,grd_z+1)

100 format (1x,e12.5,1x,e12.5,1x,e12.5,1x,e12.5)
        
        
     tim = tim + dt ! 時間ステップの積算

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! 結果出力
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ! write (10+IFN,100) &
        !      (real(II)*dx, tim, real(rho(II)), II=1,grd)
                
!        call HistoryPut('dw',dw)
                
!        call grobal(Gu, grd, total_mass)
!        write (60, 200) tim, total_mass
     
     if ( mod(IFN,INTV) == 0 ) then 
        do j = 1, grd_z
           do i = 1, grd_x
              ! du(i,j)   = dble(2)*(drho_u(i,j)+rho_u_B(i,j)) &
              !      & /(rho(i,j)+rho(i-1,j)) &
              !      & -u_B(i,j)
              du(i,j)   = dble(2)*drho_u_b2(i,j)/(rho_b2(i+1,j)+rho_b2(i,j))
              dw(i,j)   = dble(2)*rho_w_b2(i,j-1)/(rho_b2(i,j)+rho_b2(i,j-1))
              !temp(i,j) =(drho_T(i,j)+rho_T_B(j))/rho(i,j)
              temp(i,j) =(drho_T_b2(i,j)+rho_T_B(j))/rho_b2(i,j)
              dens(i,j) =rho_b2(i,j)
              pres(i,j) =p(i,j)
              theta(i,j)= ((drho_T_b2(i,j)+rho_T_B(j))/rho_b2(i,j))*((p_s/p(i,j))**(R/(c_v+R)))
           end do
        end do

        ! dens(1:grd_x,1:grd_z) = rho(1:grd_x,1:grd_z)
        ! pres(1:grd_x,1:grd_z) = p(1:grd_x,1:grd_z)

     ! write(6,*) "kekka rho_w,drho_u,rho"
     ! do i=1,grd_x
     !    write(6,*)  rho_w(i,0),drho_u(i,0),rho(i,0), p(i,0)
     ! end do
     
        call HistoryPut('du',du)
        call HistoryPut('dw',dw)
        call HistoryPut('temp',temp)
        call HistoryPut('dens',dens)
        call HistoryPut('pres',pres)
        call HistoryPut('theta',theta)
     end if
  end do

  call HistoryClose
  STOP
end program main
