!----------------------------------------------------------------------
!  Copyright (c) 2002--2026 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  SPMODEL 2D-CHANNEL ブシネスク対流モデル
!
! Solving 2-D Boussinesq fluid system (TT problem)
!     d\zeta/dt + J(\psi,\zeta) = PrRa dT/dx + Pr\nabla\zeta
!     dT/dt + J(\psi,T) - d\psi/dx = \nabla T
!     \nabla\psi = \zeta
!     psi = zeta = T = 0 at y=0,1
!
!履歴  2002/01/25  竹広真一
!      2002/04/15  竹広真一
!      2002/08/20  竹広真一
!      2004/01/26  小高正嗣: gt4f90io に対応
!      2005/03/15  竹広真一  xy -> yx に変更
!      2026/07/02  竹広真一  古い記述を update
!
program bsncnv0

  use et_module
  use gtool_history
  implicit none

 !---- 空間解像度設定 ----
  integer, parameter  :: im=32, jm=16      ! 格子点の設定(X,Y)
  integer, parameter  :: km=10, lm=10      ! 切断波数の設定(X,Y)

 !---- 座標変数など ----
  real(8), parameter  :: xmin=0, xmax=2     ! X 座標範囲
  real(8), parameter  :: ymin=0, ymax=1     ! Y 座標範囲

 !---- 変数 ----
  real(8)             :: yx_Psi(0:jm,0:im-1)    ! 格子データ(流線)
  real(8)             :: yx_Temp(0:jm,0:im-1)   ! 格子データ(温度/温度擾乱)
  real(8)             :: yx_Zeta(0:jm,0:im-1)   ! 格子データ(渦度)

  real(8)             :: et_Psi(-km:km,0:lm)    ! スペクトルデータ(流線)
  real(8)             :: et_Temp(-km:km,0:lm)   ! スペクトルデータ(温度擾乱)
  real(8)             :: et_Zeta(-km:km,0:lm)   ! スペクトルデータ(渦度)

  real(8)             :: x_TempBndry(0:im-1,2)  ! 境界値
  real(8)             :: e_TempBndry(-km:km,2)  ! 境界値

  real(8)             :: e_ZetaBndry(-km:km,2)  ! 境界値

 !---- 時間積分パラメター ----

  real(8), parameter :: dt=1e-4                  ! 時間ステップ間隔
  integer, parameter :: nt=5000, ndisp=100    ! 時間積分数, 表示ステップ

 !---- 物理パラメター ----
  character(len=2), parameter :: TempBC='DD'    ! 温度境界条件(DD/DN/ND/NN)
  real(8), parameter :: Ra=1.0e4                ! レイリー数
  real(8), parameter :: Pr=1.0                  ! プランドル数

 !---- その他 ----
  real(8), parameter  :: pi=3.1415926535897932385D0
  integer :: it = 0

 !---------------- 座標値の設定 ---------------------
  call et_Initial(im,jm,km,lm,xmin,xmax,ymin,ymax)

 !------------------- 初期値設定 ----------------------
  yx_Temp = 1 -yx_Y
  yx_Temp(jm/2,im/2) = yx_Temp(jm/2,im/2) + 0.01        ! 温度擾乱場
  yx_Psi  = 0.0
  yx_Zeta = 0.0

  et_Temp = et_yx(yx_Temp)
  et_Psi  = et_yx(yx_Psi)
  et_Zeta = et_yx(yx_Zeta)

 !------------------- 境界値設定 ----------------------

  x_TempBndry(:,1) = 0
  x_TempBndry(:,2) = 1
  e_TempBndry(:,1) = e_x(x_TempBndry(:,1))
  e_TempBndry(:,2) = e_x(x_TempBndry(:,2))

  e_ZetaBndry = 0

  call output_gtool_init                            ! ヒストリー初期化
  call output_gtool

 !------------------- 時間積分 ----------------------
  do it=1,nt                                         ! Euler 法による時間積分
     et_Temp = et_Temp + &
          dt*( -et_Jacobian_et_et(et_Psi,et_Temp) &
               + et_Lapla_et(et_Temp) )

     call et_Boundaries(et_Temp,e_TempBndry,cond=TempBC)

     et_Zeta = et_Zeta + &
          dt*( - et_Jacobian_et_et(et_Psi,et_Zeta) &
               + Pr*Ra*et_Dx_et(et_Temp) + Pr*et_Lapla_et(et_Zeta)    )

     call et_boundaries(et_Zeta,e_ZetaBndry,cond='DD')

     et_Psi = et_LaplaInv_et(et_Zeta)

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        yx_Zeta = yx_et(et_Zeta)
        call output_gtool
     endif
  enddo

  call output_gtool_close                           ! ヒストリー後処理
  stop

contains

    subroutine output_gtool_init
      call HistoryCreate( &                               ! ヒストリー作成
           file='bsncnv_0_et.nc', title='Boussinesq convection (TT)', &
           source='Sample program of gtool_history/gtool', &
           institution='GFD_Dennou Club davis/spmodel project',&
           dims=(/'x','y','t'/), dimsizes=(/im,jm+1,0/),    &
           longnames=(/'X-coordinate','Y-coordinate','time        '/),&
           units=(/'1','1','1'/),                           &
           origin=0.0, interval=real(ndisp*dt) )

      call HistoryPut('x',x_X)                            ! 変数出力
      call HistoryAddattr('x','topology','circular')      ! 周期属性
      call HistoryAddattr('x','modulo',xmax-xmin)         ! 周期属性
      call HistoryPut('y',y_Y)                            ! 変数出力

      call HistoryAddVariable( &                          ! 変数定義
           varname='psi', dims=(/'x','y','t'/), & 
           longname='stream function', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='zeta', dims=(/'x','y','t'/), & 
           longname='voticity', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='temp', dims=(/'x','y','t'/), & 
           longname='temperature', units='1', xtype='double')

    end subroutine output_gtool_init

    subroutine output_gtool
      yx_Temp = yx_et(et_Temp)
      yx_Psi  = yx_et(et_Psi)
      yx_Zeta = yx_et(et_Zeta)

      write(6,*) 'it = ',it
      call HistoryPut('t',real(it*dt))
      call HistoryPut('psi',transpose(yx_Psi))
      call HistoryPut('zeta',transpose(yx_Zeta))
      call HistoryPut('temp',transpose(yx_Temp))
    end subroutine output_gtool

    subroutine output_gtool_close
      call HistoryClose
    end subroutine output_gtool_close

end program bsncnv0