!スペクトルデータpsi, theta_s, chiから様々なグリッドデータに変換するプログラム
!
!使い方
!	*ソース内の設定を適宜変更する。(パラメータ、切り出し開始点・終了点など)
!	*コンパイルして実行ファイル(a.out)をつくる
!	*第１引数にスペクトルデータが入れられたnetcdfファイル(data3D.nc)、第２・３引数に出力ファイル名を指定して実行
!	(例) a.out data3D.nc g1.nc g2.nc 
!	この場合g1.ncに半整数グリッド上のデータ、g2.ncに整数グリッド上のデータが書き出される
!
!注意
!	*GEPフラックス = 一般化したEPフラックス
!
!このコードには多くのムダやムラがあります。お気づきの点があれば適宜修正して下さい。
!Mac OSX 10.5.6 Macbook Core Duo上でg95を用いると正常に動作しない不具合を確認しています。
!
!2009.3.27  山本 博基　(hiroki＠kugi.kyoto-u.ac.jp)
!
!-----------------------------------------------------------------

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


program uvwtheta
implicit none
call main 	!メインルーチン

end
!----------------------------------------------------------
subroutine main 
	use gt4_history
	use constants
	use wrap_sj
	implicit none
	!****************************適宜変更する
	integer, parameter :: initdays=0			!切り出し始め位置
	integer, parameter :: enddays=100			!切り出し終わり位置
	integer, parameter :: ddays= 1					!データ抽出間隔	
	!****************************
	real(8) ::  psi((NN+1)*(NN+1),KM+2), theta_s((NN+1)*(NN+1),KM+2)
	real(8) ::  chival((NN+1)*(NN+1),KM+2), chi((NN+1)*(NN+1),KM+1)
	real(8) :: z_half(KM),z_int(KM+1)
	
	real(8) :: theta_g(IM,JM,KM+2),theta_gval(IM,JM,KM)
	real(8) :: U(IM,JM,KM+2),  V(IM,JM,KM+2)
	real(8) :: Uval(IM,JM,KM),  Vval(IM,JM,KM), M(IM,JM,KM)
	real(8) :: w_g(IM,JM,KM+1), mpsi(IM,JM,KM+1), w_gval(IM,JM,KM)
	real(8) :: mpsival(IM,JM,KM)
	real(8) :: lon(IM), lat(JM)
	
	real(8) :: w_s((NN+1)*(NN+1),KM+2)
	real(8) :: chi1z((NN+1)*(NN+1),KM+2)
	
	real(8) :: theta_lambda(IM,JM,KM+2), theta_phi(IM,JM,KM+2)
	real(8) :: PV(IM, JM, KM)
	real(8) :: zeta(IM,JM,KM+2)
	
	!-------GEP計算用
	real(8) :: NF_phi(JM,KM+1), NF_z(JM,KM+1), divNF(JM,KM), NF_phi_phi(JM,KM+1),divNF_mean(JM,KM)
	real(8) :: NFP(JM,KM+1), NF_phi_phi_half(JM,KM), theta_z(IM,JM,KM+1)
	real(8) :: U_mean(JM,KM+1), theta_mean(JM,KM+1), w_half(IM,JM,KM+1), Uval_mean(JM,KM+1)
	real(8) :: V_mean(JM,KM+1), w_mean(JM,KM+1), Uval_mean_phi(JM,KM+1), U_half_mean(JM,KM+2)
	real(8) :: theta_phi_mean(JM,KM+1), theta_z_mean(JM,KM+1), theta_phi_int(IM,JM,KM+1)
	real(8) :: vtheta_dash(IM,JM,KM+1), wtheta_dash(IM,JM,KM+1), uv_dash(IM,JM,KM+1), uw_dash(IM,JM,KM+1)
	real(8) :: vtheta_dash_mean(JM,KM+1), wtheta_dash_mean(JM,KM+1), uv_dash_mean(JM,KM+1), uw_dash_mean(JM,KM+1)
	real(8) :: mpsi_mean(JM,KM+1), resmpsi(JM,KM+1)
	
	real(8) :: U_int(IM,JM,KM+1),V_int(IM,JM,KM+1),theta_int(IM,JM,KM+1)
	
	real(8) :: M_mean(JM,KM+1), Mv_bar(JM,KM+1), Mw_bar(JM,KM+1), Mv_dag(JM,KM+1), Mw_dag(JM,KM+1)
	real(8) :: w_dag(JM,KM+1),  v_dag(JM,KM+1)
	real(8) :: Mv_bar_phi(JM,KM+1), Mv_bar_phi_half(JM,KM)
	real(8) :: Mv_dag_phi(JM,KM+1), Mv_dag_phi_half(JM,KM)
	real(8) :: div_M_bar(JM,KM), div_M_dag(JM,KM)
	
	
	type(GT_HISTORY) :: hst_data, hst_EP
	
	
	character*80 filename,outname1,outname2, ctt
	character*5 :: time,cttt
	character*80 :: cuttime
	integer :: i,k, ct, j, days
	
	call SJ_init
	call philambda
	
	call getarg(1,filename)
	call getarg(2,outname1)
	call getarg(3,outname2)
	
	w_gval=0.0
	
	time='time='
	
	do ct=initdays, enddays, ddays
		write(ctt,*)  ct
		
		cttt = adjustl(ctt)
		cuttime = time//cttt
		
		print *, cuttime
		
		call HistoryGet(file=filename, varname='psi,', array=psi, range=cuttime)
		call HistoryGet(file=filename, varname='theta_s', array=theta_s,range=cuttime)
		call HistoryGet(file=filename, varname='chi', array=chival, range=cuttime)
		
		
		chi = chival(1:(NN+1)*(NN+1),1:(KM+1))
		
		do k=1, KM+1
			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 を求める。
			
			call S_S2G(theta_s(:,k), theta_g(:,:,k))
		enddo
		
		do i=1 , IM
			do j=1 , JM
				theta_g(i,j,KM+2) = theta_g(i,j,KM+1)			!境界条件
				theta_g(i,j,1) = theta_g(i,j,2)
				w_g(i,j,KM+1) = 0.0
				w_g(i,j,1) = 0.0
			enddo
		enddo
		theta_gval(:,:,:)=theta_g(:,:,2:KM+1)
		
		call S2GUV(psi,chi1z,U,V)   ! psi chi1z から  U V  を求める。
		
		
		do i=1 , IM
			do j=1 , JM
				do k=2, KM+1
					M(i,j,k-1) = U(i,j,k)*R + R*R*Omega*cphi(j)*cphi(j)
					
					Uval(i,j,k-1)=U(i,j,k)/cphi(j)
					Vval(i,j,k-1)=V(i,j,k)/cphi(j)
					
				enddo
			enddo
		enddo
		
		call makempsi(chi,mpsi)
		mpsival(:,:,:) = mpsi(:,:,1:KM)
		do k=1,KM+1
			do j=1,JM
				mpsi_mean(j,k) = sum(mpsi(:,j,k))/real(IM)
			enddo
		enddo
		
		! PVを計算----------
		
		call Grad(theta_s,theta_lambda,theta_phi)
		call ZETAplus(psi,zeta)
		
		PV = 0.0
		
		do k=2,KM+1
			do j=1, JM
				do i=1, IM
					
					PV(i,j,k-1) = -theta_lambda(i,j,k)*(V(i,j,k+1)-V(i,j,k-1))/(cphi(j)*2.0*dz*R) &
					& +theta_phi(i,j,k)*(U(i,j,k+1)-U(i,j,k-1))/(cphi(j)*2.0*dz*R) &
					& +(zeta(i,j,k)+2.0*Omega*sphi(j))*(theta_g(i,j,k+1)-theta_g(i,j,k-1))/(2.0*dz)
					
				enddo
			enddo
		enddo
		
		! GEPフラックスの計算---------------------------------------------------------------------
		!Uint, Vint, thetaint生成
		do k=1, KM+1
			do j=1, JM
				do i=1, IM
					U_int(i,j,k) = (U(i,j,k+1)+U(i,j,k))/2.0
					V_int(i,j,k) = (V(i,j,k+1)+V(i,j,k))/2.0
					theta_int(i,j,k) = (theta_g(i,j,k+1)+theta_g(i,j,k))/2.0
					theta_phi_int(i,j,k) = (theta_phi(i,j,k+1)+theta_phi(i,j,k))/2.0
				enddo
			enddo
		enddo
		!meanの生成
		do k=1, KM+1
			do j=1, JM
				V_mean(j,k) = sum(V_int(:,j,k))/real(IM)
				U_mean(j,k) = sum(U_int(:,j,k))/real(IM)
				Uval_mean(j,k) = U_mean(j,k)/cphi(j)
				theta_mean(j,k) = sum(theta_int(:,j,k))/real(IM)
				w_mean(j,k) = sum(w_g(:,j,k))/real(IM)
			enddo
		enddo
		do k=1, KM+2
			do j=1, JM
				U_half_mean(j,k) = sum(U(:,j,k))/real(IM)
			enddo
		enddo
		! dashの生成
		do k=1,KM+1
			do j=1,JM
				do i=1, IM
					vtheta_dash(i,j,k) = (theta_int(i,j,k)-theta_mean(j,k))*(V_int(i,j,k)-V_mean(j,k))/cphi(j)
					wtheta_dash(i,j,k) = (theta_int(i,j,k)-theta_mean(j,k))*(w_g(i,j,k)-w_mean(j,k))
					uv_dash(i,j,k) = (U_int(i,j,k)-U_mean(j,k))*(V_int(i,j,k)-V_mean(j,k))/(cphi(j)*cphi(j))
					uw_dash(i,j,k) = (U_int(i,j,k)-U_mean(j,k))*(w_g(i,j,k)-w_mean(j,k))/cphi(j)
				enddo
			enddo
		enddo
		! dashdasubarの生成
		do k=1, KM+1
			do j=1, JM
				vtheta_dash_mean(j,k) = sum(vtheta_dash(:,j,k))/real(IM)
				wtheta_dash_mean(j,k) = sum(wtheta_dash(:,j,k))/real(IM)
				uv_dash_mean(j,k) = sum(uv_dash(:,j,k))/real(IM)
				uw_dash_mean(j,k) = sum(uw_dash(:,j,k))/real(IM)
			enddo
		enddo
		!theta_phi_mean theta_z_meanの生成
		do k=1, KM+1
			do j=1, JM
				theta_phi_mean(j,k) = (sum(theta_phi_int(:,j,k)))/real(IM)
			enddo
		enddo
		do k=1, KM+1
			do j=1, JM
				do i=1, IM
					theta_z(i,j,k) = (theta_g(i,j,k+1)-theta_g(i,j,k))/dz
				enddo
				theta_z_mean(j,k) = sum(theta_z(:,j,k))/real(IM)
			enddo
		enddo
		
		!NFP生成
		do k=1, KM+1
			do j=1, JM
				!GEPフラックスのとき用
				NFP(j,k) = (vtheta_dash_mean(j,k)*theta_z_mean(j,k)-wtheta_dash_mean(j,k)*theta_phi_mean(j,k)/R)&
				&/(theta_phi_mean(j,k)*theta_phi_mean(j,k)/(R*R)+theta_z_mean(j,k)*theta_z_mean(j,k))
				!伝統的なEPフラックス用
				!NFP(j,k) = vtheta_dash_mean(j,k)/theta_z_mean(j,k)
				
			enddo
		enddo
		!print *, maxval(NFP), minval(NFP)
		!NF_phi, NF_z生成
		call divphi(Uval_mean, Uval_mean_phi)
		do j=1, JM
			do k=1,KM+1
				NF_phi(j,k) = R*cphi(j)*(((U_half_mean(j,k+1)-U_half_mean(j,k))/&
										  &(dz*cphi(j)))*NFP(j,k)-uv_dash_mean(j,k))
				
				NF_z(j,k) = R*cphi(j)*((2.0*Omega*sphi(j)-Uval_mean_phi(j,k)/R)&
									   &*NFP(j,k)-uw_dash_mean(j,k))
			enddo
			!NF_phi(j,1) = -R*cphi(j)*uv_dash_mean(j,1)
			!NF_z(j,1) = 0.0
			!NF_phi(j,KM+1) = -R*cphi(j)*uv_dash_mean(j,KM+1)
			!NF_z(j,KM+1) = 0.0
		enddo
		! div F の生成
		call divphi(NF_phi,NF_phi_phi)
		do k=1,KM
			do j=1,JM
				NF_phi_phi_half(j,k) = (NF_phi_phi(j,k+1)+NF_phi_phi(j,k))/2.0
			enddo
		enddo
		do k=1,KM
			do j=1, JM
				divNF(j,k) = NF_phi_phi_half(j,k)/R+(NF_z(j,k+1)-NF_z(j,k))/(dz)
			enddo
		enddo
		
		!残差循環
		resmpsi = mpsi_mean  - NFP(:,1:KM+1)
		
		!----------------------------------------------
		!M_mean作成
		do k=1, KM+1
			do j=1, JM
				M_mean(j,k) = U_mean(j,k)*R + Omega*R*R*cphi(j)*cphi(j)
			enddo
		enddo
		!w_dag作成
		call divphi(NFP,w_dag)
		w_dag = w_dag/R
		!v_dag作成
		!内点
		do k=2, KM
			do j=1, JM
				v_dag(j,k) = -(NFP(j,k+1)-NFP(j,k-1))/(2.0*dz)
			enddo
		enddo
		!上端
		do j=1, JM
			v_dag(j,KM+1) = -(NFP(j,KM+1)-NFP(j,KM))/(dz)
		enddo
		!下端
		do j=1,JM
			v_dag(j,1) = -0.5*(1+alpha)*(NFP(j,2)-NFP(j,1))/dz
		enddo
		!Mv_bar, Mw_bar, Mv_dag, Mw_dag作成
		do k=1, KM+1
			do j=1, JM
				Mv_bar(j,k) = M_mean(j,k)*V_mean(j,k)/cphi(j)
			enddo
		enddo
		Mw_bar = M_mean*w_mean
		Mv_dag = M_mean*v_dag
		Mw_dag = M_mean*w_dag
		!div_M_bar作成
		call divphi(Mv_bar,Mv_bar_phi)
		do k=1,KM
			do j=1,JM
				Mv_bar_phi_half(j,k) = (Mv_bar_phi(j,k+1)+Mv_bar_phi(j,k))/2.0
			enddo
		enddo
		do k=1,KM
			do j=1, JM
				div_M_bar(j,k) = Mv_bar_phi_half(j,k)/R+(Mw_bar(j,k+1)-Mw_bar(j,k))/(dz)
			enddo
		enddo
		!div_M_dag作成
		call divphi(Mv_dag,Mv_dag_phi)
		do k=1,KM
			do j=1,JM
				Mv_dag_phi_half(j,k) = (Mv_dag_phi(j,k+1)+Mv_dag_phi(j,k))/2.0
			enddo
		enddo
		do k=1,KM
			do j=1, JM
				div_M_dag(j,k) = Mv_dag_phi_half(j,k)/R+(Mw_dag(j,k+1)-Mw_dag(j,k))/(dz)
			enddo
		enddo
		
		
		if (ct==initdays) then
			
			call HistoryCreate(file=outname1,title='Data',source='DataConversion.f90',&
							   &institution='user',&
							   & dims=(/'lon ','lat ','z   ','time'/),dimsizes=(/IM,JM,KM,0/),&
							   &longnames=(/'lon ','lat ','z   ','time'/), &
							   & units=(/'deg ','deg ','m   ','days'/), origin=real(initdays),interval=1.0,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
				z_half(k) = real(k-1)*dz+dz/2.0
			enddo 
			
			call HistoryPut('lon',lon,hst_data)
			call HistoryPut('lat',lat,hst_data)
			call HistoryPut('z',z_half,hst_data)
			
			call HistoryAddVariable(varname='u',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='u',units='m/s',xtype='float',history=hst_data)
			call HistoryAddVariable(varname='v',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='v',units='m/s',xtype='float',history=hst_data)
			call HistoryAddVariable(varname='theta',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='theta',units='K',xtype='float',history=hst_data)
			
			call HistoryAddVariable(varname='M',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='angular_momentum',units='m/s',xtype='float',history=hst_data)
			
			call HistoryAddVariable(varname='PV',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='potential vorticity',units='K/ms',xtype='float',history=hst_data)
			
			
			call HistoryAddVariable(varname='epflx_div',dims=(/'lat ','z   ','time'/),&
									&longname='GEPflux_div',units='m s-2',xtype='float',history=hst_data)
			
			call HistoryAddVariable(varname='M_bar_div',dims=(/'lat ','z   ','time'/),&
									&longname='M_bar_div',units='m2 s-2',xtype='float',history=hst_data)
			call HistoryAddVariable(varname='M_dag_div',dims=(/'lat ','z   ','time'/),&
									&longname='M_dag_div',units='m2 s-2',xtype='float',history=hst_data)
			
			
			
			
			
			
			call HistoryCreate(file=outname2,title='Data',source='DataConversion.f90',&
							   &institution='user',&
							   & dims=(/'lon ','lat ','z   ','time'/),dimsizes=(/IM,JM,KM+1,0/),&
							   &longnames=(/'lon ','lat ','z   ','time'/), &
							   & units=(/'deg ','deg ','m   ','days'/), origin=real(initdays),interval=1.0,history=hst_EP)
			
			
			do k=1,KM+1
				z_int(k) = real(k-1)*dz
			enddo
			
			call HistoryPut('lon',lon,hst_EP)
			call HistoryPut('lat',lat,hst_EP)
			call HistoryPut('z',z_int,hst_EP)
			
			call HistoryAddVariable(varname='epflx_y',dims=(/'lat ','z   ','time'/),&
									&longname='GEPflux_y',units='m2 s-2',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='epflx_z',dims=(/'lat ','z   ','time'/),&
									&longname='GEPflux_z',units='m2 s-2',xtype='float',history=hst_EP)
			
			call HistoryAddVariable(varname='mpsi',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='meridional_psi',units='m2/s',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='rmpsi',dims=(/'lat ','z   ','time'/),&
									&longname='r_meridional_psi',units='m2/s',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='w',dims=(/'lon ','lat ','z   ','time'/),&
									&longname='w',units='m/s',xtype='float',history=hst_EP)
			
			call HistoryAddVariable(varname='uv_dash_mean',dims=(/'lat ','z   ','time'/),&
									&longname='uv_dash_mean',units='m2/s2',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='uw_dash_mean',dims=(/'lat ','z   ','time'/),&
									&longname='uw_dash_mean',units='m2/s2',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='vtheta_dash_mean',dims=(/'lat ','z   ','time'/),&
									&longname='vtheta_dash_mean',units='K m/s',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='wtheta_dash_mean',dims=(/'lat ','z   ','time'/),&
									&longname='wtheta_dash_mean',units='K m/s',xtype='float',history=hst_EP)
			
			call HistoryAddVariable(varname='Mv_bar',dims=(/'lat ','z   ','time'/),&
									&longname='Mv_bar',units='m3/s2',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='Mw_bar',dims=(/'lat ','z   ','time'/),&
									&longname='Mw_bar',units='m3/s2',xtype='float',history=hst_EP)
			
			call HistoryAddVariable(varname='Mv_dag',dims=(/'lat ','z   ','time'/),&
									&longname='Mv_dag',units='m3/s2',xtype='float',history=hst_EP)
			call HistoryAddVariable(varname='Mw_dag',dims=(/'lat ','z   ','time'/),&
									&longname='Mw_dag',units='m3/s2',xtype='float',history=hst_EP)
			
			
		endif
		
		
		
		call HistoryPut('u', Uval,hst_data)
		call HistoryPut('v', Vval,hst_data)
		call HistoryPut('w', w_g,hst_EP)
		call HistoryPut('theta',  theta_gval,hst_data)
		
		
		call HistoryPut('M',M,hst_data)
		call HistoryPut('mpsi',mpsi,hst_EP)
		call HistoryPut('rmpsi',resmpsi,hst_EP)
		
		call HistoryPut('PV',PV,hst_data)
		
		call HistoryPut('epflx_y',NF_phi,hst_EP)
		call HistoryPut('epflx_z',NF_z,hst_EP)
		call HistoryPut('epflx_div',divNF,hst_data)
		
		call HistoryPut('uv_dash_mean',uv_dash_mean,hst_EP)
		call HistoryPut('uw_dash_mean',uw_dash_mean,hst_EP)
		call HistoryPut('vtheta_dash_mean',vtheta_dash_mean,hst_EP)
		call HistoryPut('wtheta_dash_mean',wtheta_dash_mean,hst_EP)
		
		call HistoryPut('M_bar_div',div_M_bar,hst_data)
		call HistoryPut('M_dag_div',div_M_dag,hst_data)
		
		call HistoryPut('Mv_bar',Mv_bar,hst_EP)
		call HistoryPut('Mw_bar',Mw_bar,hst_EP)
		call HistoryPut('Mv_dag',Mv_dag,hst_EP)
		call HistoryPut('Mw_dag',Mw_dag,hst_EP)
		
	enddo
	
	call HistoryClose(hst_data)
	call HistoryClose(hst_EP)
	
	
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*(-0.33333*(sphi(j)*sphi(j)-1.0/3.0)+0.125*((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 zbibun1(x,x1)   ! スペクトルデータの整数グリッドの値 x (k=0 ~ KM)から半整数グリッドのdx/dz (k=1/2 ~ KM-1/2)をもとめる。
	use constants
	use wrap_sj
	implicit none
	real(8) :: x((NN+1)*(NN+1),KM+1)
	real(8) :: x1((NN+1)*(NN+1),KM+2)
	integer :: n,k
	
	do n=1, (NN+1)*(NN+1)
		do k=1, KM
			x1(n,k+1)=(x(n,k+1)-x(n,k))/dz
		enddo
	enddo
	
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 ZETAplus(psi,zeta) ! lap psi (- 2/a^2 psi) を求める
	use constants
	use wrap_sj
	implicit none
	real(8) :: psi((NN+1)*(NN+1),KM+2)
	real(8) :: lap_psi((NN+1)*(NN+1),KM+2), zeta(IM,JM,KM+2)
	
	integer :: i,k
	
	do k=1,KM+2
		do i=1,(NN+1)*(NN+1)
			
			lap_psi(i,k)=-psi(i,k)*(SJn(i))*(SJn(i)+1.0) 
			
		enddo
		
		call S_S2G(lap_psi(:,k),zeta(:,:,k))
		
	enddo
endsubroutine

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

subroutine lap2(x,y) ! y = -nabla^2 x つまり　y =　n(n+1)xをもとめる。(wを求める) 
	use constants
	use wrap_sj
	implicit none
	real(8) :: x((NN+1)*(NN+1),KM+1)
	real(8) :: y((NN+1)*(NN+1),KM+1)
	
	integer :: i,k
	
	do i=1,(NN+1)*(NN+1)
		
		do k=1,KM+1
			y(i,k) = x(i,k)*(SJn(i))*(SJn(i)+1.0)
			!print *, n , m
		enddo
	enddo
	
endsubroutine
!-----------------------------------------------------------------------------------

subroutine Grad(s,g1,g2)		! スペクトル　→　グリッド　勾配
	use constants
	use wrap_sj
	implicit none
	
	real(8) :: s((NN+1)*(NN+1),KM+2), slon((NN+1)*(NN+1),KM+2)
	real(8) ::  g1(IM,JM,KM+2), g2(IM,JM,KM+2)
	integer :: i,j,k
	
	
	!IPOW = 1
	!IFLAG = -1
	!call SJTS2G(NN,IM,IM,JM,JM,KM+2,s,g1,snIT,snT,snY,snIPK,snPK,snRK,snIA,snA,snQ,snWS,snWW,IPOW,IFLAG)
	!IFLAG = 1
	!call SJTS2G(NN,IM,IM,JM,JM,KM+2,s,g2,snIT,snT,snY,snIPK,snPK,snRK,snIA,snA,snQ,snWS,snWW,IPOW,IFLAG)
	do k=1,KM+2
		
		call S_S2G_dlat(s(:,k), g2(:,:,k))
		call S_S2S_dlon(s(:,k), slon(:,k))
		call S_S2G(slon(:,k),g1(:,:,k))
		
		do j=1, JM
			do i = 1, IM
				
				g1(i,j,k) = g1(i,j,k)*incphi(j)
				g2(i,j,k) = g2(i,j,k)
			enddo
		enddo
		
	enddo 
	
	
	
	
endsubroutine
!-----------------------------------------------------------------------------

subroutine divphi(g1,g2)		! 2次元グリッド　→　2次元グリッド　phi発散
	use constants
	use wrap_sj
	implicit none
	
	real(8) :: s((NN+1)*(NN+1),KM+1)
	real(8) ::  g1(JM,KM+1), g2(JM,KM+1)
	real(8) :: g31(IM,JM,KM+1), g32(IM,JM,KM+1)
	integer :: i, k
	
	do i=1, IM
		g31(i,:,:) = g1(:,:)
	enddo
	!print *, maxval(g31), minval(g31)
	do k =1 ,KM+1
		!IPOW = 1
		!IFLAG = 1
		!call SJTG2S(NN,IM,IM,JM,JM,KM+1,g31,s,snIT,snT,snY,isnIPK,isnPK,isnRK,snIA,snA,isnQ,isnWS,isnWW,IPOW,IFLAG)
		
		call S_G2S_div(g31(:,:,k), s(:,k))
		call S_S2G(s(:,k),g32(:,:,k))
		
	enddo 
	
	g2(:,:) = g32(1,:,:)
	!print *, maxval(g32), minval(g32)
	
endsubroutine
!-----------------------------------------------------------------------------


subroutine zbibun2(x,x1)   ! グリッドデータの半整数グリッドの値 x から整数グリッドのdx/dz をもとめる。
	use constants
	use wrap_sj
	implicit none
	real(8) :: x(IM,JM,KM+2)
	real(8) :: x1(IM,JM,KM+1)
	integer :: i,j,k
	
	do i=1, IM
		do j=1, JM
			do k=2, KM
				x1(i,j,k)=(x(i,j,k+1)-x(i,j,k))/dz
			enddo
		enddo
	enddo
endsubroutine
!-------------------------------------------------------------------
!--------------------------------------------------------------------

subroutine makempsi(chi,mpsi)
	
	use constants
	use wrap_sj
	implicit none
	
	real(8) :: chi((NN+1)*(NN+1),KM+1)
	real(8) ::  mpsi(IM,JM,KM+1)
	integer :: i,j,k
	
	real(8) :: s((NN+1)*(NN+1),KM+1),g(IM,JM,KM+1)
	
	!IFLAG = 1
	!IPOW = 1
	!call SJTS2G(NN,IM,IM,JM,JM,KM+1,chi,g,snIT,snT,snY,isnIPK,isnPK,isnRK,snIA,snA,isnQ,isnWS,isnWW,IPOW,IFLAG)
	
	do k =1, KM+1
		
		call S_S2G_dlat(chi(:,k), g(:,:,k))
		
		do i=1,IM
			do j=1,JM
				
				mpsi(i,j,k) = g(i,j,k)*R
				
			enddo
		enddo
	enddo
	
endsubroutine

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