!=====================================================================
! 共有変数の初期値設定のモジュール
! 2002/01/20
!=====================================================================
module com_var
  use prm  
  use prmphys, only : A_con, R_con
  implicit none
  real(8) , save :: x(1:IMAX)              ! x 座標値
  real(8) , save :: zk(1:KMAX)             ! ζ座標値
  real(8) , save :: zk_half(0:KMAX)             ! ζ座標値
  real(8) , save :: wgt(1:KMAX)


  real(8) , save :: delta_zk_j(0:KMAX-1)   ! Δζj
  real(8) , save :: delta_zk_j1_2(0:KMAX)  ! Δζj+1/2
  real(8) , save :: Apara(1:IMAX, 1:JMAX, 1:KMAX)

  integer        :: i,j, k             ! do 変数

  real(8), parameter :: Hequ(1:IMAX, 1:JMAX) = 0.0D0  ! アイソスタシー平衡高度

contains

!=====================================================================
! 温度依存パラメータ関数の初期値
! 2002/01/26
!=====================================================================
  subroutine  A_T(Apara_flg, Temp1, Tpm)
    implicit none
    real(8), intent(in) :: Temp1(1:IMAX, 1:JMAX, 1:KMAX)
    real(8), intent(in) :: Tpm(1:IMAX, 1:JMAX, 1:KMAX)
    integer ,intent(in) :: Apara_flg 
    real(8)             :: Aleft(1:IMAX/2, 1:JMAX, 1:KMAX)
    real(8)             :: Aright(1:IMAX/2, 1:JMAX, 1:KMAX)
    real(8)             :: Q0   ! 活性化エネルギー
    real(8)             :: A0   ! パラメータ

    if(Apara_flg == 0)then
       Apara = A_con
    else
       do i = 1, IMAX
          do j=1, JMAX
             do k = 1, KMAX
                if(Temp1(i, j, k) <=Tpm(i, j, k)-10.0D0)then
                   Q0 = 60.D3
                   A0 = 3.61D-13 ! 1.14D-5/s_a
                else
                   Q0 = 139.D3
                   A0 = 1.73D3!5.47D-10/s_a
                end if
                Apara(i, j, k) = A0 * exp(-Q0/(R_con*Temp1(i, j, k)))  
!             write(6, *)'Apara', Apara(i, j, k)*s_a, i,j ,k, Temp1(i, j, k)
             end do
          end do
       end do
       Aleft(:, :, :)= Apara(1:15, :, :)
       do i = IMAX/2+2, IMAX
          Aright(IMAX-i+1, :, :)= Apara(i, :, :)
       end do
    
       do j =1, JMAX
          do k=1, KMAX
             do i =1, IMAX/2
                if((Aleft(i, j, k)/Aright(i, j, k)  < 0.9999).or.  &
                     (Aleft(i, j, k)/Aright(i, j, k)  > 1.0001))then
                   write(6, *)i, j, k, 'Aprara =', Aleft(i, j, k), Aright(i, j, k), &
                        & 'Temp1=', Temp1(i, j, k), Temp1(IMAX+1-i, j, k)
                   
                end if
             end do
          end do
       end do
    end if


  end subroutine A_T


!=====================================================================
! 座標軸の設定
! allocate 属性を使って層を選択できるようにする
!=====================================================================
  subroutine get_frame

    x  = ((/(i, i=1, IMAX)/)-1) * dx

! 3 層
!    zk = (/0.10D0, 0.30D0, 0.70D0/)
!    zk_half= (/0.0D0, 0.20D0, 0.4D0, 1.0D0/)

! 11 層 
!    if(KMAX == 11)then
       zk = (/0.01D0, 0.03D0, 0.05D0, 0.07D0, 0.09D0, 0.175D0,0.325D0, 0.475D0, 0.625D0, 0.775D0, 0.925D0/)  

       zk_half = (/0.00D0, 0.02D0, 0.04D0, 0.06D0, 0.08D0, 0.10D0, 0.25D0, 0.40D0, 0.55D0, 0.70D0, 0.85D0, 1.0D0/)


! 0.05 刻み
!    zk = ((/(i, i=1, KMAX)/)) * 0.05D0 -0.025D0
!    zk_half = ((/(i, i=0, KMAX)/)) * 0.05D0

! 0.1 刻み
!    zk = ((/(i, i=1, KMAX)/)) * 0.1D0 -0.05D0
!    zk_half = ((/(i, i=0, KMAX)/)) * 0.1D0


    ! 各格子点の大きさ
    do k = 1, KMAX
       wgt(k) = (zk_half(k) -zk_half(k-1))
    end do


!    elseif(KMAX == 20)then
! 20 層
!       zk = (/0.0, 0.01, 0.02, 0.03, 0.05, 0.07, 0.10, 0.13, 0.17, 0.20, 0.25, &
!            & 0.32, 0.40, 0.47, 0.55, 0.63, 0.70, 0.77, 0.85, 0.93, 1.0/)
!    elseif(KMAX == 40)then
! 40 層
!       zk = (/0.0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.04, 0.05, 0.06, &
!            & 0.07, 0.085, 0.10, 0.115, 0.13, 0.15, 0.17, 0.18, 0.20, 0.22, &
!            & 0.25, 0.28, 0.32, 0.36, 0.40, 0.43, 0.47, 0.52, 0.55, 0.58,  &
!            & 0.63, 0.67, 0.70, 0.73, 0.77, 0.80, 0.85, 0.89, 0.93, 0.96, 1.0/)
!    end if
! 等間隔grid
!    zk = (/0.0D0, 0.1D0, 0.2D0, 0.3D0, 0.4D0, 0.5D0, 0.6D0, 0.7D0, 0.8D0, 0.9D0, 1.0D0/)  

  end subroutine get_frame

!=====================================================================
! ζ座標軸に関する変数達
! 格子間隔の設定
!  2002/01/29
!=====================================================================

!  subroutine get_zeta_hensu(delta_zk_j, delta_zk_j1_2)
!    real(8) , intent(out)  :: delta_zk_j(0:KMAX-1)     !Δζj
!    real(8) , intent(out)  :: delta_zk_j1_2(0:KMAX)    !Δζ_{j+1/2}
    

    ! ζ 格子間隔の設定
    ! 境界条件(上下境界で間隔が半分になっているので 2 倍)
!    delta_zk_j1_2(0) = 2.0D0*(zk(1) - zk(0))         + epsilon(0.0)
!!    delta_zk_j1_2(0) = (zk(1) - zk(0))         
!    delta_zk_j1_2(KMAX) = 2.0D0*(zk(KMAX) - zk(KMAX-1)) + epsilon(0.0) 
!!    delta_zk_j1_2(km) = (zk(km) - zk(km-1))
    

 !   delta_zk_j(0)   = zk(1) - zk(0) 
    
!    do j = 1, KMAX-1 
!       delta_zk_j(j) = zk(j+1) - zk(j) 
                                                ! 整数格子点での格子間隔
!       delta_zk_j1_2(j) = (zk(j+1) - zk(j-1))/2.0D0 
                                                ! 半格子点での格子間隔
!    end do




!  end subroutine get_zeta_hensu
end module com_var



