
!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, only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, &
       & RegZMin, RegZMax
  use fileset, only: fileset_init
  use timeset, only: timeset_init, NstepTShort, DelTShort
  use physprm, only: physprm_init
  use expname, only: expname_init
  use basicprm, only: basicprm_init
  use linlib,  only: linlib_init
  use if_display

  implicit none

  real(8), allocatable  :: ss_Exner_A(:,:)
  real(8), allocatable  :: fs_VelX_A(:,:)
  real(8), allocatable  :: sf_VelZ_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)               :: Time
  integer               :: t
!  integer, parameter    :: NStepMax = 50
!  integer               :: NStepInterval

  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 timeset_init
  
  !--- 物理量の初期化 
  call physprm_init

  !--- 実験名の初期化
  call expname_init

  !--- 基本場の値を設定
  call basicprm_init

  !--- 行列計算ライブラリの初期化. 
  !--- 引数はスカラーの z 方向の計算領域の大きさ
  call linlib_init(RegZMax - RegZMin)  

  !--- ファイル出力を行う時間間隔を決定
!  if (NStepTShort > NStepMax) then
!     NStepInterval = NStepTShort / NStepMax
!  else
!     NStepInterval = 1
!  end if
  
  !--- 出力ファイルの定義
  call OpenRegFile
  call OpenDimFile

!!!
!!! CFL 条件のチェック
!!!
  call cflchk


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

  !--- 初期擾乱の設定
  call ArareDisturb


!!!
!!!計算
!!!
  !--- とりあえず長い時間ステップは無視
  fs_Fx_B = 0.0d0
  sf_Fz_B = 0.0d0


  !--- 最初の時間
  Time = 0.0d0
  call OutputRegFile
  call OutputDimFile

  !--- 短いタイムステップでの積分
  do t = 1, NstepTShort
     Time = Time + DelTShort
     write(*,*) "Time: ", Time
     
     !--- 速度の収束
     call VelDiv(fs_VelX_B, sf_VelZ_B, ss_VelDiv_B)
     
     !--- 速度 u の計算
     call VelXHeVi(ss_Exner_B, ss_VelDiv_B, fs_Fx_B, fs_VelX_B, fs_VelX_A)
     
     !--- エクスナー関数の計算
     call ExnerHeVi(   &
          & ss_VelDiv_B, sf_Fz_B, fs_VelX_A, sf_VelZ_B, ss_Exner_B, ss_Exner_A)
     
     !--- 速度 w の計算
     call VelZHeVi(    &
          & ss_Exner_B, ss_Exner_A, ss_VelDiv_B, sf_Fz_B, sf_VelZ_B, sf_VelZ_A)
     
     !--- ループを回すための処置
     ss_Exner_B  = ss_Exner_A 
     fs_VelX_B = fs_VelX_A
     sf_VelZ_B = sf_VelZ_A
     
     !--- ファイル出力
!     if (mod(t, NStepInterval) == 0) then
        call OutputRegFile
        call OutputDimFile
!     end if
     
  end do
  
  !--- 出力ファイルのクローズ
  call CloseRegFile
  call CloseDimFile


contains
  
  subroutine ArareAlloc
    allocate(ss_Exner_A(DimXMin:DimXMax, DimZMin:DimZMax),  &
         & fs_VelX_A(DimXMin:DimXMax, DimZMin:DimZMax)  , &
         & sf_VelZ_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
    use basicprm, only: ss_ExnerBasicZ

    ss_Exner_A = ss_ExnerBasicZ
    ss_Exner_B = ss_Exner_A

    fs_VelX_A = 0.0d0
    fs_VelX_B = 0.0d0

    sf_VelZ_A = 0.0d0
    sf_VelZ_B = 0.0d0

    fs_Fx_B = 0.0d0
    sf_Fz_B = 0.0d0

    ss_VelDiv_B = 0.0d0

    Time = 0.0d0
  end subroutine ArareArrayInit


  subroutine ArareDisturb
    use gridset, only: s_X, s_Z, Xmin, Xmax, Zmin, Zmax, &
         & DimXMin, DimXMax, DimZMin, DimZMax
    use physprm, only: pi

    integer :: i, k
    real(8) :: mu
    real(8) :: sigma

    !--- エクスナー関数の擾乱
    mu = ( Zmax - Zmin ) * 1.0d0 / 5.0d0
    sigma = ( Zmax - Zmin ) / 20.0d0
    
    do k = DimZMin, DimZMax
       ss_Exner_A(:,k) = ss_Exner_A(:,k) &
            & + 10.0d0 * ss_Exner_A(:,k) &
           & * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) &
            & / (dsqrt(2.0d0 * pi) * sigma)
    end do

!    mu = ( Xmax - Xmin ) * 4.0d0  / 5.0d0
!    sigma = ( Xmax - Xmin ) / 20.0d0
!    
!    do i = DimXMin, DimXMax
!       ss_Exner_A(i,:) = ss_Exner_A(i,:) &
!            & + 100.0d0 * ss_Exner_A(i,:) &
!            & * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) &
!           & / (dsqrt(2.0d0 * pi) * sigma)
!    end do

    ss_Exner_B = ss_Exner_A

  end subroutine ArareDisturb


  !--- 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','t'/),                           &
         & dimsizes=(/NX+1, NZ, 0/),                       &
         & longnames=(/'X-coordinate',                     &
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/), 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','t'/),                           &
         & dimsizes=(/NX, NZ+1, 0/),                       &
         & longnames=(/'X-coordinate',                     & 
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/),                          &
         & 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','t'/),                         &
         & dimsizes=(/NX, NZ, 0/),                       &
         & longnames=(/'X-coordinate',                   &
         &             'Z-coordinate',                   &
         &             'Time        '/),                 &
         & units=(/'m','m','s'/), 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','t'/),        &
         & longname='nondimensional pressure', units='1',&
         & xtype='double',                               & 
         & history = ss_reg)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelX', dims=(/'x','z','t'/),             &
         & longname='zonal velocity',         &
         & units='m/s', xtype='double',                & 
         & history = fs_reg)        
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelZ', dims=(/'x','z','t'/),             &
         & longname='vertical velocity',         &
         & units='m/s', xtype='double',                & 
         & history = sf_reg)    
    
  end subroutine OpenRegFile


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

    call HistoryPut('t', Time, ss_reg)
    call HistoryPut('t', Time, fs_reg)
    call HistoryPut('t', Time, sf_reg)

    call HistoryPut('Exner', &
         & ss_Exner_A(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('VelX', &
         & fs_VelX_A(RegXMin:RegXMax, RegZMin+1:RegZMax), fs_reg)
    call HistoryPut('VelZ', &
         & sf_VelZ_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','t'/),                           &
         & dimsizes=(/N, N, 0/),                           &
         & longnames=(/'X-coordinate',                     &
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/), 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','t'/),                           &
         & dimsizes=(/N, N, 0/),                           &
         & longnames=(/'X-coordinate',                     & 
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/),                          &
         & 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','t'/),                         &
         & dimsizes=(/N, N, 0/),                         &
         & longnames=(/'X-coordinate',                   &
         &             'Z-coordinate',                   &
         &             'Time        '/),                 &
         & units=(/'m','m','s'/), 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','t'/),        &
         & longname='nondimensional pressure', units='1',&
         & xtype='double',                               & 
         & history = ss_dim)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelX', dims=(/'x','z','t'/),        &
         & longname='zonal velocity',                   &
         & units='m/s', xtype='double',                 & 
         & history = fs_dim)        
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelZ', dims=(/'x','z','t'/),        &
         & longname='vertical velocity',                &
         & units='m/s', xtype='double',                 & 
         & history = sf_dim)    
    
  end subroutine OpenDimFile


  subroutine OutputDimFile
    call HistoryPut('t', Time, ss_dim)
    call HistoryPut('t', Time, fs_dim)
    call HistoryPut('t', Time, sf_dim)

    call HistoryPut('Exner', ss_Exner_A, ss_dim)
    call HistoryPut('VelX', fs_VelX_A, fs_dim)
    call HistoryPut('VelZ', sf_VelZ_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
