!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program Arare
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: arare.f90,v 1.10 2005/02/07 10:12:29 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!非静力学モデル deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!方程式系は準圧縮系.
!
!== Future Plans
!

                                                                 !=end

program arare
                                                                 !=begin
  !== Dependency
  use gt4_history
  use dc_trace, only: BeginSub, EndSub, DbgMessage
  use dc_message
  use debugset, only: debugset_init
  use gridset,  only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, &
       &              RegZMin, RegZMax, XMin, XMax, ZMin, ZMax, &
       &              f_X, s_X, f_Z, s_Z
  use fileset,  only: fileset_init, cfgfile, basicfile, distrbfile
  use expname,  only: expname_init
  use linlib,   only: linlib_init
  use bcset,    only: bcset_init
  use arareprm, only: arareprm_init 
  use differentiate_center4

                                                                 !=end
  implicit none

  real(8), allocatable  :: ss_PotTemp_nl(:,:)
  real(8), allocatable  :: ss_Exner_nl(:,:)
  real(8), allocatable  :: fs_VelX_nl(:,:)
  real(8), allocatable  :: sf_VelZ_nl(:,:)

  real(8), allocatable  :: ss_DVelXDXAnl_nl(:,:)
  real(8), allocatable  :: ss_DVelXDX_nl(:,:)
  real(8), allocatable  :: ss_ErrorDVelXDX_nl(:,:)

  real(8), parameter :: pi=3.14159d0  
  integer  ::  i,k
  

  type(gt_history)      :: ss_reg, fs_reg, sf_reg
  type(gt_history)      :: ss_dim, fs_dim, sf_dim
  

!!!
!!!初期化
!!!
  !--- I/O ファイル名の初期化
  call fileset_init

  !--- デバグ設定の初期化
  call debugset_init(cfgfile)

  !--- 格子点情報の初期化
  call gridset_init(cfgfile)
  
  !--- 実験名の初期化
  call expname_init(cfgfile)
  
  !--- 計算の設定を読み込む
  call arareprm_init(cfgfile)

  !--- 境界条件を設定
  call bcset_init(cfgfile)

  !--- 行列計算ライブラリの初期化. 
  !--- 引数はスカラーの z 方向の計算領域の大きさ
  call linlib_init(cfgfile, RegZMax - RegZMin)
  
  !--- 出力ファイルの定義
  call ArareOpenRegFile
  call ArareOpenDimFile


!!!
!!!変数の初期化
!!!  
  !--- 配列の割り当て, 初期化
  call ArareAlloc


!!!
!!!微分テスト
!!!

  !--- 元の関数

  do k = DimZMin, DimZMax  
     fs_VelX_nl(:,k) = sin(2.0d0 * pi * f_X(:) / (XMax - XMin) ) 
  end do

  !--- 解析的に微分した関数

  do k = DimZMin, DimZMax
     ss_DVelXDXAnl_nl(:,k) =    &
          & 2.0d0 * pi / (XMax - XMin)  & 
          & * cos(2.0d0 * pi * s_X(:) / (XMax - XMin) )
  end do


! 微分計算

  ss_DVelXDX_nl = ss_dx_fs(fs_VelX_nl)

! 誤差

  ss_ErrorDVelXDX_nl = ss_DVelXDX_nl - ss_DVelXDXAnl_nl

  !--- ファイル出力
  call ArareOutputRegFile
  call ArareOutputDimFile

  !--- 微分演算テスト用ファイル出力

  call ArareOutputDiff
  
  !--- 出力ファイルのクローズ
  call ArareCloseRegFile
  call ArareCloseDimFile
  
  
contains


!!!
!!!配列の割り当て, 初期化
!!!
  subroutine ArareAlloc
    use distrbset, only: fs_VelXDistrb, sf_VelZDistrb, &
         &               ss_ExnerDistrb, ss_PotTempDistrb

    call BeginSub("ArareAlloc", &
&                 fmt="%c",     &
&                 c1="Allocate and initialize allocatable variables.")

    allocate( fs_VelX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTemp_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_DVelXDXAnl_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_DVelXDX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_ErrorDVelXDX_nl(DimXMin:DimXMax, DimZMin:DimZMax) )

    fs_VelX_nl = 0.0d0

    sf_VelZ_nl = 0.0d0

    ss_Exner_nl = 0.0d0
    
    ss_PotTemp_nl = 0.0d0

    !速度場の初期値
    fs_VelX_nl = 0.0d0
    sf_VelZ_nl = 0.0d0

    !無次元圧力の初期値
    ss_Exner_nl = 0.0d0

    !温位の初期値
    ss_PotTemp_nl = 0.0d0

    call EndSub("ArareAlloc")

  end subroutine ArareAlloc


!!!
!!!出力ファイルの定義. 計算領域のみ
!!! 
  subroutine ArareOpenRegFile
    use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax, &
         & NX, NZ, f_X, f_Z, s_X, s_Z
    use expname, only: exptitle, expsrc, expinst
    use fileset, only: fs_RegFile, sf_RegFile, ss_RegFile

    call BeginSub("ArareOpneRegFile", &
&                 fmt="%c",           &
&                 c1="Initialize output file name and variables in physical region.")
    
    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = fs_RegFile,                               &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z'/),                           &
         & dimsizes=(/NX+1, NZ/),                       &
         & longnames=(/'X-coordinate',                     &
         &             'Z-coordinate'/),                   &
         & units=(/'m','m'/), origin=0.0,              &
         & interval=0.0,                     &
         & history = fs_reg)
    
    ! 変数出力
    call HistoryPut('x', f_X(RegXMin: RegXMax), fs_reg)
    call HistoryPut('z', s_Z(RegXMin+1: RegXMax), fs_reg)
    
    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = sf_RegFile,                               &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z'/),                           &
         & dimsizes=(/NX, NZ+1/),                       &
         & longnames=(/'X-coordinate',                     & 
         &             'Z-coordinate'/),                   &
         & units=(/'m','m'/),                          &
         & origin=0.0,                                     &
         & interval=0.0,                                   &
         & history = sf_reg)
    
    ! 変数出力
    call HistoryPut('x', s_X(RegXMin+1: RegXMax), sf_reg)
    call HistoryPut('z', f_Z(RegXMin: RegXMax), sf_reg)

    ! ヒストリー作成
    call HistoryCreate(                                  &
         & file = ss_RegFile,                            &         
         & title = exptitle,                             &
         & source = expsrc,                              &
         & institution = expinst,                        &
         & dims=(/'x','z'/),                         &
         & dimsizes=(/NX, NZ/),                       &
         & longnames=(/'X-coordinate',                   &
         &             'Z-coordinate'/),                 &
         & units=(/'m','m'/), origin=0.0,            &
         & interval=0.0,                                 &
         & history = ss_reg)
    
    ! 変数出力
    call HistoryPut('x', s_X(RegXMin+1: RegXMax), ss_reg)
    call HistoryPut('z', s_Z(RegXMin+1: RegXMax), ss_reg)
    
    ! 無次元圧力
    call HistoryAddVariable(                             &
         & varname='Exner', dims=(/'x','z'/),        &
         & longname='nondimensional pressure', units='1',&
         & xtype='double',                               & 
         & history = ss_reg)    
    ! 仮温位
    call HistoryAddVariable(                            &
         & varname='PotTemp', dims=(/'x','z'/),      &
         & longname='virtual potential temperature',    &
         & units='K', xtype='double',                 & 
         & history = ss_reg)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelX', dims=(/'x','z'/),             &
         & longname='zonal velocity',         &
         & units='m/s', xtype='double',                & 
         & history = fs_reg)        
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelZ', dims=(/'x','z'/),             &
         & longname='vertical velocity',         &
         & units='m/s', xtype='double',                & 
         & history = sf_reg)    
    
    call EndSub("ArareOpneRegFile")
    
  end subroutine ArareOpenRegFile


!!!
!!!ファイル出力. 長い時間ステップで出力. 計算領域のみ
!!!
  subroutine ArareOutputRegFile
    use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax

    call BeginSub("ArareOutputRegFile", &
&                 fmt="%c",             &
&                 c1="Output variables in physical region.")

    call HistoryPut('Exner', &
         & ss_Exner_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('PotTemp', &
         & ss_PotTemp_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('VelX', &
         & fs_VelX_nl(RegXMin:RegXMax, RegZMin+1:RegZMax), fs_reg)
    call HistoryPut('VelZ', &
         & sf_VelZ_nl(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)

    call EndSub("ArareOutputRegFile")

  end subroutine ArareOutputRegFile


!!!
!!!出力ファイルのクローズ. 計算領域のみ
!!!  
  subroutine ArareCloseRegFile

    call BeginSub("ArareCloseRegFile", &
&                 fmt="%c", &
&                 c1="Close output file name for physical region.")

    call HistoryClose( ss_reg )
    call HistoryClose( fs_reg )
    call HistoryClose( sf_reg )

    call EndSub("ArareCloseRegFile")

  end subroutine ArareCloseRegFile


!!!
!!!出力ファイルの定義. 全領域
!!!
  subroutine ArareOpenDimFile
    use gridset, only: f_X, f_Z, s_X, s_Z
    use expname, only: exptitle, expsrc, expinst
    use fileset, only: fs_DimFile, sf_DimFile, ss_DimFile

    integer::   N

    call BeginSub("ArareOpneDimFile", &
&                  fmt="%c",         &
&     c1="Initialize output file name and variables in all computational region.")

    N = size(f_X, 1)

    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = fs_DimFile,                              &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z'/),                           &
         & dimsizes=(/N, N/),                           &
         & longnames=(/'X-coordinate',                     &
         &             'Z-coordinate'/),                   &
         & units=(/'m','m'/), origin=0.0,              &
         & interval=0.0, &
         & history = fs_dim)
    
    ! 変数出力
    call HistoryPut('x', f_X, fs_dim)
    call HistoryPut('z', s_Z, fs_dim)
    
    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = sf_DimFile,                              &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z'/),                           &
         & dimsizes=(/N, N/),                           &
         & longnames=(/'X-coordinate',                     & 
         &             'Z-coordinate'/),                   &
         & units=(/'m','m'/),                          &
         & origin=0.0,                                     &
         & interval=0.0,                                   &
         & history = sf_dim)
    
    ! 変数出力
    call HistoryPut('x', s_X, sf_dim)
    call HistoryPut('z', f_Z, sf_dim)

    ! ヒストリー作成
    call HistoryCreate(                                  &
         & file = ss_DimFile,                            &
         & title = exptitle,                             &
         & source = expsrc,                              &
         & institution = expinst,                        &
         & dims=(/'x','z'/),                         &
         & dimsizes=(/N, N/),                         &
         & longnames=(/'X-coordinate',                   &
         &             'Z-coordinate'/),                 &
         & units=(/'m','m'/), origin=0.0,            &
         & interval=0.0,                                 &
         & history = ss_dim)
    
    ! 変数出力
    call HistoryPut('x', s_X, ss_dim)
    call HistoryPut('z', s_Z, ss_dim)
    
    ! 無次元圧力
    call HistoryAddVariable(                             &
         & varname='Exner', dims=(/'x','z'/),        &
         & longname='nondimensional pressure', units='1',&
         & xtype='double',                               & 
         & history = ss_dim)    
    ! 仮温位
    call HistoryAddVariable(                            &
         & varname='PotTemp', dims=(/'x','z'/),      &
         & longname='virtual potential temperature',    &
         & units='K', xtype='double',                   & 
         & history = ss_dim)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelX', dims=(/'x','z'/),        &
         & longname='zonal velocity',                   &
         & units='m/s', xtype='double',                 & 
         & history = fs_dim)        
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelZ', dims=(/'x','z'/),        &
         & longname='vertical velocity',                &
         & units='m/s', xtype='double',                 & 
         & history = sf_dim)    

    call EndSub("ArareOpneDimFile")

  end subroutine ArareOpenDimFile


!!!
!!!ファイル出力. 長い時間ステップの値を出力. 全領域
!!!
  subroutine ArareOutputDimFile

    call BeginSub("ArareOutputDimFile", &
&                 fmt="%c",            &
&                 c1="Output variables in all computational region.")

    call HistoryPut('Exner', ss_Exner_nl, ss_dim)
    call HistoryPut('PotTemp', ss_PotTemp_nl, ss_dim)
    call HistoryPut('VelX',  fs_VelX_nl, fs_dim)
    call HistoryPut('VelZ',  sf_VelZ_nl, sf_dim)

    call EndSub("ArareOutputDimFile")

  end subroutine ArareOutputDimFile


!!!
!!!出力ファイルのクローズ. 全領域
!!!  
  subroutine ArareCloseDimFile

    call BeginSub("ArareCloseDimFile", &
&                 fmt="%c",            &
&                  c1="Close output file name for all computational region.")

    call HistoryClose( ss_dim )
    call HistoryClose( fs_dim )
   c all HistoryClose( sf_dim )

    call EndSub("ArareCloseDimFile")

  end subroutine ArareCloseDimFile

!!!
!!!微分演算テスト用出力サブルーチン
!!!

  subroutine ArareOutputDiff

    use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax

    !---計算領域用ヒストリに変数追加

    call HistoryAddVariable(                             &
         & varname='DVelXDX', dims=(/'x','z'/),        &
         & longname='numerically derived x-velosity', units='1/s',&
         & xtype='double',                               & 
         & history = ss_reg)    
    call HistoryAddVariable(                             &
         & varname='DVelXDXAnl', dims=(/'x','z'/),        &
         & longname='analytically derived x-velosity', units='1/s',&
         & xtype='double',                               & 
         & history = ss_reg)    
    call HistoryAddVariable(                             &
         & varname='ErrorDVelXDX', dims=(/'x','z'/),        &
         & longname='error of numerical differenciation', units='1/s',&
         & xtype='double',                               & 
         & history = ss_reg)    

    !---全領域用ヒストリに変数追加

    call HistoryAddVariable(                             &
         & varname='DVelXDX', dims=(/'x','z'/),        &
         & longname='numecally derived x-velosity', units='1',&
         & xtype='double',                               & 
         & history = ss_dim)    
    call HistoryAddVariable(                             &
         & varname='DVelXDXAnl', dims=(/'x','z'/),        &
         & longname='analytically derived x-velosity', units='1',&
         & xtype='double',                               & 
         & history = ss_dim)    
    call HistoryAddVariable(                             &
         & varname='ErrorDVelXDX', dims=(/'x','z'/),        &
         & longname='error of numerical differenciation', units='1',&
         & xtype='double',                               & 
         & history = ss_dim)    

    !--- 計算領域用ヒストリ出力

    call HistoryPut('DVelXDX',   &
         & ss_DVelXDX_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('DVelXDXAnl',    &
         & ss_DVelXDXAnl_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('ErrorDVelXDX',    &
         & ss_ErrorDVelXDX_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)

    !--- 全領域用ヒストリ出力

    call HistoryPut('DVelXDX', ss_DVelXDX_nl, ss_dim)
    call HistoryPut('DVelXDXAnl', ss_DVelXDXAnl_nl, ss_dim)
    call HistoryPut('ErrorDVelXDX', ss_ErrorDVelXDX_nl, ss_dim)

  end subroutine ArareOutputDiff 

end program arare
