!= Program ArareMPI
!
! Authors::   SUGIYAMA Ko-ichiro ($B?y;39L0lO/(B), ODAKA Masatsugu ($B>.9b@5;L(B)
! Version::   $Id: ararempi.f90,v 1.13 2011-03-01 03:04:02 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! $BHs@ENO3X%b%G%k(B deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * $BJ}Dx<07O$O=`05=L7O(B.
!
!== Future Plans
!
!

program ararempi
  !
  !$BHs@ENO3X%b%G%k(B deepconv/arare. 
  !

  !----- $B%b%8%e!<%kFI$_9~$_(B ------

  !-----   $B7?@k8@(B, $BJ8;zNs=hM}(B   ----
  use dc_types,       only : STRING
  use dc_string,      only : StoA

  !-----   $B%a%C%;!<%8=PNO(B   -----
  use dc_message,     only: MessageNotify

  !  $B%3%^%s%I%i%$%s0z?t2r<a(B
  use argset,        only : argset_init  

  !-----    $B4IM}%b%8%e!<%k(B   -----
  !  $B2=3XNL7W;;%b%8%e!<%k(B
  use ChemCalc, only: ChemCalc_init
  use chemdata, only: chemdata_init

  !  MPI $B$N=i4|2=(B
  use mpiset,        only : mpii_init, mpii_end, myrank

  !  $BF~=PNO%U%!%$%kL>4IM}%b%8%e!<%k(B
  use filesetmpi,    only : fileset_init, &
    &                       InitFile

  !  $B%G%P%C%0=PNO4IM}%b%8%e!<%k(B
  use debugset,      only : debugset_init

  !  $B;~4V4IM}%b%8%e!<%k(B
  use timeset,       only : timeset_init, &
    &                       NstepLong, NstepShort, DelTimeLong, DelTimeShort, &
    &                       NstepDisp

  !  $B3J;RE@4IM}%b%8%e!<%k(B 
  use gridset,       only : gridset_init, &
    &                       DimXMin, DimXMax, DimZMin, DimZMax, SpcNum

  !  $B4pK\>l@_Dj%b%8%e!<%k(B
  use basicset,      only : basicset_init, &
    &                       xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, &
    &                       xz_VelSoundBasicZ

  !  $B@Q;;CM4IM}%b%8%e!<%k(B
  use StorePotTemp,      only : StorePotTemp_init, StorePotTempClean, &
    &                           StorePotTempCond
  use StoreMixRt,        only : StoreMixRt_init, StoreMixRtClean, &
    &                           StoreMixRtCond, StoreMixRtFill1, StoreMixRtFill2
  use StoreBuoy,         only : StoreBuoy_init, StoreBuoyClean
  use StoreStab,         only : StoreStab_init, StoreStabClean
  use StoreMom,          only : StoreMom_Init, StoreMomClean

  !  $B<>=a%k!<%A%s@_Dj%b%8%e!<%k(B
  use moistset,      only: moistset_init

  !-----    $B2<@A$1%b%8%e!<%k(B   -----
  !  $B?tCMK`;$7W;;%b%8%e!<%k(B 
  use damping,       only : damping_init, & 
    &                        DampSponge_xz, DampSponge_pz, DampSponge_xr

  !  $B;~4V@QJ,%U%#%k%?!<%b%8%e!<%k(B
  use timefilter,    only : AsselinFilter_xz, AsselinFilter_xr, &
    &                        AsselinFilter_pz, AsselinFilter_xza

  !  $B6-3&>r7oE,MQ%b%8%e!<%k(B
  use boundarympi,    only : BoundaryXCyc_xz, BoundaryZSym_xz,     &
    &                      BoundaryXCyc_xza, BoundaryZSym_xza,   &
    &                      BoundaryXCyc_pz, BoundaryZSym_pz,     &
    &                      BoundaryXCyc_xr, BoundaryZAntiSym_xr

  !  CFL $B>r7o3NG'%b%8%e!<%k(B
  use cflcheck,      only : CFLCheckTimeShort, &
    &                       CFLCheckTimeLongVelX, CFLCheckTimeLongVelZ

  !  $BIi$N<>=aNL$NJdE67W;;%b%8%e!<%k(B
  use fillnegative,  only : FillNegative_init, xza_FillNegative_xza

  !-----    $BF~=PNO%b%8%e!<%k(B   -----
  !  $B%j%9%?!<%H%U%!%$%kF~=PNO%b%8%e!<%k(B
  use RestartFileIOmpi, only : ReStartFile_Open, ReStartFile_OutPut, &
    &                       ReStartFile_Close, ReStartFile_Get

  !  $B%R%9%H%j%U%!%$%kF~=PNO%b%8%e!<%k(B
  use HistoryFileIOmpi, only : HistoryFile_Open, HistoryFile_OutPut, &
    &                       HistoryFile_Close

  !-----       $BNO3X2aDx(B        -----
  !  $BNO3X2aDx7W;;MQ4X?t%b%8%e!<%k(B
  use DynFunc,       only : xz_AdvScalar, xz_AdvKm, xza_AdvScalar, pz_AdvVelX, &
    &                       xr_Buoy, xr_AdvVelZ, pz_GradPi
  
  !  $BNO3X2aDx1"2rK!7W;;MQ4X?t%b%8%e!<%k(B
  use DynImpFunc,    only : xz_Exner_init, xz_Exner, xr_GradPi
  
  !-----       $BJ*M}2aDx(B        -----
  !  $B?tCM3H;67W;;MQ%b%8%e!<%k(B
  use NumDiffusion,  only : NumDiffusion_Init, xz_NumDiffScalar, xz_NumDiffKm, &
    &                       xza_NumDiffScalar, pz_NumDiffVelX, xr_NumDiffVelZ

  !  $BMpN.3H;67W;;MQ%b%8%e!<%k(B
  use Turbulence,   only : Turbulence_Init, &
    &                      xz_TurbScalar, xza_TurbScalar, pz_TurbVelX, &
    &                      xr_TurbVelZ  , xz_ShearKm    , xz_DispKm,   &
    &                      xz_DispHeat,   xz_TurbKm
  
  !  $BJ|<M6/@)7W;;MQ%b%8%e!<%k(B
  use Radiation,    only : Radiation_init,  &
    &                      xz_RadHeatConst  
  
  !  $BCOI=%U%i%C%/%97W;;MQ%b%8%e!<%k(B
!  use HeatFlux,     only : xz_HeatFluxBulk, xz_MixRtFluxBulk
  use HeatFlux,     only : xz_HeatFluxDiff, xza_MixRtFluxDiff

  !  $B<>=aK0OBD4@aK!7W;;MQ%b%8%e!<%k(B
  use MoistAdjust,  only : MoistAdjustSvapPress, MoistAdjustNH4SH

  !  $B1@J*M}%Q%i%a%?%j%<!<%7%g%s(B
  use WarmRainPrm,  only : WarmRainPrm_Init, xz_Rain2GasHeat, xza_Rain2Gas, &
    &                      xza_Rain2GasNH4SH, xz_Rain2GasHeatNH4SH,         &
    &                      xza_Cloud2Rain, xza_FallRain

  !  $B<>=a5$2t$NIbNO7W;;MQ%b%8%e!<%k(B
  use MoistBuoyancy,only : MoistBuoy_Init, xz_BuoyMoistKm, xr_BuoyMolWt, &
    &                      xr_BuoyDrag

  !  $B0BDjEY$N7W;;(B
  use ECCM, only: ECCM_Stab

  !$B0EL[$N7?@k8@6X;_(B
  implicit none

  !$BFbItJQ?t(B
  character(80) :: cfgfile
  real(8), allocatable :: pz_VelXBl(:,:)
  real(8), allocatable :: pz_VelXNl(:,:)
  real(8), allocatable :: pz_VelXAl(:,:)
  real(8), allocatable :: pz_VelXNs(:,:)
  real(8), allocatable :: pz_VelXAs(:,:)
  real(8), allocatable :: xr_VelZBl(:,:)
  real(8), allocatable :: xr_VelZNl(:,:)
  real(8), allocatable :: xr_VelZAl(:,:)
  real(8), allocatable :: xr_VelZNs(:,:)
  real(8), allocatable :: xr_VelZAs(:,:)
  real(8), allocatable :: xz_ExnerBl(:,:)
  real(8), allocatable :: xz_ExnerNl(:,:)
  real(8), allocatable :: xz_ExnerAl(:,:)
  real(8), allocatable :: xz_ExnerNs(:,:)
  real(8), allocatable :: xz_ExnerAs(:,:)
  real(8), allocatable :: xz_PotTempWork(:,:)
  real(8), allocatable :: xz_PotTempAllBl(:,:)
  real(8), allocatable :: xz_PotTempAllNl(:,:)
  real(8), allocatable :: xz_PotTempBl(:,:)
  real(8), allocatable :: xz_PotTempNl(:,:)
  real(8), allocatable :: xz_PotTempAl(:,:)
  real(8), allocatable :: xz_KmBl(:,:)
  real(8), allocatable :: xz_KmNl(:,:)
  real(8), allocatable :: xz_KmAl(:,:)
  real(8), allocatable :: xz_KhBl(:,:)
  real(8), allocatable :: xz_KhNl(:,:)
  real(8), allocatable :: xz_KhAl(:,:)
  real(8), allocatable :: xza_MixRtWork(:,:,:)
  real(8), allocatable :: xza_MixRtAllBl(:,:,:)
  real(8), allocatable :: xza_MixRtAllNl(:,:,:)
  real(8), allocatable :: xza_MixRtBl(:,:,:)
  real(8), allocatable :: xza_MixRtNl(:,:,:)
  real(8), allocatable :: xza_MixRtAl(:,:,:)

  real(8), allocatable :: pz_AccelVelXNl(:,:)
  real(8), allocatable :: xr_AccelVelZNl(:,:)
  real(8), allocatable :: xza_DelMixRt(:,:,:)

  real(8)              :: Time
  real(8)              :: ReStartTime(2)
  real(8), allocatable :: DelTimeLFrog(:)
  real(8)              :: DelTimeEular
  integer, allocatable :: NStepEular(:)
  integer              :: NStepLFrog

  real(8), allocatable :: Km1(:,:), Km2(:,:), Km3(:,:), Km4(:,:), Km5(:,:), Km6(:,:), Km7(:,:)

  integer :: & 
    & t,     & 
    & tau,   & ! do $B%k!<%WJQ?t(B ; do loop variable  
    & t1,    & ! do $B%k!<%WJQ?t(B ; do loop variable  
    & t2       ! do $B%k!<%WJQ?t(B ; do loop variable 

  logical :: CloudFlag = .true.
  logical :: TurbFlag  = .true.


  !$B%3%^%s%I%i%$%s0z?t$N2r<a(B
  !
  call mpii_init

  ! $B=i4|2=<jB3$-(B Initialize procedure 
  !
  call MainInit

  !----------------------------------------------------------------------
  ! $B?tCM@QJ,(B
  !----------------------------------------------------------------------
!  do t1 = 1, NstepLFrog / NstepDisp
!    do t2 = 1, NstepDisp
  do t = 1, NstepLFrog   

      !$B;~9o$N@_Dj(B
      Time = Time + DelTimeLong
!      t = (t1 - 1) * NstepDisp + t2
      
      !----------------------------------------------------------------
      ! $B12G4@-78?t(B, $B123H;678?t$r5a$a$k(B.
      !----------------------------------------------------------------
      if (TurbFlag) then 
      xz_KmAl =                                                       &
        & xz_KmBl                                                     &
        & + DelTimeLFrog(t)                                           &
        &   * (                                                       &
        &     + xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)               &
        &     + xz_NumDiffKm(xz_KmBl)                                 &
        &     + xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) &
        &     + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)             &   
        &     + xz_TurbKm(xz_KmBl)                                    &
        &     + xz_DispKm(xz_KmBl)                                    &
        &    )

      Km1 = xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)
      call StoreDKmDtAdv(Km1)
      Km2 = xz_NumDiffKm(xz_KmBl)
      call StoreDKmDtDiff(Km2)
      Km3 = xz_BuoyKm(xz_PotTempBl)
      call StoreDKmDtBuoyT(Km3)
      Km4 = xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) 
      call StoreDKmDtBuoyM(Km4 - km3)
      Km5 = xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)
      call StoreDKmDtShear(Km5)
      Km6 = xz_TurbKm(xz_KmBl)
      call StoreDKmDtTurb(Km6)
      Km7 = xz_DispKm(xz_KmBl)
      call StoreDKmDtDisp(Km7)

      !$BCM$N>e8B2<8B$N@_Dj(B
      !  * $BCM$O@5$K$J$k$3$H$rJ]>Z$9$k(B
      !  * $BCM$N>e8B$O(B 800 $B$H$9$k(B. $BCfEg7r2p(B(1994, $B3X0LO@J8(B)$B;2>H(B
      xz_KmAl = max( 0.0d0, min( xz_KmAl, 800.0d0 ) )

      else
        xz_KmAl = 0.0d0
        xz_KmNl = 0.0d0
      end if      

      !$B6-3&>r7o(B
      call BoundaryXCyc_xz( xz_KmAl )
      call BoundaryZSym_xz( xz_KmAl )
      
      !$B123H;6Dj?t$r7h$a$k(B
      xz_KhAl = 3.0d0 * xz_KmAl
      
      !----------------------------------------------------------------
      ! $B290L$N0\N.7W;;(B.
      !----------------------------------------------------------------    
      !$B;~4V@QJ,(B
      xz_PotTempAllNl = xz_PotTempNl + xz_PotTempBasicZ
      xz_PotTempAllBl = xz_PotTempBl + xz_PotTempBasicZ
      xz_PotTempAl =                                                  &
        &   xz_PotTempBl                                              &
        & + DelTimeLFrog(t)                                           &
        &   * (                                                       &
        &     + xz_AdvScalar( xz_PotTempAllNl,  pz_VelXNl, xr_VelZNl) &
        &     + xz_TurbScalar(xz_PotTempAllBl,  xz_KhBl)              &
        &     + xz_NumDiffScalar(xz_PotTempBl)                        &
        &     + xz_DispHeat( xz_KmBl )                                &
        &     + xz_RadHeatConst( xz_ExnerBl )                         &
        &     + xz_HeatFluxDiff( xz_PotTempNl )                       &
!!      &     + xz_HeatFluxBulk( xz_PotTempNl )                     &
!!      &     + xz_NewtonCool( xz_PotTempBl )                       &
        &      )

      !$B6-3&>r7o(B
      call BoundaryXCyc_xz( xz_PotTempAl )
      call BoundaryZSym_xz( xz_PotTempAl )
      
      !----------------------------------------------------------------
      ! $B6E=L@.J,$N:.9gHf$N0\N.7W;;(B.
      !----------------------------------------------------------------
      xza_MixRtAllNl = xza_MixRtNl + xza_MixRtBasicZ
      xza_MixRtAllBl = xza_MixRtBl + xza_MixRtBasicZ
      xza_MixRtAl =                                                & 
        &   xza_MixRtBl                                            &
        & + DelTimeLFrog(t)                                        &
        &   * (                                                    &
        &    + xza_AdvScalar(xza_MixRtAllNl,  pz_VelXNl, xr_VelZNl)&
        &    + xza_TurbScalar(xza_MixRtAllBl, xz_KhBl)             &
        &    + xza_NumDiffScalar(xza_MixRtBl)                      &
        &    + xza_FallRain(xza_MixRtBl)                           &
        &    + xza_MixRtFluxDiff(xza_MixRtNl)                      &
        &   )

      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !$B0\N.$K$h$C$FIi$K$J$C$?ItJ,$rKd$a$k(B
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = xza_FillNegative_xza( xza_MixRtWork ) 
      
      !$BKd$a$?(B/$B:o$C$?NL$rJ]4I(B
      call StoreMixRtFill1( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) )

      !$B6-3&>r7o(B
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )
      
!      ! Negative values mixing ratios are corrected. 
!      xza_MixRtWork = xza_MixRtAl
!      xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtWork )
!      call StoreMixRtFill2( &
!        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
!        & )

      !$B6-3&>r7o(B
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !-------------------------------------------------------------
      ! $BCH$+$$1+$N%Q%i%a%?%j%<!<%7%g%s(B.
      !   $B1@(B<-->$B1+(B $B$NJQ49$r9T$&(B
      !-------------------------------------------------------------
      !$B1@$+$i1+$X$NJQ49J,$rDI2C$9$k(B. 
      !  xza_Cloud2Rain $B4X?t$O(B, $B0z?t$H$7$F;~4V9o$_$r<h$k$3$H$G(B, 
      !  $B;~4V@QJ,CM$r=PNO$9$k(B. 

      !$B$3$l$^$G$NCM$r:n6HG[Ns$KJ]4I(B
      xza_MixRtWork = xza_MixRtAl
      
      !$B1+$X$NJQ2=NL$r7W;;(B
      xza_MixRtAl   = xza_MixRtWork &
        &             + xza_Cloud2Rain( xza_MixRtAl, DelTimeLFrog(t) )
      
      !$B1@$+$i1+$X$NJQ49NL$rJ]4I(B
      call StoreMixRtCond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) 
      
      !$B6-3&>r7o(B
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !-------------------------------------------------------------
      ! $B<>=aK0OBD4@aK!(B
      !   $B?e>x5$(B<-->$B1@$NJQ49$r9T$&(B.
      !-------------------------------------------------------------
      !$B$3$l$^$G$NCM$r:n6HG[Ns$KJ]4I(B
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork  = xza_MixRtAl
      
      !$B<>=aD4@aK!$rE,MQ(B
      call MoistAdjustSvapPress( xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) )
      call MoistAdjustNH4SH( xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) )
      
      !$B<>=aD4@aK!$K$h$k290L$H:.9gHf$NJQ2=NL$rJ]4I(B
      call StorePotTempCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) )
      call StoreMixRtCond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) 

      !$B6-3&>r7o(B
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !-------------------------------------------------------------
      ! $BCH$+$$1+$N%Q%i%a%?%j%<!<%7%g%s(B.
      !   $B>x5$(B<-->$B1+(B $B$NJQ49$r9T$&(B
      !-------------------------------------------------------------
      !$B1+$+$i>x5$$X$NJQ49$KH<$&:.9gHf$NJQ2=$r7W;;(B
      !  xza_Rain2Gas $B4X?t$O(B, $B0z?t$H$7$F;~4V9o$_$r<h$k$3$H$G(B, 
      !  $B;~4V@QJ,CM$r=PNO$9$k(B. 
     
      !$B$3$l$^$G$NCM$r:n6HG[Ns$KJ]4I(B
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork = xza_MixRtAl
      
      !$B1+$+$i>x5$$X$N:.9gHfJQ2=$r5a$a$k(B
      !  $B290L$N7W;;$K$*$$$F(B, $B:.9gHfJQ2=$,I,MW$H$J$k$?$a(B, 
      !  $B:.9gHfJQ2=$r(B 1 $B$D$NG[Ns$H$7$FMQ0U$9$k(B.
      xza_DelMixRt = 0.0d0
      xza_DelMixRt =                                                    &
        & (                                                             &
        &   + xza_Rain2Gas(                                             &
        &        xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
        &       )                                                       &
        &   + xza_Rain2GasNH4SH(                                        &
        &        xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
        &       )                                                       &
        &  )    
     
      !$B290L$N7W;;(B. $B1+$+$i>x5$$X$NJQ49$KH<$&@xG.!&H?1~G.$rDI2C(B.
      xz_PotTempAl =                                                                      &
        & xz_PotTempWork                                                                  &
        & + (                                                                             &
        &      + xz_Rain2GasHeat(xz_PotTempAl, xz_ExnerNl, xza_DelMixRt, DelTimeLFrog(t)) &
        &      + xz_Rain2GasHeatNH4SH(xz_ExnerNl, xza_DelMixRt, DelTimeLFrog(t))          &
        &    )
     
      !$B:.9gHf$N7W;;(B. $B1+$+$i>x5$$X$NJQ49J,$rDI2C(B
      xza_MixRtAl   = xza_MixRtWork + xza_DelMixRt
      
      !$B1+$+$i>x5$$NCM$NJ]4I(B
      call StorePotTempCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) ) 
      call StoreMixRtCond( xza_DelMixRt / DelTimeLFrog(t) ) 
     
      !$B6-3&>r7o(B
      call BoundaryXCyc_xz( xz_PotTempAl )
      call BoundaryZSym_xz( xz_PotTempAl )
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !-------------------------------------------------------------
      ! $BD9$$;~4V%9%F%C%W$G$N(B, $BB.EY$N0\N.(B, $B3H;6(B, $B?tCMG4@-(B, $BIbNO(B
      !-------------------------------------------------------------
      pz_AccelVelXNl =                                    &
        & + pz_AdvVelX(pz_VelXNl, xr_VelZNl)              &
        & + pz_TurbVelX(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
        & + pz_NumDiffVelX(pz_VelXBl)
      
      xr_AccelVelZNl =                                    &
        & + xr_AdvVelZ(xr_VelZNl, pz_VelXNl)              &
        & + xr_Buoy(xz_PotTempNl)                         &    
        & + xr_BuoyMolWt(xza_MixRtNl)                     &
        & + xr_BuoyDrag(xza_MixRtNl)                      &
        & + xr_TurbVelZ(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
        & + xr_NumDiffVelZ(xr_VelZBl)                        

      !-------------------------------------------------------------
      ! $BC;$$;~4V%9%F%C%W$N=i4|CM:n@.(B
      !-------------------------------------------------------------
      pz_VelXNs  = pz_VelXBl
      xr_VelZNs  = xr_VelZBl
      xz_ExnerNs = xz_ExnerBl

      !-------------------------------------------------------------
      ! $BC;$$;~4V%9%F%C%W$N@QJ,(B. $B%*%$%i!<K!$rMxMQ(B
      !-------------------------------------------------------------
      Eular: do tau = 1, NstepEular(t)

        !-------------------------------------------------------------
        ! $BB.EY(B u $B$N7W;;(B
        !-------------------------------------------------------------
        pz_VelXAs = &
          & pz_VelXNs                                          &
          &  + DelTimeEular                                    &
          &    * (                                             &
          &     - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs)  &
          &     + pz_AccelVelXNl                               &
          &     )
        
        !$B6-3&>r7o(B
        call BoundaryXCyc_pz( pz_VelXAs )
        call BoundaryZSym_pz( pz_VelXAs )
        
        !-------------------------------------------------------------
        ! $B%(%/%9%J!<4X?t$N7W;;(B
        !-------------------------------------------------------------
        xz_ExnerAs = xz_Exner( &
          & xr_AccelVelZNl,    &
          & pz_VelXNs,         &
          & pz_VelXAs,         &
          & xr_VelZNs,         &
          & xz_ExnerNs)
        
        !$B6-3&>r7o(B
        call BoundaryXCyc_xz( xz_ExnerAs )
        call BoundaryZSym_xz( xz_ExnerAs )
        
        !-------------------------------------------------------------
        ! $BB.EY(B w $B$N7W;;(B
        !-------------------------------------------------------------
        xr_VelZAs =                                                    &
          & xr_VelZNs                                                  &
          &  + DelTimeEular                                            &
          &    * (                                                     &
          &     - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) &
          &     + xr_AccelVelZNl                                       &
          &     )
        
        !$B6-3&>r7o(B
        call BoundaryXCyc_xr( xr_VelZAs )
        call BoundaryZAntiSym_xr( xr_VelZAs )

        !-------------------------------------------------------------
        ! $BC;$$;~4V%9%F%C%W$N%k!<%W$r2s$9$?$a$N=hCV(B
        !-------------------------------------------------------------
        xz_ExnerNs  = xz_ExnerAs
        pz_VelXNs   = pz_VelXAs
        xr_VelZNs   = xr_VelZAs
        
      end do Eular

      !----------------------------------------------------------------
      ! $B:G=*E*$JC;$$;~4V%9%F%C%W$G$NCM$rD9$$;~4V%9%F%C%W$G$NCM$H$_$J$9(B
      !----------------------------------------------------------------
      xz_ExnerAl  = xz_ExnerAs
      pz_VelXAl   = pz_VelXAs
      xr_VelZAl   = xr_VelZAs

      !-------------------------------------------------------------
      ! $B%"%C%;%j%s$N;~4V%U%#%k%?(B.  Nl $B$NCM$r%U%#%k%?%j%s%0(B
      !-------------------------------------------------------------
      call AsselinFilter_xz(xz_ExnerAl, xz_ExnerNl, xz_ExnerBl)
      call AsselinFilter_pz(pz_VelXAl, pz_VelXNl, pz_VelXBl )
      call AsselinFilter_xr(xr_VelZAl, xr_VelZNl, xr_VelZBl )
      call AsselinFilter_xz(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl)
      call AsselinFilter_xz(xz_KmAl, xz_KmNl, xz_KmBl)
      call AsselinFilter_xza(xza_MixRtAl, xza_MixRtNl, xza_MixRtBl)
    
      !-------------------------------------------------------------
      ! $B%9%]%s%8AX(B
      !-------------------------------------------------------------
      call DampSponge_pz( pz_VelXAl,    pz_VelXBl,    DelTimeLFrog(t) )
      call DampSponge_xr( xr_VelZAl,    xr_VelZBl,    DelTimeLFrog(t) )
      call DampSponge_xz( xz_ExnerAl,   xz_ExnerBl,   DelTimeLFrog(t) )
      call DampSponge_xz( xz_PotTempAl, xz_PotTempBl, DelTimeLFrog(t) )
      call DampSponge_xz( xz_KmAl,      xz_KmBl,      DelTimeLFrog(t) )

      !
      ! Negative values mixing ratios are corrected. 
      !
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtWork )
      call StoreMixRtFill2( &
        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
        & )

      !--------------------------------------------------------------
      ! $B2r@OCM(B
      !---------------------------------------------------------------
      call ECCM_Stab( xz_PotTempAl, xz_ExnerAl, xza_MixRtAl )
      
      !----------------------------------------------------------------
      ! $BD9$$;~4V%9%F%C%W$N%k!<%W$r2s$9$?$a$N=hCV(B
      !----------------------------------------------------------------
      xz_PotTempBl = xz_PotTempNl
      xza_MixRtBl  = xza_MixRtNl
      xz_ExnerBl   = xz_ExnerNl
      pz_VelXBl    = pz_VelXNl
      xr_VelZBl    = xr_VelZNl
      xz_KmBl      = xz_KmNl
      xz_KhBl      = xz_KhNl
      
      xz_PotTempNl = xz_PotTempAl
      xza_MixRtNl  = xza_MixRtAl
      xz_ExnerNl   = xz_ExnerAl
      pz_VelXNl    = pz_VelXAl
      xr_VelZNl    = xr_VelZAl
      xz_KmNl      = xz_KmAl
      xz_KhNl      = xz_KhAl

    end do
    !----------------------------------------------------------------
    ! $B%U%!%$%k=PNO(B
    !----------------------------------------------------------------
    ! CFL $B>r7o$N%A%'%C%/(B
!    if (myrank == 0) then 
!      call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )
!    end if
!    call CFLCheckTimeLongVelX( pz_VelXNl )
!    call CFLCheckTimeLongVelZ( xr_VelZNl )
    
    !$B%R%9%H%j%U%!%$%k=PNO(B
    call HistoryFile_Output(  &
      &  Time,                &
      &  xz_PotTempNl,        &
      &  xz_ExnerNl,          &
      &  pz_VelXNl,           &
      &  xr_VelZNl,           &
      &  xza_MixRtNl,         &
      &  xz_KmNl,             &
      &  xz_KhNl              &
      & )

    !$B@Q;;CM$N%/%j%"(B
    call StorePotTempClean
    call StoreMixRtClean
    call StoreBuoyClean
    call StoreStabClean
    call StoreMomClean

    !----------------------------------------------------------------
    ! $B%j%9%?!<%H%U%!%$%k$N:n@.(B
    !----------------------------------------------------------------
    call ReStartFile_Open( t1 )
    call ReStartFile_OutPut(                            &
      & Time - DelTimeLong,                             &
      & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
      & xza_MixRtBl,  xz_KmBl,    xz_KhBl               )
    call ReStartFile_OutPut(                            &
      & Time,                                           &
      & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, &
      & xza_MixRtNl,  xz_KmNl,    xz_KhNl               )
    call ReStartFile_Close
    
  end do
  
  !----------------------------------------------------------------
  ! $B=PNO%U%!%$%k$N%/%m!<%:(B
  !----------------------------------------------------------------
  call HistoryFile_Close

  !----------------------------------------------------
  ! MPI END
  !----------------------------------------------------
  call mpii_end

contains

  subroutine MainInit

    implicit none

    !NAMELIST $B%U%!%$%kL>$NFI$_9~$_(B
    if (myrank == 0) then 
      call MessageNotify( "M", "main", "Namelist file is '%c'", c1=trim(cfgfile) )
    end if
    call argset_init(cfgfile)
    
    !$B%G%P%C%0@_Dj(B
    call debugset_init(cfgfile, myrank)
    
    !$BJ*<AFC@-$N=i4|2=(B
    call chemdata_init()
    
    !$B;~9o$K4X$9$k@_Dj$N=i4|2=(B
    !  NAMELIST $B$+$iI,MW$J>pJs$rFI$_<h$j(B, $B;~4V4XO"$NJQ?t$N@_Dj$r9T$&(B. 
    call timeset_init(cfgfile)

    !$B3J;RE@>pJs$N=i4|2=(B
    !  NAMELIST $B$+$i>pJs$rF@$F(B, $B3J;RE@$r7W;;$9$k(B
    call gridset_init(cfgfile)
    
    !$B2=3X7W;;%k!<%A%s$N=i4|2=(B
    call chemcalc_init()
    
    !$B4pK\>l$N>pJs$N=i4|2=(B
    !  NAMELIST $B$+$i>pJs$rF@$F(B, $B4pK\>l$r@_Dj$9$k(B.
    call basicset_init(cfgfile)
    
    !I/O $B%U%!%$%kL>$N=i4|2=(B
    !  NAMELIST $B%U%!%$%kL>$r;XDj$7(B, deepconv/arare $B$N(B
    !  $B=PNO%U%!%$%kL>$r(B NAMELIST $B$+$iF@$k(B
    call fileset_init(cfgfile)
    
    !$B<>=a%k!<%A%s$N6&M-JQ?t$N=i4|2=(B
    call moistset_init()
    
    !$B@Q;;CM$rJ]4I$9$k$?$a$N%b%8%e!<%k$N=i4|2=(B
    !  NAMELIST $B$+$i>pJs$rF@$F(B, $B4pK\>l$r@_Dj$9$k(B.
    call StorePotTemp_init( )
    call StoreMixRt_init( )
    call StoreBuoy_init( )  
    call StoreStab_init( )  
    call StoreMom_init( )  
    
    !$BFbItJQ?t$N=i4|2=(B. $B$H$j$"$($:%<%m$rF~$l$FCM$r3NDj$5$;$F$*$/(B. 
    call ArareAlloc
    
    !$BM=JsJQ?t$N=i4|2=(B
    !  ReStartFile $B$,@_Dj$5$l$F$$$k>l9g$K$O%U%!%$%k$rFI$_9~$_(B, 
    !  $B@_Dj$5$l$F$$$J$$>l9g$K$O(B, $B%G%U%)%k%H$N4pK\>l$H>qMp>l$r:n$k(B. 
    if (trim(InitFile) /= '') then    
      
      !$B4pK\>l(B, $B>qMp>l$N=i4|CM$r(B netCDF $B%U%!%$%k$+$i<hF@$9$k(B.
      call ReStartFile_Get(            &
        & ReStartTime,                 &
        & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
        & xza_MixRtBl,  xz_KmBl,    xz_KhBl,              &
        & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, &
        & xza_MixRtNl,  xz_KmNl,    xz_KhNl               )
      
    else
      !$B%G%U%)%k%H@_Dj$N4pK\>l(B, $B>qMp>l$r:n@.$9$k(B. 
      call BasicEnv()  ! $B?eJ?0lMM$J$N$G(B single CPU $BHG$HF1$8%k!<%A%s$r;H$($k(B. 
      call DisturbEnvMPI(cfgfile,                            &
        & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
        & xza_MixRtBl,  xz_KmBl,    xz_KhBl               )
      
      ! 1 $B%k!<%WL\$@$1$O(B, bl $B$H(B nl $B$NCM$rF1$8$K$7$F$*$/(B. 
      ! 1 $B%9%F%C%WL\$O%*%$%i!<K!$G2r$/I,MW$,$"$k$,(B, 
      ! 1 $B%9%F%C%WL\$H$=$l0J30$N%9%F%C%W$rJL!9$K%3!<%G%#%s%0$7$?$/$J$$0Y(B
      xz_PotTempNl = xz_PotTempBl
      xz_ExnerNl   = xz_ExnerBl
      pz_VelXNl    = pz_VelXBl
      xr_VelZNl    = xr_VelZBl
      xza_MixRtNl  = xza_MixRtBl
      xz_KmNl      = xz_KmBl
      xz_KhNl      = xz_KhBl
      
    end if
    
    !----------------------------------------------------------------------
    ! $B%Q%C%1!<%87?%b%8%e!<%k$N=i4|2=(B
    !   $B%G%U%)%k%H$NCM$+$iJQ99$9$kI,MW$N$"$k%k!<%A%s$N$_=i4|2=(B
    !----------------------------------------------------------------------
    call Damping_Init( cfgfile )      !$BGH$N8:?j78?t$N=i4|2=(B
    !  call NumDiffusion_Init3( cfgfile ) !$B?tCM3H;69`$N=i4|2=(B
    call NumDiffusion_Init()          !$B?tCM3H;69`$N=i4|2=(B
    call Turbulence_Init()            !$BMpN.7W;;$N=i4|2=(B
    call WarmRainPrm_Init( cfgfile )  !$BCH$+$$1+$N%Q%i%a%?%j%<!<%7%g%s$N=i4|2=(B
    call FillNegative_Init( xza_MixRtBasicZ, xz_DensBasicZ)
    call Radiation_Init( cfgfile )    !$BJ|<M6/@)$N=i4|2=(B
    call MoistBuoy_Init()             !$BJ,;RNL$KBP$9$kIbNO7W;;%k!<%A%s$N=i4|2=(B
    call xz_Exner_Init()              !$B1"2rK!$N=i4|2=(B  
    

    !----------------------------------------------------------------------
    ! $B;~9o$H%k!<%W2s?t$N=i4|2=(B 
    !----------------------------------------------------------------------
    NstepLFrog   = NstepLong 
    NstepEular   = NstepShort 
    DelTimeLFrog = DelTimeLong * 2.0d0 
    DelTimeEular = DelTimeShort
    if ( trim(InitFile) /= '') then    
      Time = ReStartTime(2)               !$B%j%9%?!<%H3+;O;~9o(B
    else
      Time = 0.0d0                          !$B7W;;3+;O;~9o(B
      NstepEular(1)   = NstepShort /2 ! 1 $B%k!<%WL\$@$1(B
      DelTimeLFrog(1) = DelTimeLong         ! 1 $B%k!<%WL\$@$1(B
    end if

    !----------------------------------------------------------------
    ! $B%U%!%$%kF~=PNO(B
    !----------------------------------------------------------------
    !$B%U%!%$%k%*!<%W%s(B
    call HistoryFile_Open( )
    
    if ( trim(InitFile) /= '') then    
      call HistoryFile_Output(  &
        &  ReStartTime(2),      &
        &  xz_PotTempNl,        &
        &  xz_ExnerNl,          &
        &  pz_VelXNl,           &
        &  xr_VelZNl,           &
        &  xza_MixRtNl,         &
        &  xz_KmNl,             &
        &  xz_KhNl            )
    end if
    
    !$B;~9o(B 0 $B$N>l9g$K$O=PNO(B
    if ( Time == 0 ) then
      call HistoryFile_Output(  &
        &  Time,                &
        &  xz_PotTempNl,        &
        &  xz_ExnerNl,          &
        &  pz_VelXNl,           &
        &  xr_VelZNl,           &
        &  xza_MixRtNl,         &
        &  xz_KmNl,             &
        &  xz_KhNl            )
    end if
    
    !----------------------------------------------------------------------
    ! $B@_Dj$N%A%'%C%/(B
    !----------------------------------------------------------------------
    !CFL $B>r7o$N%A%'%C%/(B
    call CFLCheckTimeShort( xz_VelSoundBasicZ )
    
  end subroutine MainInit

!--------------------------------------------------------------------------
  subroutine ArareAlloc
    !
    !$B=i4|2=$H$7$F(B, $BG[Ns$rDj5A$7(B, $BCM$H$7$F%<%m$rBeF~$9$k(B.
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BG[Ns3d$jEv$F(B
    allocate(                                                  &
      & pz_VelXBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xr_VelZBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xz_ExnerBl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax),          &
!
      & xz_PotTempWork(DimXMin:DimXMax, DimZMin:DimZMax),      &
      & xz_PotTempAllBl(DimXMin:DimXMax, DimZMin:DimZMax),     &
      & xz_PotTempAllNl(DimXMin:DimXMax, DimZMin:DimZMax),     &
      & xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempAl(DimXMin:DimXMax, DimZMin:DimZMax),        &
!
      & xz_KmBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xz_KhBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xza_MixRtWork(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
      & xza_MixRtAllBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),&
      & xza_MixRtAllNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),&
      & xza_MixRtBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
      & xza_MixRtNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
      & xza_MixRtAl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
!
      & pz_AccelVelXNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xr_AccelVelZNl(DimXMin:DimXMax, DimZMin:DimZMax),        & 
!
      & xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),  &
!
      & DelTimeLFrog(NstepLong), NStepEular(NStepLong)               ) 

    pz_VelXNs = 0.0d0;     pz_VelXAs = 0.0d0    
    xr_VelZNs = 0.0d0;     xr_VelZAs = 0.0d0
    xz_ExnerNs = 0.0d0;    xz_ExnerAs = 0.0d0

    pz_VelXBl = 0.0d0;     pz_VelXNl = 0.0d0;     pz_VelXAl = 0.0d0
    xr_VelZBl = 0.0d0;     xr_VelZNl = 0.0d0;     xr_VelZAl = 0.0d0
    xz_ExnerBl = 0.0d0;    xz_ExnerNl = 0.0d0;    xz_ExnerAl = 0.0d0
    xz_KmBl = 0.0d0;       xz_KmNl = 0.0d0;       xz_KmAl = 0.0d0
    xz_KhBl = 0.0d0;       xz_KhNl = 0.0d0;       xz_KhAl = 0.0d0
    xz_PotTempBl = 0.0d0;  xz_PotTempNl = 0.0d0;  xz_PotTempAl = 0.0d0
    xza_MixRtBl = 0.0d0;    xza_MixRtNl = 0.0d0;    xza_MixRtAl = 0.0d0

    pz_AccelVelXNl = 0.0d0
    xr_AccelVelZNl = 0.0d0

    xza_DelMixRt = 0.0d0

    DelTimeEular = 0.0d0
    DelTimeLFrog = 0.0d0; 
    NStepEular = 0.0d0
    NstepLFrog = 0.0d0
    
  end subroutine ArareAlloc
    
end program ararempi
