module gravity_wave_eq !空間差分
  use phys_const
  use CalcParameter
  use boundary_condition
  use differential
  implicit none
contains
  subroutine 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 )
      
  real(8), intent(in) :: & ! 使用する値
       & rho(-1:grd_x+2,0:grd_z+1),      &
       & drho_u(-1:grd_x+2,0:grd_z+1),   &
       & rho_w(-1:grd_x+2,0:grd_z+1),    &
       & drho_T(-1:grd_x+2,0:grd_z+1),   &
       & p(-1:grd_x+2,0:grd_z+1),        &
       & K_h(1:grd_z),                   &
       & K_vu(1:grd_z),                  &
       & K_vw(1:grd_z),                  &
       & Rf(1:grd_z)

  real(8), intent(in) :: & ! 基本場の量
       & rho_B(0:grd_z+1),    &
       & p_B(0:grd_z+1),      &
       & rho_u_B(0:grd_z+1),  &
       & rho_T_B(0:grd_z+1)

  
  real(8), intent(out) :: & !出力結果
       & 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)

  integer :: i, j
  !計算用の値(hが1/2の層, fが密度と同じ層)
  real(8) :: &
       & rho_u_hf(-1:grd_x+2,0:grd_z+1), &
       & rho_T(-1:grd_x+2,0:grd_z+1),    &
       & u_hh(-1:grd_x+2,0:grd_z+2),     &
       & w_hh(-1:grd_x+2,0:grd_z+2),     &
       & T(-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),       &
       ! 粘性項
       & drho_u_dx4(1:grd_x,1:grd_z),   &
       & drho_u_dz2(1:grd_x,1:grd_z),   &
       & drho_w_dx4(1:grd_x,1:grd_z+1), &
       & drho_w_dz2(1:grd_x,1:grd_z+1), &
       & drho_T_dx4(1:grd_x,1:grd_z),   &
       & drho_T_dz2(1:grd_x,1:grd_z)
       

  !---------------------------------------------------------------------
  do j = 1, grd_z
     do i = 1, grd_x
        rho_u_hf(i,j) = drho_u(i,j)+rho_u_B(j)
     end do
  end do
  
  ! do j = 1, grd_z+1
  !    do i = 1, grd_x
  !       rho_w_fh(i,j) = rho_w(i,j)
  !    end do
  ! end do

  !境界条件
  call uT_boundary(rho_u_hf,grd_z+1)
  
  ! call cyclic_boundary(rho_w_fh)
  ! call dw_boundary(rho_w_fh)

  do j = 1, grd_z
     do i = 1, grd_x
        u_hh(i,j) = &
             & + rho_u_hf(i,j+1)/(rho(i+1,j+1)+rho(i,j+1)) &
             & + rho_u_hf(i,j)/(rho(i+1,j)+rho(i,j))
! 
        ! w_hh(i,j) = &
        !      &  rho_w_fh(i,j)/(rho(i,j)+rho(i,j-1)) &
        !      & +rho_w_fh(i-1,j)/(rho(i-1,j)+rho(i-1,j-1))
        w_hh(i,j) = &
             & + rho_w(i+1,j)/(rho(i+1,j+1)+rho(i+1,j)) &
             & + rho_w(i,j)/(rho(i,j+1)+rho(i,j))
     end do
  end do

  do j = 1, grd_z
     do i = 1, grd_x
        rho_T(i,j) = drho_T(i,j)+rho_T_B(j)
        T(i,j)     = rho_T(i,j)/rho(i,j)
     end do   
  end do
  ! 境界条件
  call p_boundary(rho_T,grd_z+1)
  call uT_boundary(u_hh,grd_z+1)
  call dw_boundary(w_hh,grd_z+1)
  call p_boundary(T,grd_z+1)


!rho, p の差分

!   write(6,*) "p,p_B,dp,drho"

  do j = 0, grd_z+1
     do i = -1, grd_x+2
        drho(i,j) = rho(i,j)-rho_B(j)
        dp(i,j)   = p(i,j)-p_B(j)
     end do
!      write(6,*) p(1,j), p_B(j), dp(1,j), drho(1,j)
  end do


  
  ! 粘性項
  !  call four_differential(drho_u, drho_u_dx4)
  call two_differential_x(drho_u, drho_u_dx4)
  call two_differential_z(drho_u, drho_u_dz2)
  !  call four_differential(rho_w, drho_w_dx4)
  call two_differential_x(rho_w, drho_w_dx4)
  call two_differential_z(rho_w, drho_w_dz2)
  !  call four_differential(drho_T, drho_T_dx4)
  call two_differential_x(drho_T, drho_T_dx4)
  call two_differential_z(drho_T, drho_T_dz2)

  !--------------------------------------------------------------------------
100 format (1x,e12.5,1x,e12.5,1x,e12.5,1x,e12.5)
  ! write(6,*) "initial condition"
  ! write(6,*) "rho_u, drho_u, rho_u_B"
  !  do i = 0, grd_z+1
  !  ! !           write(6,100) rho(4,i), drho_u(4,i), rho_w(4,i), drho_T(4,i)
  !     write(6,100)         rho_u_hf(i,j),drho_u(i,j),rho_u_B(i,j)

  !  end do
  
  ! do i = 0, grd_x+2
  !    write(6,100) rho(i,0), rho(i,1)
  ! end do
  

  
  
  ! call boundary condition
  do j=1, grd_z
     do i=1, grd_x
        !continuity equation-----------------------------------------------
        rho_dt(i,j) =             &
             &  (                 &
             &    - (rho_u_hf(i,j) - rho_u_hf(i-1,j))/dx &
             &    - (rho_w(i,j) - rho_w(i,j-1))/dz &
             &  )
        
        !thermodynamic eq--------------------------------------------------
        drho_T_dt(i,j) =                           &
             &  (                                  &
             &    - (                      &
             &        rho_u_hf(i,j)                &
             &        *(T(i+1,j)+T(i,j))*dble(0.5) &
             &       -rho_u_hf(i-1,j)              &
             &        *(T(i,j)+T(i-1,j))*dble(0.5) &
             &      )/dx                           &
             &    - (                      &
             &        rho_w(i,j)                   &
             &        *(T(i,j+1)+T(i,j))*dble(0.5) &
             &       -rho_w(i,j-1)                 &
             &        *(T(i,j)+T(i,j-1))*dble(0.5) &
             &      )/dz                           &
             &    - (R/c_v)*rho_T(i,j)                 &
             &      * (                                &
             &          dble(2)*( &
             &            rho_u_hf(i,j)/(rho(i+1,j)+rho(i,j))   &
             &           -rho_u_hf(i-1,j)/(rho(i,j)+rho(i-1,j)) &
             &          )/dx &
             &         +dble(2)*( &
             &            rho_w(i,j)/(rho(i,j+1)+rho(i,j))   &
             &           -rho_w(i,j-1)/(rho(i,j)+rho(i,j-1)) &
             &          )/dz &
             &        )                                &
             &    + K_H(j)*drho_T_dx4(i,j)             &
             &    + K_vu(j)*drho_T_dz2(i,j)            &
             &    - Rf(j)*drho_T(i,j)                  &
             &  )
  !motion equation x-------------------------------------------------
        drho_u_dt(i,j) =                   &
             &  (                          &
             &    - (                                              &
             &        ((rho_u_hf(i+1,j)+rho_u_hf(i,j))*dble(0.5))**2 &
             &        /rho(i+1,j)                                    &
             &       -((rho_u_hf(i,j)+rho_u_hf(i-1,j))*dble(0.5))**2 &
             &        /rho(i,j)                                  &
             &      )/dx                                           &
             !
             &    - (                                           &
             ! &        ((rho_u_hf(i,j+1)+rho_u_hf(i,j))*dble(0.5)) & !rho_u_hh
             ! &        *w_hh(i,j+1)                              &
             ! &       -((rho_u_hf(i,j)+rho_u_hf(i,j-1))*dble(0.5)) & !rho_u_hh
             ! &        *w_hh(i,j)                                &
             &        ((rho_w(i+1,j)+rho_w(i,j))*dble(0.5)) & 
             &        *u_hh(i,j)                                &
             &       -((rho_w(i+1,j-1)+rho_w(i,j-1))*dble(0.5)) & 
             &        *u_hh(i,j-1)                              &
             &      )/dz                                        &
             !
             &    - (dp(i+1,j) - dp(i,j))/dx  &
             &    + K_h(j)*drho_u_dx4(i,j)  &
             &    + K_vu(j)*drho_u_dz2(i,j) &
             &    - Rf(j)*drho_u(i,j)       &
             &  )
  !motion equation z-------------------------------------------------
        drho_w_dt(i,j) =                   &
             &  (                          &
             ! &    + (rho_w_hh(i,j)*u_hh(i,j)           &
             ! &       - rho_w_hh(i-1,j)*u_hh(i-1,j))/dx &
             &    - ( &
             ! &        ((rho_w(i+1,j)+rho_w(i,j))*dble(0.5)) & !rho_w_hh
             ! &        *u_hh(i+1,j)                              &
             ! &       -((rho_w(i,j)+rho_w(i-1,j))*dble(0.5)) & !rho_w_hh
             ! &        *u_hh(i,j)                                &
             &        ((rho_u_hf(i,j+1)+rho_u_hf(i,j))*dble(0.5))     &
             &        *w_hh(i,j)                                      &
             &       -((rho_u_hf(i-1,j+1)+rho_u_hf(i-1,j))*dble(0.5)) &
             &        *w_hh(i-1,j)                                    &
             &      )/dx                                              &
             !
             &    - ( &
             &        ((rho_w(i,j+1)+rho_w(i,j))*dble(0.5))**2 &
             &        /rho(i,j+1)                              &
             &       -((rho_w(i,j)+rho_w(i,j-1))*dble(0.5))**2 &
             &        /rho(i,j)                                &
             &      )/dz                                       &
             !
             &    - (dp(i,j+1) - dp(i,j))/dz          &
             &    - g*(drho(i,j+1)+drho(i,j))/dble(2) &
             &    + K_h(j)*drho_w_dx4(i,j)  &
             &    + K_vw(j)*drho_w_dz2(i,j) &
             &    - Rf(j)*rho_w(i,j)        &
             &  )
     end do
  end do
  !

  ! do i = 1, grd_x
  !    drho_w_dt(i,grd_z+1) =  0.0
  !    drho_w_dt(i,1) =  0.0
  ! end do


  
  ! write(6,*) "rho_dt,drho_w_dt,hydro"
  ! do i = 1, grd_z
  !    write(6,100) rho_dt(1,i),drho_w_dt(1,i),hydro(1,i)
  ! end do

 
  return
  
end subroutine gravity_wave

end module gravity_wave_eq
