module constants
	implicit none
	!無次元パラメータ
	real(8), parameter :: thermal_Rossby =10.0**(3.0)			!外部熱ロスビー数
	real(8), parameter :: in_tau_Omega = 10.0**(0.0)				!1/(tau * Omega)

	!定数の定義
	
	real(8), parameter :: Rossby =sqrt(thermal_Rossby*10.0)				!外部熱ロスビー数調整用パラメータ
	!Rossby = 1.0 で地球程度(Held and Hou (1980)と同等の設定：熱ロスビー数は0.1程度)
	!Rossby = X のとき熱ロスビー数のみが(エクマン数などは変わらないように) X^2倍 になる。
	real(8), parameter :: tau_para = (10.0**2)*in_tau_Omega
	
	integer, parameter :: restart = 0      ! restart = 1ならリスタートファイル(restart.nc)を読む。restart = 2なら2Dのリスタートファイル(restart2D.nc)を読む
	integer, parameter :: initial_condition = 0 !initila_condition = 1 ならu=u_e, theta=theta_e の基準温位場x温度風の解
	
	real, parameter :: init_day = 000.0
	real, parameter :: output_interval = 1
	integer, parameter :: integration_days = 100
	
	
	integer, parameter :: NN = 85				! 切断波数
	integer, parameter :: MM = NN				! 切断波数
	integer, parameter :: KM = 32				! 鉛直層数
	integer, parameter :: NM = NN+1
	integer, parameter :: IM = NN*3+1		! 東西グリッド数
	integer, parameter :: JM = IM/2		! 南北グリッド数
	
	real(8), parameter :: pi=3.1415926535d0
	real(8), parameter :: H = 8000.0			!領域の高さ
	real(8), parameter :: dz = H / real(KM,8)			!1層の厚さ
	real(8),parameter :: indz = 1.0/dz
	real(8), parameter :: dt = 60.0*20.0		!時間刻み
	real(8), parameter ::  R = 6.4d6			!惑星半径
	real(8), parameter :: inR=1.0/R
	real(8), parameter :: delta_H = 1.0/3.0
	real(8), parameter :: delta_V = 0.0
	
	!通常粘性
	!real(8), parameter :: nu_H = 1.0d8/Rossby, kappa_H=0.0
	!4階の超粘性(T85で緩和時間2.4時間)
	!real(8), parameter :: nu_H =8.9d1, kappa_H=0.0 !水平拡散係数、水平熱拡散係数
	!4階の超粘性(T42で緩和時間2.4時間)
	real(8), parameter :: nu_H =1.45d3, kappa_H=0.0 !水平拡散係数、水平熱拡散係数
	!4階の超粘性(T21で緩和時間24時間)
	!real(8), parameter :: nu_H =2.2d3, kappa_H=0.0 !水平拡散係数、水平熱拡散係数
	
	real(8), parameter :: nu_V = 10.0/Rossby , kappa_V = 10.0/Rossby	!鉛直拡散係数、鉛直熱拡散係数
	real(8), parameter :: C = 0.005/Rossby		!摩擦係数
	!real(8), parameter :: C = 0.000
	real(8), parameter :: Omega = (7.272205419010587d-5)/Rossby		!自転角速度
	!real(8), parameter :: Omega = 0.0
	real(8), parameter :: tau = 20.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=9.8				!重力加速度
	!real(8), parameter :: ga=0.0
	real(8), parameter :: theta_0=250.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
	
	
	! SJPACKで用いられる配列
	
	real(8) :: SJP(JM/2*(NN+4)), SJR((NN+1)*(2*NM-NN-1)+1), SJT(IM*6)
	integer :: SJIT(4)
	
	
	integer :: ni((NN+1)*(NN+1)), mi((NN+1)*(NN+1))
	real(8) :: SJn((NN+1)*(NN+1)), SJC((NN+1)*(NN+1)), SJD((NN+1)*(NN+1)*2)
	real(8) :: phi(JM),cphi(JM),sphi(JM),incphi(JM)
	real(8) :: lambda(IM)
	
	contains
	subroutine SJ_init
		integer :: i
		print *, "Initializing....."
		call SJINIT(NN,NM,JM,IM, SJP,SJR,SJIT,SJT)
		call SJINIC(NN,SJC)
		call SJINID(NN,SJD)
		
		
		!$omp parallel do
		do i=1 , (NN+1)*(NN+1)
			call SJL2NM(NN,i,ni(i),mi(i))
			SJn(i) = real(ni(i),8)
		enddo
		!omp end parallel do
	endsubroutine 
	
	subroutine philambda
		integer :: i, j
		real(8) :: LTQ(JM*(NN+1)),  LTR((NN+1)*(NN+1))
		
		
		call LTINIT(NN,JM,LTQ,LTR)
		call LTOGRD(JM,phi,LTQ)
		
		
		!$omp parallel do 
		do i=1 , IM 
			lambda(i) = 2.0*pi/real(IM,8)*real(i-1,8)
			!print *, lambda(i)
		enddo
		!$omp end parallel do
		
		!$omp parallel do
		do j=1, JM
			cphi(j) = dcos(phi(j))
			sphi(j) = dsin(phi(j))
			incphi(j) = 1.0/cphi(j)
		enddo
		
		!omp end parallel do
		
	endsubroutine
	
endmodule
!-----------------------------------------------------------------
