module constants
implicit none

!無次元パラメータ
real(8), parameter :: thermal_Rossby =10.0**(2)			!外部熱ロスビー数
real(8), parameter :: in_tau_Omega = 10.0**(-3.0)				!1/(tau * Omega)
real(8), parameter :: E_H = 10.0**(1.0)
real(8), parameter :: E_V = 10.0**(-2.0)

!定数の定義
real(8), parameter :: Rossby =sqrt(thermal_Rossby*10.0**(-4.0))		!外部熱ロスビー数調整用パラメータ
															

real(8), parameter :: tau_para = (10.0)*in_tau_Omega

integer, parameter :: restart = 0      ! restart = 1ならリスタートファイル(restart2D.nc)を読む
integer, parameter :: initial_condition = 0 !initila_condition = 1 ならu=u_e, theta=theta_e の基準温位場x温度風の解

real, parameter :: init_day = 0.0
real, parameter :: output_interval = 100.0
integer, parameter :: integration_days = 50000


integer, parameter :: NN = 42				! 切断波数
integer, parameter :: MM = NN				! 切断波数
integer, parameter :: KM = 50+8				! 鉛直層数
integer, parameter :: NM = NN+1
!integer, parameter :: IM = 1		! 東西グリッド数
integer, parameter :: JM = (NN*3+2)/2		! 南北グリッド数

real(8), parameter :: pi=3.1415926535d0
real(8), parameter :: H = 5.0d4			!領域の高さ
real(8), parameter :: delta_z1 = H / real(KM-8,8)			!1層の厚さ（内部）
real(8), parameter :: delta_z2 = delta_z1/2.0			!1層の厚さ(境界付近)
real(8), parameter :: delta_z3 = delta_z1/4.0			!1層の厚さ(境界付近)

!real(8), parameter :: indz = 1.0/dz
real(8), parameter :: dt = 60.0*5.0*1.0		!時間刻み
real(8), parameter ::  R = 6.05d6			!惑星半径
real(8), parameter :: inR=1.0/R
real(8), parameter :: delta_H = 1.0/10.0
real(8), parameter :: delta_V = 0.0

!通常粘性
!real(8), parameter :: nu_H = 1.0d5, kappa_H=0.0
real(8), parameter :: nu_H = E_H*1.0d7/Rossby, kappa_H=0.0	!水平拡散係数、水平熱拡散係数
!real(8), parameter :: nu_H =0.0, kappa_H=0.0	!水平拡散係数、水平熱拡散係数

real(8), parameter :: nu_V = E_V*1000.0/Rossby , kappa_V = E_V*1000.0/Rossby	!鉛直拡散係数、鉛直熱拡散係数
!real(8), parameter :: nu_V = 0.0/Rossby , kappa_V = 0.0/Rossby	!鉛直拡散係数、鉛直熱拡散係数
!real(8), parameter :: C = 0.005/Rossby		!摩擦係数
!real(8), parameter :: C = 0.000
real(8), parameter :: Omega = (3.0d-7)/Rossby		!自転角速度
real(8), parameter :: inOmega = 1.0/Omega
!real(8), parameter :: Omega = 0.0
real(8), parameter :: tau = 25.0*24.0*60.0*60.0*Rossby/tau_para			!ニュートン加熱冷却の緩和時間
real(8), parameter :: intau = 1.0/tau
!real(8), parameter :: tau = 20.0*24.0*60.0*60.0d20
real(8), parameter :: ga=8.9				!重力加速度
!real(8), parameter :: ga=0.0
real(8), parameter :: theta_0=450.0			!平均温位
integer, parameter :: daystep = 86400/DINT(dt)		!1日計算するのに必要なステップ数
integer, parameter :: maxdaystep = daystep*integration_days		!積分日数
!real(8) :: alpha=(2.0*nu_V-C*dz)/(2.0*nu_V+C*dz)
real(8) :: alpha=-1.0


! LJPACKで用いられる配列

real(8) :: LJP(JM/2*(NN+4)), LJR((NN+1)*(2*NM-NN-1)+1)

real(8) :: LJn(NN+1), LJC((NM+1)*(NM+1))
real(8) :: phi(JM),cphi(JM),sphi(JM),incphi(JM)

real(8) :: z(KM+1),zhat(KM+2)
real(8) :: indz(KM), indzhat(KM+1), indzhat2(2:KM+1)

contains

subroutine LJ_init
integer :: j 
print *, "Initializing....."
call LJINIT(NN,NM,JM,LJP,LJR)
call SJINIC(NN,LJC)	

do j=1,NN+1
LJn(j) = real(j-1,8)
!print *, j, LJn(j)
enddo
endsubroutine

subroutine philambda
integer :: j
real(8) :: LTQ(JM*(MM+1)),  LTR((NN+1)*(NN+1))
call LTINIT(NN,JM,LTQ,LTR)
call LTOGRD(JM,phi,LTQ)
!$omp parallel do 
do j=1, JM
cphi(j) = dcos(phi(j))
sphi(j) = dsin(phi(j))
incphi(j) = 1.0/cphi(j)
!print *, j, phi(j)
enddo

!omp end parallel do

endsubroutine

subroutine definez
integer :: k

!鉛直格子生成
z(1) = 0.0
do k = 2, KM+1
!z(k) = real((k-1),8)*H/real(KM,8)
if(k < 6) then
z(k) = delta_z3 + z(k-1)
elseif(k < 8) then
z(k) = delta_z2 + z(k-1)
elseif(k < KM-4)then
z(k) = delta_z1 + z(k-1)
elseif(k<KM-2)then
z(k) = delta_z2 + z(k-1)
else
z(k) = delta_z3 + z(k-1)
endif
enddo




zhat(1) = -0.5*(z(2)+z(1))
do k = 2, KM+1
zhat(k) = 0.5*(z(k)+z(k-1)) 
enddo
zhat(KM+2) = 2.0*z(KM+1)-zhat(KM+1)
do k = 1, KM
indz(k) = 1.0/(z(k+1)-z(k))
enddo
do k = 1, KM+1
indzhat(k) = 1.0/(zhat(k+1)-zhat(k))
enddo
do k =2 ,KM+1
indzhat2(k) = 1.0/(zhat(k+1)-zhat(k-1))
enddo

do k= 1,KM+1
print *, k,zhat(k)
print *, k,z(k)
enddo
print *, KM+2, zhat(KM+2)

endsubroutine


endmodule
!-----------------------------------------------------------------
