!----------------------------------------------------------------------
!  Copyright (c) 2002--2026 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! Sample program for SPML(based on gt4f90io and ISPACK)
!
!   2002/08/19 S.Takehiro
!   2004/01/26 M.Odaka
!   2004/06/03 M.Odaka
!   2004/06/17 M.Odaka
!   2005/03/16 S.Takehiro
!   2026/07/02 S.Takehiro  esc_module => et_module
! 
! Solving 2-D Boussinesq fluid system (TT problem)
!     d\zeta/dt + J(\psi,\zeta) = PrRt dT/dx + PrRc dC/dx + Pr\nabla\zeta
!     dT/dt + J(\psi,T) - d\psi/dx = \nabla T
!     dC/dt + J(\psi,C) - d\psi/dx = Tau \nabla C
!     \nabla\psi = \zeta
!     psi = zeta = T = C = 0 at y=0,1
!
program ddfcnv0

  use et_module
  use gtool_history
  implicit none

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

  !---- 変数 ----
  real(8)            :: yx_Psi(0:jm,0:im-1)    ! 格子データ(流線)
  real(8)            :: yx_Temp(0:jm,0:im-1)   ! 格子データ(温度/温度擾乱)
  real(8)            :: yx_Comp(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_TempA(-km:km,0:lm)  ! スペクトルデータ(温度擾乱,t+1)
  real(8)            :: et_CompA(-km:km,0:lm)  ! スペクトルデータ(組成擾乱,t+1)
  real(8)            :: et_ZetaA(-km:km,0:lm)  ! スペクトルデータ(渦度,t+1)

  real(8)            :: et_TempB(-km:km,0:lm)  ! スペクトルデータ(温度擾乱,t)
  real(8)            :: et_CompB(-km:km,0:lm)  ! スペクトルデータ(組成擾乱,t)
  real(8)            :: et_ZetaB(-km:km,0:lm)  ! スペクトルデータ(渦度,t)

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

  !---- 座標変数など ----
  real(8), parameter :: xmin=0.0, xmax=2.0     ! 領域範囲(X 方向) 
  real(8), parameter :: ymin=0.0, ymax=1.0     ! 領域範囲(Y 方向) 

  !---- 時間積分パラメター ----
  real(8), parameter :: dt=5d-6                ! 時間ステップ間隔
  ! integer, parameter :: nt=5000, ndisp=250    ! 時間積分数, 表示ステップ
  integer, parameter :: nt=200000, ndisp=4000    ! 時間積分数, 表示ステップ

  !---- 物理パラメター ----
  real(8), parameter :: Ra=-1.2e4              ! 温度レイリー数
  real(8), parameter :: Rc= 1.0e4              ! 組成レイリー数
  real(8), parameter :: Pr=1.0                 ! プランドル数
  real(8), parameter :: Tau=0.1                ! 組成拡散/熱拡散の比

  integer            :: it=0                   ! DO 変数

  !---------------- 座標値の設定 ---------------------
  call et_initial(im,jm,km,lm,xmin,xmax,ymin,ymax)    ! ISPACK初期化

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

  et_Psi   = et_yx(yx_Psi)
  et_TempA = et_yx(yx_Temp) ; et_TempB = et_TempA
  et_CompA = et_yx(yx_Comp) ; et_CompB = et_CompA
  et_ZetaA = et_yx(yx_Zeta) ; et_ZetaB = et_ZetaA

  e_TempBndry = 0.0D0
  e_CompBndry = 0.0D0

  e_ZetaBndry = 0.0D0

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

  !------------------- 時間積分 ----------------------
  do it=1,nt
     et_TempA = et_TempB + &
          dt*( - et_Jacobian_et_et(et_Psi,et_TempB) &
               + et_Dx_et(et_Psi) &
               + et_Lapla_et(et_TempB) )
     et_CompA = et_CompB + &
          dt*( - et_Jacobian_et_et(et_Psi,et_CompB) &
               + et_Dx_et(et_Psi) &
               + Tau*et_Lapla_et(et_CompB) )
     et_ZetaA = et_ZetaB + &
          dt*( - et_Jacobian_et_et(et_Psi,et_ZetaB) &
               + Pr*Ra*et_Dx_et(et_TempB) &
               + Pr*Rc*et_Dx_et(et_CompB) &
               + Pr*et_Lapla_et(et_ZetaB)    )

     call et_Boundaries(et_TempA,e_TempBndry,cond='DD')
     call et_Boundaries(et_CompA,e_CompBndry,cond='DD')
     call et_boundaries(et_ZetaA,e_ZetaBndry,cond='DD')
     et_Psi   = et_LaplaInv_et(et_ZetAa)

     et_TempB = et_TempA
     et_CompB = et_CompA
     et_ZetaB = et_ZetaA

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

  call output_gtool_close
  stop

contains
  subroutine output_gtool_init
    call HistoryCreate( &                               ! ヒストリー作成
         file='ddfcnv-et-0.nc', title='Double diffusive convection', &
         source='Sample program of spmodel', &
         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')
    call HistoryAddVariable( &                          ! 変数定義
         varname='comp', dims=(/'x','y','t'/), & 
         longname='composition', units='1', xtype='double')

  end subroutine output_gtool_init

  subroutine output_gtool
    if ( Ra .ge. 0 )then
       yx_Temp = yx_et(et_TempA) + 1-yx_Y       ! 出力は全温度場
    else
       yx_Temp = yx_et(et_TempA) + yx_Y         ! 出力は全温度場
    endif
    if ( Rc .ge. 0 ) then
       yx_Comp = yx_et(et_CompA) + 1-yx_Y       ! 出力は全組成場
    else
       yx_Comp = yx_et(et_CompA) -yx_Y          ! 出力は全組成場
    endif
    yx_Psi  = yx_et(et_Psi)
    yx_Zeta = yx_et(et_ZetaA)

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

  subroutine output_gtool_close
    call HistoryClose
  end subroutine output_gtool_close

end program
