!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 :: u_avr(:,:)       ! x 方向速度
  real(8), allocatable :: u_err(:,:)       ! x 方向速度
  real(8), allocatable :: w(:,:)       ! z 方向速度
  real(8), allocatable :: w_avr(:,:)       ! x 方向速度
  real(8), allocatable :: w_err(:,:)       ! x 方向速度
  real(8), allocatable :: theta_bs(:)  ! 温位の基本場
  real(8), parameter   :: circ = 3.14159 ! 円周率
!  real(8)              :: t
  integer              :: avrloop  ! 平均繰り返し回数

  type(GT_HISTORY)  :: hist_u, hist_w, hist_u_mgn, hist_w_mgn ! ヒストリ名

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

  !--- gt4f90io 初期化
  call gt4_init

  !--- 初期値設定
  do i = imin, imax
     do k = kmin, kmax
        u(i,k) = cos(2 * x(i) * circ / (xmax - xmin) )    &
             &   * cos(2 * zs(k) * circ / (zmax - zmin) )
        w(i,k) = cos(2 * xs(i) * circ / (xmax - xmin) )   &
             &   * cos(2 * z(k) * circ / (zmax - zmin) )
     end do
  end do

  !--- 境界条件
  call boundary(3, u)
  call boundary(3, w)

  !--- 平均の計算
  avrloop = 5     ! 繰り返し平均をとる回数

  u_avr(imin:imax, kmin:kmax) = u_avr_s(s_avr_u(u(imin:imax, kmin:kmax)))
  w_avr(imin:imax, kmin:kmax) = w_avr_s(s_avr_w(w(imin:imax, kmin:kmax)))

  if (avrloop > 2) then 
     do i = 1, avrloop-1
        u_avr(imin:imax, kmin:kmax) =  &
             & u_avr_s(s_avr_u(u_avr(imin:imax, kmin:kmax)))
        w_avr(imin:imax, kmin:kmax) =  &
             & w_avr_s(s_avr_w(w_avr(imin:imax, kmin:kmax)))
     end do
  end if

  !--- 標準出力で値を書き出す
!  do i = imin, imax
!     write(*,*) 'i=',i
!     do k = kmin, kmax
!        write(*,*) u(i,k)
!     end do
!  end do

  !--- 誤差
  u_err(imin:imax, kmin:kmax) = &       
       &   u_avr(imin:imax, kmin:kmax)   &
       & - u(imin:imax, kmin:kmax)      
  w_err(imin:imax, kmin:kmax) = &       
       &   w_avr(imin:imax, kmin:kmax)   &
       & - w(imin:imax, kmin:kmax)   

  !--- 標準出力で値を書き出す
!  do i = 0, im-1
!     write(*,*) 'i=',i
!     do k = 0, km-1
!        write(*,*) u_err(i,k)
!     end do
!  end do

  !--- 出力
  call gt4_output

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

  stop

  contains
    
  subroutine gt4_init

    !--- u のヒストリ初期化

    call HistoryCreate( &
         & file = ncfile_u,                               &
         & 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_u)                   
  
    call HistoryPut('x', x(0:im-1), hist_u)
    call HistoryPut('z', zs(0:km-1), hist_u)
    
    call HistoryAddVariable( &  ! 変数 u を追加
         & varname='u', dims=(/'x', 'z', 't'/),       &
         & longname='u', units='m/s', xtype='real',   &
         & history=hist_u )
    call HistoryAddVariable( &  ! 変数 u_avr を追加
         & varname='u_avr', dims=(/'x', 'z', 't'/),          &
         & longname='u average', units='m/s', xtype='real',  &
         & history=hist_u )
    call HistoryAddVariable( &  ! 変数 u_err を追加
         & varname='u_err', dims=(/'x', 'z', 't'/),       &
         & longname='u error', units='m/s', xtype='real', &
         & history=hist_u )

    !--- w のヒストリ初期化

    call HistoryCreate( &
         & file = ncfile_w,                               &
         & 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_w)                   
  
    call HistoryPut('x', xs(0:im-1), hist_w)
    call HistoryPut('z', z(0:km-1), hist_w)
    
    call HistoryAddVariable( &  ! 変数 w を追加
         & varname='w', dims=(/'x', 'z', 't'/),       &
         & longname='w', units='m/s', xtype='real',   &
         & history=hist_w )
    call HistoryAddVariable( &  ! 変数 w_avr を追加
         & varname='w_avr', dims=(/'x', 'z', 't'/),          &
         & longname='w average', units='m/s', xtype='real',  &
         & history=hist_w )
    call HistoryAddVariable( &  ! 変数 w_err を追加
         & varname='w_err', dims=(/'x', 'z', 't'/),       &
         & longname='w error', units='m/s', xtype='real', &
         & history=hist_w )

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

    call HistoryCreate( &
         & file = ncfile_u_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_u_mgn)          
  
    call HistoryPut('x', x(imin:imax), hist_u_mgn)
    call HistoryPut('z', zs(kmin:kmax), hist_u_mgn)
    
    call HistoryAddVariable( &  ! 変数 u を追加
         & varname='u', dims=(/'x', 'z', 't'/),       &
         & longname='u', units='m/s', xtype='real',   &
         & history=hist_u_mgn )
    call HistoryAddVariable( &  ! 変数 u_avr を追加
         & varname='u_avr', dims=(/'x', 'z', 't'/),          &
         & longname='u average', units='m/s', xtype='real',  &
         & history=hist_u_mgn )
    call HistoryAddVariable( &  ! 変数 u_err を追加
         & varname='u_err', dims=(/'x', 'z', 't'/),       &
         & longname='u error', units='m/s', xtype='real', &
         & history=hist_u_mgn )

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

    call HistoryCreate( &
         & file = ncfile_w_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_w_mgn)          
  
    call HistoryPut('x', xs(imin:imax), hist_w_mgn)
    call HistoryPut('z', z(kmin:kmax), hist_w_mgn)
    
    call HistoryAddVariable( &  ! 変数 w を追加
         & varname='w', dims=(/'x', 'z', 't'/),       &
         & longname='w', units='m/s', xtype='real',   &
         & history=hist_w_mgn )
    call HistoryAddVariable( &  ! 変数 w_avr を追加
         & varname='w_avr', dims=(/'x', 'z', 't'/),          &
         & longname='w average', units='m/s', xtype='real',  &
         & history=hist_w_mgn )
    call HistoryAddVariable( &  ! 変数 w_err を追加
         & varname='w_err', dims=(/'x', 'z', 't'/),       &
         & longname='w error', units='m/s', xtype='real', &
        & history=hist_w_mgn )

  end subroutine gt4_init

!--- 変数出力

  subroutine gt4_output
    call HistoryPut('t', 0.0, hist_u)
  
    call HistoryPut('u', u(0:im-1, 0:km-1), hist_u)
    call HistoryPut('u_avr', u_avr(0:im-1, 0:km-1), hist_u)
    call HistoryPut('u_err', u_err(0:im-1, 0:km-1), hist_u)

    call HistoryPut('t', 0.0, hist_w)
  
    call HistoryPut('w', w(0:im-1, 0:km-1), hist_w)
    call HistoryPut('w_avr', w_avr(0:im-1, 0:km-1), hist_w)
    call HistoryPut('w_err', w_err(0:im-1, 0:km-1), hist_w)
    
    call HistoryPut('t', 0.0, hist_u_mgn)
  
    call HistoryPut('u', u(imin:imax, kmin:kmax), hist_u_mgn)
    call HistoryPut('u_avr', u_avr(imin:imax, kmin:kmax), hist_u_mgn)
    call HistoryPut('u_err', u_err(imin:imax, kmin:kmax), hist_u_mgn)
    call HistoryPut('t', 0.0, hist_u_mgn)
  
    call HistoryPut('t', 0.0, hist_w_mgn)

    call HistoryPut('w', w(imin:imax, kmin:kmax), hist_w_mgn)
    call HistoryPut('w_avr', w_avr(imin:imax, kmin:kmax), hist_w_mgn)
    call HistoryPut('w_err', w_err(imin:imax, kmin:kmax), hist_w_mgn)

  end subroutine gt4_output

  subroutine gt4_close
    call HistoryClose(hist_u)
    call HistoryClose(hist_w)
    call HistoryClose(hist_u_mgn)
    call HistoryClose(hist_w_mgn)
  end subroutine gt4_close  

end program arare
