!--------------------------------------------------------------------------
!  Copyright (C) 2002--2005 SPMODEL Development Group. All rights reserved.
!-------------------------------------------------------------------------
! Sample program for SPML(based on g4f90io and ISPACK)   
!
!    2002/09/19 S.Takehiro
!    2004/01/26 M.Odaka
!    2004/06/17 M.Odaka
!    2005/03/16 S.Takehiro
!    2026/07/08 S.Takehiro
! 
! Solving 2-D Boussinesq fluid system
!     d\zeta/dt + J(\psi,\zeta) = -g/\rho_0 d\rho/dx + \nu\nabla\zeta + F_\zeta
!     d\rho/dt + J(\psi,\rho) + d\rho_0/dz d\psi/dx = D\nabla\rho
!     \nabla\psi = \zeta
!     psi = zeta = d\rho/dy = 0 at y=0,1
!
program igwave1

  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_Rho(0:jm,0:im-1)   ! 格子データ(密度擾乱)
  real(8)            :: yx_Zeta(0:jm,0:im-1)  ! 格子データ(渦度)

  real(8)            :: yx_ZetaSrc(0:jm,0:im-1)  ! 格子データ(渦度源)

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

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

  real(8)            :: e_RhoBndry(-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)            :: time                  ! 時間
  real(8), parameter :: dt=5e-2               ! 時間ステップ間隔
  integer, parameter :: nt=20000, ndisp=250   ! 時間積分数, 表示ステップ

  !---- 物理パラメター ----
  real(8), parameter :: Grav=9.8d0              ! 重力加速度

  real(8), parameter :: Visc=1.0d-5             ! 動粘性係数
  real(8), parameter :: Diff=1.0d-9             ! 拡散係数

  !---- 実験設定パラメター ----
  real(8), parameter :: Rho0   = 1.0d3          ! 初期の密度成層
  real(8), parameter :: DRho0dz=-1.0d0          ! 初期の密度成層
  real(8), parameter :: A0=1.0d-5               ! 渦度源の振幅
  real(8), parameter :: Sigma=2.0d-2            ! 渦度源の大きさ
  real(8), parameter :: X0 = (xmin+xmax)/2      ! 渦度源の位置(X)
  real(8), parameter :: Y0 = (ymin+ymax)/2      ! 渦度源の位置(Y)

  ! real(8), parameter :: Omega=0.3               ! 渦度源の振動数(BV との比)
  real(8)            :: Omega=0.3               ! 渦度源の振動数(BV との比)
  real(8)            :: Brunt                   ! BV 振動数

  integer            :: it = 0                  ! DO 変数

  Brunt=sqrt(Grav*abs(Drho0Dz)/Rho0)            ! BV 振動数

  write(6,*) 'Freq. of forcing/Brunt Vaisala freq. ?'
  read(5,*) Omega

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

  !------------------- 初期値設定 ----------------------
  !yx_Rho  = Rho0 + DRho0dz * yx_Y                        ! 密度成層

  yx_Rho  = 0.0
  yx_Psi  = 0.0
  yx_Zeta = 0.0

  et_RhoA  = et_yx(yx_Rho)  ; et_RhoB  = et_RhoA
  et_ZetaA = et_yx(yx_Zeta) ; et_ZetaB = et_ZetaA
  et_Psi   = et_yx(yx_Psi) 

  !------------------- 境界値設定 ----------------------
  e_RhoBndry  = 0.0D0
  e_ZetaBndry = 0.0D0
  
  !------------------- 出力設定 ----------------------

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

  !------------------- 時間積分 ----------------------
  do it=1,nt
     time = it*dt

     yx_ZetaSrc = yx_et(et_dx_et(et_yx( &
          A0 * exp(-((yx_X-X0)**2 + (yx_Y-Y0)**2)/Sigma**2) &
          * sin(Omega*Brunt*time))))

     et_RhoA = et_RhoB + &
          dt*( - et_Jacobian_et_et(et_Psi,et_RhoB) &
          - Drho0dz * et_yx(yx_et(et_dx_et(et_Psi))) &
          + Diff * et_Lapla_et(et_RhoB) )

     call et_Boundaries(et_RhoA,e_RhoBndry,cond='NN')
     
     et_ZetaA = et_ZetaB + &
          dt*( - et_Jacobian_et_et(et_Psi,et_ZetaB) &
          - Grav/Rho0 * et_yx(yx_et(et_Dx_et(et_RhoB))) &
          + Visc * et_Lapla_et(et_ZetaB) &
          + et_yx(yx_ZetaSrc) )

     call et_boundaries(et_ZetaA,e_ZetaBndry,cond='DD')
     et_Psi = et_LaplaInv_et(et_ZetaA)

     et_RhoB  = et_RhoA
     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='igwave1.nc', title='Internal gravity waves', &
         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=(/'m','m','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='m2/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
         varname='zeta', dims=(/'x','y','t'/), & 
         longname='vorticity', units='1/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
         varname='rho', dims=(/'x','y','t'/), & 
         longname='density', units='kg/m3', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
         varname='zetasrc', dims=(/'x','y','t'/), & 
         longname='vorticity forcing', units='1/s2', xtype='double')
  end subroutine output_gtool_init

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

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

  subroutine output_gtool_close
    call HistoryClose
  end subroutine output_gtool_close

end program igwave1
