mod_atmos_phy_mp_driver.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_atmos_phy_mp_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_mp_driver_tracer_setup
  public :: atmos_phy_mp_driver_setup
  public :: atmos_phy_mp_driver_calc_tendency
  public :: atmos_phy_mp_driver_adjustment
  public :: atmos_phy_mp_driver_qhyd2qtrc

  interface abstract
     subroutine qhyd2qtrc( &
          KA, KS, KE, IA, IS, IE, JA, JS, JE, &
          QV, QHYD, &
          QTRC,     &
          QNUM      )
       use scale_precision
       use scale_atmos_hydrometeor, only: n_hyd
       use mod_atmos_phy_mp_vars, only: qa_mp
       integer, intent(in) :: ka, ks, ke
       integer, intent(in) :: ia, is, ie
       integer, intent(in) :: ja, js, je
       real(RP), intent(in) :: qv(ka,ia,ja)
       real(RP), intent(in) :: qhyd(ka,ia,ja,n_hyd)
       real(RP), intent(out) :: qtrc(ka,ia,ja,qa_mp)
       real(RP), intent(in), optional :: qnum(ka,ia,ja,n_hyd)
     end subroutine qhyd2qtrc
  end interface abstract
  procedure(qhyd2qtrc), pointer :: atmos_phy_mp_user_qhyd2qtrc => null()
  public :: atmos_phy_mp_user_qhyd2qtrc

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical,  private :: mp_do_precipitation   = .true.
  logical,  private :: mp_do_negative_fixer  = .true.
  real(RP), private :: mp_limit_negative     = 0.1_rp
  integer,  private :: mp_ntmax_sedimentation = 1
  real(RP), private :: mp_max_term_vel = 10.0_rp
  real(RP), private :: mp_cldfrac_thleshold
  integer,  private :: mp_nstep_sedimentation
  real(RP), private :: mp_rnstep_sedimentation
  real(DP), private :: mp_dtsec_sedimentation

  integer, private, allocatable :: hist_vterm_id(:)
  integer, private              :: hist_nf_rhoh_id
  integer, private              :: hist_nf_dens_id
  integer, private              :: hist_nf_engi_id
  integer, private              :: monit_nf_mass_id
  integer, private              :: monit_nf_engi_id
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_driver_tracer_setup
    use scale_prc, only: &
       prc_abort
    use scale_tracer, only: &
       tracer_regist
    use scale_atmos_hydrometeor, only: &
       atmos_hydrometeor_regist
    use scale_atmos_phy_mp_kessler, only: &
       atmos_phy_mp_kessler_ntracers,            &
       atmos_phy_mp_kessler_nwaters,             &
       atmos_phy_mp_kessler_nices,               &
       atmos_phy_mp_kessler_tracer_names,        &
       atmos_phy_mp_kessler_tracer_descriptions, &
       atmos_phy_mp_kessler_tracer_units
    use scale_atmos_phy_mp_tomita08, only: &
       atmos_phy_mp_tomita08_ntracers,            &
       atmos_phy_mp_tomita08_nwaters,             &
       atmos_phy_mp_tomita08_nices,               &
       atmos_phy_mp_tomita08_tracer_names,        &
       atmos_phy_mp_tomita08_tracer_descriptions, &
       atmos_phy_mp_tomita08_tracer_units
    use scale_atmos_phy_mp_sn14, only: &
       atmos_phy_mp_sn14_ntracers,            &
       atmos_phy_mp_sn14_nwaters,             &
       atmos_phy_mp_sn14_nices,               &
       atmos_phy_mp_sn14_tracer_names,        &
       atmos_phy_mp_sn14_tracer_descriptions, &
       atmos_phy_mp_sn14_tracer_units
    use scale_atmos_phy_mp_suzuki10, only: &
       atmos_phy_mp_suzuki10_tracer_setup,        &
       atmos_phy_mp_suzuki10_ntracers,            &
       atmos_phy_mp_suzuki10_nwaters,             &
       atmos_phy_mp_suzuki10_nices,               &
       atmos_phy_mp_suzuki10_nccn,                &
       atmos_phy_mp_suzuki10_tracer_names,        &
       atmos_phy_mp_suzuki10_tracer_descriptions, &
       atmos_phy_mp_suzuki10_tracer_units
    use mod_atmos_admin, only: &
       atmos_phy_mp_type, &
       atmos_sw_phy_mp
    use mod_atmos_phy_mp_vars, only: &
       qa_mp, &
       qs_mp, &
       qe_mp
    implicit none

    integer :: QS2
    !---------------------------------------------------------------------------

    log_newline
    log_info("ATMOS_PHY_MP_driver_tracer_setup",*) 'Setup'

    if ( atmos_sw_phy_mp ) then
       select case ( atmos_phy_mp_type )
       case ( 'KESSLER' )
          call atmos_hydrometeor_regist( &
               atmos_phy_mp_kessler_nwaters,                & ! [IN]
               atmos_phy_mp_kessler_nices,                  & ! [IN]
               atmos_phy_mp_kessler_tracer_names(:),        & ! [IN]
               atmos_phy_mp_kessler_tracer_descriptions(:), & ! [IN]
               atmos_phy_mp_kessler_tracer_units(:),        & ! [IN]
               qs_mp                                        ) ! [OUT]
          qa_mp = atmos_phy_mp_kessler_ntracers
       case ( 'TOMITA08' )
          call atmos_hydrometeor_regist( &
               atmos_phy_mp_tomita08_nwaters,                & ! [IN]
               atmos_phy_mp_tomita08_nices,                  & ! [IN]
               atmos_phy_mp_tomita08_tracer_names(:),        & ! [IN]
               atmos_phy_mp_tomita08_tracer_descriptions(:), & ! [IN]
               atmos_phy_mp_tomita08_tracer_units(:),        & ! [IN]
               qs_mp                                         ) ! [OUT]
          qa_mp = atmos_phy_mp_tomita08_ntracers
       case( 'SN14' )
          call atmos_hydrometeor_regist( &
               atmos_phy_mp_sn14_nwaters,                   & ! [IN]
               atmos_phy_mp_sn14_nices,                     & ! [IN]
               atmos_phy_mp_sn14_tracer_names(1:6),         & ! [IN]
               atmos_phy_mp_sn14_tracer_descriptions(1:6),  & ! [IN]
               atmos_phy_mp_sn14_tracer_units(1:6),         & ! [IN]
               qs_mp                                        ) ! [OUT]
          call tracer_regist( qs2,                          & ! [OUT]
               5,                                           & ! [IN]
               atmos_phy_mp_sn14_tracer_names(7:11),        & ! [IN]
               atmos_phy_mp_sn14_tracer_descriptions(7:11), & ! [IN]
               atmos_phy_mp_sn14_tracer_units(7:11)         ) ! [IN]
          qa_mp = atmos_phy_mp_sn14_ntracers
       case ( 'SUZUKI10' )
          call atmos_phy_mp_suzuki10_tracer_setup
          call atmos_hydrometeor_regist( &
               atmos_phy_mp_suzuki10_nwaters,                & ! [IN]
               atmos_phy_mp_suzuki10_nices,                  & ! [IN]
               atmos_phy_mp_suzuki10_tracer_names(:),        & ! [IN]
               atmos_phy_mp_suzuki10_tracer_descriptions(:), & ! [IN]
               atmos_phy_mp_suzuki10_tracer_units(:),        & ! [IN]
               qs_mp                                         ) ! [OUT]
          if( atmos_phy_mp_suzuki10_nccn > 0 ) then
             call tracer_regist( qs2,                                                                     & ! [OUT]
                                 atmos_phy_mp_suzuki10_nccn,                                              & ! [IN]
                                 atmos_phy_mp_suzuki10_tracer_names( atmos_phy_mp_suzuki10_nwaters &
                                                                          + atmos_phy_mp_suzuki10_nices   &
                                                                          + 2 : ),                        & ! [IN]
                                 atmos_phy_mp_suzuki10_tracer_descriptions( atmos_phy_mp_suzuki10_nwaters &
                                                                          + atmos_phy_mp_suzuki10_nices   &
                                                                          + 2 : ),                        & ! [IN]
                                 atmos_phy_mp_suzuki10_tracer_units( atmos_phy_mp_suzuki10_nwaters &
                                                                          + atmos_phy_mp_suzuki10_nices   &
                                                                          + 2 : )                         ) ! [IN]
          end if
          qa_mp = atmos_phy_mp_suzuki10_ntracers
       case default
          log_error("ATMOS_PHY_MP_driver_tracer_setup",*) 'ATMOS_PHY_MP_TYPE is invalud: ', trim(atmos_phy_mp_type)
          call prc_abort
       end select

       qe_mp = qs_mp + qa_mp - 1

    else
       qa_mp = 0
       qs_mp = -1
       qe_mp = -2
    end if

    return
  end subroutine atmos_phy_mp_driver_tracer_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_driver_setup
    use scale_prc, only: &
       prc_abort
    use scale_atmos_grid_cartesc, only: &
       cdz => atmos_grid_cartesc_cdz
    use scale_const, only: &
       eps => const_eps
    use scale_time, only: &
       time_dtsec_atmos_phy_mp
    use scale_atmos_phy_mp_kessler, only: &
       atmos_phy_mp_kessler_setup
    use scale_atmos_phy_mp_tomita08, only: &
       atmos_phy_mp_tomita08_setup
    use scale_atmos_phy_mp_sn14, only: &
       atmos_phy_mp_sn14_setup
    use scale_atmos_phy_mp_suzuki10, only: &
       atmos_phy_mp_suzuki10_setup
    use scale_file_history, only: &
       file_history_reg
    use scale_monitor, only: &
       monitor_reg, &
       monitor_put
    use mod_atmos_admin, only: &
       atmos_phy_mp_type, &
       atmos_sw_phy_mp
    use mod_atmos_phy_mp_vars, only: &
       atmos_phy_mp_cldfrac_thleshold, &
       sflx_rain => atmos_phy_mp_sflx_rain, &
       sflx_snow => atmos_phy_mp_sflx_snow
    use mod_atmos_phy_mp_vars, only: &
       qs_mp, &
       qe_mp
    implicit none

    namelist / param_atmos_phy_mp / &
       mp_do_precipitation,     &
       mp_do_negative_fixer,    &
       mp_limit_negative,       &
       mp_ntmax_sedimentation,  &
       mp_max_term_vel,         &
       mp_cldfrac_thleshold

    real(RP) :: ZERO(ka,ia,ja)

    integer :: nstep_max

    integer :: iq
    integer :: ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("ATMOS_PHY_MP_driver_setup",*) 'Setup'

    if ( atmos_sw_phy_mp ) then

       mp_cldfrac_thleshold = eps

       !--- read namelist
       rewind(io_fid_conf)
       read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
       if( ierr < 0 ) then !--- missing
          log_info("ATMOS_PHY_MP_driver_setup",*) 'Not found namelist. Default used.'
       elseif( ierr > 0 ) then !--- fatal error
          log_error("ATMOS_PHY_MP_driver_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
          call prc_abort
       endif
       log_nml(param_atmos_phy_mp)

       nstep_max = ceiling( ( time_dtsec_atmos_phy_mp * mp_max_term_vel ) / minval( cdz(:) ) )
       mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )

       mp_nstep_sedimentation  = mp_ntmax_sedimentation
       mp_rnstep_sedimentation = 1.0_rp / real(mp_ntmax_sedimentation,kind=rp)
       mp_dtsec_sedimentation  = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation

       atmos_phy_mp_cldfrac_thleshold = mp_cldfrac_thleshold

       log_newline
       log_info("ATMOS_PHY_MP_driver_setup",*) 'Enable negative fixer?                    : ', mp_do_negative_fixer
       log_info("ATMOS_PHY_MP_driver_setup",*) 'Value limit of negative fixer (abs)       : ', abs(mp_limit_negative)
       log_info("ATMOS_PHY_MP_driver_setup",*) 'Enable sedimentation (precipitation)?     : ', mp_do_precipitation
       log_info("ATMOS_PHY_MP_driver_setup",*) 'Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, 'step'
       log_info("ATMOS_PHY_MP_driver_setup",*) 'DT of sedimentation                       : ', mp_dtsec_sedimentation, '[s]'

       select case ( atmos_phy_mp_type )
       case ( 'KESSLER' )
          call atmos_phy_mp_kessler_setup
       case ( 'TOMITA08' )
          call atmos_phy_mp_tomita08_setup( &
               ka, ks, ke, ia, is, ie, ja, js, je )
       case ( 'SN14' )
          call atmos_phy_mp_sn14_setup( &
               ka, ia, ja )
       case ( 'SUZUKI10' )
          call atmos_phy_mp_suzuki10_setup( &
               ka, ia, ja )
       end select

       ! history putput
       if ( mp_do_precipitation ) then
          allocate( hist_vterm_id(qs_mp+1:qe_mp) )
          do iq = qs_mp+1, qe_mp
             call file_history_reg( 'Vterm_'//trim(tracer_name(iq)), 'terminal velocity of '//trim(tracer_name(iq)), 'm/s', hist_vterm_id(iq) )
          end do
       end if

       ! monitor
       if ( mp_do_negative_fixer ) then
          call file_history_reg( "RHOH_MP_NF",   "sensible heat by the negative fixer",          "J/m3/s", & ! [IN]
                                 hist_nf_rhoh_id                                                            ) ! [OUT]
          call file_history_reg( "DENS_t_MP_NF", "vapor supply by the negative fixer",           "kg/m3/s", & ! [IN]
                                 hist_nf_dens_id                                                            ) ! [OUT]
          call file_history_reg( "ENGI_t_MP_NF", "internal energy supply by the negative fixer", "J/m3/s",  & ! [IN]
                                 hist_nf_engi_id                                                            ) ! [OUT]
          call monitor_reg( "QTOTTND_NF", "vapor supply by the negative fixer", "kg",           & ! [IN]
                            monit_nf_mass_id,                                                   & ! [OUT]
                            isflux=.true.                                                       ) ! [IN]
          call monitor_reg( "ENGITND_NF", "internal energy supply by the negative fixer", "J",  & ! [IN]
                            monit_nf_engi_id,                                                   & ! [OUT]
                            isflux=.true.                                                       ) ! [IN]
          zero(:,:,:) = 0.0_rp
          call monitor_put( monit_nf_mass_id, zero(:,:,:) )
          call monitor_put( monit_nf_engi_id, zero(:,:,:) )
       end if

    else

       log_info("ATMOS_PHY_MP_driver_setup",*) 'this component is never called.'
       log_info("ATMOS_PHY_MP_driver_setup",*) 'SFLX_rain and SFLX_snow is set to zero.'
       sflx_rain(:,:) = 0.0_rp
       sflx_snow(:,:) = 0.0_rp

    endif


    return
  end subroutine atmos_phy_mp_driver_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_driver_adjustment
    use scale_const, only: &
       pre00 => const_pre00
    use scale_atmos_hydrometeor, only: &
       atmos_hydrometeor_dry, &
       i_qv, &
       qla, &
       qia, &
       qhs, &
       qhe
    use scale_atmos_phy_mp_common, only: &
       atmos_phy_mp_negative_fixer
    use mod_atmos_vars, only: &
       temp,  &
       cvtot, &
       cptot, &
       dens,  &
       rhot,  &
       qtrc
    use mod_atmos_phy_mp_vars, only: &
       qa_mp
    use scale_time, only: &
       dt => time_dtsec
    use scale_file_history, only: &
       file_history_query, &
       file_history_put
    use scale_monitor, only: &
       monitor_put

    real(RP) :: rhoh  (ka,ia,ja)
    real(RP) :: DENS_d(ka,ia,ja)
    real(RP) :: ENGI_d(ka,ia,ja)
    real(RP) :: Rtot

    logical :: do_put_rhoh
    logical :: do_put_dens
    logical :: do_put_engi

    integer :: k, i, j, iq

    if ( mp_do_negative_fixer .and. (.not. atmos_hydrometeor_dry) ) then

       call file_history_query( hist_nf_rhoh_id, do_put_rhoh )
       call file_history_query( hist_nf_dens_id, do_put_dens )
       call file_history_query( hist_nf_engi_id, do_put_engi )

       if ( monit_nf_mass_id > 0 .or. monit_nf_engi_id > 0 .or. &
            do_put_rhoh .or. do_put_dens .or. do_put_engi ) then
          call atmos_phy_mp_negative_fixer( &
               ka, ks, ke, ia, 1, ia, ja, 1, ja, qla, qia, &
               mp_limit_negative,                     & ! [IN]
               dens(:,:,:), temp(:,:,:),              & ! [INOUT]
               cvtot(:,:,:), cptot(:,:,:),            & ! [INOUT]
               qtrc(:,:,:,i_qv), qtrc(:,:,:,qhs:qhe), & ! [INOUT]
               rhoh = rhoh,                           & ! [OUT, optional]
               dens_diff = dens_d, engi_diff = engi_d ) ! [OUT, optional]
       else
          call atmos_phy_mp_negative_fixer( &
               ka, ks, ke, ia, 1, ia, ja, 1, ja, qla, qia, &
               mp_limit_negative,                    & ! [IN]
               dens(:,:,:), temp(:,:,:),             & ! [INOUT]
               cvtot(:,:,:), cptot(:,:,:),           & ! [INOUT]
               qtrc(:,:,:,i_qv), qtrc(:,:,:,qhs:qhe) ) ! [INOUT]
       end if

       !$omp parallel private(Rtot)

       !$omp do
       do j = 1, ja
       do i = 1, ia
       do k = ks, ke
          rtot = cptot(k,i,j) - cvtot(k,i,j)
          rhot(k,i,j) = pre00 / rtot * ( dens(k,i,j) * temp(k,i,j) * rtot / pre00 )**( cvtot(k,i,j) / cptot(k,i,j) )
       end do
       end do
       end do
       !$omp end do nowait

       ! for non-mass tracers, such as number density
       do iq = qhe+1, qa_mp
       !$omp do
       do j = 1, ja
       do i = 1, ia
       do k = ks, ke
          qtrc(k,i,j,iq) = max( qtrc(k,i,j,iq), 0.0_rp )
       end do
       end do
       end do
       end do

       !$omp end parallel

       if ( do_put_rhoh ) then
          !$omp parallel do
          do j = js, je
          do i = is, ie
          do k = ks, ke
             rhoh(k,i,j) = rhoh(k,i,j) / dt
          end do
          end do
          end do
          call file_history_put( hist_nf_rhoh_id, rhoh(:,:,:) )
       end if
       if ( monit_nf_mass_id > 0 .or. do_put_dens ) then
          !$omp parallel do
          do j = js, je
          do i = is, ie
          do k = ks, ke
             dens_d(k,i,j) = dens_d(k,i,j) / dt
          end do
          end do
          end do
          call file_history_put( hist_nf_dens_id, dens_d(:,:,:) )
          call monitor_put( monit_nf_mass_id, dens_d(:,:,:) )
       end if
       if ( monit_nf_engi_id > 0 .or. do_put_engi ) then
          !$omp parallel do
          do j = js, je
          do i = is, ie
          do k = ks, ke
             engi_d(k,i,j) = engi_d(k,i,j) / dt
          end do
          end do
          end do
          call file_history_put( hist_nf_engi_id, engi_d(:,:,:) )
          call monitor_put( monit_nf_engi_id, engi_d(:,:,:) )
       end if

    end if

    return
  end subroutine atmos_phy_mp_driver_adjustment

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_driver_calc_tendency( update_flag )
    use scale_const, only: &
       pre00 => const_pre00
    use scale_time, only: &
       dt_mp => time_dtsec_atmos_phy_mp
    use scale_atmos_grid_cartesc_real, only: &
       real_cz => atmos_grid_cartesc_real_cz, &
       real_fz => atmos_grid_cartesc_real_fz, &
       atmos_grid_cartesc_real_vol,       &
       atmos_grid_cartesc_real_totvol,    &
       atmos_grid_cartesc_real_volwxy,    &
       atmos_grid_cartesc_real_totvolwxy, &
       atmos_grid_cartesc_real_volzuy,    &
       atmos_grid_cartesc_real_totvolzuy, &
       atmos_grid_cartesc_real_volzxv,    &
       atmos_grid_cartesc_real_totvolzxv
    use scale_statistics, only: &
       statistics_checktotal, &
       statistics_total
    use scale_file_history, only: &
       file_history_in
    use scale_atmos_hydrometeor, only: &
       i_qv,  &
       n_hyd, &
       qha,   &
       qhs,   &
       qhe,   &
       qla,   &
       qia
    use scale_atmos_phy_mp_common, only: &
       atmos_phy_mp_precipitation_upwind,  &
       atmos_phy_mp_precipitation_semilag, &
       atmos_phy_mp_precipitation_momentum
    use scale_atmos_refstate, only: &
       refstate_dens => atmos_refstate_dens
    use scale_atmos_phy_mp_kessler, only: &
       atmos_phy_mp_kessler_adjustment, &
       atmos_phy_mp_kessler_terminal_velocity
    use scale_atmos_phy_mp_tomita08, only: &
       atmos_phy_mp_tomita08_adjustment, &
       atmos_phy_mp_tomita08_terminal_velocity
    use scale_atmos_phy_mp_sn14, only: &
       atmos_phy_mp_sn14_tendency, &
       atmos_phy_mp_sn14_terminal_velocity
    use scale_atmos_phy_mp_suzuki10, only: &
       atmos_phy_mp_suzuki10_tendency, &
       atmos_phy_mp_suzuki10_terminal_velocity
    use scale_file_history, only: &
       file_history_query, &
       file_history_put
    use mod_atmos_admin, only: &
       atmos_phy_mp_type,    &
       atmos_phy_precip_type
    use mod_atmos_vars, only: &
       dens   => dens_av, &
       momz   => momz_av, &
       u, &
       v, &
       w, &
       pott, &
       qtrc   => qtrc_av, &
       dens_t => dens_tp, &
       momz_t => momz_tp, &
       rhou_t => rhou_tp, &
       rhov_t => rhov_tp, &
!       RHOT_t => RHOT_tp, &
       rhoq_t => rhoq_tp, &
       rhoh   => rhoh_p, &
       temp, &
       pres, &
       qdry, &
       cvtot, &
       cptot, &
       exner, &
       rhot   => rhot_av
    use mod_atmos_phy_mp_vars, only: &
       qa_mp, &
       qs_mp, &
       qe_mp, &
       dens_t_mp => atmos_phy_mp_dens_t,    &
       momz_t_mp => atmos_phy_mp_momz_t,    &
!       RHOT_t_MP => ATMOS_PHY_MP_RHOT_t,    &
       rhou_t_mp => atmos_phy_mp_rhou_t,    &
       rhov_t_mp => atmos_phy_mp_rhov_t,    &
       rhoq_t_mp => atmos_phy_mp_rhoq_t,    &
       rhoh_mp   => atmos_phy_mp_rhoh,      &
       evaporate => atmos_phy_mp_evaporate, &
       sflx_rain => atmos_phy_mp_sflx_rain, &
       sflx_snow => atmos_phy_mp_sflx_snow
    use mod_atmos_phy_ae_vars, only: &
       ccn_t => atmos_phy_ae_ccn_t
    implicit none

    logical, intent(in) :: update_flag

    real(RP) :: RHOE_t(ka,ia,ja)
    real(RP) :: TEMP1 (ka,ia,ja)
    real(RP) :: CPtot1(ka,ia,ja)
    real(RP) :: CVtot1(ka,ia,ja)
    real(RP) :: CCN   (ka,ia,ja)
    real(RP) :: vterm (ka,qs_mp+1:qe_mp)
    real(RP), target :: QTRC1(ka,ia,ja,qs_mp:qe_mp)

    real(RP) :: FLX_hydro(ka)
    real(RP) :: DENS2    (ka)
    real(RP) :: TEMP2    (ka)
    real(RP) :: PRES2    (ka)
    real(RP) :: CPtot2   (ka)
    real(RP) :: CVtot2   (ka)
    real(RP) :: RHOE     (ka)
    real(RP) :: RHOE2    (ka)
    real(RP) :: RHOQ2    (ka,qs_mp+1:qe_mp)
    real(RP) :: mflux    (ka)
    real(RP) :: sflux    (2)

    real(RP) :: FZ  (ka)
    real(RP) :: FDZ (ka)
    real(RP) :: RFDZ(ka)
    real(RP) :: RCDZ(ka)

    real(RP) :: CPtot_t(ka,ia,ja), CVtot_t(ka,ia,ja)
    real(RP) :: CP_t, CV_t

    real(RP) :: precip   (ia,ja)

    ! for history output
    real(RP), allocatable :: vterm_hist(:,:,:,:)
    integer :: hist_vterm_idx(qs_mp+1:qe_mp)
    logical :: flag
    integer :: ih

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

    if ( update_flag ) then

       ccn(:,:,:) = ccn_t(:,:,:) * dt_mp

       select case ( atmos_phy_mp_type )
       case ( 'KESSLER' )
!OCL XFILL
          temp1(:,:,:) = temp(:,:,:)
!OCL XFILL
          qtrc1(:,:,:,qs_mp:qe_mp) = qtrc(:,:,:,qs_mp:qe_mp)
!OCL XFILL
          cvtot1(:,:,:) = cvtot(:,:,:)
!OCL XFILL
          cptot1(:,:,:) = cptot(:,:,:)

          call atmos_phy_mp_kessler_adjustment( &
               ka, ks, ke, ia, is, ie, ja, js, je, &
               dens(:,:,:), pres(:,:,:), dt_mp,                                      & ! [IN]
               temp1(:,:,:), qtrc1(:,:,:,qs_mp:qe_mp), cptot1(:,:,:), cvtot1(:,:,:), & ! [INOUT]
               rhoe_t(:,:,:), evaporate(:,:,:)                                       ) ! [OUT]

          do iq = qs_mp, qe_mp
          do j = js, je
          do i = is, ie
          do k = ks, ke
             rhoq_t_mp(k,i,j,iq) = ( qtrc1(k,i,j,iq) - qtrc(k,i,j,iq) ) * dens(k,i,j) / dt_mp
          enddo
          enddo
          enddo
          enddo

          do j = js, je
          do i = is, ie
          do k = ks, ke
             cptot_t(k,i,j) = ( cptot1(k,i,j) - cptot(k,i,j) ) / dt_mp
             cvtot_t(k,i,j) = ( cvtot1(k,i,j) - cvtot(k,i,j) ) / dt_mp
          end do
          end do
          end do

       case ( 'TOMITA08' )
!OCL XFILL
          do j = js, je
          do i = is, ie
          do k = ks, ke
             temp1(k,i,j) = temp(k,i,j)
          end do
          end do
          end do
!OCL XFILL
          do iq = qs_mp, qe_mp
          do j = js, je
          do i = is, ie
          do k = ks, ke
             qtrc1(k,i,j,iq) = qtrc(k,i,j,iq)
          end do
          end do
          end do
          end do
!OCL XFILL
          do j = js, je
          do i = is, ie
          do k = ks, ke
             cvtot1(k,i,j) = cvtot(k,i,j)
          end do
          end do
          end do
!OCL XFILL
          do j = js, je
          do i = is, ie
          do k = ks, ke
             cptot1(k,i,j) = cptot(k,i,j)
          end do
          end do
          end do

          call atmos_phy_mp_tomita08_adjustment( &
               ka, ks, ke, ia, is, ie, ja, js, je, &
               dens(:,:,:), pres(:,:,:), ccn(:,:,:), dt_mp,                          & ! [IN]
               temp1(:,:,:), qtrc1(:,:,:,qs_mp:qe_mp), cptot1(:,:,:), cvtot1(:,:,:), & ! [INOUT]
               rhoe_t(:,:,:), evaporate(:,:,:)                                       ) ! [OUT]

          do iq = qs_mp, qe_mp
          do j = js, je
          do i = is, ie
          do k = ks, ke
             rhoq_t_mp(k,i,j,iq) = ( qtrc1(k,i,j,iq) - qtrc(k,i,j,iq) ) * dens(k,i,j) / dt_mp
          enddo
          enddo
          enddo
          enddo

          do j = js, je
          do i = is, ie
          do k = ks, ke
             cptot_t(k,i,j) = ( cptot1(k,i,j) - cptot(k,i,j) ) / dt_mp
             cvtot_t(k,i,j) = ( cvtot1(k,i,j) - cvtot(k,i,j) ) / dt_mp
          end do
          end do
          end do

       case ( 'SN14' )

          call atmos_phy_mp_sn14_tendency( &
               ka, ks, ke, ia, is, ie, ja, js, je, &
               dens(:,:,:), w(:,:,:), qtrc(:,:,:,qs_mp:qe_mp), pres(:,:,:), temp(:,:,:),                     & ! [IN]
               qdry(:,:,:), cptot(:,:,:), cvtot(:,:,:), ccn(:,:,:), dt_mp, real_cz(:,:,:), real_fz(:,:,:),   & ! [IN]
               rhoq_t_mp(:,:,:,qs_mp:qe_mp), rhoe_t(:,:,:), cptot_t(:,:,:), cvtot_t(:,:,:), evaporate(:,:,:) ) ! [OUT]

       case ( 'SUZUKI10' )

          call atmos_phy_mp_suzuki10_tendency( ka, ks,  ke, ia, is, ie, ja, js, je, kijmax, &
                                               dt_mp,                                  & ! [IN]
                                               dens(:,:,:),  pres(:,:,:), temp(:,:,:), & ! [IN]
                                               qtrc(:,:,:,qs_mp:qe_mp), qdry(:,:,:),   & ! [IN]
                                               cptot(:,:,:), cvtot(:,:,:),             & ! [IN]
                                               ccn(:,:,:),                             & ! [IN]
                                               rhoq_t_mp(:,:,:,qs_mp:qe_mp),           & ! [OUT]
                                               rhoe_t(:,:,:),                          & ! [OUT]
                                               cptot_t(:,:,:), cvtot_t(:,:,:),         & ! [OUT]
                                               evaporate(:,:,:)                        ) ! [OUT]

       end select


       do j = js, je
       do i = is, ie
       do k = ks, ke
          rhoh_mp(k,i,j) = rhoe_t(k,i,j) &
                  - ( cptot_t(k,i,j) + log( pres(k,i,j) / pre00 ) * ( cvtot(k,i,j) / cptot(k,i,j) * cptot_t(k,i,j) - cvtot_t(k,i,j) ) ) &
                  * dens(k,i,j) * temp(k,i,j)
!          RHOT_t_MP(k,i,j) = RHOE_t(k,i,j) / ( EXNER(k,i,j) * CPtot(k,i,j) ) &
!                  - RHOT(k,i,j) * CPtot_t(k,i,j) / CPtot(k,i,j) &
!                  - RHOT(k,i,j) * CVtot(k,i,j) / ( CPtot(k,i,j) - CVtot(k,i,j) ) &
!                  * log( EXNER(k,i,j) ) * ( CPtot_t(k,i,j) / CPtot(k,i,j) - CVtot_t(k,i,j) / CVtot(k,i,j) )
       end do
       end do
       end do

       if ( mp_do_precipitation ) then

          call prof_rapstart('MP_Precipitation', 2)

          ! prepare for history output
          hist_vterm_idx(:) = -1
          ih = 0
          do iq = qs_mp+1, qe_mp
             call file_history_query( hist_vterm_id(iq), flag )
             if ( flag ) then
                ih = ih + 1
                hist_vterm_idx(iq) = ih
             end if
          end do
          if ( ih > 0 ) then
             allocate( vterm_hist(ka,ia,ja,ih) )
             vterm_hist(:,:,:,:) = 0.0_rp
          end if

          !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
          !$omp shared (KA,KS,KE,IS,IE,JS,JE,QS_MP,QE_MP,QHA,QLA,QIA, &
          !$omp         PRE00, &
          !$omp         ATMOS_PHY_MP_TYPE, ATMOS_PHY_PRECIP_TYPE, &
          !$omp         dt_MP,MP_NSTEP_SEDIMENTATION,MP_DTSEC_SEDIMENTATION,MP_RNSTEP_SEDIMENTATION, &
          !$omp         REAL_CZ,REAL_FZ, &
          !$omp         DENS,MOMZ,U,V,RHOT,TEMP,PRES,QTRC,CPtot,CVtot,EXNER, &
          !$omp         DENS_t_MP,MOMZ_t_MP,RHOU_t_MP,RHOV_t_MP,RHOQ_t_MP,RHOH_MP, &
          !$omp         SFLX_rain,SFLX_snow, &
          !$omp         REFSTATE_dens, &
          !$omp         vterm_hist,hist_vterm_idx) &
          !$omp private(i,j,k,iq,step, &
          !$omp         FZ,FDZ,RFDZ,RCDZ, &
          !$omp         DENS2,TEMP2,PRES2,CPtot2,CVtot2,RHOE,RHOE2,RHOQ2, &
          !$omp         vterm,mflux,sflux,FLX_hydro,CP_t,CV_t)
          do j = js, je
          do i = is, ie

             fz(1:ka) = real_fz(1:ka,i,j)

             fdz(ks-1) = real_cz(ks,i,j) - real_fz(ks-1,i,j)
             rfdz(ks-1) = 1.0_rp / fdz(ks-1)
             do k = ks, ke
                fdz(k) = real_cz(k+1,i,j) - real_cz(k  ,i,j)
                rfdz(k) = 1.0_rp / fdz(k)
                rcdz(k) = 1.0_rp / ( real_fz(k  ,i,j) - real_fz(k-1,i,j) )
             enddo

             do k = ks, ke
                dens2(k)  = dens(k,i,j)
                temp2(k)  = temp(k,i,j)
                pres2(k)  = pres(k,i,j)
                cptot2(k) = cptot(k,i,j)
                cvtot2(k) = cvtot(k,i,j)
                rhoe(k)   = temp(k,i,j) * cvtot(k,i,j) * dens2(k)
                rhoe2(k)  = rhoe(k)
             end do
             do iq = qs_mp+1, qe_mp
             do k = ks, ke
                rhoq2(k,iq) = dens2(k) * qtrc(k,i,j,iq)
             end do
             end do

             sflx_rain(i,j) = 0.0_rp
             sflx_snow(i,j) = 0.0_rp
             flx_hydro(:) = 0.0_rp
             do step = 1, mp_nstep_sedimentation

                select case ( atmos_phy_mp_type )
                case ( 'KESSLER' )
                   call atmos_phy_mp_kessler_terminal_velocity( &
                        ka, ks, ke, &
                        dens2(:), rhoq2(:,:), & ! [IN]
                        refstate_dens(:,i,j), & ! [IN]
                        vterm(:,:)            ) ! [OUT]
                case ( 'TOMITA08' )
                   call atmos_phy_mp_tomita08_terminal_velocity( &
                        ka, ks, ke, &
                        dens2(:), temp2(:), rhoq2(:,:), & ! [IN]
                        vterm(:,:)                      ) ! [OUT]
                case ( 'SN14' )
                   call atmos_phy_mp_sn14_terminal_velocity( &
                        ka, ks, ke, &
                        dens2(:), temp2(:), rhoq2(:,:), pres2(:), & ! [IN]
                        vterm(:,:)                                ) ! [OUT]
                case ( 'SUZUKI10' )
                   call atmos_phy_mp_suzuki10_terminal_velocity( &
                        ka,        & ! [IN]
                        vterm(:,:) ) ! [OUT]
                case default
                   vterm(:,:) = 0.0_rp ! tentative
                end select

                ! store to history output
                do iq = qs_mp+1, qe_mp
                   if ( hist_vterm_idx(iq) > 0 ) then
                      do k = ks, ke
                         vterm_hist(k,i,j,hist_vterm_idx(iq)) = vterm_hist(k,i,j,hist_vterm_idx(iq)) &
                                                              + vterm(k,iq) * mp_rnstep_sedimentation
                      end do
                   end if
                end do

                select case ( atmos_phy_precip_type )
                case ( 'Upwind-Euler' )
                   call atmos_phy_mp_precipitation_upwind( &
                        ka, ks, ke, qe_mp-qs_mp, qla, qia, &
                        temp2(:), vterm(:,:),   & ! [IN]
                        fdz(:), rcdz(:),        & ! [IN]
                        mp_dtsec_sedimentation, & ! [IN]
                        i, j,                   & ! [IN]
                        dens2(:), rhoq2(:,:),   & ! [INOUT]
                        cptot2(:), cvtot2(:),   & ! [INOUT]
                        rhoe2(:),               & ! [INOUT]
                        mflux(:), sflux(:)      ) ! [OUT]
                case ( 'Semilag' )
                   call atmos_phy_mp_precipitation_semilag( &
                        ka, ks, ke, qe_mp-qs_mp, qla, qia, &
                        temp2(:), vterm(:,:),   & ! [IN]
                        fz(:), fdz(:), rcdz(:), & ! [IN]
                        mp_dtsec_sedimentation, & ! [IN]
                        i, j,                   & ! [IN]
                        dens2(:), rhoq2(:,:),   & ! [INOUT]
                        cptot2(:), cvtot2(:),   & ! [INOUT]
                        rhoe2(:),               & ! [INOUT]
                        mflux(:), sflux(:)      ) ! [OUT]
                case default
                   mflux(:) = 0.0_rp
                   sflux(:) = 0.0_rp
                end select

                do k = ks, ke
                   temp2(k) = rhoe2(k) / ( dens2(k) * cvtot2(k) )
                end do

                do k = ks-1, ke-1
                   flx_hydro(k) = flx_hydro(k) + mflux(k) * mp_rnstep_sedimentation
                enddo

                sflx_rain(i,j) = sflx_rain(i,j) - sflux(1) * mp_rnstep_sedimentation
                sflx_snow(i,j) = sflx_snow(i,j) - sflux(2) * mp_rnstep_sedimentation

             enddo

!OCL XFILL
             do k = ks, ke
                dens_t_mp(k,i,j) = ( dens2(k) - dens(k,i,j) ) / dt_mp
             end do

             do k = ks, ke
                cp_t = ( cptot2(k) - cptot(k,i,j) ) / dt_mp
                cv_t = ( cvtot2(k) - cvtot(k,i,j) ) / dt_mp
                rhoh_mp(k,i,j) = rhoh_mp(k,i,j) &
                     + ( rhoe2(k) - rhoe(k) ) / dt_mp &
                     - ( cp_t + log( pres(k,i,j) / pre00 ) * ( cvtot(k,i,j) / cptot(k,i,j) * cp_t - cv_t ) ) &
                     * dens(k,i,j) * temp(k,i,j)
!                RHOT_t_MP(k,i,j) = RHOT_t_MP(k,i,j) &
!                     + ( RHOE2(k) - RHOE(k) ) / ( dt_MP * EXNER(k,i,j) * CPtot(k,i,j) ) &
!                     - RHOT(k,i,j) * CP_t / CPtot(k,i,j) &
!                     - RHOT(k,i,j) * CVtot(k,i,j) / ( CPtot(k,i,j) - CVtot(k,i,j) ) &
!                     * log( EXNER(k,i,j) ) * ( CP_t / CPtot(k,i,j) - CV_t / CVtot(k,i,j) )
             end do

             do iq = qs_mp+1, qe_mp
             do k  = ks, ke
                rhoq_t_mp(k,i,j,iq) = rhoq_t_mp(k,i,j,iq) &
                     + ( rhoq2(k,iq) - dens(k,i,j) * qtrc(k,i,j,iq) ) / dt_mp
             enddo
             enddo

             call atmos_phy_mp_precipitation_momentum( &
                     ka, ks, ke, &
                     dens(:,i,j), momz(:,i,j), u(:,i,j), v(:,i,j),        & ! [IN]
                     flx_hydro(:),                                        & ! [IN]
                     rcdz(:), rfdz(:),                                    & ! [IN]
                     momz_t_mp(:,i,j), rhou_t_mp(:,i,j), rhov_t_mp(:,i,j) ) ! [OUT]

          enddo
          enddo

          ! history output
          do iq = qs_mp+1, qe_mp
             if ( hist_vterm_idx(iq) > 0 ) &
                call file_history_put( hist_vterm_id(iq), vterm_hist(:,:,:,hist_vterm_idx(iq)) )
          end do
          if ( allocated( vterm_hist ) ) deallocate( vterm_hist )

          call prof_rapend  ('MP_Precipitation', 2)

       end if

!OCL XFILL
       do j = js, je
       do i = is, ie
          precip(i,j) = sflx_rain(i,j) + sflx_snow(i,j)
       end do
       end do

       call file_history_in( sflx_rain(:,:),   'RAIN_MP',   'surface rain rate by MP',          'kg/m2/s',  fill_halo=.true. )
       call file_history_in( sflx_snow(:,:),   'SNOW_MP',   'surface snow rate by MP',          'kg/m2/s',  fill_halo=.true. )
       call file_history_in( precip(:,:),   'PREC_MP',   'surface precipitation rate by MP', 'kg/m2/s',  fill_halo=.true. )
       call file_history_in( evaporate(:,:,:), 'EVAPORATE', 'evaporated cloud number',          'num/m3/s', fill_halo=.true. )

       call file_history_in( dens_t_mp(:,:,:), 'DENS_t_MP', 'tendency DENS in MP',         'kg/m3/s'  , fill_halo=.true. )
       call file_history_in( momz_t_mp(:,:,:), 'MOMZ_t_MP', 'tendency MOMZ in MP',         'kg/m2/s2' , fill_halo=.true. )
       call file_history_in( rhou_t_mp(:,:,:), 'RHOU_t_MP', 'tendency RHOU in MP',         'kg/m2/s2' , fill_halo=.true. )
       call file_history_in( rhov_t_mp(:,:,:), 'RHOV_t_MP', 'tendency RHOV in MP',         'kg/m2/s2' , fill_halo=.true. )
!       call FILE_HISTORY_in( RHOT_t_MP(:,:,:), 'RHOT_t_MP', 'tendency RHOT in MP',         'K*kg/m3/s', fill_halo=.true. )
       call file_history_in( rhoh_mp(:,:,:), 'RHOH_MP',   'diabatic heating rate in MP', 'J/m3/s',    fill_halo=.true. )
       do iq = qs_mp, qe_mp
          call file_history_in( rhoq_t_mp(:,:,:,iq), trim(tracer_name(iq))//'_t_MP', &
                        'tendency rho*'//trim(tracer_name(iq))//' in MP', 'kg/m3/s', fill_halo=.true. )
       enddo

    endif

    !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
    !$omp shared(KS,KE,IS,IE,JS,JE, &
    !$omp        DENS_t_MP,MOMZ_t_MP,RHOU_t_MP,RHOV_t_MP,RHOH_MP, &
    !$omp        DENS_t,MOMZ_t,RHOU_t,RHOV_t,RHOH)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       dens_t(k,i,j) = dens_t(k,i,j) + dens_t_mp(k,i,j)
       momz_t(k,i,j) = momz_t(k,i,j) + momz_t_mp(k,i,j)
       rhou_t(k,i,j) = rhou_t(k,i,j) + rhou_t_mp(k,i,j)
       rhov_t(k,i,j) = rhov_t(k,i,j) + rhov_t_mp(k,i,j)
!       RHOT_t(k,i,j) = RHOT_t(k,i,j) + RHOT_t_MP(k,i,j)
       rhoh(k,i,j) = rhoh(k,i,j) + rhoh_mp(k,i,j)
    enddo
    enddo
    enddo

    do iq = qs_mp, qe_mp
    !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ &
    !$omp shared(JS,JE,IS,IE,KS,KE,RHOQ_t,iq,RHOQ_t_MP)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       rhoq_t(k,i,j,iq) = rhoq_t(k,i,j,iq) + rhoq_t_mp(k,i,j,iq)
    enddo
    enddo
    enddo
    enddo

    if ( statistics_checktotal ) then
       call statistics_total( ka, ks, ke, ia, is, ie, ja, js, je, &
                              dens_t_mp(:,:,:), 'DENS_t_MP',       &
                              atmos_grid_cartesc_real_vol(:,:,:),  &
                              atmos_grid_cartesc_real_totvol       )
       call statistics_total( ka, ks, ke, ia, is, ie, ja, js, je, &
                              momz_t_mp(:,:,:), 'MOMZ_t_MP',         &
                              atmos_grid_cartesc_real_volwxy(:,:,:), &
                              atmos_grid_cartesc_real_totvolwxy      )
       call statistics_total( ka, ks, ke, ia, is, ie, ja, js, je, &
                              rhoh_mp(:,:,:),   'RHOH_MP',      &
                              atmos_grid_cartesc_real_vol(:,:,:), &
                              atmos_grid_cartesc_real_totvol      )
!       call STATISTICS_total( KA, KS, KE, IA, IS, IE, JA, JS, JE, &
!                              RHOT_t_MP(:,:,:), 'RHOT_t_MP',      &
!                              ATMOS_GRID_CARTESC_REAL_VOL(:,:,:), &
!                              ATMOS_GRID_CARTESC_REAL_TOTVOL      )

       do iq = qs_mp, qe_mp
          call statistics_total( ka, ks, ke, ia, is, ie, ja, js, je, &
                                 rhoq_t_mp(:,:,:,iq), trim(tracer_name(iq))//'_t_MP', &
                                 atmos_grid_cartesc_real_vol(:,:,:),                  &
                                 atmos_grid_cartesc_real_totvol                       )
       enddo
    endif

    return
  end subroutine atmos_phy_mp_driver_calc_tendency

  subroutine atmos_phy_mp_driver_qhyd2qtrc( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       QV, QHYD, &
       QTRC,     &
       QNUM      )
    use scale_prc, only: &
       prc_abort
    use scale_atmos_hydrometeor, only: &
       n_hyd
    use mod_atmos_phy_mp_vars, only: &
       qa_mp
    use mod_atmos_admin, only: &
       atmos_phy_mp_type
    use scale_atmos_phy_mp_kessler, only: &
       atmos_phy_mp_kessler_qhyd2qtrc
    use scale_atmos_phy_mp_tomita08, only: &
       atmos_phy_mp_tomita08_qhyd2qtrc
    use scale_atmos_phy_mp_sn14, only: &
       atmos_phy_mp_sn14_qhyd2qtrc
    use scale_atmos_phy_mp_suzuki10, only: &
       atmos_phy_mp_suzuki10_qhyd2qtrc
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in) :: QV  (ka,ia,ja)
    real(RP), intent(in) :: QHYD(ka,ia,ja,n_hyd)

    real(RP), intent(out) :: QTRC(ka,ia,ja,qa_mp)

    real(RP), intent(in), optional :: QNUM(ka,ia,ja,n_hyd)

    integer :: k, i, j

    select case( atmos_phy_mp_type )
    case ( "NONE" )
       if ( associated( atmos_phy_mp_user_qhyd2qtrc ) ) then
          call atmos_phy_mp_user_qhyd2qtrc( ka, ks, ke, ia, is, ie, ja, js, je, &
                                            qv(:,:,:), qhyd(:,:,:,:), & ! [IN]
                                            qtrc(:,:,:,:),            & ! [OUT]
                                            qnum=qnum                 ) ! [IN]
       end if
    case ( "KESSLER" )
       !$omp parallel do OMP_SCHEDULE_
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,1) = qv(k,i,j)
       end do
       end do
       end do
       call atmos_phy_mp_kessler_qhyd2qtrc( ka, ks, ke, ia, is, ie, ja, js, je, &
                                            qhyd(:,:,:,:),  & ! [IN]
                                            qtrc(:,:,:,2:)  ) ! [OUT]
    case ( "TOMITA08" )
       !$omp parallel do OMP_SCHEDULE_
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,1) = qv(k,i,j)
       end do
       end do
       end do
       call atmos_phy_mp_tomita08_qhyd2qtrc( ka, ks, ke, ia, is, ie, ja, js, je, &
                                             qhyd(:,:,:,:),  & ! [IN]
                                             qtrc(:,:,:,2:)  ) ! [OUT]
    case ( "SN14" )
       !$omp parallel do OMP_SCHEDULE_
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,1) = qv(k,i,j)
       end do
       end do
       end do
       call atmos_phy_mp_sn14_qhyd2qtrc( ka, ks, ke, ia, is, ie, ja, js, je, &
                                         qhyd(:,:,:,:),  & ! [IN]
                                         qtrc(:,:,:,2:), & ! [OUT]
                                         qnum=qnum       ) ! [IN]
    case ( "SUZUKI10" )
       !$omp parallel do OMP_SCHEDULE_
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,1) = qv(k,i,j)
       end do
       end do
       end do
       call atmos_phy_mp_suzuki10_qhyd2qtrc( ka, ks, ke, ia, is, ie, ja, js, je, &
                                             qhyd(:,:,:,:),  & ! [IN]
                                             qtrc(:,:,:,2:), & ! [OUT]
                                             qnum=qnum       ) ! [IN]
    case default
       log_error("ATMOS_PHY_MP_driver_qhyd2qtrc",*) 'ATMOS_PHY_MP_TYPE (', trim(atmos_phy_mp_type), ') is not supported'
       call prc_abort
    end select

    return
  end subroutine atmos_phy_mp_driver_qhyd2qtrc

end module mod_atmos_phy_mp_driver