!-----------------------------------------------------------------
include "module_BPM2D.f90"
include "wrap_lj.f90"


program busi1
use constants
implicit none

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

contains

!----------------------------------------------------------
subroutine main		!メインルーチン
	use gt4_history
	use constants
	use wrap_lj
	
	implicit none
	! スペクトルデータ
	real(8) :: psi((NN+1),KM+2) , oldpsi((NN+1),KM+2) 
	real(8) :: kpsi((NN+1),KM+2,4)			!ルンゲクッタ用
	real(8) :: theta_s((NN+1),KM+2), oldtheta_s((NN+1),KM+2)
	real(8) :: ktheta_s((NN+1),KM+2,4)		!ルンゲクッタ用
	real(8) :: chi2z((NN+1),KM+1) , oldchi2z((NN+1),KM+1)
	real(8) :: kchi2z((NN+1),KM+1,4)			!ルンゲクッタ用
	
	real(8) :: chi((NN+1),KM+1)
	real(8) :: zeta_s((NN+1),KM+2)
	real(8) :: E_s((NN+1),KM+2)
	real(8) :: wtheta_s((NN+1),KM+2),Vtheta_s((NN+1),KM+2),Utheta_s((NN+1),KM+2)
	real(8) :: chi1z((NN+1),KM+2)
	real(8) :: w_s((NN+1),KM+1)
	real(8) :: theta_es((NN+1),KM+2)
	
	! グリッドデータ
	real(8) :: psigrid(JM,KM+2)
	real(8) :: chigrid(JM,KM+1)
	
	real(8) :: theta_g(JM,KM+2)
	real(8) :: U(JM,KM+2),  V(JM,KM+2)
	real(8) :: zeta_g(JM,KM+2)
	real(8) :: w_g(JM,KM+1)
	real(8) :: A_g(JM,KM+2)
	real(8) :: A1z(JM,KM+1)
	real(8) :: B_g(JM,KM+2)
	real(8) :: B1z(JM,KM+1)
	real(8) :: E_g(JM,KM+2)
	real(8) :: Utheta_g(JM,KM+2)
	real(8) :: Utheta1z(JM,KM+1)
	real(8) :: Vtheta_g(JM,KM+2)
	real(8) :: Vtheta1z(JM,KM+1)
	real(8) :: wtheta_g(JM,KM+2)
	real(8) :: rdm
	real(8) :: totalMold, totalMnew
	real(8) :: spct(NN+1)
	real(8) :: chi2zval((NN+1),KM+2)
	
	integer :: i,j,k,n,tstep,	time	! tstep =  ルンゲクッタのカウンター
	
	type(GT_HISTORY) :: hst_data1, hst_data2, hst_restart
	character(len=80) :: filename, filename2
	
	totalMold = 0.0
	
	
	call LJ_init
	call philambda
	
	
	
	call defineTE(theta_es)
	do i=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) = - Omega*((1+((real(k-1,8)-0.5)*dz)*(2.0*ga*delta_H)/&
										 &(R*R*Omega*Omega))**0.5-1)*sphi(j)
			enddo
		enddo
		psigrid(:,KM+2) = psigrid(:,KM+1)
		psigrid(:,1) = psigrid(:,2)
		do k =1, KM+2
			call L_G2S(psigrid(:,k), psi(:,k))
		enddo
		
		do k=1, KM+2
			call L_S2G(theta_es(:,k), theta_g(:,k))
			!do j =1, JM
			!call RANDOM_NUMBER(rdm)
			!theta_g(j,k) = theta_g(j,k)*(1+0.01*(rdm-0.5))
			!enddo
			call L_G2S(theta_g(:,k), theta_s(:,k))
		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 j=1,JM
				call RANDOM_NUMBER(rdm)
				!theta_g(j,k) = rdm*0.1
				theta_g(j,k) = 0.0
			enddo
			call L_G2S(theta_g(:,k),theta_s(:,k))
		enddo
		!omp end prallel do
		print *, "Initial condition was set with constant PT in rest."
	endif
	
	
	
	!print *, "Initial condition was set."
	
	!-------------------リスタートファイルの読み込み：この部分を有効にするとリスタートファイルを読む
	
	
	if(restart == 1) then
		filename = 'restart2D.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
	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(j,k) = sqrt(3.0)*cos(phi(j))*cos(lambda(i))/R
	!enddo
	!enddo
	!enddo
	
	!do k=1,KM+2
	!call L_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)
						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)
						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)
						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 L_S2S_lap(psi(:,k),zeta_s(:,k))		! psi から zeta_s を求める。
				
				call L_S2G_dual(theta_s(:,k),zeta_s(:,k),theta_g(:,k),zeta_g(:,k))   ! theta_s から theta_g ,  zeta_s から zeta_g  を求める（内点のみ）
				
				call L_S2S_lap(chi(:,k),w_s(:,k))	     ! chi から w_s を求める。
				w_s(:,k) = -w_s(:,k)
				
				call L_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)		          !境界条件
				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
				theta_g(j,KM+2) = theta_g(j,KM+1)			!境界条件
				theta_g(j,1) = theta_g(j,2)
				zeta_g(j,KM+2) = zeta_g(j,KM+1)
				zeta_g(j,1) = alpha*zeta_g(j,2)
				w_g(j,KM+1) = 0.0
				w_g(j,1) = 0.0
			enddo
			!$omp end parallel do
			
			! 非線形項の計算。（グリッド）--------------------------------------
			
			!$omp parallel do
			do k=2, KM+1
				do j=1,JM
					A_g(j,k) = (zeta_g(j,k)+2.0*Omega*sphi(j))*U(j,k) + (w_g(j,k)+w_g(j,k-1))*0.5*(V(j,k+1)-V(j,k-1))*0.5*indz
					B_g(j,k) = (zeta_g(j,k)+2.0*Omega*sphi(j))*V(j,k) - (w_g(j,k)+w_g(j,k-1))*0.5*(U(j,k+1)-U(j,k-1))*0.5*indz
				enddo
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do k=1, KM+2
				do j=1,JM
					E_g(j,k) = (U(j,k)*U(j,k)+V(j,k)*V(j,k))*0.5*incphi(j)*incphi(j)
					Utheta_g(j,k) = U(j,k)*theta_g(j,k)
					Vtheta_g(j,k) = V(j,k)*theta_g(j,k)
				enddo
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do k=1, KM+1
				do j=1,JM
					wtheta_g(j,k) = w_g(j,k)*(theta_g(j,k+1)+theta_g(j,k))*0.5
				enddo
			enddo
			!$omp end parallel do
			
			
			
			! 非線形項のスペクトル変換
			!$omp parallel do
			do k =1, KM+1
				call L_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
					A1z(j,k)=(A_g(j,k+1)-A_g(j,k))*indz
					B1z(j,k)=(B_g(j,k+1)-B_g(j,k))*indz
					Utheta1z(j,k)=(Utheta_g(j,k+1)-Utheta_g(j,k))*indz
					Vtheta1z(j,k)=(Vtheta_g(j,k+1)-Vtheta_g(j,k))*indz
				enddo
			enddo
			!$omp end parallel do
			
			!$omp parallel do
			do i=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,B_g,kpsi,tstep)
				call nextchi2z(chi2z,A1z,theta_s,E_s,kchi2z,tstep)
				call nexttheta(theta_s,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,w_g,chi,theta_s,psi,chi2z,hst_data1,hst_data2,hst_restart,TotalMnew)
					
					!open(8,FILE='psiwave.txt')	!数値粘性が不十分でないかの確認用に出力する
					!k=KM+1
					!do i=1,(NN+1)
					!spct(i) = spct(i) + psi(i,k)*psi(i,k)
					!enddo
					!do i=1, (NN)
					!write(8,*)  i, spct(i)!*(i*(i+1))
					!enddo
					!close(8)
					if(abs(totalMold - totalMnew) < 0.0001)then		!定常状態の判定。
						call HistoryClose(hst_data1)
						call HistoryClose(hst_data2)
						call HistoryClose(hst_restart)
						stop
						
						elseif(totalMnew == totalMold)then
						call HistoryClose(hst_data1)
						call HistoryClose(hst_data2)
						call HistoryClose(hst_restart)
						stop
						
						else 
						totalMold = totalMnew 
					endif
					
					
					
				endif
			endif
		enddo
		
		!$omp parallel do
		do i=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_data1)
	call HistoryClose(hst_data2)
	call HistoryClose(hst_restart)
	
endsubroutine

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

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

!---------------------------------------------------------------------------------
subroutine useLAPACK(x2,x)   !三重対角行列で表される連立1次方程式を解く
	use constants
	implicit none
	real(8) :: x2((NN+1),KM+1)
	real(8) :: x((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) 
		
		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_lj
	implicit none
	real(8) :: x((NN+1),KM+2)
	real(8) :: y((NN+1),KM+2)
	real(8) :: a(JM,KM+2)
	real(8) :: b(JM,KM+2)
	real(8) ::  g1a(JM,KM+2) , g2a(JM,KM+2)
	integer :: i , k ,j ,n
	
	
	!$omp parallel do 
	do k = 2, KM+1
		
		call L_S2G_dlat_dual(x(:,k),y(:,k),g1a(:,k),g2a(:,k))
		
		!call LJCS2Y(NN,s1a,sy1a,LJC)
		!call LJCS2Y(NN,s2a,sy2a,LJC)
		!call LJMS2G(NN,NM,NN+1,JM,sy1a,sy2a,g1a,g2a,LJIT,LJT,LJP,Q,LJR,WS1,WS2,WG,W1,W2,1)
		
		!call LJCS2X(NN,s1b,sx1b) 
		!call LJCS2X(NN,s2b,sx2b)
		!call LJMS2G(NN,NM,NN,JM,sx1b,sx2b,g1b,g2b,LJIT,LJT,LJP,Q,LJR,WS1,WS2,WG,W1,W2,0)
		
		do j=1 , JM
			a(j,k) = -g1a(j,k)*R*cphi(j)     !  Uの生成
			b(j,k) = g2a(j,k)*R*cphi(j)    !  Vの生成
		enddo
		
	enddo
	!$omp end parallel do 
	
	!$omp parallel do 
	do j=1,JM
		a(j,1)=alpha*a(j,2)			!境界条件
		a(j,KM+2)=a(j,KM+1)
		b(j,1) = alpha*b(j,2)				!境界条件
		b(j,KM+2) = b(j,KM+1)
	enddo
	!$omp end parallel do 
	
endsubroutine
!-------------------------------------------------------------------------------------
subroutine nextpsi(psi,B,kpsi,tstep)
	use constants
	use wrap_lj
	implicit none
	real(8) :: psi((NN+1),KM+2),kpsi((NN+1),KM+2,4)
	real(8) :: B(JM,KM+2)
	real(8) :: lB(JM,KM+2)
	
	integer :: tstep
	
	
	real(8) :: s2((NN+1),KM+2)
	integer :: i , k ,j 
	
	!$omp parallel do 
	do k=2 , KM+1
		do j =1, JM
			lB(j,k) = B(j,k)*incphi(j)
		enddo
		
		call L_G2S_div(lB(:,k),s2(:,k))
		
		
		do i=2 , (NN+1)
			!通常の粘性
			kpsi(i,k,tstep+1) =  ((s2(i,k))/(R*LJn(i)*(LJn(i)+1.0)) - &
								  &nu_H*(LJn(i)*(LJn(i)+1.0)-2.0)*psi(i,k)*inR*inR +&
								  &nu_V*(psi(i,k+1)-2.*psi(i,k)+psi(i,k-1))*indz*indz)*dt
			!4階の超粘性
			!kpsi(i,k,tstep+1) =  ((s2(i,k))/(R*LJn(i)*(LJn(i)+1.)) - &
			!&nu_H*(LJn(i)*(LJn(i)+1.0)*LJn(i)*(LJn(i)+1.0)-2.0*LJn(i)*(LJn(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,theta_s,E_s,kchi2z,tstep)
	use constants
	use wrap_lj
	implicit none
	real(8) :: theta_s((NN+1),KM+2),chi2z((NN+1),KM+1),E_s((NN+1),KM+2)
	real(8) :: A1z(JM,KM+1),kchi2z((NN+1),KM+1,4)
	integer :: tstep
	real(8) :: lA(JM,KM+1)
	
	real(8) :: s1((NN+1),KM+1)
	integer :: i , k ,j 
	
	
	
	!$omp parallel do
	do k=2 , KM
		
		do j =1, JM
			lA(j,k) = A1z(j,k)*incphi(j)
		enddo
		
		call L_G2S_div(lA(:,k), s1(:,k))
		
		do i=2 , (NN+1)
			!通常の粘性
			kchi2z(i,k,tstep+1) =  ((s1(i,k))/(R*LJn(i)*(LJn(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*LJn(i)*(LJn(i)+1.)-2.0)*chi2z(i,k)*inR*inR &
									&+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*LJn(i)*(LJn(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*LJn(i)*(LJn(i)+1.)*LJn(i)*(LJn(i)+1.)-2.0*LJn(i)*(LJn(i)+1.))&
			!&*chi2z(i,k)*inR*inR +&
			!&nu_V*(chi2z(i,k+1)-2.*chi2z(i,k)+chi2z(i,k-1))*indz*indz)*dt
			
		enddo
	enddo
	!$omp end parallel do
	
endsubroutine

!---------------------------------------------------------------------------------
subroutine nexttheta(theta_s,Vtheta_g,wtheta_s,theta_es,ktheta_s,tstep)
	use constants
	use wrap_lj
	implicit none
	real(8) :: wtheta_s((NN+1),KM+1)
	real(8) :: theta_s((NN+1),KM+2)
	real(8) :: ktheta_s((NN+1),KM+2,4)
	real(8) :: Vtheta_g(JM,KM+2)
	real(8) :: lV(JM,KM+2)
	
	real(8) :: theta_es((NN+1),KM+2)
	integer :: tstep
	
	real(8) :: s2((NN+1),KM+2)
	integer :: i , k ,j 
	
	
	!$omp parallel do
	do k=2 , KM+1
		do j =1, JM
			lV(j,k) = Vtheta_g(j,k)*incphi(j)
		enddo
		
		call L_G2S_div(lV(:,k), s2(:,k))
		
		do i=1 , (NN+1)
			ktheta_s(i,k,tstep+1) =  (-(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*LJn(i)*(LJn(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,w_g,chi,theta_s,psi,chi2z,hst_data1, hst_data2,hst_restart,TotalM)  !出力ルーチン
	use constants
	use gt4_history
	use wrap_lj
	
	implicit none
	real(8) :: U(JM,KM+2), lat(JM), hight(KM+2), Uval(JM,KM+2)
	real(8) :: V(JM,KM+2), Vval(JM,KM+2), M(JM,KM+2)
	real(8) :: theta_g(JM,KM+2), w_g(JM, KM+1), PV(JM,KM+2)
	real(8) :: theta_s((NN+1),KM+2), chi(NN+1, KM+1)
	real(8) :: totalPT,totalM,totalDM,mpsi(JM,KM+1)
	real(8) :: psi((NN+1),KM+2), chi2z((NN+1),KM+1), chi2zval((NN+1),KM+2)
	
	integer :: time, j,k
	real(8) :: z_half(KM+2), z_int(KM+1)
	!character :: name,fname
	type(GT_HISTORY) :: hst_data1, hst_restart, hst_data2
	
	if (time==daystep*INT(output_interval)) then
		!netcdfに書き込まれるメタ情報。適当に書き換える。
		!g95を用いるときの注意：文字列配列の要素文字列の長さをそろえる必要がある。
		call HistoryCreate(file='data2D-g1.nc',title='BPM2D',source='BPM2D.f90',institution='user',&
						   & dims=(/'lat ','z   ','time'/),dimsizes=(/JM,KM+2,0/),longnames=(/'lat ','z   ','time'/), &
						   & units=(/'deg ','m   ','days'/), origin=init_day,interval=output_interval, history=hst_data1)
		
		
		do j=1, JM 
			lat(j) = phi(j)*180.0/pi
		enddo
		do k=1 ,KM+2
			z_half(k) = dz*real(k-1,8) - 0.5*dz
		enddo
		
		call HistoryPut('lat ',lat,hst_data1)
		
		call HistoryPut('z   ',z_half,hst_data1)
		
		call HistoryAddVariable(varname='u',dims=(/'lat ','z   ','time'/),&
								&longname='u',units='m/s',xtype='float',history=hst_data1)
		call HistoryAddVariable(varname='v',dims=(/'lat ','z   ','time'/),&
								&longname='v',units='m/s',xtype='float',history=hst_data1)
		call HistoryAddVariable(varname='theta',dims=(/'lat ','z   ','time'/),&
								&longname='theta',units='K',xtype='float',history=hst_data1)
		
		call HistoryAddVariable(varname='M',dims=(/'lat ','z   ','time'/),&
								&longname='absolute_angular_momentum',units='m/s',xtype='float',history=hst_data1)
		
		call HistoryAddVariable(varname='PV',dims=(/'lat ','z   ','time'/),&
								&longname='potential vorticity',units='K/ms',xtype='float',history=hst_data1)
		
		
		call HistoryAddVariable(varname='totalPT',dims=(/'time'/),&
								&longname='totalPT',units='K',xtype='float',history=hst_data1)
		call HistoryAddVariable(varname='totalM',dims=(/'time'/),&
								&longname='totalM',units='m/s',xtype='float',history=hst_data1)
		call HistoryAddVariable(varname='totalDM',dims=(/'time'/),&
								&longname='totalDM',units='MKS',xtype='float',history=hst_data1)
		
		
		call HistoryCreate(file='data2D-g2.nc',title='BPM2D',source='BPM2D.f90',institution='user',&
						   & dims=(/'lat ','z   ','time'/),dimsizes=(/JM,KM+1,0/),longnames=(/'lat ','z   ','time'/), &
						   & units=(/'deg ','m   ','days'/), origin=init_day,interval=output_interval, history=hst_data2)
		
		do k=1,KM+1
			z_int(k) = real(k-1)*dz
		enddo
		
		call HistoryPut('lat',lat,hst_data2)
		call HistoryPut('z',z_int,hst_data2)
		
		call HistoryAddVariable(varname='w',dims=(/'lat ','z   ','time'/),&
								&longname='w',units='m/s',xtype='float',history=hst_data2)
		call HistoryAddVariable(varname='mpsi',dims=(/'lat ','z   ','time'/),&
								&longname='meridional_psi',units='m2/s',xtype='float',history=hst_data2)
		
		
		
	endif
	
	call conserves(theta_s,U,totalPT,totalM,totalDM)
	call makempsi(chi, mpsi)
	
	call HistoryPut('totalPT',totalPT,hst_data1)
	call HistoryPut('totalM',totalM,hst_data1)
	call HistoryPut('totalDM',totalDM,hst_data1)
	
	do j = 1, JM 
		Uval(j,:) = U(j,:)*incphi(j)
		M(j,:) = U(j,:)*R + R*Omega*Omega*cphi(j)+cphi(j)
		Vval(j,:) = V(j,:)*incphi(j)
	enddo
	call makePV(Uval, theta_g, theta_s, psi, PV)
	
	
	call HistoryPut('u', Uval,hst_data1)
	call HistoryPut('v', Vval,hst_data1)
	call HistoryPut('w', w_g,hst_data2)
	call HistoryPut('theta',  theta_g,hst_data1)
	call HistoryPut('M',M,hst_data1)
	call HistoryPut('mpsi',mpsi,hst_data2)
	
	call HistoryPut('PV',PV,hst_data1)
	
	!-----------------------------------------
	!リスタートファイル
	if (time==daystep*INT(output_interval)) then
		
		call HistoryCreate(file='restart2D.nc',title='BPM2D',source='BPM2D-s.f90',institution='user',&
						   & dims=(/'wn','z '/),dimsizes=(/(NN+1),KM+2/),longnames=(/'wn','z '/), &
						   & units=(/'1/m','m  '/),history=hst_restart)
		
		call HistoryPut('wn',LJn,hst_restart)
		call HistoryPut('z ',z_half,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)
		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_lj
	implicit none
	real(8) :: theta_s((NN+1),KM+2)
	real(8) :: U(JM,KM+2),U_s((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 L_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 L_G2S(U(:,1),U_s(:,1))
	totalDM = -2.0*C*(U_s(1,2)-U_s(1,1))*indz
	
	print *, totalPT,totalM,totalDM
	
endsubroutine

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

subroutine makempsi(chi,mpsi)		!子午面流線関数の計算
	use constants
	use wrap_lj
	implicit none
	
	real(8) :: chi((MM+1),KM+1)
	real(8) ::  mpsi(JM ,KM+1)
	integer :: i,j,k
	
	do k=1,KM+1
		
		call L_S2G_dlat(chi(:,k), mpsi(:,k))
		
		do j=1,JM
			mpsi(j,k) = mpsi(j,k)*R
			
		enddo
	enddo
	
endsubroutine

!-------------------------------------------------------------------------------
subroutine makePV(Uval, theta_g,theta_s,psi,PV)    ! PVの計算
	use constants
	use wrap_lj
	implicit none
	real(8) :: theta_g(JM, KM+2), theta_s(NN+1, KM+2), psi(NN+1, KM+2), PV(JM, KM+2)
	real(8) :: theta_phi(JM,KM+2), zeta_g(JM,KM+2), zeta_s(NN+1, KM+2), Uval(JM,KM+1)
	integer :: j, k
	
	
	PV = 0.0
	do k =2 , KM+1
		call L_S2G_dlat(theta_s(:,k), theta_phi(:,k))
		call L_S2S_lap(psi(:,k), zeta_s(:,k))
		call L_S2G(zeta_s(:,k), zeta_g(:,k))
		
		do j=1, JM
			PV(j,k-1) = theta_phi(j,k)*(Uval(j,k+1)-Uval(j,k-1))*0.5*indz*inR  &
			& +(zeta_g(j,k)+2.0*Omega*sphi(j))*(theta_g(j,k+1)-theta_g(j,k-1))*0.5*indz
		enddo
		
	enddo
	
	do j =1, JM
		PV(j,1) = PV(j,2)
		PV(j,KM+2) = PV(j,KM+1)
	enddo
	
endsubroutine


endprogram