!c Description:
!c  雲密度の計算. 単位は [g/l] (昔の計算に合わせるため)
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-06-02    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_clouddensity
  interface

     subroutine CloudDensity(temp, pres, property, molfr, cldns)
       use prm_type, only: type_property
       real(8), intent(in)                :: temp(:)    
       real(8), intent(in)                :: pres(:)
       real(8), intent(in)                :: molfr(:,:)
       real(8), intent(out), allocatable  :: cldns(:,:)
       type(type_property), intent(in)    :: property(:)
     end subroutine CloudDensity

  end interface
end module interface_clouddensity



subroutine CloudDensity(temp, pres, property, molfr, cldns)
  use prm_phys, only: gasr
  use prm_type, only: type_property

  implicit none
  real(8), intent(in)               :: temp(:)    
  real(8), intent(in)               :: pres(:)
  real(8), intent(in)               :: molfr(:,:)
  real(8), intent(out), allocatable :: cldns(:,:)
  type(type_property), intent(in)   :: property(:)

  real(8), allocatable :: dx(:,:) 
  real(8), allocatable :: dp(:)
  integer              :: step
  integer              :: smax
  integer              :: smax_gas
  integer              :: snum
  integer              :: num
  
  
  !配列のサイズの取得
  step = size(molfr, 1)
  smax = size(molfr, 2)
  smax_gas = count(property%phase == 1)
  
  !配列の定義
  allocate(cldns(step, smax), dp(step), dx(step, smax))
  
  !モル比の変化を計算
  do num = 2, step - 1
     dx(num,:) = (molfr(num+1,:) - molfr(num-1,:)) / 2.0D0
  end do
  dx(1,:) = dx(2,:)
  dx(step,:) = dx(step - 1,:)
  
  !圧力差の計算
  do num = 2, step - 1
     dp(num) = (pres(num + 1) - pres(num - 1)) / 2.0d0
  end do
  dp(1) = dp(2)
  dp(step) = dp(step - 1)
  
  !雲量計算
  do snum = 1, smax_gas  
     cldns(:,snum) = abs((property(snum)%molwt * dx(:, snum) * pres**2.0D0) &
          & / (gasr * temp * dp)) * 1.0D-3
  end do
  
  !配列の解放
  deallocate(dp, dx)
  
end subroutine CloudDensity

  
