program sample

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

 !---- 空間解像度設定 ----
  integer,parameter  :: jm=192, km=48  ! 格子点の設定(南北,鉛直)
  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(0:km) :: x          ! Gauss-Lobatto 選点 (-1≦x≦1)
  real(8), dimension(0:km) :: z_Rad       ! x の座標変換 (0≦z_Rad≦H)
  real(8), dimension(0:km) :: theta          ! x=\cos\theta (\pi≦\theta≦0)

  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       ! vT
  real(8), dimension(1:jm,0:lm) :: yt_wT       ! 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       ! cos(緯度)

  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  ! 1ステップ前の右辺行列
  real(8), dimension(0:nm,0:lm) :: F_Xi  ! 右辺行列
  real(8), dimension(0:nm,0:lm) :: F_XiOld  ! 1ステップ前の右辺行列
  real(8), dimension(0:nm,0:lm) :: F_Temp  ! 右辺行列
  real(8), dimension(0:nm,0:lm) :: F_TempOld  ! 1ステップ前の右辺行列

 !---- その他 ----
  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 :: mu0max=0.20 ! 加熱中心の最大緯度
  real(8) :: 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):上端ノイマン,下端ディリクレ

 !---- 初期値ファイル指定----
  character(len=100) :: initial_file='init.nc' ! 初期値データファイル名
  real(8)            :: restart_time=1000.0         ! 時刻指定 (日)

 !---- 出力ファイル指定----
  character(len=100) :: output_file='test.nc'

 !---------------- 座標値の設定 ---------------------
  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 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

  if ( initial_file /= '' ) then
     call initial_read_gtool(initial_file,restart_time, &
                      ya_VelLon,ya_VelLat,ya_VelRad,ya_Temp)
     lt_Vor=-lt_la(la_DivLat_ya(ya_VelLon))*(1.0/ri)
     lt_Psi=la_LaplaInv_la(lt_Vor)*(ri**2)
     lt_Xi=-la_LaplaInv_la(lt_ya(ya_VelRad))*(ri**2)
     lt_Temp = lt_ya(ya_Temp)
     call MK_TempE(time,lt_TempE)
  else
     call set_initial_values(ya_Psi,ya_Xi,ya_Temp)
     lt_Psi = lt_ya(ya_Psi)
     lt_Xi  = lt_ya(ya_Xi)
     lt_Temp = lt_ya(ya_Temp)
     call MK_TempE(time,lt_TempE)
  endif

  call output_gtool_init
  call output_gtool                                 !初期値出力

 !------------------- 時間積分(Crank-Nicolson/Adams-Bashforth法) --------------------
 !---1ステップ目はCrank-Nicolson&Eulerで積分---
    time = initial_time + dt
   call MK_TempE(time,lt_TempE)      ! 放射平衡温位場を与える
   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                                 !1ステップ目出力
 !---2ステップ目以降はCrank-Nicolson & AdamsBashforth(CNAB2)で積分---
  do it=2,nt
    time = initial_time + it * dt
   call MK_TempE(time,lt_TempE)      ! 放射平衡温位場を与える
 !--右辺係数行列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 ! 鉛直拡散係数! 鉛直拡散係数 ((2/H)**2 を掛けられた)
    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
    !---トロイダルポテンシャルの係数行列---
    A_Psi=MatrixI-((dt*nu)/2.0)*MatrixD_sq
      !---トロイダルポテンシャルの境界条件(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
    !---ポロイダルポテンシャルの境界条件(lm-1,lm行目)---
    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行目)---
      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 ! 鉛直拡散係数 ((2/H)**2 を掛けられた)
  character(len=4), intent(in), optional :: var ! (var1: トロイダル,温位),(var2: ポロイダル)
  real(8), dimension(0:nm,0:lm)  :: lt_d2f ! 二階微分した展開係数
 !--初期化---
  B_f=0.0
 !--場合分け---
  if (var .eq. 'var1') then      ! トロイダル,温位用
   B_f=lt_f+(dt*nu/2.0)*lt_DRad2_lt(lt_DRad2_lt(lt_f))
  else if (var .eq. 'var2') then   ! ポロイダル用
   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 ! (var1: トロイダル,温位),(var2: ポロイダル)
  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

  !----非線形項の時間変化計算(d*/dt)----
  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))))
   !---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) 
   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)
    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  ! 温位
     ya_Psi = 0.d0
     ya_Xi  = 0.d0
     ya_Temp = T0
  end subroutine set_initial_values

 !------------------- 初期値設定(ファイルから) ----------------------
  subroutine initial_read_gtool(filename,time, &
                      ya_VelLon,ya_VelLat,ya_VelRad,ya_Temp)
    character(len=*), intent(in)   :: filename
    real(8), intent(in)            :: time
    real(8), dimension(1:jm,0:km), intent(out) :: ya_VelLon !東西流
    real(8), dimension(1:jm,0:km), intent(out) :: ya_VelLat !南北流
    real(8), dimension(1:jm,0:km), intent(out) :: ya_VelRad !鉛直流
    real(8), dimension(1:jm,0:km), intent(out) :: ya_Temp   !温位

    call HistoryGet(filename,'u',ya_VelLon,range='t=1000')
    call HistoryGet(filename,'v',ya_VelLat,range='t=1000')
    call HistoryGet(filename,'w',ya_VelRad,range='t=1000')
    call HistoryGet(filename,'T',ya_Temp,range='t=1000')
  end subroutine initial_read_gtool

 !------------------- 放射平衡温位場----------------------
  subroutine MK_TempE(t,lt_TE)
  real(8), intent(in) :: t                  ! 時刻
  real(8), dimension(0:nm,0:lm), intent(out) :: lt_TE ! 放射平衡温位場
  mu0 = mu0max*sin((2.0*pi*(t/86400.0))/(360.0))
     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
  lt_TE = lt_ya(ya_TempE)
  end subroutine MK_TempE

  !----子午面流線関数を求める----
  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

 !-------------------スペクトル変換----------------------
  !-----実空間での値(緯度,高度)から展開係数(波数,波数)への変換-----
  function lt_ya(ya_f)
    real(8), dimension(1:jm,0:km), intent(in) :: ya_f     ! 実空間での値(緯度,高度)
    real(8), dimension(0:nm,0:lm)  :: lt_ya     ! 展開係数(波数,波数)
   !---水平方向に球面調和関数展開する---
   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     ! 展開係数(波数,波数)
    real(8), dimension(1:jm,0:km)  :: ya_lt     ! 実空間での値(緯度,高度)
   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  ! 実空間での値(緯度,高度) 
    real(8), dimension(1:jm,0:lm)   :: yt_ya     ! 展開係数(緯度,波数)
   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     ! 展開係数(波数,高度)
    real(8), dimension(0:nm,0:lm)  :: lt_la     ! 展開係数(波数,波数)
   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  ! 展開係数(波数,波数) 
    real(8), dimension(0:nm,0:km)   :: la_lt     ! 展開係数(波数,高度) 
   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

 !------------------- 鉛直微分----------------------
  !-----微分した関数のスペクトルデータを求める((2/H)**2 あり)-----
  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

  !-----微分した関数のスペクトルデータを求める((2/H)**2 なし)-----
  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=output_file, &
           title='Fang Tung', &
           source='test.f90', &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'lat   ','height','t     '/), dimsizes=(/jm,km+1,0/),&
           longnames=(/'Latitude','Height  ','time    '/), &
           units=(/'degree','km    ','day   '/),   &
           origin=real(time), interval=real(ndisp*dt), &
           xtypes=(/'double','double','double'/))

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

    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='vertical 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', 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 HistoryPut('u',ya_VelLon)
    call HistoryPut('v',ya_VelLat)
    call HistoryPut('w',ya_VelRad)
    call HistoryPut('T',ya_Temp)
    call HistoryPut('msf',ya_msf)
  end subroutine output_gtool

  subroutine output_gtool_close
    call HistoryClose
  end subroutine output_gtool_close

end program sample
