program ArareInitData
  !
  !非静力学モデル deepconv/arare. 
  !

  !----- モジュール読み込み ------

  !-----   型宣言, 文字列処理   ----
  use dc_types,       only : STRING, DP
  use dc_string,      only : StoA
  use gtool_history,  only : HistoryPut

  !-----   メッセージ出力   -----
  use dc_message,     only: MessageNotify

  !  コマンドライン引数解釈
  use argset,        only : argset_init  

  !  MPI の初期化
  !  Initialize MPI wrapper
  !
  use mpi_wrapper,   only : MPIWrapperInit, MPIWrapperFinalize

  !-----    管理モジュール   -----
  !  入出力ファイル名管理モジュール
  use fileset,    only : fileset_init

  !  格子点管理モジュール 
  use gridset, only : gridset_init, imin, imax, jmin, jmax, kmin, kmax, ncmax

  !  基本場設定モジュール
  use basicset,      only : basicset_init

  use axesset,       only : axesset_init

  use constants,     only : constants_init, molwtdry, pressbasis, gasrdry, cpdry, cvdry

  !  湿潤ルーチン設定モジュール
  use composition,      only: composition_init, molwtwet

  ! 化学計算
  use ChemCalc, only: ChemCalc_init

  !  境界条件適用モジュール
  use setmargin

  use initialdata_basic, only: initialdata_basic_nml, initialdata_basic_dry, initialdata_basic_moist, initialdata_basic_strat, initialdata_basic_isothermal, ID_FlaginitdataBasic, ID_FlaginitDataBasic_dry, ID_FlaginitDataBasic_moist , ID_FlaginitDataBasic_isothermal
  use initialdata_disturb, only: initialdata_disturb_nml, initialdata_disturb_random, initialdata_disturb_gaussXZ, initialdata_disturb_gaussXY, initialdata_disturb_gaussXYZ, initialdata_disturb_dryregion, initialdata_disturb_moist, ID_FlaginitdataPTemp, ID_FlaginitdataExner, ID_FlaginitdataQMix, ID_Flaginitdata_GaussXY, ID_Flaginitdata_GaussXZ, ID_Flaginitdata_GaussXYZ, ID_Flaginitdata_random, ID_FlaginitdataQMix_dryregion, ID_FlaginitdataQMix_moist
  use initialdata_yamasaki1983


  !-----    入出力モジュール   -----
  !  リスタートファイル入出力モジュール
  use RestartFileIO, only : ReStartFileIO_Init, ReStartFileIO_Finalize, rstat
  use Arare4FileIO,  only : Arare4FileIO_Init, Arare4FileIO_Var_Get, Arare4FileIO_BZ_Get, Arare4FileIO_MMC_Var_Get, Arare4FileIO_MMC_BZ_Get

  !暗黙の型宣言禁止
  implicit none

  !内部変数
  real(DP), allocatable :: pyz_VelX(:,:,:)
  real(DP), allocatable :: xqz_VelY(:,:,:)
  real(DP), allocatable :: xyr_VelZ(:,:,:)
  real(DP), allocatable :: xyz_Exner(:,:,:)
  real(DP), allocatable :: xyz_PTemp(:,:,:)
  real(DP), allocatable :: xyz_Km(:,:,:)
  real(DP), allocatable :: xyz_Kh(:,:,:)
  real(DP), allocatable :: xyz_CDens(:,:,:)
  real(DP), allocatable :: xyzf_QMix(:,:,:,:)
  real(DP), allocatable :: xyz_DensBZ(:,:,:)
  real(DP), allocatable :: xyz_PressBZ(:,:,:)
  real(DP), allocatable :: xyz_ExnerBZ(:,:,:)
  real(DP), allocatable :: xyz_TempBZ(:,:,:)
  real(DP), allocatable :: xyz_PTempBZ(:,:,:)
  real(DP), allocatable :: xyz_VelSoundBZ(:,:,:)
  real(DP), allocatable :: xyzf_QMixBZ(:,:,:,:)
  real(DP), allocatable :: xyz_EffMolWtBZ(:,:,:)
  real(DP), allocatable :: z_TempBZ(:)
  real(DP), allocatable :: z_PressBZ(:)
  real(DP), allocatable :: xyzf_MolFr(:,:,:,:)
  real(DP), allocatable :: zf_MolFr(:,:)
  real(DP), allocatable :: xyzf_QMixDivMolWt(:,:,:,:)
  real(DP)              :: Time
  integer               :: i, j, s
  character(STRING)     :: cfgfile = ""
  integer            :: ID_FlaginitData = 0
  integer, parameter :: ID_FlaginitData_default = 1
  integer, parameter :: ID_FlaginitData_Yamasaki1983 = 2
  integer, parameter :: ID_FlaginitData_Takemi2007 = 3
  integer, parameter :: ID_FlaginitData_arare4 = 4
  integer, parameter :: ID_FlaginitData_arare4_mmc = 5
  
  ! MPI
  !
  call MPIWrapperInit

  !NAMELIST ファイル名の読み込み
  !
  call argset_init(cfgfile) !(out)

  ! nml から情報を取り出す
  !
  call deepconvInitialData(cfgfile)
  
  !格子点情報の初期化
  !  NAMELIST から情報を得て, 格子点を計算する
  !
  call gridset_init(cfgfile)
  
  ! 化学計算ルーチンの初期化
  ! Initialization of chemical routines.
  !
  call chemcalc_init( )
  
  ! 軸の情報の初期化
  ! Initialization of axis variables.
  !
  call axesset_init( cfgfile ) !(in)

  ! 定数の情報の初期化
  ! Initialization of constant variables.
  !
  call constants_init( cfgfile ) !(in)

  ! 湿潤過程共有変数の初期化
  ! Initialization of common variables for moist process.
  !
  call composition_init( cfgfile )
  
  ! I/O ファイル名の初期化
  ! Initialization of output file name. 
  !
  call fileset_init( cfgfile ) ! (in)    

  !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. 
  !
  call ArareAlloc

  !---------------------------------------------------------------  
  ! 基本場, 擾乱場を作成する. 
  !
  Time = 0.0d0

  call MessageNotify( "M", "main", "Making Initial data...." )

  select case(ID_FlaginitData)      
    
  !----------------------------------    
  ! deepconv/arare4 の初期値を読み込む
  !
  case (ID_FlaginitData_arare4)
    call MessageNotify( "M", "main", "Making Initial data using arare4 files ...." )
    call arare4fileio_init(cfgfile)
    call Arare4fileio_Var_Get( pyz_VelX, xqz_VelY, xyr_VelZ, xyz_PTemp, xyz_Exner, xyzf_QMix, xyz_Km, xyz_Kh, xyz_CDens   )   ! (out)
    call Arare4fileio_BZ_Get( xyz_PressBZ, xyz_ExnerBZ, xyz_TempBZ, xyz_PTempBZ, xyz_DensBZ, xyz_VelSoundBZ, xyzf_QMixBZ, xyz_EffMolWtBZ )

  !----------------------------------    
  ! deepconv/arare4 火星版の初期値を読み込む
  !
  case (ID_FlaginitData_arare4_mmc)
    call MessageNotify( "M", "main", "Making Initial data using arare4 files ...." )
    call arare4fileio_init(cfgfile)
    call Arare4fileio_MMC_Var_Get( pyz_VelX, xqz_VelY, xyr_VelZ, xyz_PTemp, xyz_Exner, xyzf_QMix, xyz_Km, xyz_Kh, xyz_CDens   )   ! (out)
    call Arare4fileio_MMC_BZ_Get( xyz_PressBZ, xyz_ExnerBZ, xyz_TempBZ, xyz_PTempBZ, xyz_DensBZ, xyz_VelSoundBZ, xyzf_QMixBZ, xyz_EffMolWtBZ )

  !----------------------------------    
  ! デフォルト(?)な初期値
  !
  case (ID_FlaginitData_default)
    call MessageNotify( "M", "main", "Making Initial data by default setting...." )
    call initialdata_basic_nml(cfgfile)
    call initialdata_disturb_nml(cfgfile)

    ! 基本場の温度・圧力・物質存在量を決める
    !
    select case(ID_FlaginitdataBasic)
      
    ! 等温/重力考慮しない場合
    !
    case(ID_FlaginitDataBasic_isothermal) 
      call MessageNotify( "M", "main", "Making Initial data (basic) named IsoThermal..." )
      call initialdata_basic_isothermal( z_TempBZ, z_PressBZ, zf_MolFr ) !(out)

    ! 乾燥断熱的な初期場
    !
    case(ID_FlaginitDataBasic_dry) 
      call MessageNotify( "M", "main", "Making Initial data (basic) named DRY..." )
      call initialdata_basic_dry( z_TempBZ, z_PressBZ, zf_MolFr ) !(out)
!      call initialdata_basic_strat( z_TempBZ, z_PressBZ ) !(inout)
      
    ! 湿潤断熱的な初期場
    !
    case(ID_FlaginitDataBasic_moist) 
      call MessageNotify( "M", "main", "Making Initial data (basic) named MOIST..." )
      call initialdata_basic_moist( z_TempBZ, z_PressBZ, zf_MolFr ) !(out)
!      call initialdata_basic_strat( z_TempBZ, z_PressBZ ) !(inout)
      
    end select
    
    ! 基本場の残りの値を決める
    !
    call DetermineBZ      


    ! 温位: 擾乱場の値を決める
    !    
    select case(ID_FlaginitdataPTemp)

    ! 温位: ガウシアンな分布を与える (XY)
    !          
    case(ID_Flaginitdata_GaussXY)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXY..." )
      call initialdata_disturb_gaussXY(xyz_PTemp)

    ! 温位: ガウシアンな分布を与える (XZ)
    !                
    case(ID_Flaginitdata_GaussXZ)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXZ..." )
      call initialdata_disturb_gaussXZ(xyz_PTemp)

    ! 温位: ガウシアンな分布を与える (XYZ)
    !                
    case(ID_Flaginitdata_GaussXYZ)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXYZ..." )
      call initialdata_disturb_gaussXYZ(xyz_PTemp)

    ! 温位: 指定された高度にランダムな分布を与える. 
    !                      
    case(ID_Flaginitdata_random)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named random..." )
      call initialdata_disturb_random(xyz_PTemp)
      
    end select

    ! 温位: マージンの値を確定
    !                
    call SetMargin_xyz( xyz_PTemp )

    ! エクスナー関数: 擾乱場の値を決める    
    !                
    select case(ID_FlaginitdataExner)

    ! エクスナー関数: ガウシアンな分布を与える (XY)
    !                      
    case(ID_Flaginitdata_GaussXY)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXY..." )
      call initialdata_disturb_gaussXY(xyz_Exner)

    ! エクスナー関数: ガウシアンな分布を与える (XZ)
    !                            
    case(ID_Flaginitdata_GaussXZ)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXZ..." )
      call initialdata_disturb_gaussXZ(xyz_Exner)
      
    ! エクスナー関数: ガウシアンな分布を与える (XYZ)
    !                      
    case(ID_Flaginitdata_GaussXYZ)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXYZ..." )
      call initialdata_disturb_gaussXYZ(xyz_Exner)
    end select

    ! エクスナー関数: マージンの値を確定
    !                
    call SetMargin_xyz( xyz_Exner )

    ! 混合比: 擾乱場の値を決める
    !                    
    select case(ID_FlaginitdataQMix)

    ! 混合比: ガウシアンな分布を与える (XY)
    !                            
    case(ID_Flaginitdata_GaussXY)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXY..." )
      call initialdata_disturb_gaussXY(xyzf_QMix(:,:,:,1))

    ! 混合比: ガウシアンな分布を与える (XZ)
    !                                  
    case(ID_Flaginitdata_GaussXZ)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXZ..." )
      call initialdata_disturb_gaussXZ(xyzf_QMix(:,:,:,1))

    ! 混合比: ガウシアンな分布を与える (XYZ)
    !                                  
    case(ID_Flaginitdata_GaussXYZ)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named GaussXYZ..." )
      call initialdata_disturb_gaussXYZ(xyzf_QMix(:,:,:,1))

    ! 混合比: 乾燥領域を作る
    !                            
    case(ID_FlaginitdataQMix_dryregion)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named DryRegion..." )
      call initialdata_disturb_dryregion(xyzf_QMix)

    ! 混合比: 湿潤断熱的な分布を与える 
    !                                  
    case(ID_FlaginitdataQMix_moist)
      call MessageNotify( "M", "main", "Making Initial data (disturb) named MOIST..." )
      call initialdata_disturb_moist(xyzf_QMix)
      
    end select

    call SetMargin_xyzf( xyzf_QMix )
   
  !----------------------------------------------------- 
  ! Yamasaki(1983)の温度と相対湿度の観測値を使用する場合 
  !
  case(ID_Flaginitdata_Yamasaki1983) 
    call MessageNotify( "M", "main", "Making Initial data named Yamasaki1983..." )
!    call initialdata_yamasaki1983_nml
    
    ! 基本場
    !
    call initialdata_yamasaki1983_set( z_TempBZ, z_PressBZ, zf_MolFr )
    
    ! 残りの基本場の値を決める
    !
    call DetermineBZ            
    
    ! ランダムな温位分布を置く
    !
    call initialdata_disturb_random(xyz_PTemp)
    call SetMargin_xyz( xyz_PTemp )

  !----------------------------------------------------- 
  ! Takemi(2007)の基本場を使用する場合
  !
  case(ID_Flaginitdata_Takemi2007)
    call MessageNotify( "M", "main", "Making Initial data named Takemi2007..." )
!      call initialdata_takemi2007_nml
!
!      ! 基本場
!      call initialdata_takemi2007_basic( z_TempBZ, z_PressBZ, zf_MolFr )
!
!      ! 残りの基本場の値を決める
!      call DetermineBZ      
!
!      ! 風速分布
!      call initialdata_takemi2007_wind( z_TempBZ, z_PressBZ, zf_MolFr ) 
!
!      ! ランダムな温位分布を置く
!      call initialdata_disturb_random()
  end select
  
  !のり代の値を決める  
  !
  call SetMargin_xyz( xyz_PTemp )
  
  !のり代の値を決める
  !
  call SetMargin_xyzf( xyzf_QMix )
  call MessageNotify( "M", "main", "Output variables into netCDF file..." )
  
  !---------------------------------
  ! リスタートファイル作成. 基本場と擾乱場を出力. 
  !
  call ReStartFileIO_Init(cfgfile)

  call HistoryPut( 't',     Time,      rstat)
  call HistoryPut( 'VelX',  pyz_VelX,  rstat)
  call HistoryPut( 'VelY',  xqz_VelY,  rstat)
  call HistoryPut( 'VelZ',  xyr_VelZ,  rstat)
  call HistoryPut( 'Exner', xyz_Exner, rstat)
  call HistoryPut( 'PTemp', xyz_PTemp, rstat)
  call HistoryPut( 'Km',    xyz_Km,    rstat)
  call HistoryPut( 'Kh',    xyz_Kh,    rstat)
  call HistoryPut( 'QMix',  xyzf_QMix, rstat)    
  call HistoryPut( 'CDens', xyz_CDens, rstat)    

  !-------------------------------------------------------------
  ! 基本場のファイル出力
  !
  call HistoryPut( 'DensBZ',   xyz_DensBZ     , rstat)
  call HistoryPut( 'ExnerBZ',  xyz_ExnerBZ    , rstat)
  call HistoryPut( 'PTempBZ',  xyz_PTempBZ  , rstat)
  call HistoryPut( 'VelSoundBZ', xyz_VelSoundBZ , rstat)
  call HistoryPut( 'TempBZ',    xyz_TempBZ     , rstat)
  call HistoryPut( 'PressBZ',   xyz_PressBZ    , rstat)
  call HistoryPut( 'QMixBZ',    xyzf_QMixBZ  , rstat)
  call HistoryPut( 'EffMolWtBZ', xyz_EffMolWtBZ, rstat)
  
  call ReStartFileIO_Finalize

  !----------------------------------------------------
  ! MPI END
  !
  call MPIWrapperFinalize
  
contains

  subroutine ArareAlloc
    !
    !初期化として, 配列を定義し, 値としてゼロを代入する.
    !

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate( pyz_VelX(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xqz_VelY(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyr_VelZ(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_Exner(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_PTemp(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_Km(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_Kh(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_CDens(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyzf_QMix(imin:imax, jmin:jmax, kmin:kmax, ncmax) )

    allocate( xyz_DensBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( xyz_PressBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( xyz_ExnerBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( xyz_TempBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( xyz_PTempBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( xyz_VelSoundBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( xyzf_QMixBZ(imin:imax,jmin:jmax,kmin:kmax,ncmax))
    allocate( xyz_EffMolWtBZ(imin:imax,jmin:jmax,kmin:kmax))
    allocate( z_TempBZ(kmin:kmax))
    allocate( z_PressBZ(kmin:kmax))
    allocate( xyzf_MolFr(imin:imax, jmin:jmax, kmin:kmax, ncmax))
    allocate( zf_MolFr(kmin:kmax, ncmax))
    allocate( xyzf_QMixDivMolWt(imin:imax,jmin:jmax,kmin:kmax, ncmax))

    pyz_VelX = 0.0d0
    xqz_VelY = 0.0d0
    xyr_VelZ = 0.0d0
    xyz_PTemp = 0.0d0
    xyz_Exner = 0.0d0
    xyz_Km = 0.0d0
    xyz_Kh = 0.0d0
    xyz_CDens = 0.0d0
    xyzf_QMix = 0.0d0

    xyz_DensBZ = 0.0d0
    xyz_PressBZ = 0.0d0
    xyz_ExnerBZ = 0.0d0
    xyz_TempBZ = 0.0d0
    xyz_PTempBZ = 0.0d0
    xyz_VelSoundBZ = 0.0d0
    xyzf_QMixBZ = 0.0d0
    xyz_EffMolWtBZ  = 0.0d0
    z_TempBZ  = 0.0d0
    z_PressBZ = 0.0d0
    xyzf_MolFr = 0.0d0
    zf_MolFr = 0.0d0 
    xyzf_QMixDivMolWt = 0.0d0

  end subroutine ArareAlloc
    

  subroutine DetermineBZ

    ! 3 次元配列に格納
    do j = jmin, jmax
      do i = imin, imax
        xyz_TempBZ(i,j,:)  = z_TempBZ(:)  
        xyz_PressBZ(i,j,:) = z_PressBZ(:)
      end do
    end do
    
    !のり代の値を決める
    !
    call SetMargin_xyz( xyz_TempBZ )
    call SetMargin_xyz( xyz_PressBZ)

    !---------------------------------------------------------------
    ! 混合比
    !---------------------------------------------------------------
    !水平方向には一様
    do i = imin, imax
      do j = jmin, jmax
        xyzf_MolFr(i,j,:,:) = zf_MolFr
      end do
    end do
 
    !気相のモル比を混合比に変換
    do s = 1, ncmax
      xyzf_QMixBZ(:,:,:,s) = xyzf_MolFr(:,:,:,s) * MolWtWet(s) / MolWtDry
    end do
    
    !  !値が小さくなりすぎないように最低値を与える
    !  where (xyzf_QMixBZ <= 1.0d-20 )
    !    xyzf_QMixBZ = 1.0d-20
    !  end where
    
    !のり代の値を決める
    !
    call SetMargin_xyzf( xyzf_QMixBZ )
    
    !---------------------------------------------------------------
    ! 分子量の効果
    !---------------------------------------------------------------
    do s = 1, ncmax
      xyzf_QMixDivMolWt(:,:,:,s) = xyzf_QMixBZ(:,:,:,s) / MolWtWet(s)
    end do
    
    xyz_EffMolWtBZ = (1.0d0 + sum(xyzf_QMixBZ,4) ) / ( MolWtDry * ((1.0d0 / MolWtDry) + sum(xyzf_QMixDivMolWt,4)) )
    
    !のり代の値を決める  
    !
    call SetMargin_xyz( xyz_EffMolWtBZ )
    
    !---------------------------------------------------------------    
    ! 温位
    !---------------------------------------------------------------
    xyz_PTempBZ = xyz_TempBZ * (PressBasis / xyz_PressBZ) ** (GasRDry / CpDry) 

    !のり代の値を決める  
    !
    call SetMargin_xyz( xyz_PTempBZ )
    
    !---------------------------------------------------------------    
    ! エクスナー関数
    !---------------------------------------------------------------
    xyz_ExnerBZ = xyz_TempBZ / xyz_PTempBZ    
    
    !のり代の値を決める
    !
    call SetMargin_xyz( xyz_ExnerBZ )

    !---------------------------------------------------------------    
    ! 密度
    !---------------------------------------------------------------
    xyz_DensBZ = PressBasis * (xyz_ExnerBZ ** (CvDry / GasRDry)) / (GasRDry * xyz_PTempBZ / xyz_EffMolWtBZ)
    
    !のり代の値を決める
    !
    call SetMargin_xyz( xyz_DensBZ )
    
    !---------------------------------------------------------------    
    ! 音速
    !---------------------------------------------------------------
    xyz_VelSoundBZ = sqrt ( CpDry * GasRDry * xyz_ExnerBZ * xyz_PTempBZ / (CvDry * xyz_EffMolWtBZ) )
    
    !のり代の値を決める
    !
    call SetMargin_xyz( xyz_VelSoundBZ )
   
    !----------------------------------------------------------
    ! BasicSet モジュールに値を設定
    !----------------------------------------------------------
    call BasicSet_Init( xyz_PressBZ, xyz_ExnerBZ, xyz_TempBZ, xyz_PTempBZ, xyz_DensBZ,  xyz_VelSoundBZ, xyzf_QMixBZ, xyz_EffMolWtBZ )
    
  end subroutine DetermineBZ


  subroutine deepconvInitialdata( cfgfile )

    use dc_iounit,     only : FileOpen    

    implicit none

    character(STRING), intent(in) :: cfgfile
    character(STRING)             :: FlaginitData = ""
    integer                       :: unit     !装置番号

    NAMELIST /initialdata_nml/ FlaginitData

    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=initialdata_nml)
    close(unit)

    if (FlaginitData == "default") then 
      ID_FlaginitData = ID_FlaginitData_default
    elseif (FlaginitData == "Yamasaki1983") then
      ID_FlaginitData = ID_FlaginitData_Yamasaki1983
    elseif (FlaginitData == "Takemi2007") then
      ID_FlaginitData = ID_FlaginitData_Takemi2007
    elseif (FlaginitData == "arare4") then
      ID_FlaginitData = ID_FlaginitData_arare4
    elseif (FlaginitData == "arare4_mmc") then
      ID_FlaginitData = ID_FlaginitData_arare4_mmc
    end if

  end subroutine deepconvInitialdata

end program ArareInitData