module constants
implicit none

!定数の定義

real(8), parameter :: Rossby = 23.7			!外部熱ロスビー数調整用パラメータ
											!Rossby = 1.0 で地球程度(Held and Hou (1980)と同等の設定：熱ロスビー数は0.1程度)
											!Rossby = X のとき熱ロスビー数のみが(エクマン数などは変わらないように) X^2倍 になる。
real(8), parameter :: tau_para = 178.0

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 = 1
integer, parameter :: integration_days = 50


integer, parameter :: NN = 21				! 切断波数
integer, parameter :: MM = NN				! 切断波数
integer, parameter :: KM = 10				! 鉛直層数
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 = 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 :: nu_H = 1.0d5, kappa_H=0.0
real(8), parameter :: nu_H =1.0d9/Rossby, kappa_H=0.0	!水平拡散係数、水平熱拡散係数
!real(8), parameter :: nu_H =0.0, kappa_H=0.0	!水平拡散係数、水平熱拡散係数

real(8), parameter :: nu_V = 10.0/Rossby , kappa_V = 10.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 = (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


! 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)

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

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