scale_atmos_sub_hydrometeor.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_atmos_hydrometeor
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_hydrometeor_setup
  public :: atmos_hydrometeor_regist
  public :: atmos_hydrometeor_lhv
  public :: atmos_hydrometeor_lhs
  public :: atmos_hydrometeor_lhf
  public :: atmos_hydrometeor_entr
  public :: atmos_hydrometeor_diagnose_number_concentration

  interface atmos_hydrometeor_lhv
     module procedure atmos_hydrometeor_lhv_0d
     module procedure atmos_hydrometeor_lhv_1d
     module procedure atmos_hydrometeor_lhv_2d
     module procedure atmos_hydrometeor_lhv_3d
  end interface atmos_hydrometeor_lhv

  interface atmos_hydrometeor_lhs
     module procedure atmos_hydrometeor_lhs_0d
     module procedure atmos_hydrometeor_lhs_1d
     module procedure atmos_hydrometeor_lhs_2d
     module procedure atmos_hydrometeor_lhs_3d
  end interface atmos_hydrometeor_lhs

  interface atmos_hydrometeor_lhf
     module procedure atmos_hydrometeor_lhf_0d
     module procedure atmos_hydrometeor_lhf_1d
     module procedure atmos_hydrometeor_lhf_2d
     module procedure atmos_hydrometeor_lhf_3d
  end interface atmos_hydrometeor_lhf

  interface atmos_hydrometeor_entr
     module procedure atmos_hydrometeor_entr_0d
     module procedure atmos_hydrometeor_entr_2d
     module procedure atmos_hydrometeor_entr_3d
  end interface atmos_hydrometeor_entr

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public            :: i_qv = -1

  integer, public, parameter :: n_hyd = 6

  integer, public, parameter :: i_hc =  1
  integer, public, parameter :: i_hr =  2
  integer, public, parameter :: i_hi =  3
  integer, public, parameter :: i_hs =  4
  integer, public, parameter :: i_hg =  5
  integer, public, parameter :: i_hh =  6

  integer, public            :: i_qc = -1
  integer, public            :: i_qr = -1
  integer, public            :: i_qi = -1
  integer, public            :: i_qs = -1
  integer, public            :: i_qg = -1
  integer, public            :: i_qh = -1

  integer, public            :: i_nc = -1
  integer, public            :: i_nr = -1
  integer, public            :: i_ni = -1
  integer, public            :: i_ns = -1
  integer, public            :: i_ng = -1
  integer, public            :: i_nh = -1

  ! hydrometeor (water + ice)
  integer, public            :: qhs = -1
  integer, public            :: qhe = -2
  ! water
  integer, public            :: qls = -1
  integer, public            :: qle = -2
  ! ice
  integer, public            :: qis = -1
  integer, public            :: qie = -2

  real(RP), public           :: lhv
  real(RP), public           :: lhs
  real(RP), public           :: lhf

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private :: cv_vapor
  real(RP), private :: cp_vapor
  real(RP), private :: cv_water
  real(RP), private :: cp_water
  real(RP), private :: cv_ice
  real(RP), private :: cp_ice

  real(RP), private :: thermodyn_emask = 1.0_rp

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_setup
    use scale_const, only: &
       cpvap          => const_cpvap,          &
       cvvap          => const_cvvap,          &
       cl             => const_cl,             &
       ci             => const_ci,             &
       lhv00          => const_lhv00,          &
       lhs00          => const_lhs00,          &
       lhf00          => const_lhf00,          &
       lhv0           => const_lhv0,           &
       lhs0           => const_lhs0,           &
       lhf0           => const_lhf0,           &
       thermodyn_type => const_thermodyn_type
    use scale_process, only: &
       prc_mpistop
    implicit none
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[HYDEROMETER] / Categ[ATMOS SHARE] / Origin[SCALElib]'

    if ( thermodyn_type == 'EXACT' ) then

       cv_vapor = cvvap
       cp_vapor = cpvap
       cv_water = cl
       cp_water = cv_water
       cv_ice   = ci
       cp_ice   = cv_ice

       lhv = lhv00
       lhs = lhs00
       lhf = lhf00
       thermodyn_emask = 1.0_rp

    elseif( thermodyn_type == 'SIMPLE' ) then

       cv_vapor = cvvap
       cp_vapor = cpvap
       cv_water = cvvap
       cp_water = cv_water
       cv_ice   = cvvap
       cp_ice   = cv_ice

       lhv = lhv0
       lhs = lhs0
       lhf = lhf0
       thermodyn_emask = 0.0_rp

    else
       write(*,*) 'xxx Not appropriate ATMOS_THERMODYN_ENERGY_TYPE. Check!', trim(thermodyn_type)
       call prc_mpistop
    endif

    return
  end subroutine atmos_hydrometeor_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_regist( &
       Q0,   &
       NV,   &
       NL,   &
       NI,   &
       NAME, &
       DESC, &
       UNIT, &
       ADVC  )
    use scale_process, only: &
      prc_mpistop
    use scale_tracer, only: &
      tracer_regist
    use scale_const, only: &
      rvap => const_rvap
    implicit none

    integer,          intent(out) :: q0
    integer,          intent(in)  :: nv
    integer,          intent(in)  :: nl
    integer,          intent(in)  :: ni
    character(len=*), intent(in)  :: name(nv+nl+ni)
    character(len=*), intent(in)  :: desc(nv+nl+ni)
    character(len=*), intent(in)  :: unit(nv+nl+ni)

    logical,          intent(in), optional :: advc(nv+nl+ni)

    real(RP) :: cv   (nv+nl+ni)
    real(RP) :: cp   (nv+nl+ni)
    real(RP) :: r    (nv+nl+ni)
    logical  :: mass (nv+nl+ni)
    logical  :: advc_(nv+nl+ni)

    integer  :: nq
    integer  :: n
    !---------------------------------------------------------------------------

    if ( i_qv > 0 ) then
       write(*,*) 'xxx tracer for hydrometeor is already registerd'
       call prc_mpistop
    endif

    if ( nv /= 1 ) then
       write(*,*) 'xxx number of vapor must be 1 at this moment'
       call prc_mpistop
    endif

    nq = 0

    do n = 1, nv
       nq = nq + 1
       cv(nq) = cv_vapor
       cp(nq) = cp_vapor
       r(nq) = rvap
    end do

    do n = 1, nl
       nq = nq + 1
       cv(nq) = cv_water
       cp(nq) = cp_water
       r(nq) = 0.0_rp
    end do

    do n = 1, ni
       nq = nq + 1
       cv(nq) = cv_ice
       cp(nq) = cp_ice
       r(nq) = 0.0_rp
    end do

    ! NQ = 1 + NL + NI,   vapor + liqid + ice

    if ( present(advc) ) then
       advc_(:) = advc(:)
    else
       advc_(:) = .true.
    endif

    do n = 1, nq
       mass(n) = .true.
    end do

    call tracer_regist( q0,         & ! [OUT]
                        nq,         & ! [IN]
                        name,       & ! [IN]
                        desc,       & ! [IN]
                        unit,       & ! [IN]
                        cv, cp, r,  & ! [IN], optional
                        advc_, mass ) ! [IN], optional

    i_qv = q0

    if ( nq > 1 ) then
       qhs = i_qv + 1
       qhe = qhs + nl + ni - 1
    endif

    if ( nl > 0 ) then
       qls = i_qv + 1
       qle = qls + nl - 1
    endif

    if ( ni > 0 ) then
       qis = qle + 1
       qie = qis + ni - 1
    endif

    return
  end subroutine atmos_hydrometeor_regist

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_0d( &
       lhv, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhv0  => const_lhv0
    implicit none

    real(RP), intent(out) :: lhv
    real(RP), intent(in)  :: temp
    !---------------------------------------------------------------------------

    lhv = lhv0 + ( cp_vapor - cp_water ) * ( temp - tem00 ) * thermodyn_emask

    return
  end subroutine atmos_hydrometeor_lhv_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_1d( &
       lhv, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhv0  => const_lhv0
    implicit none

    real(RP), intent(out) :: lhv (KA)
    real(RP), intent(in)  :: temp(KA)

    integer :: k
    !---------------------------------------------------------------------------

    do k = ks, ke
       lhv(k) = lhv0 + ( cp_vapor - cp_water ) * ( temp(k) - tem00 ) * thermodyn_emask
    enddo

    return
  end subroutine atmos_hydrometeor_lhv_1d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_2d( &
       lhv, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhv0  => const_lhv0
    implicit none

    real(RP), intent(out) :: lhv (IA,JA)
    real(RP), intent(in)  :: temp(IA,JA)

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
       lhv(i,j) = lhv0 + ( cp_vapor - cp_water ) * ( temp(i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_lhv_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhv_3d( &
       lhv, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhv0  => const_lhv0
    implicit none

    real(RP), intent(out) :: lhv (KA,IA,JA)
    real(RP), intent(in)  :: temp(KA,IA,JA)

    integer :: k, i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhv(k,i,j) = lhv0 + ( cp_vapor - cp_water ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_lhv_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_0d( &
       lhs, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhs0  => const_lhs0
    implicit none

    real(RP), intent(out) :: lhs
    real(RP), intent(in)  :: temp
    !---------------------------------------------------------------------------

    lhs = lhs0 + ( cp_vapor - cp_ice ) * ( temp - tem00 ) * thermodyn_emask

    return
  end subroutine atmos_hydrometeor_lhs_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_1d( &
       lhs, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhs0  => const_lhs0
    implicit none

    real(RP), intent(out) :: lhs (KA)
    real(RP), intent(in)  :: temp(KA)

    integer :: k
    !---------------------------------------------------------------------------

    do k = ks, ke
       lhs(k) = lhs0 + ( cp_vapor - cp_ice ) * ( temp(k) - tem00 ) * thermodyn_emask
    enddo

    return
  end subroutine atmos_hydrometeor_lhs_1d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_2d( &
       lhs, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhs0  => const_lhs0
    implicit none

    real(RP), intent(out) :: lhs (IA,JA)
    real(RP), intent(in)  :: temp(IA,JA)

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
       lhs(i,j) = lhs0 + ( cp_vapor - cp_ice ) * ( temp(i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_lhs_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhs_3d( &
       lhs, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhs0  => const_lhs0
    implicit none

    real(RP), intent(out) :: lhs (KA,IA,JA)
    real(RP), intent(in)  :: temp(KA,IA,JA)

    integer :: k, i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhs(k,i,j) = lhs0 + ( cp_vapor - cp_ice ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_lhs_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_0d( &
       lhf, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhf0  => const_lhf0
    implicit none

    real(RP), intent(out) :: lhf
    real(RP), intent(in)  :: temp
    !---------------------------------------------------------------------------

    lhf = lhf0 + ( cp_water - cp_ice ) * ( temp - tem00 ) * thermodyn_emask

    return
  end subroutine atmos_hydrometeor_lhf_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_1d( &
       lhf, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhf0  => const_lhf0
    implicit none

    real(RP), intent(out) :: lhf (KA)
    real(RP), intent(in)  :: temp(KA)

    integer :: k
    !---------------------------------------------------------------------------

    do k = ks, ke
       lhf(k) = lhf0 + ( cp_water - cp_ice ) * ( temp(k) - tem00 ) * thermodyn_emask
    enddo

    return
  end subroutine atmos_hydrometeor_lhf_1d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_2d( &
       lhf, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhf0  => const_lhf0
    implicit none

    real(RP), intent(out) :: lhf (IA,JA)
    real(RP), intent(in)  :: temp(IA,JA)

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
       lhf(i,j) = lhf0 + ( cp_water - cp_ice ) * ( temp(i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_lhf_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_lhf_3d( &
       lhf, &
       temp )
    use scale_const, only: &
       tem00 => const_tem00, &
       lhf0  => const_lhf0
    implicit none

    real(RP), intent(out) :: lhf (KA,IA,JA)
    real(RP), intent(in)  :: temp(KA,IA,JA)

    integer :: k, i, j
    !---------------------------------------------------------------------------

    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke
       lhf(k,i,j) = lhf0 + ( cp_water - cp_ice ) * ( temp(k,i,j) - tem00 ) * thermodyn_emask
    enddo
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_lhf_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr_0d( &
       entr, &
       temp, &
       pres, &
       q,    &
       Rq    )
    use scale_const, only: &
       eps   => const_eps,   &
       pre00 => const_pre00, &
       tem00 => const_tem00, &
       rdry  => const_rdry,  &
       cpdry => const_cpdry, &
       rvap  => const_rvap,  &
       lhv0  => const_lhv0,  &
       lhf0  => const_lhf0,  &
       psat0 => const_psat0
    use scale_tracer, only: &
       qa
    implicit none

    real(RP), intent(out) :: entr
    real(RP), intent(in)  :: temp
    real(RP), intent(in)  :: pres
    real(RP), intent(in)  :: q(QA)
    real(RP), intent(in)  :: Rq(QA)

    real(RP) :: qdry, Rtot
    real(RP) :: logT_T0, pres_dry, pres_vap

    integer  :: iqw
    !---------------------------------------------------------------------------

    logt_t0 = log( temp / tem00 )

    qdry = 1.0_rp
    rtot = 0.0_rp
    do iqw = 1, qa
       qdry = qdry - q(iqw)
       rtot = rtot + q(iqw) * rq(iqw)
    enddo
    rtot = rtot + rdry * qdry

    ! dry air + vapor
    pres_dry = max( pres * qdry * rdry / rtot, eps )
    entr = qdry * cpdry * logt_t0 &
         - qdry * rdry  * log( pres_dry / pre00 )

    if ( i_qv > 0 ) then
       pres_vap = max( pres * q(i_qv) * rvap / rtot, eps )
       entr = entr + q(i_qv) * cp_vapor * logt_t0 &
                   - q(i_qv) * rvap  * log( pres_vap / psat0 ) &
                   + q(i_qv) * lhv0 / tem00
    endif

    ! liquid water
    if ( qls > 0 ) then
       do iqw = qls, qle
          entr = entr + q(iqw) * cp_water * logt_t0
       enddo
    endif

    ! ice
    if ( qis > 0 ) then
       do iqw = qis, qie
          entr = entr + q(iqw) * cp_ice * logt_t0 &
                      - q(iqw) * lhf0 / tem00
       enddo
    endif

    return
  end subroutine atmos_hydrometeor_entr_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr_2d( &
       entr, &
       temp, &
       pres, &
       q,    &
       Rq    )
    use scale_const, only: &
       eps   => const_eps,   &
       pre00 => const_pre00, &
       tem00 => const_tem00, &
       cpdry => const_cpdry, &
       rdry  => const_rdry,  &
       rvap  => const_rvap,  &
       lhv0  => const_lhv0,  &
       lhf0  => const_lhf0,  &
       psat0 => const_psat0
    use scale_tracer, only: &
       qa
    implicit none

    real(RP), intent(out) :: entr(IA,JA)
    real(RP), intent(in)  :: temp(IA,JA)
    real(RP), intent(in)  :: pres(IA,JA)
    real(RP), intent(in)  :: q   (IA,JA,QA)
    real(RP), intent(in)  :: Rq  (QA)

    real(RP) :: qdry, Rtot
    real(RP) :: logT_T0, pres_dry, pres_vap

    integer  :: i, j, iqw
    !---------------------------------------------------------------------------

    ! dry air + vapor
    do j = jsb, jeb
    do i = isb, ieb

       logt_t0 = log( temp(i,j) / tem00 )

       qdry = 1.0_rp
       rtot = 0.0_rp
       do iqw = 1, qa
          qdry = qdry - q(i,j,iqw)
          rtot = rtot + q(i,j,iqw) * rq(iqw)
       enddo
       rtot = rtot + rdry * qdry

       ! dry air + vapor
       pres_dry = max( pres(i,j) * qdry    * rdry / rtot, eps )
       entr(i,j) = qdry    * cpdry * logt_t0 &
                 - qdry    * rdry  * log( pres_dry / pre00 )

       if ( i_qv > 0 ) then
          pres_vap = max( pres(i,j) * q(i,j,i_qv) * rvap / rtot, eps )
          entr(i,j) = entr(i,j) + q(i,j,i_qv) * cp_vapor * logt_t0 &
                                - q(i,j,i_qv) * rvap  * log( pres_vap / psat0 ) &
                                + q(i,j,i_qv) * lhv0 / tem00
       endif

       ! liquid water
       if ( qls > 0 ) then
          do iqw = qls, qle
             entr(i,j) = entr(i,j) + q(i,j,iqw) * cp_water * logt_t0
          enddo
       endif

       ! ice
       if ( qis > 0 ) then
          do iqw = qis, qie
             entr(i,j) = entr(i,j) + q(i,j,iqw) * cp_ice * logt_t0 &
                                   - q(i,j,iqw) * lhf0 / tem00
          enddo
       endif

    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_entr_2d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_entr_3d( &
       entr, &
       temp, &
       pres, &
       q,    &
       Rq    )
    use scale_const, only: &
       eps   => const_eps,   &
       pre00 => const_pre00, &
       tem00 => const_tem00, &
       cpdry => const_cpdry, &
       rdry  => const_rdry,  &
       rvap  => const_rvap,  &
       lhv0  => const_lhv0,  &
       lhf0  => const_lhf0,  &
       psat0 => const_psat0
    use scale_tracer, only: &
       qa
    implicit none

    real(RP), intent(out) :: entr(KA,IA,JA)
    real(RP), intent(in)  :: temp(KA,IA,JA)
    real(RP), intent(in)  :: pres(KA,IA,JA)
    real(RP), intent(in)  :: q   (KA,IA,JA,QA)
    real(RP), intent(in)  :: Rq  (QA)

    real(RP) :: qdry, Rtot
    real(RP) :: logT_T0, pres_dry, pres_vap

    integer  :: k, i, j, iqw
    !---------------------------------------------------------------------------

    ! dry air + vapor
    do j = jsb, jeb
    do i = isb, ieb
    do k = ks, ke

       logt_t0 = log( temp(k,i,j) / tem00 )

       qdry = 1.0_rp
       rtot = 0.0_rp
       do iqw = 1, qa
          qdry = qdry - q(k,i,j,iqw)
          rtot = rtot + q(k,i,j,iqw) * rq(iqw)
       enddo
       rtot = rtot + rdry * qdry

       ! dry air + vapor
       pres_dry = max( pres(k,i,j) * qdry    * rdry / rtot, eps )
       entr(k,i,j) = qdry    * cpdry * logt_t0 &
                   - qdry    * rdry  * log( pres_dry / pre00 )

       if ( i_qv > 0 ) then
          pres_vap = max( pres(k,i,j) * q(k,i,j,i_qv) * rvap / rtot, eps )
          entr(k,i,j) = entr(k,i,j) + q(k,i,j,i_qv) * cp_vapor * logt_t0 &
                                    - q(k,i,j,i_qv) * rvap  * log( pres_vap / psat0 ) &
                                    + q(k,i,j,i_qv) * lhv0 / tem00
       endif

       ! liquid water
       if ( qls > 0 ) then
          do iqw = qls, qle
             entr(k,i,j) = entr(k,i,j) + q(k,i,j,iqw) * cp_water * logt_t0
          enddo
       endif

       ! ice
       if ( qis > 0 ) then
          do iqw = qis, qie
             entr(k,i,j) = entr(k,i,j) + q(k,i,j,iqw) * cp_ice * logt_t0 &
                                       - q(k,i,j,iqw) * lhf0 / tem00
          enddo
       endif

    enddo
    enddo
    enddo

    return
  end subroutine atmos_hydrometeor_entr_3d

  !-----------------------------------------------------------------------------
  subroutine atmos_hydrometeor_diagnose_number_concentration( &
       QTRC )
    use scale_const, only: &
       pi => const_pi
    implicit none

    real(RP), intent(inout) :: qtrc(:,:,:,:)

    real(RP), parameter :: dc   =  20.e-6_rp ! typical particle diameter for cloud  [m]
    real(RP), parameter :: dr   = 200.e-6_rp ! typical particle diameter for rain   [m]
    real(RP), parameter :: di   =  80.e-6_rp ! typical particle diameter for ice    [m]
    real(RP), parameter :: ds   =  80.e-6_rp ! typical particle diameter for snow   [m]
    real(RP), parameter :: dg   = 200.e-6_rp ! typical particle diameter for grapel [m]
    real(RP), parameter :: rhow =  1000.0_rp ! typical density for warm particles   [kg/m3]
    real(RP), parameter :: rhof =   100.0_rp ! typical density for frozen particles [kg/m3]
    real(RP), parameter :: rhog =   400.0_rp ! typical density for grapel particles [kg/m3]
    real(RP), parameter :: b    =     3.0_rp ! assume spherical form

    real(RP) :: piov6
    !---------------------------------------------------------------------------

    piov6 = pi / 6.0_rp

    if ( i_nc > 0 ) qtrc(:,:,:,i_nc) = qtrc(:,:,:,i_qc) / ( (piov6*rhow) * dc**b )
    if ( i_nr > 0 ) qtrc(:,:,:,i_nr) = qtrc(:,:,:,i_qr) / ( (piov6*rhow) * dr**b )
    if ( i_ni > 0 ) qtrc(:,:,:,i_ni) = qtrc(:,:,:,i_qi) / ( (piov6*rhof) * di**b )
    if ( i_ns > 0 ) qtrc(:,:,:,i_ns) = qtrc(:,:,:,i_qs) / ( (piov6*rhof) * ds**b )
    if ( i_ng > 0 ) qtrc(:,:,:,i_ng) = qtrc(:,:,:,i_qg) / ( (piov6*rhog) * dg**b )

    return
  end subroutine atmos_hydrometeor_diagnose_number_concentration

end module scale_atmos_hydrometeor