!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gtool_history/gtool4 and ISPACK   2001/12/26 S.Takehiro
! 
! Solving 2-D Boussinesq fluid system (TT problem)
!     d\zeta/dt + J(\psi,\zeta) + \beta d\psi/dx 
!          = -r \zeta + \nabla^2 \zeta + \curl \tau,
!     \nabla^2\psi = \zeta,
!     \curl \tau = -\nu dU/dy,
!     psi = \zeta = 0 at y=0,1
!
  program windocean

    use c2pack
    use gtool_history
    implicit none

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

 !---- 変数 ----
    real(8)            :: gg_psi(0:jm,0:im-1)    ! 格子データ(流線)
    real(8)            :: gg_zeta(0:jm,0:im-1)   ! 格子データ(渦度)
    real(8)            :: gg_u(0:jm,0:im-1)      ! 格子データ(風)

    real(8)            :: es_psi(-km:km,lm)      ! スペクトルデータ(流線)
    real(8)            :: es_zeta(-km:km,lm)     ! スペクトルデータ(渦度)
    real(8)            :: es_cltau(-km:km,lm)    ! スペクトルデータ(応力回転)

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

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

 !---- 物理パラメター ----
    real(8), parameter :: beta=1.0               ! βパラメター
    real(8), parameter :: nu=1.0e-4              ! 水平粘性
    real(8), parameter :: edump=1.0e-2           ! エクマン収束による渦減衰
    real(8), parameter :: nlin=0.0               ! 非線形スイッチ(0 or 1)

    real(8), parameter :: pi=3.141592653589793
    integer            :: i, j, it               ! DO 変数

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

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

 !------------------- 強制設定 ----------------------
    gg_u     = -cos(2*pi*gg_x)
    es_cltau = -es_dx_es(es_gg(gg_u))

 !------------------- 初期値設定 ----------------------
    gg_psi  = 0.0
    gg_zeta = 0.0

    es_psi  = es_gg(gg_psi) ; es_zeta = es_gg(gg_zeta)

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

 !------------------- 時間積分 ----------------------
    do it=1,nt
       es_zeta = es_zeta + &
            dt*( - nlin * es_jacobian_es_es(es_psi,es_zeta) &
                 - beta * es_gg(gg_ec(ec_dy_es(es_psi)))  &
                 - edump * es_zeta + nu * es_lapla_es(es_zeta) &
                 + es_cltau  )
       es_psi = es_laplainv_es(es_zeta)

       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='windocean.nc', title='Wind driven circulation', &
           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',gg_x(1,0:im-1))                   ! 変数出力
      call HistoryPut('y',gg_y(0:jm,1))                     ! 変数出力

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

    subroutine output_gtool4
      gg_psi  = gg_es(es_psi)
      gg_zeta = gg_es(es_zeta)

      write(6,*) 'it = ',it
      call HistoryPut('psi',gg_psi)
      call HistoryPut('zeta',gg_zeta)
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program windocean
