!-----------------------------------------------------------------

include "module_BPM3D_omp.f90"
include "wrap_sj.f90"

program busi1
use constants
implicit none

call main	!メインルーチン呼び出し

contains

!----------------------------------------------------------
subroutine main		!メインルーチン
	use gtool_history
	use constants
	use wrap_sj
	
	implicit none
	! スペクトルデータ
	real(8) :: psi((NN+1)*(NN+1),KM+2) , oldpsi((NN+1)*(NN+1),KM+2) 
	real(8) :: kpsi((NN+1)*(NN+1),KM+2,4)			!ルンゲクッタ用
	real(8) :: theta_s((NN+1)*(NN+1),KM+2), oldtheta_s((NN+1)*(NN+1),KM+2)
	real(8) :: ktheta_s((NN+1)*(NN+1),KM+2,4)		!ルンゲクッタ用
	real(8) :: chi2z((NN+1)*(NN+1),KM+1) , oldchi2z((NN+1)*(NN+1),KM+1)
	real(8) :: kchi2z((NN+1)*(NN+1),KM+1,4)			!ルンゲクッタ用
	
	real(8) :: chi((NN+1)*(NN+1),KM+1)
	real(8) :: zeta_s((NN+1)*(NN+1),KM+2)
	real(8) :: E_s((NN+1)*(NN+1),KM+2)
	real(8) :: wtheta_s((NN+1)*(NN+1),KM+2),Vtheta_s((NN+1)*(NN+1),KM+2),Utheta_s((NN+1)*(NN+1),KM+2)
	real(8) :: chi1z((NN+1)*(NN+1),KM+2)
	real(8) :: w_s((NN+1)*(NN+1),KM+1)
	real(8) :: theta_es((NN+1)*(NN+1),KM+2)
	
	! グリッドデータ
	real(8) :: psigrid(IM,JM,KM+2)
	real(8) :: chigrid(IM,JM,KM+1)
	
	real(8) :: theta_g(IM,JM,KM+2)
	real(8) :: U(IM,JM,KM+2),  V(IM,JM,KM+2)
	real(8) :: zeta_g(IM,JM,KM+2)
	real(8) :: w_g(IM,JM,KM+1)
	real(8) :: A_g(IM,JM,KM+2)
	real(8) :: A1z(IM,JM,KM+1)
	real(8) :: B_g(IM,JM,KM+2)
	real(8) :: B1z(IM,JM,KM+1)
	real(8) :: E_g(IM,JM,KM+2)
	real(8) :: Utheta_g(IM,JM,KM+2)
	real(8) :: Utheta1z(IM,JM,KM+1)
	real(8) :: Vtheta_g(IM,JM,KM+2)
	real(8) :: Vtheta1z(IM,JM,KM+1)
	real(8) :: wtheta_g(IM,JM,KM+2)
	real(8) :: rdm
	
	
	real(8) :: spct(NN+1)
	real(8) :: chi2zval((NN+1)*(NN+1),KM+2)
	
	integer :: i,j,k,n,tstep,	time	! tstep =  ルンゲクッタのカウンター
	
	type(GT_HISTORY) :: hst_data, hst_restart
	character(len=80) :: filename, filename2
	
	real(8) :: psi2D(NN+1,KM+2), theta_s2D(NN+1,KM+2), chi2z2D(NN+1,KM+2)
	
	
	call SJ_init
	call philambda
	
	
	
	call defineTE(theta_es)
	do i=1,(NN+1)*(NN+1)
		theta_es(i,KM+2) = theta_es(i,KM+1)
		theta_es(i,1) = theta_es(i,2)
	enddo
	print *, "Theta_e was defined"
	
	
	
	call RANDOM_SEED
	if(initial_condition==1)then
		do k=2, KM+1
			do j=1, JM
				psigrid(:,j,k) = - R*Omega*(dsqrt(1+((real(k-1,8)-0.5)*dz/H)*(2.0*ga*delta_H)/(R*R*Omega*Omega))-1)*sphi(j)
			enddo
		enddo
		psigrid(:,:,KM+2) = psigrid(:,:,KM+1)
		psigrid(:,:,1) = psigrid(:,:,2)
		
		do k=1, KM+2
			call S_S2G(theta_es(:,k), theta_g(:,:,k))
			do j =1, JM
				do i = 1, IM
					call RANDOM_NUMBER(rdm)
					theta_g(i,j,k) = theta_g(i,j,k)*(1+0.01*(rdm-0.5))
				enddo
			enddo
		enddo
		print *, "Initial condition was set with thermal wind in theta_e "
		chi2z = 0.0
		
		else
		
		!omp parallel do private(rdm)
		do k=1,KM+2
			do i=1,IM
				do j=1,JM
					call RANDOM_NUMBER(rdm)
					theta_g(i,j,k) = rdm*0.1
					!theta_g = theta_0*0.0
				enddo
			enddo
			call S_G2S(theta_g(:,:,k),theta_s(:,k))
		enddo
		!omp end prallel do
		print *, "Initial condition was set with constant PT in rest."
	endif
	
	
	
	!-------------------リスタートファイルの読み込み：この部分を有効にするとリスタートファイルを読む
	
	
	if(restart == 1) then
		filename = 'restart.nc'
		call HistoryGet(file=filename, varname='psi', array=psi)
		call HistoryGet(file=filename, varname='theta_s', array=theta_s)
		call HistoryGet(file=filename, varname='chi2z', array=chi2zval)
		
		do k=1,KM+1
			do n=1,MM+1
				chi2z(n,k)=chi2zval(n,k)
			enddo
		enddo
		
		elseif(restart==2)then
		
		filename = 'restart2D.nc'
		call HistoryGet(file=filename, varname='psi', array=psi2D)
		call HistoryGet(file=filename, varname='theta_s', array=theta_s2D)
		call HistoryGet(file=filename, varname='chi2z', array=chi2z2D)
		
		psi = 0.0
		theta_s = 0.0
		chi2z = 0.0
		
		psi(1:NN+1,:) = psi2D
		theta_s(1:NN+1,:) = theta_s2D
		chi2z(1:NN+1,:) = chi2z2D(:,1:KM+1)
		
		do k=2,KM+1
			call S_S2G(theta_s(:,k), theta_g(:,:,k))
			do j=1,JM
				do i=1,IM
					call RANDOM_NUMBER(rdm)
					theta_g(i,j,k) = theta_g(i,j,k)*(1+0.01*(rdm-0.5))
				enddo
			enddo
			call S_G2S(theta_g(:,:,k), theta_s(:,k))
		enddo
		
		theta_s(:,1) = theta_s(:,2)
		theta_s(:,KM+2) = theta_s(:,KM+1)
		
		
	endif
	!------------------------------------------
	
	
	!!!!!!!!!!!!!!!!!!!!!!!!!水平の力学チェック    P^1_1(mu)=root3 cos(phi)  Y^1_1=root3 cos(phi)cos(lambda)
	!do i=1,IM
	!do j=1,JM
	!do k=1, KM+2
	!psigrid(i,j,k) = sqrt(3.0)*cos(phi(j))*cos(lambda(i))/R
	!enddo
	!enddo
	!enddo
	
	!do k=1,KM+2
	!call S_G2S(psigrid(:,:,k),psi(:,k)) 
	!enddo
	
	do time=1,maxdaystep			!------------time step !--------
		
		oldpsi = psi
		oldtheta_s = theta_s
		oldchi2z = chi2z 
		
		do tstep = 0, 3
			
			if (tstep == 1) then
				!$omp parallel do 
				do k=2,KM+1
					do i=1,(NN+1)*(NN+1)
						psi(i,k) = oldpsi(i,k) + kpsi(i,k,1)*0.5
						theta_s(i,k) = oldtheta_s(i,k) + ktheta_s(i,k,1)*0.5
						chi2z(i,k) = oldchi2z(i,k) + kchi2z(i,k,1)*0.5
					enddo
				enddo
				!$omp end parallel do 
				
				else if (tstep == 2) then
				!$omp parallel do 
				do k=2,KM+1
					do i=1,(NN+1)*(NN+1)
						psi(i,k) = oldpsi(i,k) + kpsi(i,k,2)*0.5
						theta_s(i,k) = oldtheta_s(i,k) + ktheta_s(i,k,2)*0.5
						chi2z(i,k) = oldchi2z(i,k) + kchi2z(i,k,2)*0.5
					enddo
				enddo
				!$omp end parallel do 
				
				else if (tstep == 3) then
				!$omp parallel do 
				do k=2,KM+1
					do i=1,(NN+1)*(NN+1)
						psi(i,k) = oldpsi(i,k) + kpsi(i,k,3)
						theta_s(i,k) = oldtheta_s(i,k) + ktheta_s(i,k,3)
						chi2z(i,k) = oldchi2z(i,k) + kchi2z(i,k,3)
					enddo
				enddo
				!$omp end parallel do 
			endif
			
			
			
			
			call useLAPACK(chi2z,chi)  ! chi2z から chi  を求める。
			
			
			
			!$omp parallel do
			do k=2, KM+1
				call S_S2S_lap(psi(:,k),zeta_s(:,k))		! psi から zeta_s を求める。
				
				call S_S2G_dual(theta_s(:,k),zeta_s(:,k),theta_g(:,:,k),zeta_g(:,:,k))   ! theta_s から theta_g ,  zeta_s から zeta_g  を求める（内点のみ）
				
				call S_S2S_lap(chi(:,k),w_s(:,k))	     ! chi から w_s を求める。
				w_s(:,k) = -w_s(:,k)
				
				call S_S2G(w_s(:,k),w_g(:,:,k))		! w_s から w_g  を求める（内点のみ）
				
				chi1z(:,k)=(chi(:,k)-chi(:,k-1))*indz    ! chi から chi1z を求める。
			enddo
			!$omp end parallel do
			
			call S2GUV(psi,chi1z,U,V)   ! psi chi1z から  U V  を求める。
			
			!$omp parallel do
			do i = 1, (NN+1)*(NN+1)		          !境界条件
				zeta_s(i,KM+2) = zeta_s(i,KM+1)
				psi(i,KM+2) = psi(i,KM+1)
				zeta_s(i,1) = alpha*zeta_s(i,2)
				psi(i,1) = alpha*psi(i,2)
				theta_s(i,KM+2) = theta_s(i,KM+1)
				theta_s(i,1) = theta_s(i,2)
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do j=1 , JM
				do i=1 , IM
					theta_g(i,j,KM+2) = theta_g(i,j,KM+1)			!境界条件
					theta_g(i,j,1) = theta_g(i,j,2)
					zeta_g(i,j,KM+2) = zeta_g(i,j,KM+1)
					zeta_g(i,j,1) = alpha*zeta_g(i,j,2)
					w_g(i,j,KM+1) = 0.0
					w_g(i,j,1) = 0.0
				enddo
			enddo
			!$omp end parallel do
			
			! 非線形項の計算。（グリッド）--------------------------------------
			
			!$omp parallel do
			do k=2, KM+1
				do j=1,JM
					do i=1,IM 
						A_g(i,j,k) = (zeta_g(i,j,k)+2.0*Omega*sphi(j))*U(i,j,k) + (w_g(i,j,k)+w_g(i,j,k-1))*0.5*(V(i,j,k+1)-V(i,j,k-1))*0.5*indz
						B_g(i,j,k) = (zeta_g(i,j,k)+2.0*Omega*sphi(j))*V(i,j,k) - (w_g(i,j,k)+w_g(i,j,k-1))*0.5*(U(i,j,k+1)-U(i,j,k-1))*0.5*indz
					enddo
				enddo
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do k=1, KM+2
				do j=1,JM
					do i=1,IM 
						E_g(i,j,k) = (U(i,j,k)*U(i,j,k)+V(i,j,k)*V(i,j,k))*0.5*incphi(j)*incphi(j)
						Utheta_g(i,j,k) = U(i,j,k)*theta_g(i,j,k)
						Vtheta_g(i,j,k) = V(i,j,k)*theta_g(i,j,k)
					enddo
				enddo
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do k=1, KM+1
				do j=1,JM
					do i=1,IM 
						wtheta_g(i,j,k) = w_g(i,j,k)*(theta_g(i,j,k+1)+theta_g(i,j,k))*0.5
					enddo
				enddo
			enddo
			!$omp end parallel do
			
			
			
			! 非線形項のスペクトル変換
			!$omp parallel do
			do k =1, KM+1
				call S_G2S_dual(E_g(:,:,k),wtheta_g(:,:,k),E_s(:,k),wtheta_s(:,k))
			enddo
			!$omp end parallel do
			
			
			! グリッドデータの半整数グリッドの値 x から整数グリッドのdx/dz をもとめる。
			!$omp parallel do
			do k=2, KM
				do j=1, JM
					do i=1, IM
						A1z(i,j,k)=(A_g(i,j,k+1)-A_g(i,j,k))*indz
						B1z(i,j,k)=(B_g(i,j,k+1)-B_g(i,j,k))*indz
						Utheta1z(i,j,k)=(Utheta_g(i,j,k+1)-Utheta_g(i,j,k))*indz
						Vtheta1z(i,j,k)=(Vtheta_g(i,j,k+1)-Vtheta_g(i,j,k))*indz
					enddo
				enddo
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do i=1, (NN+1)*(NN+1)       !境界条件
				chi2z(i,KM+1) = 0.0
				chi2z(i,1) = chi(i,2)*(1.0-alpha)*indz*indz
			enddo
			!$omp end parallel do
			
			if (tstep < 4)then
				call nextpsi(psi,A_g,B_g,kpsi,tstep)
				call nextchi2z(chi2z,A1z,B1z,theta_s,E_s,kchi2z,tstep)
				call nexttheta(theta_s,Utheta_g,Vtheta_g,wtheta_s,theta_es,ktheta_s,tstep)
			endif
			
			
			
			!---------------ouput--------------
			if(mod(time,daystep*INT(output_interval))==0)then
				if(tstep==0)then
					print *, time/daystep+INT(init_day), 'days'
					!print *, "u,v", maxval(U), maxval(V)
					call output(time,U,V,theta_g,chi,theta_s,psi,chi2z,hst_data,hst_restart)
					open(8,FILE='psiwave.txt')	!数値粘性が不十分でないかの確認用に出力する
					k=KM+1
					do i=2,(NN+1)*(NN+1)
						spct(ni(i)) = spct(ni(i)) + psi(i,k)*psi(i,k)
					enddo
					do i=1, (NN)
						write(8,*)  i, spct(i)!*(i*(i+1))
					enddo
					close(8)
					
				endif
			endif
		enddo
		
		!$omp parallel do
		do i=1,(NN+1)*(NN+1)
			do k=2,KM+1
				psi(i,k) = oldpsi(i,k) + (kpsi(i,k,1) +kpsi(i,k,2)*2.0 + kpsi(i,k,3)*2.0 + kpsi(i,k,4))/6.0
				theta_s(i,k) = oldtheta_s(i,k) + (ktheta_s(i,k,1) + ktheta_s(i,k,2)*2.0 + ktheta_s(i,k,3)*2.0 + ktheta_s(i,k,4))/6.0
				chi2z(i,k) = oldchi2z(i,k) + (kchi2z(i,k,1) + kchi2z(i,k,2)*2.0 +  kchi2z(i,k,3)*2.0 + kchi2z(i,k,4))/6.0 
			enddo
		enddo
		!$omp end parallel do
		
	enddo
	
	call HistoryClose(hst_data)
	call HistoryClose(hst_restart)
	
endsubroutine

!---------------------------------------------------------------------------------

subroutine defineTE(theta_es)	!ニュートン加熱冷却の基準温位場生成
	use constants
	use wrap_sj
	implicit none
	real(8) :: theta_es((NN+1)*(NN+1),KM+2)
	real(8) :: TE(IM,JM,KM+2)
	
	integer :: i,j,k
	
	
	do k=1, KM+2
		do i=1,IM
			do j=1,JM
				TE(i,j,k)=theta_0*(-delta_H*(sphi(j)*sphi(j)-1.0/3.0)+delta_V*((real(k-1,8)-0.5)*dz/H-0.5))
				!TE(i,j,k)=theta_0*0.0
			enddo
		enddo
		call S_G2S(TE(:,:,k), theta_es(:,k))
	enddo
	
endsubroutine

!---------------------------------------------------------------------------------
subroutine useLAPACK(x2,x)   !三重対角行列で表される連立1次方程式を解く
	use constants
	implicit none
	real(8) :: x2((NN+1)*(NN+1),KM+1)
	real(8) :: x((NN+1)*(NN+1),KM+1)
	
	real(8) :: D(KM-1)
	real(8) :: DU(KM-2),DL(KM-2),DU2(KM-3)
	real(8) :: B(KM-1,1)
	integer :: info , ii , k, IPIV(KM-1),num
	character :: TRANS
	
	
	!初期化
	D = -2.0
	DU = 1.0
	DL = 1.0
	TRANS = 'N'
	! ここでLAPACKの三重対格行列を解くサブルーチンをコールする。
	call DGTTRF(KM-1,DL,D,DU,DU2,IPIV,info)
	
	!$omp parallel do firstprivate(D,DU,DL,DU2,B,IPIV,num,info,TRANS)
	do ii=1, (NN+1)*(NN+1) 
		
		do k=2,KM
			B(k-1,1) = x2(ii,k)*dz*dz
		enddo
		num =1
		
		call DGTTRS(TRANS,KM-1,num,DL,D,DU,DU2,IPIV,B,KM-1,info)
		
		if(info /= 0)then
			print *,'eroor at lapack'
		endif
		do k=2,KM
			x(ii,k) = B(k-1,1) 
		enddo
		x(ii,1)=0.0			!境界条件  chi(1)=chi(KM+1)=0
		x(ii,KM+1)=0.0
	enddo
	!$omp end parallel do 
endsubroutine

!-------------------------------------------------------------------
subroutine S2GUV(x,y,a,b)		!スペクトルデータからU,Vのグリッドデータを計算
	use constants
	use wrap_sj
	implicit none
	real(8) :: x((NN+1)*(NN+1),KM+2)
	real(8) :: y((NN+1)*(NN+1),KM+2)
	real(8) :: a(IM,JM,KM+2)
	real(8) :: b(IM,JM,KM+2)
	real(8) ::  sx1b((NN+1)*(NN+1), KM+2)
	real(8) ::  sx2b((NN+1)*(NN+1), KM+2)
	real(8) ::  g1a(IM,JM,KM+2) , g2a(IM,JM,KM+2)
	real(8) ::  g1b(IM,JM,KM+2) , g2b(IM,JM,KM+2)
	integer :: i , k ,j ,n
	
	!real(8) :: Q(JM/2*11), WS1(2*(NN+2)),WS2(2*(NN+2)),W1((JM+1)*IM),W2((JM+1)*IM),WG((IM+2)*JM)
	
	!$omp parallel do 
	do k = 2, KM+1
		
		call S_S2G_dlat_dual(x(:,k),y(:,k),g1a(:,:,k),g2a(:,:,k))
		call S_S2S_dlon(x(:,k),sx1b(:,k))
		call S_S2S_dlon(y(:,k),sx2b(:,k))
		call S_S2G_dual(sx1b(:,k),sx2b(:,k),g1b(:,:,k),g2b(:,:,k))
		
		!call SJCS2Y(NN,s1a,sy1a,SJC)
		!call SJCS2Y(NN,s2a,sy2a,SJC)
		!call SJMS2G(NN,NM,NN+1,IM,JM,sy1a,sy2a,g1a,g2a,SJIT,SJT,SJP,Q,SJR,WS1,WS2,WG,W1,W2,1)
		
		!call SJCS2X(NN,s1b,sx1b) 
		!call SJCS2X(NN,s2b,sx2b)
		!call SJMS2G(NN,NM,NN,IM,JM,sx1b,sx2b,g1b,g2b,SJIT,SJT,SJP,Q,SJR,WS1,WS2,WG,W1,W2,0)
		
		do i=1 , IM
			do j=1 , JM
				a(i,j,k) = -g1a(i,j,k)*R*cphi(j) +g2b(i,j,k)*R     !  Uの生成
				b(i,j,k) = g1b(i,j,k)*R +g2a(i,j,k)*R*cphi(j)    !  Vの生成
			enddo
		enddo
		
	enddo
	!$omp end parallel do 
	
	!$omp parallel do 
	do i=1,IM
		do j=1,JM
			a(i,j,1)=alpha*a(i,j,2)			!境界条件
			a(i,j,KM+2)=a(i,j,KM+1)
			b(i,j,1) = alpha*b(i,j,2)				!境界条件
			b(i,j,KM+2) = b(i,j,KM+1)
		enddo
	enddo
	!$omp end parallel do 
	
endsubroutine
!-------------------------------------------------------------------------------------
subroutine nextpsi(psi,A,B,kpsi,tstep)
	use constants
	use wrap_sj
	implicit none
	real(8) :: psi((NN+1)*(NN+1),KM+2),kpsi((NN+1)*(NN+1),KM+2,4)
	real(8) :: A(IM,JM,KM+2),B(IM,JM,KM+2)
	real(8) :: lA(IM,JM,KM+2),lB(IM,JM,KM+2)
	
	integer :: tstep
	
	
	real(8) :: s1((NN+1)*(NN+1),KM+2) , s2((NN+1)*(NN+1),KM+2)
	real(8) :: sx1((NN+1)*(NN+1),KM+2)
	
	integer :: i , k ,j 
	
	!$omp parallel do 
	do k=2 , KM+1
		do j =1, JM
			do i = 1, IM
				lA(i,j,k) = A(i,j,k)*incphi(j)*incphi(j)
				lB(i,j,k) = B(i,j,k)*incphi(j)
			enddo
		enddo
		
		call S_G2S_div(lB(:,:,k),s2(:,k))
		call S_G2S(lA(:,:,k), sx1(:,k))
		call S_S2S_dlon(sx1(:,k), s1(:,k))
		
		
		do i=2 , (NN+1)*(NN+1)
			!通常の粘性
			!kpsi(i,k,tstep+1) =  ((ls1(i,k)+ls2(i,k))/(R*SJn(i)*(SJn(i)+1.)) - &
			!&nu_H*(SJn(i)*(SJn(i)+1.0)-2.0)*psi(i,k)/(R*R) +&
			!&nu_V*(psi(i,k+1)-2.*psi(i,k)+psi(i,k-1))*indz*indz)*dt
			!4階の超粘性
			kpsi(i,k,tstep+1) =  ((s1(i,k)+s2(i,k))/(R*SJn(i)*(SJn(i)+1.)) - &
								  &nu_H*(SJn(i)*(SJn(i)+1.0)*SJn(i)*(SJn(i)+1.0)-2.0*SJn(i)*(SJn(i)+1.0))*psi(i,k)*inR*inR +&
								  &nu_V*(psi(i,k+1)-2.*psi(i,k)+psi(i,k-1))*indz*indz)*dt
		enddo
		
	enddo
	!$omp end parallel do 
	
endsubroutine


!---------------------------------------------------------------------------------
subroutine nextchi2z(chi2z,A1z,B1z,theta_s,E_s,kchi2z,tstep)
	use constants
	use wrap_sj
	implicit none
	real(8) :: theta_s((NN+1)*(NN+1),KM+2),chi2z((NN+1)*(NN+1),KM+1),E_s((NN+1)*(NN+1),KM+2)
	real(8) :: A1z(IM,JM,KM+1),B1z(IM,JM,KM+1),kchi2z((NN+1)*(NN+1),KM+1,4)
	integer :: tstep
	real(8) :: lA(IM,JM,KM+1),lB(IM,JM,KM+1)
	
	real(8) :: s1((NN+1)*(NN+1),KM+1) , s2((NN+1)*(NN+1),KM+1)
	real(8) :: sx2((NN+1)*(NN+1),KM+1)
	
	integer :: i , k ,j 
	
	!real(8) :: WS(2*(NN+2)), Q(JM/2*7),WG((IM+2)*JM),WWS(2*(NN+1))
	
	
	!$omp parallel do
	do k=2 , KM
		
		do j =1, JM
			do i = 1, IM
				lA(i,j,k) = A1z(i,j,k)*incphi(j)
				lB(i,j,k) = B1z(i,j,k)*incphi(j)*incphi(j)
			enddo
		enddo
		
		call S_G2S_div(lA(:,:,k), s1(:,k))
		call S_G2S(lB(:,:,k), sx2(:,k))
		call S_S2S_dlon(sx2(:,k), s2(:,k))
		
		!call SJTG2S(NN,NM,NN+1,IM,JM,lsy1,lA,SJIT,SJT,SJP,Q,SJR,WS,WG,W,1)
		
		!call SJTG2S(NN,NM,NN,IM,JM,lsx2,lB,SJIT,SJT,SJP,Q,SJR,WWS,WG,W,0)
		
		!call SJCY2S(NN,lsy1,ls1,SJC)
		!call SJCS2X(NN, lsx2, ls2)
		
		
		
		do i=2 , (NN+1)*(NN+1)
			!通常の粘性
			!kchi2z(i,k,tstep+1) =  ((-s1(i,k)+s2(i,k))/(R*SJn(i)*(SJn(i)+1.)) &
			!&         -((ga/theta_0)*(theta_s(i,k+1)+theta_s(i,k))*0.5+(E_s(i,k+1)-E_s(i,k))*indz)/(R*R)   &
			!&            - nu_H*(2.0*SJn(i)*(SJn(i)+1.)-2.0)*chi2z(i,k)/(R*R) +&
			!&nu_V*(chi2z(i,k+1)-2.*chi2z(i,k)+chi2z(i,k-1))*indz*indz)*dt
			
			!4階の超粘性
			kchi2z(i,k,tstep+1) =  ((s1(i,k)-s2(i,k))/(R*SJn(i)*(SJn(i)+1.)) &
									&         -((ga/theta_0)*(theta_s(i,k+1)+theta_s(i,k))*0.5+(E_s(i,k+1)-E_s(i,k))*indz)*inR*inR   &
									&            - nu_H*(2.0*SJn(i)*(SJn(i)+1.)*SJn(i)*(SJn(i)+1.)-2.0*SJn(i)*(SJn(i)+1.))&
									&*chi2z(i,k)*inR*inR +&
									&nu_V*(chi2z(i,k+1)-2.*chi2z(i,k)+chi2z(i,k-1))*indz*indz)*dt
			
			!if(kchi2z(i,j,tstep+1)>0.0)then
			!print *, i,k,tstep+1,kchi2z(i,j,tstep+1)
			!endif
			
		enddo
	enddo
	!$omp end parallel do
	
endsubroutine

!---------------------------------------------------------------------------------
subroutine nexttheta(theta_s,Utheta_g,Vtheta_g,wtheta_s,theta_es,ktheta_s,tstep)
	use constants
	use wrap_sj
	implicit none
	real(8) :: wtheta_s((NN+1)*(NN+1),KM+1)
	real(8) :: theta_s((NN+1)*(NN+1),KM+2)
	real(8) :: ktheta_s((NN+1)*(NN+1),KM+2,4)
	real(8) :: Utheta_g(IM,JM,KM+2),Vtheta_g(IM,JM,KM+2)
	real(8) :: lU(IM,JM,KM+2),lV(IM,JM,KM+2)
	
	real(8) :: theta_es((NN+1)*(NN+1),KM+2)
	integer :: tstep
	
	real(8) :: s1((NN+1)*(NN+1),KM+2) , s2((NN+1)*(NN+1),KM+2)
	real(8) :: sx1((NN+1)*(NN+1),KM+2)
	integer :: i , k ,j 
	
	!real(8) :: WS(2*(NN+2)), Q(JM/2*7),WG((IM+2)*JM),W((JM+1)*IM), wws(2*(NN+1))
	
	!$omp parallel do
	do k=2 , KM+1
		do j =1, JM
			do i = 1, IM
				lU(i,j,k) = Utheta_g(i,j,k)*incphi(j)*incphi(j)
				lV(i,j,k) = Vtheta_g(i,j,k)*incphi(j)
			enddo
		enddo
		
		call S_G2S_div(lV(:,:,k), s2(:,k))
		call S_G2S(lU(:,:,k), sx1(:,k))
		call S_S2S_dlon(sx1(:,k), s1(:,k))
		
		!call SJTG2S(NN,NM,NN+1,IM,JM,lsy2,lV,SJIT,SJT,SJP,Q,SJR,WS,WG,W,1)
		!call SJTG2S(NN,NM,NN,IM,JM,lsx1,lU,SJIT,SJT,SJP,Q,SJR,wws,WG,W,0)
		
		!call SJCY2S(NN,lsy2,ls2,SJC)
		!call SJCS2X(NN, lsx1, ls1)
		
		
		
		
		do i=1 , (NN+1)*(NN+1)
			ktheta_s(i,k,tstep+1) =  (-(s1(i,k)+s2(i,k))*inR  &
									  &                 -((wtheta_s(i,k)-wtheta_s(i,k-1))*indz)-(theta_s(i,k)-theta_es(i,k))*intau  &
									  &                - kappa_H*SJn(i)*(SJn(i)+1.)*theta_s(i,k)*inR*inR&
									  & +kappa_V*(theta_s(i,k+1)-2.*theta_s(i,k)+theta_s(i,k-1))*indz*indz)*dt
		enddo
		
	enddo
	!$omp end parallel do
	
endsubroutine
!---------------------------------------------------------------------------------
subroutine output(time,u,v,theta_g,chi,theta_s,psi,chi2z,hst_data,hst_restart)  !出力ルーチン
	use constants
	use wrap_sj
	use gtool_history
	
	implicit none
	real(8) :: u(IM,JM,KM+2), lon(IM),lat(JM),hight(KM+2),uval(IM,JM,KM+2)
	real(8) :: v(IM,JM,KM+2), vval(IM,JM,KM+2),m(IM,JM,KM+2)
	real(8) :: theta_g(IM,JM,KM+2),theta_s((NN+1)*(NN+1),KM+2)
	real(8) :: totalPT,totalM,totalDM,mpsi(IM,JM,KM+2),chi((NN+1)*(NN+1),KM+1), chival((NN+1)*(NN+1),KM+2)
	real(8) :: psi((NN+1)*(NN+1),KM+2), chi2z((NN+1)*(NN+1),KM+1), chi2zval((NN+1)*(NN+1),KM+2)
	
	integer :: time, i,j,k
	real(8) :: wn_isp((NN+1)*(NN+1))
	!character :: name,fname
	type(GT_HISTORY) :: hst_data, hst_restart
	
	if (time==daystep*INT(output_interval)) then
		!netcdfに書き込まれるメタ情報。適当に書き換える。
		!g95を用いるときの注意：文字列配列の要素文字列の長さをそろえる必要がある。
		call HistoryCreate(file='data3D.nc',title='BPM3D',source='BPM3D-s.f90',institution='user',&
						   & dims=(/'wn  ','z   ','time'/),dimsizes=(/(NN+1)*(NN+1),KM+2,0/),longnames=(/'wn  ','z   ','time'/), &
						   & units=(/'1/m ','m   ','days'/), origin=init_day,interval=output_interval, history=hst_data)
		
		
		do i=1,IM
			lon(i) = lambda(i)*180.0/pi
		enddo
		do j=1, JM 
			lat(j) = phi(j)*180.0/pi
		enddo
		do k=1 ,KM+2
			hight(k) = dz*real(k-1,8)
		enddo
		do i=1, (NN+1)*(NN+1)
			wn_isp(i) = i
		enddo
		
		call HistoryPut('wn  ',wn_isp,hst_data)
		
		call HistoryPut('z   ',hight,hst_data)
		
		call HistoryAddVariable(varname='psi',dims=(/'wn  ','z   ','time'/),&
								&longname='psi',units='m2/s',xtype='float',history=hst_data)
		
		call HistoryAddVariable(varname='chi',dims=(/'wn  ','z   ','time'/),&
								&longname='chi',units='m3/s',xtype='float',history=hst_data)
		
		call HistoryAddVariable(varname='theta_s',dims=(/'wn  ','z   ','time'/),&
								&longname='theta_s',units='K',xtype='float',history=hst_data)
		
		
		call HistoryAddVariable(varname='totalPT',dims=(/'time'/),&
								&longname='totalPT',units='K',xtype='float',history=hst_data)
		call HistoryAddVariable(varname='totalM',dims=(/'time'/),&
								&longname='totalM',units='m/s',xtype='float',history=hst_data)
		call HistoryAddVariable(varname='totalDM',dims=(/'time'/),&
								&longname='totalDM',units='MKS',xtype='float',history=hst_data)
		
	endif
	
	call conserves(theta_s,u,totalPT,totalM,totalDM)
	
	do j=1,(NN+1)*(NN+1)
		do k=1,KM+1
			chival(j,k)= chi(j,k)
		enddo
		chival(j,KM+2)=chi(j,KM+1)
	enddo
	
	
	
	
	
	call HistoryPut('psi', psi, hst_data)
	call HistoryPut('chi', chival, hst_data)
	call HistoryPut('theta_s',  theta_s, hst_data)
	
	
	call HistoryPut('totalPT',totalPT,hst_data)
	call HistoryPut('totalM',totalM,hst_data)
	call HistoryPut('totalDM',totalDM,hst_data)
	
	if(totalM /= totalM)then		!NaNの判定
		call HistoryClose(hst_data)
		call HistoryClose(hst_restart)
		print *, "numerical instability!!"
		stop
	endif
	
	
	!-----------------------------------------
	!リスタートファイル
	if (time==daystep*INT(output_interval)) then
		
		call HistoryCreate(file='restart.nc',title='BPM2D',source='BPM2D-s.f90',institution='user',&
						   & dims=(/'wn','z '/),dimsizes=(/(NN+1)*(NN+1),KM+2/),longnames=(/'wn','z '/), &
						   & units=(/'1/m','m  '/),history=hst_restart)
		
		call HistoryPut('wn',SJn,hst_restart)
		call HistoryPut('z ',hight,hst_restart)
		
		call HistoryAddVariable(varname='psi',dims=(/'wn','z '/),&
								&longname='psi',units='MKS',xtype='double',history=hst_restart)
		call HistoryAddVariable(varname='chi2z',dims=(/'wn','z '/),&
								&longname='theta_s',units='MKS',xtype='double',history=hst_restart)
		call HistoryAddVariable(varname='theta_s',dims=(/'wn','z '/),&
								&longname='chi2z',units='MKS',xtype='double',history=hst_restart)
	endif
	
	
	do j=1,(NN+1)*(NN+1)
		do k=1,KM+1
			chi2zval(j,k)= chi2z(j,k)
		enddo
		chi2zval(j,KM+2)=chi2z(j,KM+1)
	enddo
	
	call HistoryPut('psi',psi,hst_restart)
	call HistoryPut('chi2z',chi2zval,hst_restart)
	call HistoryPut('theta_s',theta_s,hst_restart)
	
endsubroutine 


!--------------------------------------------------------------------

subroutine conserves(theta_s,U,totalPT,totalM,totalDM)
	use constants
	use wrap_sj
	
	real(8) :: theta_s((NN+1)*(NN+1),KM+2)
	real(8) :: U(IM,JM,KM+2),U_s((NN+1)*(NN+1),KM+2)
	real(8) :: totalPT,totalM,totalDM
	
	integer :: k
	
	totalPT = 0.0
	totalM = 0.0
	
	
	!$omp parallel do
	do k=2,KM+1
		call S_G2S(U(:,:,k), U_s(:,k))
	enddo
	!omp end parallel do 
	do k=2, KM+1
		totalPT = totalPT + 2.0*theta_s(1,k)*dz
		totalM = totalM + 2.0*U_s(1,k)*dz
	enddo
	
	call S_G2S(U(:,:,1),U_s(:,1))
!	totalDM = -2.0*C*(U_s(1,2)-U_s(1,1))*indz
	totalDM = nu_V*(U_s(1,3)-2.0*U_s(1,2)+U_s(1,1))*indz*indz
	print *, totalPT,totalM,totalDM
	
	
endsubroutine

!-------------------------------------------------------------------------------------
endprogram
