!c Description: 
!c   2 次元非静力学モデル
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-11-11    高橋こう子 作成
!c   1.1        2003-11-12    杉山耕一朗 修正
!c   1.1        2003-11-17    杉山耕一朗 
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

program arare
  use fileset
  use gridset
  use if_avr
  use if_diff
  use if_boundary
  use gt4_history

  !--- 暗黙の型宣言禁止
  implicit none

  !--- ループを回すのに使う変数
  integer  :: i, k

  !--- 変数
  real(8), allocatable :: pi(:,:)      ! エクスナー関数
  real(8), allocatable :: u(:,:)       ! x 方向速度
  real(8), allocatable :: dpi_dz_cul(:,:)       ! x 方向速度
  real(8), allocatable :: dpi_dz_anl(:,:)       ! x 方向速度
  real(8), allocatable :: err(:,:)       ! x 方向速度
  real(8), allocatable :: w(:,:)       ! z 方向速度
  real(8), allocatable :: theta_bs(:)  ! 温位の基本場
  real(8), parameter   :: circ = 3.14159 ! 円周率
!  real(8)              :: t

  type(GT_HISTORY)  :: hist_s, hist_s_mgn ! ヒストリ名

!!!
!!!初期化
!!!
  !--- モジュールの初期化
  call fileset_init
  call gridset_init
  
  !--- 擾乱場の計算
  allocate(pi(imin:imax, kmin:kmax), &
       & u(imin:imax, kmin:kmax),     &
       & dpi_dz_cul(imin:imax, kmin:kmax), &
       & dpi_dz_anl(imin:imax, kmin:kmax), &
       & err(imin:imax, kmin:kmax), &
       & w(imin:imax, kmin:kmax),     &
       & theta_bs(kmin:kmax))

  !--- gt4f90io 初期化
  call gt4_init

  !--- 初期値設定
  do k = kmin, kmax
     pi(imin:imax, k) = sin(2 * zs(k) * circ / (zmax - zmin) ) 
  end do

  !--- 境界条件
  call boundary(2, pi)

  !--- 微分計算
  dpi_dz_cul(imin:imax, kmin:kmax) = s_dz_s( pi(imin:imax, kmin:kmax) ) 

  !-- 解析的に微分した時の関数
  do k = kmin, kmax
     dpi_dz_anl(imin:imax, k) = & 
          & 2 * circ / (zmax - zmin) * cos(2 * zs(k) * circ / (zmax - zmin) ) 
  end do

  call boundary(2, dpi_dz_anl)

  !--- 誤差
  err(imin:imax, kmin:kmax) = &       
       &   dpi_dz_cul(imin:imax, kmin:kmax)   &
       & - dpi_dz_anl(imin:imax, kmin:kmax)      

  !--- 出力
  call gt4_output

  !--- ヒストリを閉じる
  call gt4_close

  stop

  contains
    
  subroutine gt4_init

    !--- ヒストリ初期化

    call HistoryCreate( &
         & file = ncfile_s,                               &
         & title ='error check',                          &
         & source = 'arare --2D cumulus model',           &
         & institution = 'kitamo@ep.sci.hokudai.ac.jp',   &
         & dims=(/'x', 'z', 't'/),                        &
         & dimsizes=(/im, km, 0/),                      &
         & longnames=(/'x-coordinate', 'z-coordinate', 'time        '/),&
         & units=(/'m', 'm', 's'/),                       &
         & origin=0.0, interval=0.0,                      &
         & history=hist_s)                   
  
    call HistoryPut('x', xs(0:im-1), hist_s)
    call HistoryPut('z', zs(0:km-1), hist_s)
    
    call HistoryAddVariable( &  ! 変数 u を追加
         & varname='pi', dims=(/'x', 'z', 't'/),       &
         & longname='original function', units='Pa', xtype='real',   &
         & history=hist_s )
    call HistoryAddVariable( &  ! 変数 dpi_dz_cul を追加
         & varname='dpi_dz_cul', dims=(/'x', 'z', 't'/),          &
         & longname='culculated difference', units='Pa/m', xtype='real',  &
         & history=hist_s )
    call HistoryAddVariable( &  ! 変数 dpi_dz_anl を追加
         & varname='dpi_dz_anl', dims=(/'x', 'z', 't'/),          &
         & longname='analytical difference', units='Pa/m', xtype='real',  &
         & history=hist_s )
    call HistoryAddVariable( &  ! 変数 err を追加
         & varname='err', dims=(/'x', 'z', 't'/),       &
         & longname='error', units='Pa/m', xtype='real', &
         & history=hist_s )

    !--- ヒストリ初期化 (マージン含)

    call HistoryCreate( &
         & file = ncfile_s_mgn,                           &
         & title ='error check',                          &
         & source = 'arare --2D cumulus model',           &
         & institution = 'kitamo@ep.sci.hokudai.ac.jp',   &
         & dims=(/'x', 'z', 't'/),                        &
         & dimsizes=(/im+2*bm, km+2*bm, 0/),              &
         & longnames=(/'x-coordinate', 'z-coordinate', 'time        '/),&
         & units=(/'m', 'm', 's'/),                       &
         & origin=0.0, interval=0.0,                      &
         & history=hist_s_mgn)                   
  
    call HistoryPut('x', xs(imin:imax), hist_s_mgn)
    call HistoryPut('z', zs(imin:kmax), hist_s_mgn)
    
    call HistoryAddVariable( &  ! 変数 u を追加
         & varname='pi', dims=(/'x', 'z', 't'/),       &
         & longname='original function', units='Pa', xtype='real',   &
         & history=hist_s_mgn )
    call HistoryAddVariable( &  ! 変数 dpi_dz_cul を追加
         & varname='dpi_dz_cul', dims=(/'x', 'z', 't'/),          &
         & longname='culculated difference', units='Pa/m', xtype='real',  &
         & history=hist_s_mgn )
    call HistoryAddVariable( &  ! 変数 dpi_dz_anl を追加
         & varname='dpi_dz_anl', dims=(/'x', 'z', 't'/),          &
         & longname='analytical difference', units='Pa/m', xtype='real',  &
         & history=hist_s_mgn )
    call HistoryAddVariable( &  ! 変数 err を追加
         & varname='err', dims=(/'x', 'z', 't'/),       &
         & longname='error', units='Pa/m', xtype='real', &
         & history=hist_s_mgn )

  end subroutine gt4_init

!--- 変数出力

  subroutine gt4_output

    !--- マージンを含まない

    call HistoryPut('t', 0.0, hist_s)
  
    call HistoryPut('pi', pi(0:im-1, 0:km-1), hist_s)
    call HistoryPut('dpi_dz_cul', dpi_dz_cul(0:im-1, 0:km-1), hist_s)
    call HistoryPut('dpi_dz_anl', dpi_dz_anl(0:im-1, 0:km-1), hist_s)
    call HistoryPut('err', err(0:im-1, 0:km-1), hist_s)

    !--- マージンを含む

    call HistoryPut('t', 0.0, hist_s_mgn)
  
    call HistoryPut('pi', pi(imin:imax, kmin:kmax), hist_s_mgn)
    call HistoryPut('dpi_dz_cul', dpi_dz_cul(imin:imax, kmin:kmax), hist_s_mgn)
    call HistoryPut('dpi_dz_anl', dpi_dz_anl(imin:imax, kmin:kmax), hist_s_mgn)
    call HistoryPut('err', err(imin:imax, kmin:kmax), hist_s_mgn)
    call HistoryPut('t', 0.0, hist_s_mgn)
  
  end subroutine gt4_output

!--- ヒストリを閉じる

  subroutine gt4_close
    call HistoryClose(hist_s)
    call HistoryClose(hist_s_mgn)

  end subroutine gt4_close  

end program arare
