!--------------------------------------------------------------------------
!  Copyright (C) 2002--2026 SPMODEL Development Group. All rights reserved.
!--------------------------------------------------------------------------
! Sample program for SPML(based on gt4f90io and ISPACK)
!
!    2002/09/24 S.Takehiro
!    2004/01/26 M.Odaka
!    2004/06/17 M.Odaka
!    2005/03/16 S.Takehiro
!    2026/07/03 S.Takehiro
! 
! Solving 2-D Boussinesq fluid system 
!     d\zeta/dt 
!        + Ud\zeta/dx - dU^2/dy^2 dpsi/dx + J(\psi,\zeta) 
!        + g/\rho_0 d\rho/dx - nu\lapla\zeta = 0,
!
!     \nabla\psi = \zeta
!
!     d\rho/dt + Ud\rho/dx + J(\psi,\rho) = D\nabla\rho
!
!     psi = zeta = 0 at y=0,1
!
! where U(y) is basic flow.
!
program kh1

  use et_module
  use gtool_history
  implicit none

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

  !---- 変数 ----
  real(8)            :: yx_Psi(0:jm,0:im-1)     ! 格子データ(流線)
  real(8)            :: yx_Zeta(0:jm,0:im-1)    ! 格子データ(渦度)
  real(8)            :: yx_Rho(0:jm,0:im-1)     ! 格子データ(密度)
  real(8)            :: yx_U(0:jm,0:im-1)       ! 格子データ(基本流)
  real(8)            :: yx_dU2dy2(0:jm,0:im-1)  ! 格子データ(基本流)

  real(8)            :: et_PsiA(-km:km,0:lm)    ! スペクトルデータ(流線,t+1)
  real(8)            :: et_ZetaA(-km:km,0:lm)   ! スペクトルデータ(渦度,t+1)
  real(8)            :: et_RhoA(-km:km,0:lm)    ! スペクトルデータ(密度,t+1)

  real(8)            :: et_PsiB(-km:km,0:lm)    ! スペクトルデータ(流線,t)
  real(8)            :: et_ZetaB(-km:km,0:lm)   ! スペクトルデータ(渦度,t)
  real(8)            :: et_RhoB(-km:km,0:lm)    ! スペクトルデータ(密度,t)

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

  !---- 時間積分パラメター ----
  real(8), parameter :: dt=0.25d-4               ! 時間ステップ間隔
  integer, parameter :: nt=50000, ndisp=1000     ! 時間積分数, 表示ステップ

  !---- 物理パラメター ----
  real(8), parameter :: Grav=9.8d2               ! 重力加速度
  real(8), parameter :: Visc=1.0d-1              ! 動粘性係数
  real(8), parameter :: Diff=1.0d-5              ! 拡散係数

  real(8), parameter :: U0 = 3.0d1               ! 基本流の分布(速度差)
  real(8), parameter :: A0 = 0.2                 ! 基本流の分布(シアー層の幅)
  real(8), parameter :: EPS= 1e-6                ! 初期擾乱の速度振幅

  real(8), parameter :: deltaRho = 0.05          ! 密度差
  real(8), parameter :: Rho0     = 1.0d0         ! 下層の密度

  !--- 超粘性 ---
  real(8), parameter :: nu=1e-10                 ! 粘性係数のfactor
  integer, parameter :: nv=5                     ! 粘性の階数
  real(8)            :: et_sVisc(-km:km,0:lm)    ! スペクトルデータ(粘性係数)
  real(8)            :: et_sDiff(-km:km,0:lm)    ! スペクトルデータ(粘性係数)

  real(8), parameter :: pi = 3.141592653589793   ! 円周率

  integer            :: i, j, k, l, it=0         ! DO 変数
  real               :: rnum

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

  !------------------- 初期値設定 ----------------------

  ! 流速分布
  yx_U      = U0 * tanh((yx_Y-(ymin+ymax)/2)/A0)
  !yx_dU2dY2 = yx_et(et_dy_et(et_dy_et(et_yx(yx_U))))     ! 境界で dU0/dy=0
  yx_dU2dY2 = -2*U0/A0**2 * sinh((yx_Y-(ymin+ymax)/2)/A0) &
             / cosh((yx_Y-(ymin+ymax)/2)/A0)**3 

  !yx_Psi = EPS * cos(2*pi/(xmax-xmin)*yx_X) &
  !             * cos(2*pi/(ymax-ymin)*yx_Y)
  !yx_Psi = 0.0 ; yx_Psi(jm/2,im/2) = EPS
  do j=0,jm
     do i=0,im-1
        call RANDOM_NUMBER(rnum)
        yx_Psi(j,i)=EPS * rnum
     enddo
  enddo

  et_PsiA  = et_yx(yx_Psi)        ;  yx_Psi = yx_et(et_PsiA) 
  et_ZetaA = et_Lapla_et(et_PsiA) ;  yx_Zeta  = yx_et(et_ZetaA)


  ! 密度分布
  yx_Rho    = Rho0 - deltaRho * (tanh((yx_Y-(ymin+ymax)/2)/A0) + 1)/2
  et_RhoA    = et_yx(yx_Rho)

  et_ZetaB = et_ZetaA
  et_PsiB  = et_PsiA
  et_RhoB  = et_RhoA

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

  !---------------------- 超粘性設定 ----------------------
  do l=1,lm
     do k=-km,km
        et_sVisc(k,l)= Visc * nu &
             * ((2*pi*k/(xmax-xmin))**2+(2*pi*l/(ymax-ymin))**2)**nv
     enddo
  enddo

  do l=0,lm
     do k=-km,km
        et_sDIff(k,l)= Diff * nu &
             * ((2*pi*k/(xmax-xmin))**2+(2*pi*l/(ymax-ymin))**2)**nv
     enddo
  enddo

  !---------------------- 時間積分 ----------------------
  do it=1,nt
     et_RhoA = ( et_RhoB + &
          dt*( - et_Jacobian_et_et(et_PsiB,et_RhoB)       &
               - et_yx(yx_U*yx_et(et_dx_et(et_RhoB)))    &
               + Diff * et_Lapla_et(et_RhoB) )           &
              )/(1+et_sDiff*dt)

     et_ZetaA = ( et_ZetaB +                                 &
          dt*( - et_Jacobian_et_et(et_PsiB,et_ZetaB)       &
               - et_yx(yx_U*yx_et(et_dx_et(et_ZetaB)))    &
               + et_yx(yx_dU2dy2*yx_et(et_dx_et(et_PsiB)))&
               - Grav/Rho0 * et_yx(yx_et(et_Dx_et(et_RhoB))) &
               + Visc * et_lapla_et(et_ZetaB) )           &
          )/(1+et_sVisc*dt)

     et_PsiA = et_LaplaInv_et(et_ZetaA)

     et_ZetaB = et_ZetaA
     et_PsiB  = et_PsiA
     et_RhoB  = et_RhoA

     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='kh1-et.nc', title='K-H instability', &
         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=(/'cm','cm','s '/),                           &
         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='vorticity', units='1', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
         varname='rho', dims=(/'x','y','t'/), & 
         longname='density', units='g/cm3', xtype='double')

    call HistoryAddVariable( &                          ! 変数定義
         varname='u0', dims=(/'y'/), & 
         longname='basic flow', units='1', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
         varname='d2u0dy2', dims=(/'y'/), & 
         longname='d2U0/dy', units='1', xtype='double')
    call HistoryPut('u0',yx_U(:,1))
    call HistoryPut('d2u0dy2',yx_dU2dy2(:,1))
  end subroutine output_gtool_init

  subroutine output_gtool
    yx_Psi  = yx_et(et_PsiA)
    yx_Zeta = yx_et(et_ZetaA)
    yx_Rho = yx_et(et_RhoA)

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

  subroutine output_gtool_close
    call HistoryClose
  end subroutine output_gtool_close

end program
