!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine 
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!$BCH$+$$1+$N%P%k%/K!(B
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!  * $B;~4V$K4X$7$F$OA0?J:9J,$;$M$P$J$i$J$$$+$J(B. 
!
subroutine WarmRainPrm( &
  & DelTime, &
  & fs_VelX_nl, sf_VelZ_nl, ss_Exner_nl, &
  & ss_PotTemp_nl, ss_PotTemp_bl, &
  & ss_MixRtV_nl, ss_MixRtV_bl, &
  & ss_MixRtC_nl, ss_MixRtC_bl, &
  & ss_MixRtR_nl,  ss_MixRtR_bl, &
  & ss_Kh_bl, &
  & 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, &
  & ss_MixRtSat_nl    )

  !== Dependency
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, &
    &                 RegXMin, RegXMax, RegZMin, RegZMax
  use bcset,    only: ss_BC
  use average,  only: ss_avr_fs, ss_avr_sf, fs_avr_ss, sf_avr_ss
  use basicset, only: ss_ExnerBasicZ, &
       &              ss_DensBasicZ,  &
       &              ss_PotTempBasicZ, &
       &              ss_MixRtVBasicZ,  &
       &              ss_PressBasicZ  
  use physset,  only: PressSfc, GasRUniv, GasRDry, MolWtDry, MolWtWet, CpDry, epsilon
  use arareset, only: Min
  use differentiate_center4, only: sf_dz_ss

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

  !== Input
  real(8), intent(in) :: DelTime !$B@QJ,;~4V(B
  real(8), intent(in) :: ss_Exner_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B%(%/%9%J!<4X?t(B
  real(8), intent(in) :: ss_PotTemp_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B290L(B
  real(8), intent(in) :: ss_PotTemp_bl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B290L(B
  real(8), intent(in) :: ss_MixRtV_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B?e>x5$NL(B
  real(8), intent(in) :: ss_MixRtV_bl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B?e>x5$NL(B
  real(8), intent(in) :: ss_MixRtC_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B1@N3:.9gHf(B
  real(8), intent(in) :: ss_MixRtC_bl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !$B1@N3:.9gHf(B
  real(8), intent(in) :: ss_MixRtR_nl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                                          !$B1+N3:.9gHf(B
  real(8), intent(in) :: ss_MixRtR_bl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                                          !$B1+N3:.9gHf(B
  real(8), intent(in) :: fs_VelX_nl(DimXMin:DimXMax, DimZMin:DimZMax)   
                                                          !$BB.EY(B u
  real(8), intent(in) :: sf_VelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)   
                                                          !$BB.EY(B w
  real(8), intent(in) :: ss_Kh_bl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B12G4@-78?t(B

  !== Out
  real(8), intent(Out)  :: ss_PotTemp_al(DimXMin:DimXMax, DimZMin:DimZMax) 
                                            !$B290L(B
  real(8), intent(Out)  :: ss_MixRtV_al(DimXMin:DimXMax, DimZMin:DimZMax) 
                                            !$B?e>x5$:.9gHf(B
  real(8), intent(Out)  :: ss_MixRtC_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B1@N3:.9gHf(B
  real(8), intent(Out)  :: ss_MixRtR_al(DimXMin:DimXMax, DimZMin:DimZMax) 
                                            !$B1+N3:.9gHf(B
  real(8), intent(Out)  :: ss_EVrv_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B>xH/$K$h$k1+?e$+$i?e>x5$$X$NJQ49(B
  real(8), intent(Out)  :: ss_CNcr_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                            !$BJ;9g@.D9$K$h$k1@?e$+$i1+?e$X$NJQ49(B
  real(8), intent(Out)  :: ss_CLcr_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B6E7k$K$h$k?e>x5$$+$i1@?e$X$NJQ49(B
  real(8), intent(Out)  :: ss_VelZRain_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B1+N3=*C<B.EY(B
  real(8), intent(Out)  :: ss_PRr_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                            !$B6E7k$K$h$k?e>x5$$+$i1@?e$X$NJQ49(B. 
  real(8), intent(Out)  :: ss_MixRtSat_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                            !$BK0OB>x5$05(B

  !== Work       
  real(8) :: ss_Gamma_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                             ! $B@xG.(B / ( $BHfG.(B * $B%(%/%9%J!<4X?t(B )
  real(8) :: ss_PressAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)   
                                             ! $B05NO(B($B>qMp(B+$BJ?6Q(B)
  real(8) :: ss_ExnerAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)   
                                             ! $BL5<!8505NO(B($B>qMp(B+$BJ?6Q(B)
  real(8) :: ss_TempAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)    
                                             ! $B29EY(B($B>qMp(B+$BJ?6Q(B)
  real(8) :: ss_PotTempAll_nl(DimXMin:DimXMax, DimZMin:DimZMax) 
                                             ! $B290L(B($B>qMp(B+$BJ?6Q(B)
  real(8) :: ss_MixRtVAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                             ! $B?e>x5$:.9gHf(B($B>qMp(B+$BJ?6Q(B)
  real(8) :: ss_MixRtVSat_nl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                             ! $BK0OB>x5$05:.9gHf(B
  real(8) :: ss_SvapPress_nl(DimXMin:DimXMax, DimZMin:DimZMax)  
                                             ! $BK0OB>x5$05%b%kHf(B
  real(8) :: ss_PotTempAll_al(DimXMin:DimXMax, DimZMin:DimZMax) 
                                             ! $BA4290L(B
  real(8) :: ss_MixRtVAll_al(DimXMin:DimXMax, DimZMin:DimZMax)  
                                             ! $BA4?e>x5$(B

  real(8) :: MixRtV_tmp
  real(8) :: MixRtC_tmp
  real(8) :: PotTemp_tmp

  real(8) :: PressAll_adj
  real(8) :: ExnerAll_adj
  real(8) :: DelMixRtV_adj

  real(8) :: TempAll_tmp
  real(8) :: PotTempAll_tmp
  real(8) :: MixRtVAll_tmp

  real(8) :: SvapPress_tmp
  real(8) :: MixRtVSat_tmp
  real(8) :: Gamma_tmp

  real(8) :: PotTemp_adj
  real(8) :: MixRtV_adj
  real(8) :: MixRtC_adj

  real(8) :: DMixRtSatDPotTemp_tmp
  real(8) :: DMixRtSatDSvapPress_tmp 

  real(8) :: DensGround         ! $BCOI=LLL)EY(B
  integer :: i, k, num          ! $BE:$(;z(B  
  real(8), parameter :: K1 = 1.0d-3  !$B%1%9%i!<$N78?t(B
  real(8), parameter :: K2 = 2.2d0   !$B%1%9%i!<$N78?t(B
  real(8), parameter :: A  = 1.0d-3  !$B1@?eNL$NogCM(B
  real(8), parameter :: N0 = 5.0d7
  real(8), parameter :: D0 = 3.66d-1

  
  !=== $B%5%V%k!<%A%s$N=i4|2=(B
  call BeginSub("Scalar", &
    &           fmt="%c",              &
    &           c1="Time integration of scalar component.")
  
  !=== $BJQ?t$KNm$rBeF~(B
  call warmrainprm_init
  
  !=== $B29EY(B, $B05NO(B, $B:.9gHf$NA4NL$r5a$a$k(B
  !==== $B>qMp@.J,$HJ?6Q@.J,$NB-$7;;(B
  ss_PotTempAll_nl = ss_PotTemp_nl  + ss_PotTempBasicZ
  ss_ExnerAll_nl   = ss_Exner_nl    + ss_ExnerBasicZ
  ss_MixRtVAll_nl  = ss_MixRtV_nl   + ss_MixRtVBasicZ
  ss_TempAll_nl    = ss_ExnerAll_nl * ss_PotTempAll_nl
  
  !=== $BK0OB>x5$05$KH<$&J*M}NL$N7W;;(B 
  !==== $BK0OB>x5$05$N%b%kHf(B
  ss_SvapPress_nl = ss_svap_ss(ss_TempAll_nl)
  
  !==== $BK0OB>x5$05$+$iK0OB:.9gHf$r7W;;(B
  ss_MixRtVSat_nl = &
    & ((MolWtWet/MolWtDry) * ss_SvapPress_nl) &
    & / (PressSfc * (ss_ExnerAll_nl **(CpDry / GasRDry)) - ss_SvapPress_nl)
 
  !=== $BJQ?t$N=i4|2=(B
  ss_Gamma_nl   = &
    & ss_LatentHeat_ss(ss_TempAll_nl) / (ss_ExnerAll_nl * CpDry) 
  DensGround = ss_DensBasicZ(RegXMin, 0)   ! X $B:BI8$O$I$3$G$bF1$8(B
    

!!! 
!!! $B7W;;$KI,MW$H$J$k(B, $B?e>x5$(B--$B1@N3(B--$B1+N34V$NJQ4978?t$r5a$a$k(B
!!!
  !=== $BJ;9g@.D9$K$h$k1@?e$+$i1+?e$X$NJQ49(B
  !==== Berry (1968) $B$N%Q%i%a%?%j%<!<%7%g%s(B
!  where (ss_MixRtC_nl > A)
!    ss_CNcr_nl = K1 * (ss_MixRtC_nl - A)
!  end where
  where (ss_MixRtC_nl > Min)
    ss_CNcr_nl =  &
      & ss_DensBasicZ * (ss_MixRtC_nl ** 3.0d0) * 1.0d6 &
      & / (60.0d0 * 2.0d0 * ss_MixRtC_nl                &
      &    + 60.0d0 * 2.66d-8 * N0 / (ss_DensBasicZ * D0)) 
  end where


  !=== $B6E7k$K$h$k?e>x5$$+$i1@?e$X$NJQ49(B
  !==== Kessler (1969) $B$N%Q%i%a%?%j%<!<%7%g%s(B
  where (ss_MixRtC_nl > Min .AND. ss_MixRtR_nl > Min)
    ss_CLcr_nl = K2 * ss_MixRtC_nl * (ss_MixRtR_nl ** 0.875d0)
  end where


  !=== $B>xH/$K$h$k1+?e$+$i?e>x5$$X$NJQ49(B
  !====== $BCfEg7r2p(B (1994) $B$GMxMQ$7$?Dj<0(B. 
  where (ss_MixRtVSat_nl > ss_MixRtVAll_nl  .AND. ss_MixRtR_nl > Min )
    ss_EVrv_nl = &
      &  4.85d-2                                        &
      &  * (ss_MixRtVSat_nl - ss_MixRtVAll_nl )         &
      &  * ((ss_DensBasicZ * ss_MixRtR_nl) ** 0.65d0 )
  end where

  
  !=== $B6E7k$K$h$k?e>x5$$+$i1@?e$X$NJQ49(B. $B=*C<B.EY$r7W;;$7$F$+$iI>2A(B.
  !==== $B1+N3=*C<B.EY(B
  where (ss_MixRtR_nl > Min)
    ss_VelZRain_nl = 12.2d0 * (ss_MixRtR_nl ** 0.125d0) 
  end where

  !==== $B6E7k$K$h$k?e>x5$$+$i1@?e$X$NJQ49(B. 
  !===== Dens $B$N(B avr $B$r<h$C$F$+$i3d$k$H(B, $B%<%m3d$,@8$8$k$N$GCm0U(B
  where (ss_MixRtR_nl > Min)
    ss_PRr_nl &
      & = ss_avr_sf(  &
      &     sf_dz_ss(ss_DensBasicZ * ss_VelZRain_nl * ss_MixRtR_nl) &
      &   ) / ss_DensBasicZ 
  end where

  
!!!
!!!$BAjJQ2=$rI=8=$9$kItJ,$r=|$$$F290L(B, $B:.9gHf$N<0$r2r$/(B
!!!
  !=== $B;~4V@QJ,(B. $B;~4VJ}8~$K%j!<%W%U%m%C%0(B, $B6u4VJ}8~$KCf?4:9J,(B
  !==== $B$3$3$G$O1@!&9_?e2aDx$O=|$/(B

  !=== $B290L(B
  ss_PotTemp_al =                                                     &
    &   ss_PotTemp_bl                                                 &
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_PotTemp_nl, fs_VelX_nl)                  &
    &     - ss_AdvZ_ss_sf(ss_PotTemp_nl, sf_VelZ_nl)                  &
    &     - ss_AdvZ_ss_sf(ss_PotTempBasicZ, sf_VelZ_nl)               &
    &     - ss_Gamma_nl * ss_EVrv_nl                                  &
    &     + ss_TurbXZ_ss_ss(ss_PotTemp_bl, ss_Kh_bl)                  &
    &     + ss_TurbZ_ss_ss(ss_PotTempBasicZ, ss_Kh_bl)                &
    &     + ss_NumDiffXZ_ss(ss_PotTemp_bl)                            &
    &    )
  
  !=== $B?e>x5$:.9gHf(B
  ss_MixRtV_al =                                                      &
    &   ss_MixRtV_bl                                                  &  
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_MixRtV_nl, fs_VelX_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_MixRtV_nl, sf_VelZ_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_MixRtVBasicZ, sf_VelZ_nl)                &
    &     + ss_TurbXZ_ss_ss(ss_MixRtV_bl, ss_Kh_bl)                   & 
    &     + ss_TurbZ_ss_ss(ss_MixRtVBasicZ, ss_Kh_bl)                 & 
    &     + ss_NumDiffXZ_ss(ss_MixRtV_bl)                             &
    &    )

  !=== $B1@N3:.9gHf(B
  ss_MixRtC_al =                                                      &     
    &   ss_MixRtC_bl                                                  &
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_MixRtC_nl, fs_VelX_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_DensBasicZ * ss_MixRtC_nl, sf_VelZ_nl)   &
    &       / ss_DensBasicZ                                           &
    &     - ss_AdvXZ_ss_fs_sf(ss_MixRtC_nl, fs_VelX_nl, sf_VelZ_nl)   & 
    &     + ss_TurbXZ_ss_ss  (ss_MixRtC_bl, ss_Kh_bl)                 &
    &     + ss_NumDiffXZ_ss  (ss_MixRtC_bl)                           &
    &    )
  
  !=== $B1+N3:.9gHf(B
  ss_MixRtR_al  =                                                     &
    &   ss_MixRtR_bl                                                  &
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_MixRtR_nl, fs_VelX_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_DensBasicZ * ss_MixRtR_nl, sf_VelZ_nl)   &
    &       / ss_DensBasicZ                                           &
    &      + ss_TurbXZ_ss_ss(ss_MixRtR_bl, ss_Kh_bl)                  &
    &      + ss_NumDiffXZ_ss( ss_MixRtR_bl )                          &
    &     )
  

  !==== $B3NG'(B
  write(*,*) "MinRtR", minval(ss_MixRtR_al), &
    & maxval(ss_MixRtR_al), &
    & sum(ss_DensBasicZ(RegXMin:RegXMax, RegZMin:RegZMax) &
    &     *  ss_MixRtR_al(RegXMin:RegXMax, RegZMin:RegZMax)) 
  write(*,*) "MinRtC", minval(ss_MixRtC_al), &
    & maxval(ss_MixRtC_al), &
    & sum(ss_DensBasicZ(RegXMin:RegXMax, RegZMin:RegZMax) &
    &     * ss_MixRtC_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtV", minval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & maxval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & sum(ss_DensBasicZ(RegXMin:RegXMax, RegZMin:RegZMax) &
    &     * (ss_MixRtV_al(RegXMin:RegXMax, RegZMin:RegZMax) &
    &        + ss_MixRtVBasicZ(RegXMin:RegXMax, RegZMin:RegZMax)) )


!!!
!!! $B1@!&9_?e2aDx$N9MN8(B. $BCM$,Ii$K$J$i$J$$$h$&$K5$$r$D$1$k(B
!!! 
  !=== $BDL>o=hM}(B
  where (ss_MixRtC_al - DelTime * (ss_CNcr_nl + ss_CLcr_nl) >= Min) 
    !=== $B1+N3:.9gHf(B
    ss_MixRtR_al = ss_MixRtR_al + DelTime * (ss_CNcr_nl + ss_CLcr_nl)

    !=== $B1@N3:.9gHf(B
    ss_MixRtC_al = ss_MixRtC_al - DelTime * (ss_CNcr_nl + ss_CLcr_nl)  
  end where

  !=== $BNc30=hM}(B: $B1@N3:.9gHf0J>e$K6E=L$,@8$8$F$7$^$&>l9g(B
  where (ss_MixRtC_al - DelTime * (ss_CNcr_nl + ss_CLcr_nl) < Min) 
    !=== $B1+N3:.9gHf(B
    ss_MixRtR_al  = ss_MixRtR_al + ss_MixRtC_al  

    !=== $B1@N3:.9gHf(B
    ss_MixRtC_al = 0.0d0 
  end where


  !=== $BDL>o=hM}(B
  where (ss_MixRtR_al - DelTime * (ss_EVrv_nl - ss_PRr_nl) >= Min) 
    !=== $B290L(B
    ss_PotTemp_al = ss_PotTemp_al - DelTime * ss_Gamma_nl * ss_EVrv_nl   
  
    !=== $B?e>x5$:.9gHf(B
    ss_MixRtV_al = ss_MixRtV_al + DelTime * ss_EVrv_nl   
    
    !=== $B1+N3:.9gHf(B
    ss_MixRtR_al  = ss_MixRtR_al - DelTime * (ss_EVrv_nl - ss_PRr_nl)
  end where
  
  !=== $BNc30=hM}(B: $B1+N3:.9gHf0J>e$K>xH/$,@8$8$k>l9g(B
  where (ss_MixRtR_al - DelTime * (ss_EVrv_nl - ss_PRr_nl) < Min)
    !=== $B290L(B
    ss_PotTemp_al = ss_PotTemp_al - ss_Gamma_nl * ss_MixRtR_al 
  
    !=== $B?e>x5$:.9gHf(B
    ss_MixRtV_al = ss_MixRtV_al + ss_MixRtR_al 
    
    !=== $B1+N3:.9gHf(B
    ss_MixRtR_al  = 0.0d0
  end where

  !==== $B3NG'(B
  write(*,*) "MinRtR", minval(ss_MixRtR_al), &
    & maxval(ss_MixRtR_al), sum(ss_MixRtR_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtC", minval(ss_MixRtC_al), &
    & maxval(ss_MixRtC_al), sum(ss_MixRtC_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtV", minval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & maxval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & sum(ss_MixRtV_al(RegXMin:RegXMax, RegZMin:RegZMax))



!!!
!!! $B<>=aK0OBD4@aK!$rMQ$$$??e>x5$(B--$B1@N34V$NJQ49(B
!!!
!  MoistAjastX: do i = DimXMin, DimXMax
!    MoistAjastZ: do k = DimZMin, DimZMax
!      
!      !=== $B:#$^$G$KF@$i$l$?>qMp@.J,$NCM$r;CDjNL$H$_$J$9(B. $BE:$(;z$rJQ99(B
!      MixRtV_tmp   = ss_MixRtV_al(i,k)
!      MixRtC_tmp   = ss_MixRtC_al(i,k)
!      PotTemp_tmp  = ss_PotTemp_al(i,k)
!      
!      !=== $B<>=aK0OBK!$G$O05NO$OJQ2=$5$;$J$$(B. 
!      !==== $B<>=aK0OBD4@a$K$h$C$F?e>x5$%b%kHf$,JQ2=$9$l$P(B, 
!      !==== $BHfG.$H5$BNDj?t$NCM$,JQ2=$7(B, $B%(%/%9%J!<4X?t$NCM$OJQ2=$9$k$O$:(B. 
!      ExnerAll_adj = ss_ExnerAll_nl(i,k)
!      PressAll_adj = ss_PressAll_nl(i,k)
!      
!      !=== $B%k!<%W2s?t$N=i4|2=(B
!      num = 0
!      
!      Adjasting: do 
!        !=== $B>qMp@.J,$HJ?6Q@.J,$NOB$r5a$a$k(B        
!        PotTempAll_tmp = PotTemp_tmp    + ss_PotTempBasicZ(i,k) 
!        MixRtVAll_tmp  = MixRtV_tmp     + ss_MixRtVBasicZ(i,k) 
!        TempAll_tmp    = PotTempAll_tmp * ExnerAll_adj
!        
!        !==== $BK0OB>x5$05(B
!        SvapPress_tmp = s_svap_s(TempAll_tmp)
!  
!        !==== $BK0OB>x5$05$+$iK0OB:.9gHf$r7W;;(B
!        MixRtVSat_tmp = &
!          & ((MolWtWet/MolWtDry) * SvapPress_tmp) &
!          & / (PressSfc * (ExnerAll_adj**(CpDry / GasRDry)) - SvapPress_tmp)
!        
!        !=== $BJQ?t$N=i4|2=(B
!        Gamma_tmp = s_LatentHeat_s(TempAll_tmp) / (ExnerAll_adj * CpDry)
!        
!        !=== $B?e>x5$J,05$,K0OB>x5$05$h$jB?$$(B or $B1@N3:.9gHf$,@5(B
!        if (MixRtVAll_tmp - MixRtVSat_tmp > Min  .OR.        &
!          & MixRtC_tmp > Min) then 
!
!          write(*,*) i,k, "Now Adjasting.... "
!          write(*,*) i,k, MixRtVAll_tmp - MixRtVSat_tmp
!          write(*,*) i,k, MixRtVAll_tmp, MixRtC_tmp
!          
!          DMixRtSatDSvapPress_tmp =  &
!            & PressSfc                              &
!            & * (MolWtWet / MolWtDry)               &
!            & * (ExnerAll_adj ** (CpDry / GasRDry)) &
!            & / ((PressSfc * (ExnerAll_adj ** (CpDry / GasRDry)) &
!            &     - SvapPress_tmp) ** 2.0d0)
!          
!          DMixRtSatDPotTemp_tmp = &
!            & DMixRtSatDSvapPress_tmp * s_DSvapDT_s(TempAll_tmp) * ExnerAll_adj
!          
!          DelMixRtV_adj =  &
!            & (MixRtVAll_tmp - MixRtVSat_tmp) &
!            &  / ( 1.0d0 + Gamma_tmp * DMixRtSatDPotTemp_tmp )          
!
!          write(*,*) i,k,"Gamma, DelQv", Gamma_tmp, DelMixRtV_adj 
!
!          PotTemp_adj = PotTemp_tmp + Gamma_tmp * DelMixRtV_adj 
!          MixRtV_adj  = MixRtV_tmp - DelMixRtV_adj             
!          MixRtC_adj  = MixRtC_tmp + DelMixRtV_adj            

!          !=== $B1@N3:.9gHf$,Ii$K$J$C$?$H$-$O(B, $B1@N3:.9gHf$r@5$K$9$k(B. 
!          if ( MixRtC_adj <= Min) then 
!            write(*,*) "MixRtC_adj <= Min"
!            ss_PotTemp_al(i,k) = PotTemp_tmp - Gamma_tmp * MixRtC_tmp
!            ss_MixRtV_al(i,k)  = MixRtV_tmp + MixRtC_tmp 
!            ss_MixRtC_al(i,k)  = 0.0d0
!            exit Adjasting
!          end if
!
!          !=== $B1@N3:.9gHf$,<}B+$7$F$$$?$i7+$jJV$7$rCf;_$9$k(B
!          if (abs((MixRtC_adj - MixRtC_tmp) / MixRtC_adj) <= Min .OR. &
!            & abs(MixRtC_adj - MixRtC_tmp) <= Min) then 
!            ss_PotTemp_al(i,k) = PotTemp_adj                 
!            ss_MixRtV_al(i,k)  = MixRtV_adj
!            ss_MixRtC_al(i,k)  = MixRtC_adj
!            exit Adjasting
!          end if
!          
          !=== $B%k!<%W2s?t$,2s$j$9$.$?$i%(%i!<$rEG$$$FB3$1$k(B
!         if ( num > 20 ) then 
!           write(*,*) "Looping is too much, num > 20"
!            write(*,*) "ss_PotTemp_al(i,k) =     ", PotTemp_adj
!            write(*,*) "ss_MixRtV_al(i,k) = ", MixRtV_adj
!            write(*,*) "ss_MixRtC_al(i,k) = ", MixRtC_adj
!            ss_PotTemp_al(i,k) = PotTemp_adj     
!            ss_MixRtV_al(i,k)  = MixRtV_adj
!            ss_MixRtC_al(i,k)  = MixRtC_adj
!            exit Adjasting
!          end if
!
!        else 
!
!          !=== $B<}B+$7$?$H$7$F=*N;(B
!          ss_PotTemp_al(i,k) = PotTemp_tmp - Gamma_tmp * MixRtC_tmp
!          ss_MixRtV_al(i,k) = MixRtV_tmp + MixRtC_tmp
!          ss_MixRtC_al(i,k) = 0.0d0          
!          exit Adjasting
!        end if

!         !---- $B7+$jJV$7$N$?$a$NJQ?tDj5A(B
!        PotTemp_tmp = PotTemp_adj
!        MixRtV_tmp  = MixRtV_adj
!        MixRtC_tmp  = MixRtC_adj
!        num         = num + 1
        
!      end do Adjasting
!      
!    end do MoistAjastZ
!  end do MoistAjastX
  

  !==== $B3NG'(B
  write(*,*) "MinRtR", minval(ss_MixRtR_al), &
    & maxval(ss_MixRtR_al), sum(ss_MixRtR_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtC", minval(ss_MixRtC_al), &
    & maxval(ss_MixRtC_al), sum(ss_MixRtC_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtV", minval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & maxval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & sum(ss_MixRtV_al(RegXMin:RegXMax, RegZMin:RegZMax))


  !=== $BCM$,Ii$K$J$C$?$H$-$N=hM}(B
  !==== $BA4NL$r5a$a$k(B
!  ss_PotTempAll_al = ss_PotTemp_al + ss_PotTempBasicZ
!  ss_MixRtVAll_al  = ss_MixRtV_al  + ss_MixRtVBasicZ

  !==== $BIi$K$J$C$?ItJ,$N7jKd$a$r9T$&(B
!  call FillNegative( ss_PotTempAll_al )
!  call FillNegative( ss_MixRtVAll_al )

!  where ( ss_PotTempAll_al <= Min )
!    ss_PotTempAll_al = Min
!  end where
!  where ( ss_MixRtVAll_al <= Min )
!    ss_MixRtVAll_al = Min
!  end where
  
  !==== $B>qMp$NCM$KLa$9(B
!  ss_PotTemp_al = ss_PotTempAll_al - ss_PotTempBasicZ
!  ss_MixRtV_al  = ss_MixRtVAll_al  - ss_MixRtVBasicZ

  !$B3NG'(B
  !==== $B3NG'(B
  write(*,*) "MinRtR", minval(ss_MixRtR_al), &
    & maxval(ss_MixRtR_al), sum(ss_MixRtR_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtC", minval(ss_MixRtC_al), &
    & maxval(ss_MixRtC_al), sum(ss_MixRtC_al(RegXMin:RegXMax, RegZMin:RegZMax))
  write(*,*) "MinRtV", minval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & maxval(ss_MixRtV_al + ss_MixRtVBasicZ), &
    & sum(ss_MixRtV_al(RegXMin:RegXMax, RegZMin:RegZMax))
  

  !=== $B6-3&>r7o(B
  call boundary(ss_BC, ss_PotTemp_al)
  call boundary(ss_BC, ss_MixRtV_al)
  call boundary(ss_BC, ss_MixRtC_al)
  call boundary(ss_BC, ss_MixRtR_al)  


  !=== $BK0OB:.9gHf$N=PNO(B
  ss_SvapPress_nl = ss_svap_ss((ss_PotTemp_al + ss_PotTempBasicZ) * ss_ExnerAll_nl)
  
  ss_MixRtSat_nl = &
    & ((MolWtWet/MolWtDry) * ss_SvapPress_nl) &
    & / (PressSfc * (ss_ExnerAll_nl **(CpDry / GasRDry)) - ss_SvapPress_nl)
  
  call EndSub("Scalar")

  
contains

  function ss_AdvXZ_ss_fs_sf(ss_Scalar, fs_VelX, sf_VelZ)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: ss_avr_fs, ss_avr_sf
    use differentiate_center4, only: fs_dx_ss, sf_dz_ss
    implicit none
    
    real(8), intent(in) :: fs_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: sf_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_AdvXZ_ss_fs_sf(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_AdvXZ_ss_fs_sf = 0.0d0  !$B=i4|2=(B
    ss_AdvXZ_ss_fs_sf = &
      &   ss_avr_fs(fs_VelX * fs_dx_ss(ss_Scalar)) &
      & + ss_avr_sf(sf_VelZ * sf_dz_ss(ss_Scalar))
    
  end function ss_AdvXZ_ss_fs_sf


  function ss_AdvZ_ss_sf(ss_Scalar, sf_VelZ)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: ss_avr_fs, ss_avr_sf
    use differentiate_center4, only: sf_dz_ss
    implicit none
    
    real(8), intent(in) :: sf_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_AdvZ_ss_sf(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_AdvZ_ss_sf = 0.0d0  !$B=i4|2=(B
    ss_AdvZ_ss_sf = ss_avr_sf(sf_VelZ * sf_dz_ss(ss_Scalar))
    
  end function ss_AdvZ_ss_sf
  
  
  function ss_AdvX_ss_fs(ss_Scalar, fs_VelX)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: ss_avr_fs, ss_avr_fs
    use differentiate_center4, only: fs_dx_ss
    implicit none
    
    real(8), intent(in) :: fs_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_AdvX_ss_fs(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_AdvX_ss_fs = 0.0d0  !$B=i4|2=(B
    ss_AdvX_ss_fs = ss_avr_fs(fs_VelX * fs_dx_ss(ss_Scalar))
    
  end function ss_AdvX_ss_fs
  
  
  function ss_TurbXZ_ss_ss(ss_Scalar, ss_Kh)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: fs_avr_ss, sf_avr_ss
    use differentiate_center4, only: ss_dx_fs, fs_dx_ss, ss_dz_sf, sf_dz_ss
    implicit none
    
    real(8), intent(in) :: ss_Kh(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_TurbXZ_ss_ss(DimXMin:DimXMax, DimZMin:DimZMax)

    ss_TurbXZ_ss_ss = 0.0d0  !$B=i4|2=(B
    ss_TurbXZ_ss_ss = &
      &   ss_dx_fs(fs_avr_ss(ss_Kh) * fs_dx_ss(ss_Scalar)) &
      & + ss_dz_sf(sf_avr_ss(ss_Kh) * sf_dz_ss(ss_Scalar))
    
  end function ss_TurbXZ_ss_ss


  function ss_TurbZ_ss_ss(ss_Scalar, ss_Kh)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: fs_avr_ss, sf_avr_ss
    use differentiate_center4, only: ss_dx_fs, fs_dx_ss, ss_dz_sf, sf_dz_ss
    implicit none

    real(8), intent(in) :: ss_Kh(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_TurbZ_ss_ss(DimXMin:DimXMax, DimZMin:DimZMax)

    ss_TurbZ_ss_ss = 0.0d0  !$B=i4|2=(B
    ss_TurbZ_ss_ss =  ss_dz_sf(sf_avr_ss(ss_Kh) * sf_dz_ss(ss_Scalar))
    
  end function ss_TurbZ_ss_ss


  function ss_NumDiffXZ_ss(ss_Scalar)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use differentiate_center4, only: ss_dx_fs, fs_dx_ss, ss_dz_sf, sf_dz_ss
    use arareset, only: NuH, NuV
    implicit none
    
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_NumDiffXZ_ss(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_NumDiffXZ_ss = 0.0d0  !$B=i4|2=(B
    ss_NumDiffXZ_ss = &
      &   NuH * (ss_dx_fs(fs_dx_ss( ss_Scalar ))) &
      & + NuV * (ss_dz_sf(sf_dz_ss( ss_Scalar ))) 
    
  end function ss_NumDiffXZ_ss


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

    ! $BF~=PNOJQ?t(B
    real(8)             :: s_DSvapDT_s          !$BK0OB>x5$05(B [Pa]
    real(8), intent(in) :: s_temp               !$B29EY(B [K]
    
    ! Antoine $B$N<0$rMxMQ(B. H2O $B$N78?t$O2=3X9)3XJXMw;2>H(B. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    s_DSvapDT_s = dexp( A - B / (C + s_temp) ) * B / ((C + s_temp) ** 2.0d0)
    
  end function s_DSvapDT_s

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

    ! $BF~=PNOJQ?t(B
    real(8)             :: s_svap_s             !$BK0OB>x5$05(B [Pa]
    real(8), intent(in) :: s_temp               !$B29EY(B [K]
    
    ! Antoine $B$N<0$rMxMQ(B. H2O $B$N78?t$O2=3X9)3XJXMw;2>H(B. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    s_svap_s = dexp( A - B / (C + s_temp) ) !Antoine $B$N<0(B
    
  end function s_svap_s

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

    ! $BF~=PNOJQ?t(B
    real(8)             :: ss_svap_ss(DimXMin:DimXMax, DimZMin:DimZMax)
                                                        !$BK0OB>x5$05(B [Pa]
    real(8), intent(in) :: ss_temp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                               !$B29EY(B [K]
    
    ! Antoine $B$N<0$rMxMQ(B. H2O $B$N78?t$O2=3X9)3XJXMw;2>H(B. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    ss_svap_ss = dexp( A - B / (C + ss_temp) ) !Antoine $B$N<0(B
    
  end function ss_svap_ss

  
  function s_LatentHeat_s(temp)
    use physset, only: GasRWet
    implicit none

    real(8), intent(in)    :: temp         !$B29EY(B
    real(8)                :: s_LatentHeat_s   !$B@xG.(B
    real(8)                :: svap_log_dt

    ! Antoine $B$N<0$rMxMQ(B. H2O $B$N78?t$O2=3X9)3XJXMw;2>H(B. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0

    svap_log_dt = B / (C + temp ** 2.0d0)
    s_LatentHeat_s  = svap_log_dt * GasRWet * ( temp ** 2.0d0 )
    
  end function s_LatentHeat_s
  
  function ss_LatentHeat_ss(ss_temp)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use physset, only: GasRWet
    implicit none
    
    real(8), intent(in) :: ss_temp(DimXMin:DimXMax, DimZMin:DimZMax) !$B29EY(B
    real(8)             :: ss_LatentHeat_ss(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: svap_log_dt(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ! Antoine $B$N<0$rMxMQ(B. H2O $B$N78?t$O2=3X9)3XJXMw;2>H(B. 
!    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0

    svap_log_dt = B / (C + ss_temp ** 2.0d0)
    ss_LatentHeat_ss  = svap_log_dt * GasRWet * ( ss_Temp ** 2.0d0 ) 
    
  end function ss_LatentHeat_ss


  
    
  subroutine warmrainprm_init
    !=== $B=PNOJQ?t(B
    ss_PotTemp_al = 0.0d0
    ss_MixRtV_al = 0.0d0
    ss_MixRtC_al = 0.0d0
    ss_MixRtR_al = 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_MixRtSat_nl = 0.0d0

    !=== $BFbItJQ?t(B
    ss_Gamma_nl = 0.0d0
    ss_PressAll_nl = 0.0d0
    ss_ExnerAll_nl = 0.0d0
    ss_TempAll_nl = 0.0d0
    ss_PotTempAll_nl = 0.0d0
    ss_MixRtVAll_nl = 0.0d0
    ss_MixRtVSat_nl = 0.0d0
    ss_PotTempAll_al = 0.0d0
    ss_SvapPress_nl = 0.0d0

    MixRtV_tmp = 0.0d0
    MixRtC_tmp = 0.0d0
    PotTemp_tmp = 0.0d0
    
    PressAll_adj = 0.0d0
    ExnerAll_adj = 0.0d0
    
    TempAll_tmp = 0.0d0
    PotTempAll_tmp = 0.0d0
    MixRtVAll_tmp = 0.0d0
    
    SvapPress_tmp = 0.0d0
    MixRtVSat_tmp = 0.0d0
    Gamma_tmp = 0.0d0
    
    PotTemp_adj = 0.0d0
    MixRtV_adj = 0.0d0
    MixRtC_adj = 0.0d0
    DensGround = 0.0d0

  end subroutine warmrainprm_init
  

end subroutine WarmRainPrm
