program sample

  use la_module
  use dc_types, only: DP
  use gtool_history
  use lumatrix
  implicit none

 !---- 空間解像度設定 ----
  ! T21 N20
  integer,parameter  :: jm=192, km=48  ! 格子点の設定(経度, 緯度, 動径: km+1)
  integer,parameter  :: nm=126, lm=48         ! 切断波数の設定(水平, 動径)

  real(8)            :: ri=6.4D6                 ! 内半径
  real(8)            :: ro=6.415D6                ! 外半径

 !---- 変数名の決まり----
  !1, ya(1:jm,0:km): 緯度格子点, 鉛直格子点データ
  !2, yt(1:jm,0:lm): 緯度格子点, 鉛直スペクトルデータ
  !3, la(0:nm,0:km): 緯度スペクトル, 鉛直格子点データ
  !4, lt(0:nm,0:lm): 緯度スペクトル, 鉛直スペクトルデータ

 !---- 変数 ----
  real(8), dimension(0:nm,0:km) :: la_Vor  ! 東西速度
  real(8), dimension(0:nm,0:lm) :: lt_Vor   ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_Psi   ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_Psi   ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_PsiOld       ! 鉛直渦度
  real(8), dimension(0:nm,0:km) :: la_Div   ! 水平発散
  real(8), dimension(1:jm,0:km) :: ya_Xi   ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_Xi   ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_XiOld       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_Temp  ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_Temp   ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_TempOld       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_TempE  ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_TempE   ! 鉛直渦度

  real(8), dimension(1:jm,0:km) :: ya_VelLon       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_VelLat       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_VelRad       ! 鉛直渦度

  real(8), dimension(1:nm,0:km) :: la_ENG       ! 運動エネルギーの水平スペクトル
  real(8), dimension(0:nm,0:lm) :: lt_ENG       ! エネルギー
  real(8), dimension(1:jm,0:km) :: ya_Coli       ! 鉛直渦度

  real(8), dimension(1:nm) :: l_lat      ! 波数0を除いた水平波数
  real(8), dimension(0:km) :: x          ! Gauss-Lobatto 選点 (解析解)
  real(8), dimension(0:km) :: z_Rad       ! z_Rad を変数変換して -1:1 にした座標
  real(8), dimension(0:km) :: theta          ! x=\cos\theta の\theta

  real(8), dimension(1:jm,0:lm) :: yt_A       ! 鉛直渦度
  real(8), dimension(1:jm,0:lm) :: yt_B       ! 鉛直渦度
  real(8), dimension(1:jm,0:lm) :: yt_vT       ! 鉛直渦度
  real(8), dimension(1:jm,0:lm) :: yt_wT       ! 鉛直渦度

 !---- サブルーチン内で使う変数----
  real(8), dimension(1:jm,0:km) :: ya_f       ! 鉛直渦度
  real(8), dimension(1:jm,0:lm) :: yt_f       ! 鉛直渦度
  real(8), dimension(0:nm,0:km) :: la_f       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_f       ! 鉛直渦度
  real(8), dimension(1:jm,0:lm) :: yt_fx       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_fx       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_fxx       ! 鉛直渦度
  real(8), dimension(1:jm)      :: y_data       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_data       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: lt_data       ! 鉛直渦度

  real(8), dimension(1:jm,0:km) :: ya_msf       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_mu       ! 鉛直渦度
  real(8), dimension(1:jm,0:km) :: ya_cos       ! 鉛直渦度
  real(8)                       :: Int_ram
  real(8), dimension(1:jm)      :: Int_vT       ! 鉛直渦度

  real(8), dimension(0:nm,0:lm,0:lm) :: A_Psi       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm,0:lm) :: A_Xi       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm,0:lm) :: A_Temp       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: B_Psi  ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: B_Xi ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: B_Temp ! 鉛直渦度
  integer, dimension(0:nm,0:lm) :: kp_Psi             ! 鉛直渦度
  integer, dimension(0:nm,0:lm) :: kp_Xi              ! 鉛直渦度
  integer, dimension(0:nm,0:lm) :: kp_Temp             ! 鉛直渦度

  real(8), dimension(0:nm,0:lm,0:lm) :: MatrixD_sq       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm,0:lm) :: MatrixD_sqsq       ! 鉛直渦度
  real(8), dimension(0:nm,0:lm,0:lm) :: MatrixI       ! 鉛直渦度

  real(8), dimension(0:nm,0:lm) :: F_Psi     ! 非線形項
  real(8), dimension(0:nm,0:lm) :: F_PsiOld  ! 非線形項
  real(8), dimension(0:nm,0:lm) :: F_Xi  ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: F_XiOld  ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: F_Temp  ! 鉛直渦度
  real(8), dimension(0:nm,0:lm) :: F_TempOld  ! 鉛直渦度

 !---- その他 ----
  real(8), parameter :: PI = acos(-1.0d0)
  real(8), parameter :: g=9.8
  real(8), parameter :: omega=(2.d0*PI)/(8.64D4)
  real(8), parameter :: delta_h=1.d0/6.d0
  real(8), parameter :: delta_v=3.8
  real(8) :: nu_V=3.5
  real(8), parameter :: mu0=0.0
  real(8), parameter :: T0=300.0
  real(8), parameter :: const=5.d-3
  real(8), parameter :: tau=60.0*60.0*24.0*20.0

  real(8) :: initial_time=0.d0
  real(8) :: dt=120.0                          ! 時間積分刻み
  integer :: nt=720000                        ! 時間積分ステップ数
  integer :: total_day=1000                       ! 総積分時間
  integer :: ndisp=720                           ! 出力間隔ステップ数
  integer :: nprog=720                           ! 出力間隔ステップ数

  integer :: i,j,k, it=0
  real(8) :: time
  real(8) :: xl
  real(8) :: Lat
  real(8) :: z
  character(len=4), parameter :: Psi='var1'     ! 渦度:上端ノイマン,下端ディリクレ
  character(len=4), parameter :: Xi='var2'      ! 水平発散(D):上端ノイマン,下端ディリクレ
  character(len=4), parameter :: Temp='var1'      ! 水平発散(D):上端ノイマン,下端ディリクレ

 !---------------- 座標値の設定 ---------------------
  call la_Initial(nm,jm)
  xl = ro-ri
  nu_V = nu_V*((2.0/xl)**2)
  do j=0,km
    x(j) = cos(((km-j)*PI/km))
    z_Rad(j) = ((x(j)+1.0)*xl)/2.d0
    theta(j) = acos(x(j))
  end do
  do j=1,nm
    l_lat(j)=j
  end do

 !---------------- 物理係数の設定 ---------------------
  do k=0,km
  ya_Coli(:,k) = 2*omega*sin(y_Lat(:))
  ya_cos(:,k) = cos(y_Lat(:))
  ya_mu(:,k) = sin(y_Lat(:))
  end do

  call MK_Matrix(MatrixD_sq,MatrixD_sqsq,MatrixI)
  call CNABMatrixA(A_Psi,A_Xi,A_Temp, &
       & dt,nu_V,MatrixD_sq,MatrixD_sqsq,MatrixI)
  call LuDeComp(A_Psi,kp_Psi)
  call LuDeComp(A_Xi,kp_Xi)
  call LuDeComp(A_Temp,kp_Temp)

 !------------------- 初期値設定 ----------------------
  time = initial_time
  call set_initial_values(ya_Psi,ya_Xi,ya_Temp,ya_TempE)

  lt_Psi = lt_ya(ya_Psi)
  lt_Xi  = lt_ya(ya_Xi)
  lt_Temp = lt_ya(ya_Temp)
  lt_TempE = lt_ya(ya_TempE)

  call output_gtool_init
  call output_gtool                                 !初期値出力

 !------------------- 時間積分(Crank-Nicolson/Adams-Bashforth法) --------------------
 !---1ステップ目はCrank-Nicolson&Eulerで積分---
    time = initial_time + dt
   lt_PsiOld = lt_Psi
   lt_XiOld  = lt_Xi
   lt_TempOld = lt_Temp
   !--線形(拡散)項計算---
    call DiffMatrixB(lt_Psi,B_Psi,dt,nu_V,Psi)
    call DiffMatrixB(lt_Xi,B_Xi,dt,nu_V,Xi)
    call DiffMatrixB(lt_Temp,B_Temp,dt,nu_V,Temp)
   !--非線形項計算---
    call NonLinear_AB(lt_Psi,lt_Xi,lt_Temp, &
       & F_Psi,F_Xi,F_Temp)
    F_PsiOld=F_Psi
    F_XiOld=F_Xi
    F_TempOld=F_Temp
  !--右辺行列作成---
    B_Psi=B_Psi+F_Psi
    B_Xi=B_Xi+F_Xi
    B_Temp=B_Temp+F_Temp
  !--境界条件適用---
    call apply_bc(B_Psi,Psi)
    call apply_bc(B_Xi,Xi)
    call apply_bc(B_Temp,Temp)
   !---値の更新---
    lt_Psi   = LUSolve(A_Psi,kp_Psi,B_Psi)
    lt_Xi    = LUSolve(A_Xi,kp_Xi,B_Xi)
    lt_Temp  = LUSolve(A_Temp,kp_Temp,B_Temp)
  call output_gtool                                 !初期値出力
 !---2ステップ意向はCrank-Nicolson & AdamsBashforth(CNAB2)で積分---
  do it=2,nt
    time = initial_time + it * dt
 !--右辺係数行列Bを作る---
   !--線形(拡散)項計算---
    call DiffMatrixB(lt_Psi,B_Psi,dt,nu_V,Psi)
    call DiffMatrixB(lt_Xi,B_Xi,dt,nu_V,Xi)
    call DiffMatrixB(lt_Temp,B_Temp,dt,nu_V,Temp)
   !--非線形項計算---
    call NonLinear_AB(lt_Psi,lt_Xi,lt_Temp, &
       & F_Psi,F_Xi,F_Temp)
   !--右辺行列作成---
    B_Psi=B_Psi+(1.0/2.0)*(3.0*F_Psi-F_PsiOld)
    B_Xi=B_Xi+(1.0/2.0)*(3.0*F_Xi-F_XiOld)
    B_Temp=B_Temp+(1.0/2.0)*(3.0*F_Temp-F_TempOld)
   !--境界条件適用---
    call apply_bc(B_Psi,Psi)
    call apply_bc(B_Xi,Xi)
    call apply_bc(B_Temp,Temp)
 !---1ステップ前の値を保存---
    lt_PsiOld = lt_Psi
    lt_XiOld = lt_Xi
    lt_TempOld = lt_Temp
    F_PsiOld=F_Psi
    F_XiOld=F_Xi
    F_TempOld=F_Temp
 !---値の更新---
    lt_Psi   = LUSolve(A_Psi,kp_Psi,B_Psi)
    lt_Xi    = LUSolve(A_Xi,kp_Xi,B_Xi)
    lt_Temp  = LUSolve(A_Temp,kp_Temp,B_Temp)

    if(mod(it,ndisp) .eq. 0)then                  ! 出力
      call output_gtool
    endif
    if(mod(it,nprog) .eq. 0)then                  ! 出力
      write (6,*) ((it*1.0)/(nt*1.0))*total_day
    endif
  enddo

  call output_gtool_close                           ! ヒストリー後処理

contains

! 以上 メインプログラム
!--------------------------------------------------------------------------------
! 以下 サブルーチン
  subroutine CNABMatrixA(A_Psi,A_Xi,A_Temp, &
       & dt,nu,MatrixD_sq,MatrixD_sqsq,MatrixI)
   !---- 変数(入力) ----
    real(DP), dimension(0:nm,0:lm,0:lm), intent(out) :: A_Psi
    real(DP), dimension(0:nm,0:lm,0:lm), intent(out) :: A_Xi
    real(DP), dimension(0:nm,0:lm,0:lm), intent(out) :: A_Temp
    real(8), intent(in) :: dt
    real(8), intent(in) :: nu
    real(DP), dimension(0:nm,0:lm,0:lm), intent(in) :: MatrixD_sq
    real(DP), dimension(0:nm,0:lm,0:lm), intent(in) :: MatrixD_sqsq
    real(DP), dimension(0:nm,0:lm,0:lm), intent(in) :: MatrixI
    real(DP), dimension(0:lm)  :: data1
    real(DP), dimension(0:lm)  :: data2
    real(DP), dimension(0:lm)  :: data3
    !---初期化---
    A_Psi=0.0; A_Xi=0.0; A_Temp=0.0
    !---Psiの係数行列---
    A_Psi=MatrixI-((dt*nu)/2.0)*MatrixD_sq
      !---Psiの境界条件(lm-1,lm行目)---
      A_Psi(:,lm-1,:) = 0.0    ! 初期化
      A_Psi(:,lm-1,lm) = (lm**2)/2.0
        do j=1,lm-1,1
          A_Psi(:,lm-1,j) = j**2
        end do
      A_Psi(:,lm,:) = 0.0    ! 初期化
      A_Psi(:,lm,0) = 1.0/2.0    ! case1
      A_Psi(:,lm,lm) = (1.0/2.0)+(nu/const)*(xl/2.0)*((lm**2)/2.0) ! case1
        do j=1,lm-1,1
          A_Psi(:,lm,j) = ((-1.0)**j)*(1.0+(nu/const)*(xl/2.0)*(j**2))   ! case1
        end do
    !---Xiの係数行列---
    A_Xi=MatrixD_sq-((dt*nu)/2.0)*MatrixD_sqsq
      !---行の並び替え(4行下げる)---
      do k=lm,4,-1
        A_Xi(:,k,:)=A_Xi(:,k-4,:)
      end do
      !---Xiの境界条件(0行目)---
      A_Xi(:,0,:)=0.0  ! 初期化
      A_Xi(:,0,0)=1.0; A_Xi(:,0,lm)=1.0
      do j=2,lm-2,2
      A_Xi(:,0,j)=2.0
      end do
      !---Xiの境界条件(1行目)---
      A_Xi(:,1,:)=0.0  ! 初期化
      do j=1,lm-1,2
       A_Xi(:,1,j)=1.0
      end do
      !---Xiの境界条件(2,3行目)---
      data1=0.0; data2=0.0; data3=0.0
      do j=2,lm,2
       data1(j)=data1(j-2)+((j-1)**2)
      end do
      do j=3,lm-1,2
       data2(j)=data2(j-2)+((j-1)**2)
      end do
      do j=2,lm-1,1
       data1(j)=data1(j)*(2.0*j)
       data2(j)=data2(j)*(2.0*j)
      end do
       data1(lm)=data1(lm)*lm
      do j=1,lm-1
       data3(j)=((-1.0)**(j+1))*(j**2)
      end do
       data3(lm)=-((lm**2)/2.0)
      do j=0,lm
       A_Xi(:,2,j)=data3(j)-((nu/const)*(xl/2.0))*data1(j)
       A_Xi(:,3,j)=data3(j)+((2.0*nu/const)*(xl/2.0))*data2(j)
      end do
    !---Tempの係数行列---
    A_Temp=MatrixI-((dt*nu)/2.0)*MatrixD_sq
      !---Tempの境界条件(lm-1,lm行目)---
      A_Temp(:,lm-1,:)=0.0  ! 初期化
      do j=1,lm-1,2
        A_Temp(:,lm-1,j) = (j**2)
      end do
      A_Temp(:,lm,:)=0.0  ! 初期化
      do j=2,lm-2,2
        A_Temp(:,lm,j) = j**2
      end do
      A_Temp(:,lm,lm) = (lm**2)/2.0
   end subroutine

  subroutine DiffMatrixB(lt_f,B_f,dt,nu,var)
  real(8), dimension(0:nm,0:lm), intent(in) :: lt_f
  real(8), dimension(0:nm,0:lm), intent(out) :: B_f ! 係数行列
  real(8), intent(in) :: dt
  real(8), intent(in) :: nu
  character(len=4), intent(in), optional :: var
  real(8), dimension(0:nm,0:lm)  :: lt_d2f ! 係数行列
 !--初期化---
  B_f=0.0
 !--初期化---
  if (var .eq. 'var1') then      ! 渦度,Dz,温位用
   B_f=lt_f+(dt*nu/2.0)*lt_DRad2_lt(lt_DRad2_lt(lt_f))
  else if (var .eq. 'var2') then   ! \chi用
   lt_d2f=lt_DRad2_lt(lt_DRad2_lt(lt_f))
   B_f=lt_d2f+(dt*nu/2.0)*lt_DRad2_lt(lt_DRad2_lt(lt_d2f))
  else
    write (6,*) "wrong coding"
  end if
 end subroutine

  subroutine apply_bc(B_f,var)
  real(8), dimension(0:nm,0:lm), intent(inout) :: B_f ! 係数行列
  character(len=4), intent(in), optional :: var
  if (var .eq. 'var1') then
    B_f(:,lm-1)=0.0; B_f(:,lm)=0.0
  else if (var .eq. 'var2') then
   do k=lm,4,-1
    B_f(:,k)=B_f(:,k-4)
   end do
    B_f(:,0)=0.0; B_f(:,1)=0.0; B_f(:,2)=0.0; B_f(:,3)=0.0
  else
    write (6,*) "wrong coding"
  end if
  end subroutine

  subroutine MK_Matrix(MatrixD_sq,MatrixD_sqsq,MatrixI)
    real(8), dimension(0:nm,0:lm,0:lm)  :: MatrixD ! 係数行列
    real(8), dimension(0:nm,0:lm,0:lm), intent(out) :: MatrixD_sq ! 係数行列
    real(8), dimension(0:nm,0:lm,0:lm), intent(out) :: MatrixD_sqsq ! 係数行列
    real(8), dimension(0:nm,0:lm,0:lm), intent(out) :: MatrixI ! 係数行列
  !---微分演算行列---
    MatrixD=0.0
    do k=0,lm-2,1
      do j=k+1,lm-1,2
        MatrixD(:,k,j)=2.0*j
      end do
    end do
    do k=1,lm-1,2
      MatrixD(:,k,lm)=lm
    end do
  !---二階微分演算行列---
    call Square(MatrixD,MatrixD_sq)
  !---四階微分演算行列---
    call Square(MatrixD_sq,MatrixD_sqsq)
  !---単位行列---
    MatrixI=0.0
    do k=0,lm,1
        MatrixI(:,k,k)=1.0
    end do
  end subroutine

  subroutine Square(MatrixL,MatrixL_sq)
    real(8), dimension(0:nm,0:lm,0:lm), intent(in) :: MatrixL ! 係数行列
    real(8), dimension(0:nm,0:lm,0:lm), intent(out) :: MatrixL_sq ! 係数行列
    MatrixL_sq=0.0
    do k=0,lm,1
      do j=0,lm,1
        do i=0,lm,1
          MatrixL_sq(:,k,j)=MatrixL_sq(:,k,j)+MatrixL(:,k,i)*MatrixL(:,i,j)
        end do
      end do
    end do
  end subroutine

  !----時間変化項計算----
  subroutine NonLinear_AB(lt_Psi,lt_Xi,lt_Temp, &
              & lt_DPsiDt,lt_DXiDt,lt_DTDt)
   !---- 変数(入力) ----
    real(DP), dimension(0:nm,0:lm), intent(in) :: lt_Psi
    real(DP), dimension(0:nm,0:lm), intent(in) :: lt_Xi
    real(DP), dimension(0:nm,0:lm), intent(in) :: lt_Temp
   !---- 変数(出力) ----
    real(DP), dimension(0:nm,0:lm), intent(out) :: lt_DPsiDt
    real(DP), dimension(0:nm,0:lm), intent(out) :: lt_DXiDt
    real(DP), dimension(0:nm,0:lm), intent(out) :: lt_DTDt
   !----流線関数, 速度成分, 運動エネルギーなどを計算
   lt_Vor = la_Lapla_la(lt_Psi)*(1.0/ri**2)
   !---速度場,温位場を求める---
   call Pot2VelTemp(lt_Psi,lt_Xi,lt_Temp, &
              & ya_VelLon,ya_VelLat,ya_VelRad,ya_Temp)
   ya_VelLon=ya_VelLon*ya_cos        ! u をU=u\cos\phi で定義
   ya_VelLat=ya_VelLat*ya_cos        ! v をV=v\cos\phi で定義
   !---A, B を求める---
   yt_A = yt_ya((ya_lt(lt_Vor)+ya_Coli)*ya_VelLon &
         & +ya_VelRad*ya_lt(lt_DRad_lt(lt_ya(ya_VelLat))))
   yt_B = yt_ya((ya_lt(lt_Vor)+ya_Coli)*ya_VelLat &
         & -ya_VelRad*ya_lt(lt_DRad_lt(lt_ya(ya_VelLon))))
   !---エネルギーの計算---
   lt_ENG = (g/T0)*lt_Temp &
    & + lt_DRad_lt(lt_ya((ya_VelLon**2+ya_VelLat**2)/(2.0*(1.0-ya_mu**2))))
   !---uT,vT,wT の計算---
   yt_vT = yt_ya(ya_VelLat*ya_Temp)
   yt_wT = yt_ya(ya_VelRad*ya_Temp)
   !----時間変化項の見積もり----
   lt_DPsiDt = &
   & dt*(-la_LaplaInv_la(la_DivMu_ya(yt_B))*ri) 
   !& dt*(-la_LaplaInv_la(la_DivMu_ya(yt_B))*ri &
   !& +nu_h*la_Lapla_la(lt_Psi)*(1.0/ri**2))
   lt_DXiDt = &
   & dt*((xl/2.0)**2)*(-la_LaplaInv_la(lt_DRad_lt(la_DivMu_ya(yt_A)))*(ri)  &
   & -lt_ENG)
   lt_DTDt = &
   & dt*(-la_DivMu_ya(yt_vT)*(1.0/ri) &
   & -lt_DRad_lt(la_ya(yt_wT)) &
   & -(lt_Temp-lt_TempE)/tau)
  end subroutine NonLinear_AB

  !----速度場の計算----
  subroutine Pot2VelTemp(lt_Psi,lt_Xi,lt_Temp, &
             & ya_VelLon,ya_VelLat,ya_VelRad,ya_Temp)
   !---- 変数(入力) ----
    real(DP), dimension(0:nm,0:lm), intent(in) :: lt_Psi
    real(DP), dimension(0:nm,0:lm), intent(in) :: lt_Xi
    real(DP), dimension(0:nm,0:lm), intent(in) :: lt_Temp
   !---- 変数(出力) ----
    real(DP), dimension(1:jm,0:km), intent(out):: ya_VelLon
    real(DP), dimension(1:jm,0:km), intent(out):: ya_VelLat
    real(DP), dimension(1:jm,0:km), intent(out):: ya_VelRad
    real(DP), dimension(1:jm,0:km), intent(out):: ya_Temp
   !---速度成分を求める---
   ya_VelLon = -ya_GradLat_la(la_lt(lt_Psi))*(1.0/ri)
   ya_VelLat = ya_GradLat_la(la_lt(lt_DRad_lt(lt_Xi)))*(1.0/ri)
   ya_VelRad = -ya_lt(la_Lapla_la(lt_Xi))*(1.0/ri**2)
   ya_Temp = ya_lt(lt_Temp)
  end subroutine Pot2VelTemp

  !------------------- 初期値設定 ----------------------
  subroutine set_initial_values(ya_Psi,ya_Xi,ya_Temp,ya_TempE)
    real(8), dimension(1:jm,0:km), intent(out) :: ya_Psi ! 渦度
    real(8), dimension(1:jm,0:km), intent(out) :: ya_Xi  ! 水平発散
    real(8), dimension(1:jm,0:km), intent(out) :: ya_Temp  ! 温位
    real(8), dimension(1:jm,0:km), intent(out) :: ya_TempE  ! 放射平衡温位
     ya_Psi = 0.d0
     ya_Xi  = 0.d0
     ya_Temp = T0
     do j=1,jm
     do k=0,km
     ya_TempE(j,k) = T0-(T0*delta_h*((sin(y_Lat(j))-mu0)**2)) &
            & +delta_v*(z_Rad(k)/1000.0)
     end do
     end do
  end subroutine set_initial_values

   !----子午面流線関数を求める----
  subroutine calc_strfc(lt_f,ya_msf)
  real(8), dimension(0:nm,0:lm), intent(in) :: lt_f ! ポロイダルポテンシャル
  real(8), dimension(1:jm,0:km), intent(out) :: ya_msf ! 子午面流線関数
   ya_msf = ya_GradLat_la(la_lt(lt_f))
   ya_msf = -(ya_cos/ri)*ya_msf
  end subroutine

   !----相対角運動量の全球積分----
  subroutine IntLatZ_ram(ya_f,Int_ram)
  real(8), dimension(1:jm,0:km), intent(in) :: ya_f ! 1階微分
  real(8),                       intent(out) :: Int_ram !
   lt_f = lt_ya(ya_f*ya_cos)
   call VertInt_Phi(lt_f,lt_data)
   ya_data=ya_lt(lt_data)
   y_data=ya_data(:,km)-ya_data(:,0)
   Int_ram=IntLat_y(y_data)
  end subroutine

   !----南北熱フラックスの鉛直積分----
  subroutine IntZ_vT(ya_f,ya_g,Int_vT)
  real(8), dimension(1:jm,0:km), intent(in) :: ya_f ! 1階微分
  real(8), dimension(1:jm,0:km), intent(in) :: ya_g ! 1階微分
  real(8), dimension(1:jm), intent(out) :: Int_vT   ! 1階微分
  real(8), dimension(0:nm,0:lm)  :: lt_fg ! 1階微分
   !----熱フラックス----
   lt_fg=lt_ya(ya_f*ya_g)
   call VertInt(lt_fg,lt_data)
   ya_data=ya_lt(lt_data)
   Int_vT=(1.0/(xl*T0))*(ya_data(:,km)-ya_data(:,0))
  end subroutine

   !----運動エネルギースペクトルの計算----
  subroutine KE_spectrum(la_zeta,la_div,la_ENG)
  real(8), dimension(0:nm,0:km), intent(in) :: la_zeta    ! 鉛直渦度
  real(8), dimension(0:nm,0:km), intent(in) :: la_div    ! 水平発散
  real(8), dimension(1:nm,0:km), intent(out) :: la_ENG   ! エネルギースペクトル
    do j=1,nm
      la_ENG(j,:)=(1.d0/4.d0)*((ri**2)/(j*(j+1))) &
       & *((la_zeta(j,:)**2)+(la_div(j,:)**2))
    end do
  end subroutine

  subroutine VertInt_Phi(lt_fx,lt_f)
  real(8), dimension(0:nm,0:lm), intent(in) :: lt_fx ! 1階微分
  real(8), dimension(0:nm,0:lm), intent(out) :: lt_f   !
   lt_f(:,lm) = (1.0/lm)*lt_fx(:,lm-1)
   lt_f(:,lm-1) = (lt_fx(:,lm-2)-(1.0/2.0)*lt_fx(:,lm))/(2.d0*(lm-1))
    do k=lm-2,1,-1
     lt_f(:,k) = (lt_fx(:,k-1)-lt_fx(:,k+1))/(2.d0*k)
    end do
   lt_f(:,0) = 0.0
    do k=1,lm-2,1
      lt_f(:,0) = lt_f(:,0)-((lt_fx(:,k-1)-lt_fx(:,k+1))/k)*((-1)**k)
    end do
   lt_f(:,0) = lt_f(:,0)+(lt_fx(:,lm-2)-(1.0/2.0)*lt_fx(:,lm))/(lm-1) &
         & -lt_fx(:,lm-1)/lm
   lt_f = lt_f*(xl/2.0)
  end subroutine

 !------------------- 変換----------------------
  !-----三次元格子点データをスペクトルデータにする-----
  function lt_ya(ya_f)
    real(8), dimension(1:jm,0:km), intent(in) :: ya_f     ! 実空間上でのfxx
    real(8), dimension(0:nm,0:lm)  :: lt_ya     ! fxx の展開係数
   !---水平方向に球面調和関数展開する---
   la_f = la_ya(ya_f)
   do k=0,lm
    lt_ya(:,k) = (1.d0/2.d0)*(la_f(:,0)*cos(k*theta(0)) &
     & +la_f(:,km)*cos(k*theta(km)))
    do j=1,km-1
      lt_ya(:,k) = lt_ya(:,k) &
      & +(la_f(:,j)*cos(k*theta(j)))
    end do
   end do
    lt_ya = lt_ya*(2.d0/lm)
 end function lt_ya

  !-----スペクトルデータを三次元格子点データにする-----
  function ya_lt(lt_f)
    real(8), dimension(0:nm,0:lm), intent(in) :: lt_f     ! fxx の展開係数
    real(8), dimension(1:jm,0:km)  :: ya_lt     ! 実空間上でのfxx
   do j=0,km
    la_f(:,j) = (1.d0/2.d0)*(lt_f(:,0)*cos(0*theta(j)) &
     & +lt_f(:,lm)*cos(lm*theta(j)))
    do k=1,lm-1
    la_f(:,j) = la_f(:,j) &
      & +(lt_f(:,k)*cos(k*theta(j)))
    end do
   end do
   !---水平方向に球面調和関数展開する---
   ya_lt = ya_la(la_f)
 end function ya_lt

  !-----三次元格子点データを鉛直方向スペクトルデータにする-----
  function yt_ya(ya_f)
    real(8), dimension(1:jm,0:km), intent(in) :: ya_f     ! 実空間上でのfxx
    real(8), dimension(1:jm,0:lm)   :: yt_ya     ! 実空間上でのfxx
   do k=0,lm
    yt_ya(:,k) = (1.d0/2.d0)*(ya_f(:,0)*cos(k*theta(0)) &
     & +ya_f(:,km)*cos(k*theta(km)))
    do j=1,km-1
    yt_ya(:,k) = yt_ya(:,k) &
      & +(ya_f(:,j)*cos(k*theta(j)))
    end do
   end do
    yt_ya = yt_ya*(2.d0/lm)
 end function yt_ya

  function lt_la(la_f)
    real(8), dimension(0:nm,0:km), intent(in) :: la_f     ! 実空間上でのfxx
    real(8), dimension(0:nm,0:lm)  :: lt_la     ! fxx の展開係数
   do k=0,lm
    lt_la(:,k) = (1.d0/2.d0)*(la_f(:,0)*cos(k*theta(0)) &
     & +la_f(:,km)*cos(k*theta(km)))
    do j=1,km-1
    lt_la(:,k) = lt_la(:,k) &
      & +(la_f(:,j)*cos(k*theta(j)))
    end do
   end do
    lt_la = lt_la*(2.d0/lm)
 end function lt_la

  !-----スペクトルデータを三次元格子点データにする-----
  function la_lt(lt_f)
    real(8), dimension(0:nm,0:lm), intent(in) :: lt_f     ! fxx の展開係数
    real(8), dimension(0:nm,0:km)   :: la_lt     ! fxx の展開係数
   do j=0,km
    la_lt(:,j) = (1.d0/2.d0)*(lt_f(:,0)*cos(0*theta(j)) &
     & +lt_f(:,lm)*cos(lm*theta(j)))
    do k=1,lm-1
    la_lt(:,j) = la_lt(:,j) &
      & +(lt_f(:,k)*cos(k*theta(j)))
    end do
   end do
 end function la_lt

 !------------------- 鉛直微分----------------------
  !-----微分した関数のスペクトルデータを求める-----
  function lt_DRad_lt(lt_f)
    real(8), dimension(0:nm,0:lm), intent(in) :: lt_f     ! f の展開係数
    real(8), dimension(0:nm,0:lm)           :: lt_DRad_lt ! fx の展開係数
   lt_DRad_lt(:,lm) = 0.d0
   lt_DRad_lt(:,lm-1) = lm*lt_f(:,lm)
   lt_DRad_lt(:,lm-2) = 2.d0*(lm-1)*lt_f(:,lm-1)
     do k=lm-3,0,-1
      lt_DRad_lt(:,k) = lt_DRad_lt(:,k+2) &
          & +2.d0*(k+1)*lt_f(:,k+1)
     end do
   lt_DRad_lt = lt_DRad_lt*(2.d0/xl)
  end function lt_DRad_lt

  !-----微分した関数のスペクトルデータを求める-----
  function lt_DRad2_lt(lt_f)
    real(8), dimension(0:nm,0:lm), intent(in) :: lt_f     ! f の展開係数
    real(8), dimension(0:nm,0:lm)           :: lt_DRad2_lt ! fx の展開係数
   lt_DRad2_lt(:,lm) = 0.d0
   lt_DRad2_lt(:,lm-1) = lm*lt_f(:,lm)
   lt_DRad2_lt(:,lm-2) = 2.d0*(lm-1)*lt_f(:,lm-1)
     do k=lm-3,0,-1
      lt_DRad2_lt(:,k) = lt_DRad2_lt(:,k+2) &
          & +2.d0*(k+1)*lt_f(:,k+1)
     end do
  end function lt_DRad2_lt
  subroutine output_gtool_init
   !---- ヒストリーファイル作成 ----
    call HistoryCreate( &
           file='ft.nc', &
           title='Fang and Tung', &
           source='test.f90', &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'lat   ','height','t     ','n     '/), dimsizes=(/jm,km+1,0,nm/),&
           longnames=(/'Latitude             ','Height               ', &
           'time                 ','horizontal wavenumber'/), &
           units=(/'degree','km    ','day   ','1     '/),   &
           origin=real(time), interval=real(ndisp*dt), &
           xtypes=(/'double','double','double','double'/))

   !---- 座標変数定義, 出力 ----
    call HistoryPut('lat   ',y_Lat/pi*180)            ! 変数出力
    call HistoryPut('height',z_Rad/1000.0)    ! 変数出力
    call HistoryPut('n     ',l_lat)    ! 変数出力

    call HistoryAddVariable( &                                ! 変数定義
           varname='u', dims=(/'lat   ','height','t     '/), &
           longname='zonal velocity', units='m/s', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='v', dims=(/'lat   ','height','t     '/), &
           longname='meridional velocity', units='m/s', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='w', dims=(/'lat   ','height','t     '/), &
           longname='vertival velocity', units='m/s', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='T', dims=(/'lat   ','height','t     '/), &
           longname='potential temperature', units='K', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='msf', dims=(/'lat   ','height','t     '/), &
           longname='Meridional StreamFunction', units='m^2/s^2', xtype='double')
    !call HistoryAddVariable( &                                ! 変数定義
    !       varname='ram', dims=(/'t'/), &
    !       longname='RelativeAngularMomentum', units='m^2/s', xtype='double')
    !call HistoryAddVariable( &                                ! 変数定義
    !       varname='engspm', dims=(/'n     ','height','t     '/), &
    !       longname='KE spectrum', units='m^2/s^2', xtype='double')
   end subroutine output_gtool_init

  subroutine output_gtool
    call HistoryPut('t',real(time/86400.0))
   !---- 出力変数の計算----
   call Pot2VelTemp(lt_Psi,lt_Xi,lt_Temp, &
                & ya_VelLon,ya_VelLat,ya_VelRad,ya_Temp)
   !----子午面流線関数 ----
   call calc_strfc(lt_Xi,ya_msf)
   !----相対角運動量の全球積分----
   !call IntLatZ_ram(ya_VelLon,Int_ram)
   !----熱フラックス----
   !----運動エネルギースペクトル----
   !la_Vor=la_lt(la_Lapla_la(lt_Psi))*(1.0/ri**2)
   !la_Div=la_lt(lt_DRad_lt(la_Lapla_la(lt_Xi)))*(1.0/ri**2)
   !call KE_spectrum(la_Vor,la_Div,la_ENG)
   !---- 物理変数出力 ----
    call HistoryPut('u',ya_VelLon)
    call HistoryPut('v',ya_VelLat)
    call HistoryPut('w',ya_VelRad)
    call HistoryPut('T',ya_Temp)
    call HistoryPut('msf',ya_msf)
   ! call HistoryPut('ram',Int_ram)
   ! call HistoryPut('engspm',la_ENG)
  end subroutine output_gtool

  subroutine output_gtool_close
    call HistoryClose
  end subroutine output_gtool_close

end program sample
