program main

! 凝結量のチェック. 
! physics/masscondense_threshold-masstest2.f90 のテスト. 
!
! 木星計算に倣い, 雲密度計算を移流拡散部と凝結部に分けた. 
! masscondense_threshold-masstest2.f90 は雲密度計算の凝結部を計算する
! プログラムである. 
! 本プログラムは masscondense_threshold-masstest2.f90 を換骨奪胎したもので, 
! 与えられた飽和比と雲密度に対して, 凝結量が期待通りに出力されるかを確かめる
! ことを目的として作ったものである. 
! 雲密度に閾値を導入した場合, 9 通りに場合分けされるので, それぞれの場合に
! ついての値を出力する. 

  implicit none ! 暗黙の型宣言禁止

  integer i             ! ループ変数
  real(8) DelTime       ! 時間ステップ
  real(8) SatRatioCr    ! 臨界飽和比
  real(8) DensCloudThreshold ! 雲密度の閾値
  real(8) SatRatio(9) ! 飽和比
  real(8) DensCloud(9) ! 雲密度
  real(8) MassCond(9) ! 凝結量

  DelTime = 1.0d0
  SatRatioCr = 1.2d0
  DensCloudThreshold = 0.1d0

  SatRatio(1) = 0.9d0
!  SatRatio(1) = 0.1d0
  SatRatio(2) = 1.1d0
  SatRatio(3) = 1.3d0
  SatRatio(4) = 0.9d0
!  SatRatio(4) = 0.99d0
  SatRatio(5) = 1.1d0
  SatRatio(6) = 1.3d0
  SatRatio(7) = 0.9d0
  SatRatio(8) = 1.1d0
  SatRatio(9) = 1.3d0

  DensCloud(1) = 0.5d0
  DensCloud(2) = 0.5d0
  DensCloud(3) = 0.5d0
  DensCloud(4) = 5.0d-2
  DensCloud(5) = 5.0d-2
  DensCloud(6) = 5.0d-2
  DensCloud(7) = 0.0d0
  DensCloud(8) = 0.0d0
  DensCloud(9) = 0.0d0


do i = 1, 9
   if (SatRatio(i) - SatRatioCr > epsilon(0.0d0)) then

! 飽和比が定めた臨界飽和比よりも大きい場合

      MassCond(i) = &
           & max( - DensCloud(i) / DelTime, &
           &   SatRatio(i) - 1.0d0 )

   else if (DensCloud(i) > DensCloudThreshold ) then

! 臨界飽和比を超えていないが, 雲密度が閾値以上である場合
! 飽和比が 1 以上のときに凝結, 1 以下のときに蒸発.

      MassCond(i) = &
           & max( - DensCloud(i) / DelTime, &
           &   SatRatio(i) - 1.0d0 )
   else

! 臨界飽和比を超えておらず, 雲密度が閾値未満の場合
! 雲密度がゼロではなく, 飽和比 1.0 未満のときに限って蒸発.
! それ以外の場合には何も起きない.

      MassCond(i) = &
           & max( - DensCloud(i) / DelTime, &
           & ( 0.5d0 * &
           &   ( 1.0d0 + sign(1.0d0, DensCloud(i) - epsilon(0.0d0)) ) &
           & ) * &
           & ( 0.5d0 * &
           &   ( 1.0d0 - sign(1.0d0, SatRatio(i) - 1.0d0) ) ) &
           & * (SatRatio(i) - 1.0d0) )
   end if
end do

  write(*,*) "***** input *****"
  write(*,*) "DelTime=", DelTime
  write(*,*) "rhos(threshold)=", DensCloudThreshold
  write(*,*) "Scr=", SatRatioCr
  write(*,*) "***** output *****"
  write(*,*) "Mcond(S=0.9,rhos=0.5)", MassCond(1)
  write(*,*) "Mcond(S=1.1,rhos=0.5)", MassCond(2)
  write(*,*) "Mcond(S=1.3,rhos=0.5)", MassCond(3)
  write(*,*) "Mcond(S=0.9,rhos=0.05)", MassCond(4)
  write(*,*) "Mcond(S=1.1,rhos=0.05)", MassCond(5)
  write(*,*) "Mcond(S=1.3,rhos=0.05)", MassCond(6)
  write(*,*) "Mcond(S=0.9,rhos=0.0)", MassCond(7)
  write(*,*) "Mcond(S=1.1,rhos=0.0)", MassCond(8)
  write(*,*) "Mcond(S=1.3,rhos=0.0)", MassCond(9)

end program main