!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program ArareInit
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: init.f90,v 1.2.4.1 2006/01/13 06:55:04 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!非静力学モデル deepconv/arare の初期値を作成する. 
!初期値ファイルには 2 ステップ分のデータが必要なので, 
!最初の 1 ステップをオイラー法で計算する
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
! *方程式系は準圧縮系.
!
!== Future Plans
!
                                                                 !=end

program ArareInit
                                                                 !=begin
  !== Dependency
  use dc_trace,   only: BeginSub, EndSub, DbgMessage
  use dc_message
  use debugset,   only: debugset_init, DebugOn
  use gridset,    only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, &
       &                RegZMin, RegZMax
  use fileset,    only: fileset_init, ReStartFile
  use timeset,    only: timeset_init, NstepShort, DelTimeLong, DelTimeShort
  use nameset,    only: nameset_init
  use physset,    only: physset_init  
  use linlib,     only: linlib_init
  use bcset,      only: bcset_init
  use arareset,   only: arareset_init
  use basicset,   only: basicset_calc_init
  use disturbset, only: disturbset_calc_init, &
       &                ss_ExnerDisturb, ss_PotTempDisturb, &
       &                fs_VelXDisturb,  sf_VelZDisturb,    &
       &                ss_MixRtVDisturb,              &
       &                ss_MixRtCDisturb,              &
       &                ss_MixRtRDisturb
                                                                 !=end
  implicit none

  real(8), allocatable  :: ss_Exner_as(:,:)
  real(8), allocatable  :: fs_VelX_as(:,:)
  real(8), allocatable  :: sf_VelZ_as(:,:)
  real(8), allocatable  :: ss_Exner_ns(:,:)
  real(8), allocatable  :: fs_VelX_ns(:,:)
  real(8), allocatable  :: sf_VelZ_ns(:,:)

  real(8), allocatable  :: fs_AdvDiffX_nl(:,:)
  real(8), allocatable  :: sf_AdvDiffZ_nl(:,:)
  
  real(8), allocatable  :: fs_VelX_nl(:,:)
  real(8), allocatable  :: sf_VelZ_nl(:,:)
  real(8), allocatable  :: fs_VelX_al(:,:)
  real(8), allocatable  :: sf_VelZ_al(:,:)

  real(8), allocatable  :: ss_PotTemp_nl(:,:)
  real(8), allocatable  :: ss_Exner_nl(:,:)
  real(8), allocatable  :: ss_PotTemp_al(:,:)
  real(8), allocatable  :: ss_Exner_al(:,:)

  real(8), allocatable  :: ss_Km_al(:,:)
  real(8), allocatable  :: ss_Kh_al(:,:)
  real(8), allocatable  :: ss_Km_nl(:,:)
  real(8), allocatable  :: ss_Kh_nl(:,:)

  real(8), allocatable  :: ss_MixRtSat_nl(:,:)   !飽和水蒸気混合比 [t=nΔt]
  real(8), allocatable  :: ss_MixRtV_al(:,:) !水蒸気混合比 [t=(n-1)Δt]
  real(8), allocatable  :: ss_MixRtV_nl(:,:) !水蒸気混合比 [t=nΔt]
  real(8), allocatable  :: ss_MixRtC_al(:,:) !水蒸気混合比 [t=(n-1)Δt]
  real(8), allocatable  :: ss_MixRtC_nl(:,:) !水蒸気混合比 [t=nΔt]
  real(8), allocatable  :: ss_MixRtR_al(:,:)  !水蒸気混合比 [t=(n-1)Δt]
  real(8), allocatable  :: ss_MixRtR_nl(:,:)  !水蒸気混合比 [t=nΔt]

  real(8), allocatable  :: sf_Buoy_nl(:,:) ! 鉛直方向の運動方程式の浮力項 [t=nΔt]
  !--- 蒸発による雨水から水蒸気への変換
  real(8), allocatable  :: &
       & ss_EVrv_nl(:,:), & 
       !--- 併合成長による雲水から雨水への変換
       & ss_CNcr_nl(:,:),&
       !--- 凝結による水蒸気から雲水への変換
       & ss_CLcr_nl(:,:),& 
       !--- 雨粒終端速度
       & ss_VelZRain_nl(:,:),& 
       !--- 凝結による水蒸気から雲水への変換. 
       & ss_PRr_nl(:,:) 

  real(8)               :: ReStartTime(2)
  real(8)               :: Time
  integer               :: tau
  character(30)         :: cfgfile            ! NAMELIST ファイル


!!!
!!!初期化
!!!
  call getarg( 1, cfgfile )
  write(*,*) "Input NAMELIST file: ", cfgfile

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

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

  !--- 物理量の初期化
  call physset_init(cfgfile)
  
  !--- 格子点情報の初期化
  call gridset_init(cfgfile)
  
  !--- 時間積分の設定の初期化
  call timeset_init(cfgfile)
  
  !--- 実験名の初期化
  call nameset_init(cfgfile)
  
  !--- 計算の設定を読み込む
  call arareset_init(cfgfile)

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

  !--- 実験設定を読み込む
  call basicset_calc_init(cfgfile)

  !--- 実験設定を読み込む
  call disturbset_calc_init(cfgfile)
  
  !--- 行列計算ライブラリの初期化. 
  !--- 引数はスカラーの z 方向の計算領域の大きさ
  call linlib_init(cfgfile, RegZMax - RegZMin)
  

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

  !--- 擾乱場の初期化
  ss_Exner_nl   = ss_ExnerDisturb
  ss_PotTemp_nl = ss_PotTempDisturb
  fs_VelX_nl    = fs_VelXDisturb
  sf_VelZ_nl    = sf_VelZDisturb
  ss_MixRtV_nl = ss_MixRtVDisturb 
  ss_MixRtC_nl = ss_MixRtCDisturb 
  ss_MixRtR_nl = ss_MixRtRDisturb 

!!!
!!!数値積分. 1 ステップ分だけ回す. 
!!!
  !--- 渦粘性係数, 渦拡散係数の計算
  call EddyViscosity(                                     &
    & ss_Km_nl,                                        & !初期値はゼロ
    & DelTimeLong,                                     & 
    & fs_VelX_nl, sf_VelZ_nl, ss_Km_nl,                & !移流
    & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl,           & !生成, 拡散, 消散
    & ss_Km_nl, ss_Exner_nl,                           & !生成, 拡散, 消散
    & ss_MixRtV_nl, ss_MixRtC_nl,            & !生成, 拡散, 消散
    & ss_MixRtR_nl,                                &
    & ss_Km_al, ss_Kh_al )
  
  !=== 暖かい雨のパラメタリゼーション
  call WarmRainPrm( &
    & DelTimeLong, & 
    & fs_VelX_nl, sf_VelZ_nl, ss_Exner_nl, &
    & ss_PotTemp_nl, ss_PotTemp_nl, &
    & ss_MixRtV_nl, ss_MixRtV_nl, &
    & ss_MixRtC_nl, ss_MixRtC_nl, &
    & ss_MixRtR_nl,  ss_MixRtR_nl, &
    & ss_Kh_nl, &
    & ss_PotTemp_al, &
    & ss_MixRtV_al, &
    & ss_MixRtC_al, &
    & ss_MixRtR_al, &
    & ss_EVrv_nl, & 
    & ss_CNcr_nl,&
    & ss_CLcr_nl,& 
    & ss_VelZRain_nl,& 
    & ss_PRr_nl )

!  call Scalar(                              & 
!    & ss_PotTemp_nl,                            &
!    & DelTimeLong,                                 &
!    & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, & !移流
!    & ss_PotTemp_nl, ss_Kh_nl,                & !生成, 乱流拡散, 数値粘性
!    & ss_PotTemp_al)
  
  !--- 長い時間ステップでの速度 u の移流, 拡散, 粘性
  call AdvectDiffuseX(                     &
    & fs_VelX_nl, sf_VelZ_nl,           & !移流
    & fs_VelX_nl, sf_VelZ_nl, ss_Km_nl, & !乱流拡散, 数値粘性
    & fs_AdvDiffX_nl )
  
  !--- 長い時間ステップでの速度 w の移流, 拡散, 粘性
  call AdvectDiffuseZ(                          &
    & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, & !移流, 浮力
    & fs_VelX_nl, sf_VelZ_nl, ss_Km_nl,      & !乱流拡散, 数値粘性
    & ss_MixRtV_nl, ss_MixRtC_nl, ss_MixRtR_nl, &
    & sf_Buoy_nl, &
    & sf_AdvDiffZ_nl )
  
  !--- 短いタイムステップの初期値を作る
  fs_VelX_ns  = fs_VelX_nl
  sf_VelZ_ns  = sf_VelZ_nl
  ss_Exner_ns = ss_Exner_nl 
  
  !--- 短いタイムステップでの積分
  do tau = 1, NstepShort / 2
    !--- 速度 u の計算
    call VelocityX(                                             &
      & fs_VelX_ns,                                          &
      & DelTimeShort,                                        &
      & fs_AdvDiffX_nl, fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
      & fs_VelX_as )
    
    !--- エクスナー関数の計算
    call Exner( &
      & sf_AdvDiffZ_nl, fs_VelX_ns, fs_VelX_as, sf_VelZ_ns, ss_Exner_ns, &
      & ss_Exner_as )
    
    !--- 速度 w の計算
    call VelocityZ(  &
      & sf_VelZ_nl, &
      & DelTimeShort,    &
      & sf_AdvDiffZ_nl, fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, ss_Exner_as, &
      & sf_VelZ_as )
    
    !--- ループを回すための処置
    ss_Exner_ns  = ss_Exner_as 
    fs_VelX_ns   = fs_VelX_as
    sf_VelZ_ns   = sf_VelZ_as
  end do
  
  !--- 最終的な速度を長い時間ステップでのものとみなす
  ss_Exner_al  = ss_Exner_as 
  fs_VelX_al   = fs_VelX_as
  sf_VelZ_al   = sf_VelZ_as     

  !--- ファイル出力
  ReStartTime(1) = 0.0d0
  ReStartTime(2) = DelTimeLong
  call MakeReStartFile( ReStartFile, ReStartTime,  &
    & fs_VelX_nl,    sf_VelZ_nl,  ss_Exner_nl,  &
    & ss_PotTemp_nl, ss_Km_nl,    ss_Kh_nl,     &
    & ss_MixRtSat_nl, ss_MixRtV_nl,       &
    & ss_MixRtC_nl, ss_MixRtR_nl,      &
    & fs_VelX_al,    sf_VelZ_al,  ss_Exner_al,  &
    & ss_PotTemp_al, ss_Km_al,    ss_Kh_al,     &
    & ss_MixRtV_al,                        &
    & ss_MixRtC_al, ss_MixRtR_al )
  
contains

  
!!!
!!!配列の割り当て, 初期化
!!!
  subroutine ArareAlloc

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

    allocate(ss_Exner_as(DimXMin:DimXMax, DimZMin:DimZMax),  &
         & fs_VelX_as(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_VelZ_as(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_VelX_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_VelZ_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Exner_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTemp_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & 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), &
         & fs_AdvDiffX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_AdvDiffZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtV_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtV_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtC_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtC_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtR_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtR_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_MixRtSat_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_Buoy_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_EVrv_nl(DimXMin:DimXMax, DimZMin:DimZMax), & 
         & ss_CNcr_nl(DimXMin:DimXMax, DimZMin:DimZMax),&
         & ss_CLcr_nl(DimXMin:DimXMax, DimZMin:DimZMax),& 
         & ss_VelZRain_nl(DimXMin:DimXMax, DimZMin:DimZMax),& 
         & ss_PRr_nl(DimXMin:DimXMax, DimZMin:DimZMax), &         
         & ss_Km_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Kh_al(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Km_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Kh_nl(DimXMin:DimXMax, DimZMin:DimZMax) )

    fs_VelX_al = 0.0d0
    fs_VelX_nl = 0.0d0

    sf_VelZ_al = 0.0d0
    sf_VelZ_nl = 0.0d0

    ss_Exner_al = 0.0d0
    ss_Exner_nl = 0.0d0

    ss_PotTemp_al = 0.0d0
    ss_PotTemp_nl = 0.0d0

    fs_VelX_as = 0.0d0
    fs_VelX_ns = 0.0d0

    sf_VelZ_as = 0.0d0
    sf_VelZ_ns = 0.0d0

    ss_Exner_as = 0.0d0
    ss_Exner_ns = 0.0d0

    fs_AdvDiffX_nl = 0.0d0
    sf_AdvDiffZ_nl = 0.0d0

    ss_MixRtSat_nl = 0.0d0
    ss_MixRtV_al = 0.0d0 
    ss_MixRtV_nl = 0.0d0
    ss_MixRtC_al = 0.0d0
    ss_MixRtC_nl = 0.0d0
    ss_MixRtR_al  = 0.0d0
    ss_MixRtR_nl  = 0.0d0

    sf_Buoy_nl     = 0.0d0
    ss_EVrv_nl     = 0.0d0
    ss_CNcr_nl     = 0.0d0
    ss_CLcr_nl     = 0.0d0
    ss_VelZRain_nl = 0.0d0
    ss_PRr_nl      = 0.0d0
    
    ss_Km_nl = 0.0d0
    ss_Kh_nl = 0.0d0

    ss_Km_al = 0.0d0
    ss_Kh_al = 0.0d0

    Time = 0.0d0

    call EndSub("ArareAlloc")

  end subroutine ArareAlloc


end program ArareInit
