! Sample program for gtool_history/gtool4 and ISPACK   2001/11/29 S.Takehiro
! 
! Solving 2-D Boussinesq fluid system (FF 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 = dT/dy = 0 at y=0,1
!
  program bsncnv2

    use c2pack
    use gtool_history
    implicit none

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

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

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

 !---- 座標変数など ----
    real(8)            :: gx(0:jm,0:im-1)       ! 格子データ(X)
    real(8)            :: gy(0:jm,0:im-1)       ! 格子データ(Y)
    real(8), parameter :: xmin=0.0, xmax=8.0    ! 領域範囲(X 方向) 
    real(8), parameter :: ymin=0.0, ymax=1.0    ! 領域範囲(Y 方向) 

 !---- 時間積分パラメター ----
    real(8), parameter :: dt=1e-3               ! 時間ステップ間隔
    integer, parameter :: nt=10000, ndisp=500  ! 時間積分数, 表示ステップ

 !---- 物理パラメター ----
    real(8), parameter :: Ra=1.0e4              ! レイリー数
    real(8), parameter :: Pr=1.0                ! プランドル数

    integer            :: i, j, it              ! DO 変数

 !---------------- 座標値の設定 ---------------------
    do i=0,im-1
       do j=0,jm
          gx(j,i)= (xmax-xmin)/im*i + xmin
          gy(j,i)= (ymax-ymin)/jm*j + ymin
       enddo
    enddo

    call c2initial(im,jm,km,lm,xmax-xmin,ymax-ymin)      ! ISPACK初期化

 !------------------- 初期値設定 ----------------------
    gtemp = 0.0 ;  gtemp(jm/2,im/2) = 0.01               ! 温度擾乱場
    gpsi  = 0.0
    gzeta = 0.0

    ctemp = c_g(gtemp) ; spsi  = s_g(gpsi) ; szeta = s_g(gzeta)

    call output_gtool4_init                              ! ヒストリー初期化
    call output_gtool4                                   ! 初期値出力

 !------------------- 時間積分 ----------------------
    do it=1,nt
       ctemp = ctemp + &
            dt*( -c_jac_sc(spsi,ctemp) + c_g(g_s(s_dx_s(spsi))) + c_lapla_c(ctemp) )
       szeta = szeta + &
            dt*( - s_jac_ss(spsi,szeta) &
                 + Pr*Ra*s_g(g_c(c_dx_c(ctemp))) + Pr*s_lapla_s(szeta)    )
       spsi = s_laplainv_s(szeta)

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

    call output_gtool4_close
    stop

  contains
    subroutine output_gtool4_init
      call HistoryCreate( &                               ! ヒストリー作成
           file='bsncnv_2.nc', title='Boussinesq convection (FF)', &
           source='Sample program of gtool_history/gtool4', &
           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',gx(1,0:im-1))                   ! 変数出力
      call HistoryPut('y',gy(0:jm,1))                     ! 変数出力

      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_gtool4_init

    subroutine output_gtool4
      gtemp = g_c(ctemp) + 1 - gy                         ! 出力は全温度場
      gpsi  = g_s(spsi)
      gzeta = g_s(szeta)

      write(6,*) 'it = ',it
      call HistoryPut('psi',transpose(gpsi))
      call HistoryPut('zeta',transpose(gzeta))
      call HistoryPut('temp',transpose(gtemp))
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program bsncnv2
