!c Description: 
!c   非静力学モデル deepconv/arare
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

program arare
  use gt4_history
  use gridset
  use fileset, only: fileset_init
  use expname, only: expname_init
  use if_diff
  implicit none

  real(8), allocatable  :: ss_Exner_A(:,:)
  real(8), allocatable  :: sf_DExnerDX_A(:,:)
  real(8), allocatable  :: sf_DExnerDXAnl_A(:,:)
  real(8), allocatable  :: sf_ExnerError_A(:,:)
  real(8), allocatable  :: fs_VelX_A(:,:)
  real(8), allocatable  :: sf_DVelXDX_A(:,:)
  real(8), allocatable  :: sf_DVelXDXAnl_A(:,:)
  real(8), allocatable  :: sf_VelXError_A(:,:)
  real(8), allocatable  :: sf_VelZ_A(:,:)
  real(8), allocatable  :: sf_DVelZDX_A(:,:)
  real(8), allocatable  :: sf_DVelZDXAnl_A(:,:)
  real(8), allocatable  :: sf_VelZError_A(:,:)
  real(8), allocatable  :: ss_Exner_B(:,:)
  real(8), allocatable  :: fs_VelX_B(:,:)
  real(8), allocatable  :: sf_VelZ_B(:,:)
  real(8), allocatable  :: fs_Fx_B(:,:)
  real(8), allocatable  :: sf_Fz_B(:,:)
  real(8), allocatable  :: ss_VelDiv_B(:,:)

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

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

!!!
!!!初期化
!!!
  !--- I/O ファイル名の初期化
  call fileset_init
  
  !--- 格子点情報の初期化
  call gridset_init
  
  !--- 出力ファイルの定義
  call OpenRegFile
  call OpenDimFile


!!!
!!!配列の確保, 初期化
!!!  
  call ArareAlloc
  call ArareArrayInit


!!!
!!!計算
!!!

  sf_DExnerDX_A = sf_dx_ss(ss_Exner_A)
  sf_ExnerError_A = sf_DExnerDX_A - sf_DExnerDXAnl_A 

  !--- ファイル出力
  call OutputRegFile
  call OutputDimFile
     
  
  !--- 出力ファイルのクローズ
  call CloseRegFile
  call CloseDimFile

  
contains
  
  subroutine ArareAlloc
    allocate(ss_Exner_A(DimXMin:DimXMax, DimZMin:DimZMax),  &
         & sf_DExnerDX_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_DExnerDXAnl_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_ExnerError_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & fs_VelX_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_DVelXDX_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_DVelXDXAnl_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_VelXError_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_VelZ_A(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_DVelZDX_A(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_DVelZDXAnl_A(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZError_A(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_B(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_B(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_B(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_Fx_B(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_Fz_B(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_VelDiv_B(DimXMin:DimXMax, DimZMin:DimZMax) )
  end subroutine ArareAlloc

  subroutine ArareArrayInit
    do i = DimXMin, DimXMax
       ss_Exner_A(i, DimZMin:DimZMax) = &
            & cos( 2 * pi * s_X(i) / (XMax - XMin) )
    end do
    ss_Exner_B = 0.0d0
    sf_DExnerDX_A = 0.0d0
    do i = DimXMin, DimXMax
       sf_DExnerDXAnl_A(i, DimZMin:DimZMax) = &
            & - 2 * pi / (XMax - XMin) * sin( 2 * pi * s_X(i) / (XMax - XMin) )
    end do
    sf_ExnerError_A = 0.0d0

    do i = DimXMin, DimXMax
       fs_VelX_A(i, DimZMin:DimZMax) =  & 
            & sin( 2 * pi * f_X(i) / (XMax - XMin) )
    end do
    fs_VelX_B = 0.0d0
    sf_DVelXDX_A = 0.0d0
    do i = DimXMin, DimXMax
       sf_DVelXDXAnl_A(i, DimZMin:DimZmax) =  &
            & - 2 * pi / (XMax - XMin) * cos( 2 * pi * s_X(i) / (XMax - XMin) )
    end do
    sf_VelXError_A = 0.0d0

    sf_VelZ_A = 0.0d0
    sf_VelZ_B = 0.0d0
    sf_DVelZDX_A = 0.0d0
    sf_DVelZDXAnl_A = 0.0d0
    sf_VelZError_A = 0.0d0

    fs_Fx_B = 0.0d0
    sf_Fz_B = 0.0d0

    ss_VelDiv_B = 0.0d0

  end subroutine ArareArrayInit

  !--- gtool4 出力関連
  subroutine OpenRegFile
    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 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='DExnerDX', dims=(/'x','z'/),        &
         & longname='nondimensional pressure diffentinated by x', units='1',&
         & xtype='double',                               & 
         & history = sf_reg)    
    call HistoryAddVariable(                             &
         & varname='DExnerDXAnl', dims=(/'x','z'/),        &
         & longname='nondimensional pressure diffentinated analitically', &
         & units='1',&
         & xtype='double',                               & 
         & history = sf_reg)    
    call HistoryAddVariable(                             &
         & varname='ExnerError', dims=(/'x','z'/),        &
         & longname='error', units='1',&
         & xtype='double',                               & 
         & history = sf_reg)    
    ! 速度
!    call HistoryAddVariable(                            &
!         & varname='VelX', dims=(/'x','z'/),             &
!         & longname='zonal velocity',         &
!         & units='m/s', xtype='double',                & 
!         & history = fs_reg)        
!    call HistoryAddVariable(                            &
!         & varname='DVelXDX', dims=(/'x','z'/),             &
!         & longname='zonal velocity diffrentinated by X',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)        
!    call HistoryAddVariable(                            &
!         & varname='DVelXDXAnl', dims=(/'x','z'/),             &
!         & longname='zonal velocity diffrentinated analitically',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)        
!    call HistoryAddVariable(                            &
!         & varname='VelXError', dims=(/'x','z'/),             &
!         & longname='error',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)        
!    ! 速度
!    call HistoryAddVariable(                            &
!         & varname='VelZ', dims=(/'x','z'/),             &
!         & longname='vertical velocity',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)    
!    call HistoryAddVariable(                            &
!         & varname='DVelZDX', dims=(/'x','z'/),             &
!         & longname='vertical velocity diffrentinated by X',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)    
!    call HistoryAddVariable(                            &
!         & varname='DVelZDXAnl', dims=(/'x','z'/),             &
!         & longname='vertical velocity diffrentinated analitically',    &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)    
!    call HistoryAddVariable(                            &
!         & varname='VelZError', dims=(/'x','z'/),             &
!         & longname='error',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_reg)     
   
  end subroutine OpenRegFile


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

    !--- エクスナー関数
    call HistoryPut('Exner', &
         & ss_Exner_A(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('DExnerDX', &
         & sf_DExnerDX_A(RegXMin+1:RegXMax, RegZMin+1:RegZMax), sf_reg)
    call HistoryPut('DExnerDXAnl', &
         & sf_DExnerDXAnl_A(RegXMin+1:RegXMax, RegZMin+1:RegZMax), sf_reg)
    call HistoryPut('ExnerError', &
         & sf_ExnerError_A(RegXMin+1:RegXMax, RegZMin+1:RegZMax), sf_reg)
    !--- 速度
!    call HistoryPut('VelX', &
!         & fs_VelX_A(RegXMin:RegXMax, RegZMin+1:RegZMax), fs_reg)
!    call HistoryPut('DVelXDX', &
!         & sf_DVelXDX_A(RegXMin:RegXMax, RegZMin+1:RegZMax), sf_reg)
!    call HistoryPut('DVelXDXAnl', &
!         & sf_DVelXDXAnl_A(RegXMin:RegXMax, RegZMin+1:RegZMax), sf_reg)
!    call HistoryPut('VelXError', &
!         & sf_VelXError_A(RegXMin:RegXMax, RegZMin+1:RegZMax), sf_reg)
    !--- 速度
!    call HistoryPut('VelZ', &
!         & sf_VelZ_A(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)
!    call HistoryPut('DVelZDX', &
!         & sf_DVelZDX_A(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)
!    call HistoryPut('DVelZDXAnl', &
!         & sf_DVelZDXAnl_A(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)
!    call HistoryPut('VelZError', &
!         & sf_VelZError_A(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)
  end subroutine OutputRegFile
  
  subroutine CloseRegFile
    call HistoryClose( ss_reg )
    call HistoryClose( fs_reg )
    call HistoryClose( sf_reg )
  end subroutine CloseRegFile


  !--- gtool4 出力関連
  subroutine OpenDimFile
    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
    
    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='DExnerDX', dims=(/'x','z'/),        &
         & longname='nondimensional pressure diffentinated by x', units='1',&
         & xtype='double',                               & 
         & history = sf_dim)    
    call HistoryAddVariable(                             &
         & varname='DExnerDXAnl', dims=(/'x','z'/),        &
         & longname='nondimensional pressure diffentinated by x', units='1',&
         & xtype='double',                               & 
         & history = sf_dim)    
    call HistoryAddVariable(                             &
         & varname='ExnerError', dims=(/'x','z'/),        &
         & longname='error', units='1',&
         & xtype='double',                               & 
         & history = sf_dim)    
    ! 速度
!    call HistoryAddVariable(                            &
!         & varname='VelX', dims=(/'x','z'/),        &
!         & longname='zonal velocity',                   &
!         & units='m/s', xtype='double',                 & 
!         & history = fs_dim)        
!    call HistoryAddVariable(                            &
!         & varname='DVelXDX', dims=(/'x','z'/),             &
!         & longname='zonal velocity diffrentinated by X',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_dim)        
!    call HistoryAddVariable(                            &
!         & varname='DVelXDXAnl', dims=(/'x','z'/),             &
!         & longname='zonal velocity diffrentinated analitically',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_dim)        
!    call HistoryAddVariable(                            &
!         & varname='VelXError', dims=(/'x','z'/),             &
!         & longname='error',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_dim)            ! 速度
    !--- 速度
!    call HistoryAddVariable(                            &
!         & varname='VelZ', dims=(/'x','z'/),        &
!         & longname='vertical velocity',                &
!         & units='m/s', xtype='double',                 & 
!         & history = sf_dim)    
!    call HistoryAddVariable(                            &
!         & varname='DVelZDX', dims=(/'x','z'/),             &
!         & longname='vertical velocity diffrentinated by X',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_dim)    
!    call HistoryAddVariable(                            &
!         & varname='DVelZDXAnl', dims=(/'x','z'/),             &
!         & longname='vertical velocity diffrentinated analitically',    &
!         & units='m/s', xtype='double',                & 
!         & history = sf_dim)    
!    call HistoryAddVariable(                            &
!         & varname='VelZError', dims=(/'x','z'/),             &
!         & longname='error',         &
!         & units='m/s', xtype='double',                & 
!         & history = sf_dim)     
    
  end subroutine OpenDimFile


  subroutine OutputDimFile
    !--- エクスナー関数
    call HistoryPut('Exner', ss_Exner_A, ss_dim)
    call HistoryPut('DExnerDX', sf_DExnerDX_A, sf_dim)
    call HistoryPut('DExnerDXAnl', sf_DExnerDXAnl_A, sf_dim)
    call HistoryPut('ExnerError', sf_ExnerError_A, sf_dim)
    !--- 速度
!    call HistoryPut('VelX', fs_VelX_A, fs_dim)
!    call HistoryPut('DVelXDX', sf_DVelXDX_A, sf_dim)
!    call HistoryPut('DVelXDXAnl', sf_DVelXDXAnl_A, sf_dim)
!    call HistoryPut('VelXError', sf_VelXError_A, sf_dim)
    !--- 速度
!    call HistoryPut('VelZ', sf_VelZ_A, sf_dim)
!    call HistoryPut('DVelZDX', sf_DVelZDX_A, sf_dim)
!    call HistoryPut('DVelZDXAnl', sf_DVelZDXAnl_A, sf_dim)
!    call HistoryPut('VelZError', sf_VelZError_A, sf_dim)
  end subroutine OutputDimFile
  
  subroutine CloseDimFile
    call HistoryClose( ss_dim )
    call HistoryClose( fs_dim )
    call HistoryClose( sf_dim )
  end subroutine CloseDimFile


end program arare
