scale_atmos_phy_mp_common.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_mp_common
  !-----------------------------------------------------------------------------
  !
  !++ 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_negative_fixer
  public :: atmos_phy_mp_saturation_adjustment
  public :: atmos_phy_mp_precipitation_upwind
  public :: atmos_phy_mp_precipitation_semilag
  public :: atmos_phy_mp_precipitation_momentum

  interface atmos_phy_mp_saturation_adjustment
     module procedure atmos_phy_mp_saturation_adjustment_3d
  end interface atmos_phy_mp_saturation_adjustment

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_negative_fixer( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, QLA, QIA, &
       limit_negative, &
       DENS, TEMP,     &
       CVtot, CPtot,   &
       QV, QTRC,       &
       RHOH,           &
       DENS_diff,      &
       ENGI_diff       )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       cvdry => const_cvdry, &
       cpdry => const_cpdry
    use scale_atmos_hydrometeor, only: &
       atmos_hydrometeor_lhv, &
       atmos_hydrometeor_lhs, &
       cv_vapor, &
       cp_vapor, &
       cv_water, &
       cp_water, &
       cv_ice,   &
       cp_ice
    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
    integer, intent(in) :: QLA, QIA

    real(RP), intent(in)    :: limit_negative

    real(RP), intent(inout) :: DENS (ka,ia,ja)
    real(RP), intent(inout) :: TEMP (ka,ia,ja)
    real(RP), intent(inout) :: CVtot(ka,ia,ja)
    real(RP), intent(inout) :: CPtot(ka,ia,ja)
    real(RP), intent(inout) :: QV   (ka,ia,ja)
    real(RP), intent(inout) :: QTRC (ka,ia,ja,qla+qia)

    real(RP), intent(out), optional :: RHOH     (ka,ia,ja)
    real(RP), intent(out), optional :: DENS_diff(ka,ia,ja)
    real(RP), intent(out), optional :: ENGI_diff(ka,ia,ja)

    real(RP) :: LHV(ka)
    real(RP) :: LHS(ka)
    real(RP) :: eng  (ka)
    real(RP) :: eng0 (ka)
    real(RP) :: dens0(ka)
    real(RP) :: diffq(ka)
    real(RP) :: diffq_min
    real(RP) :: work

    logical :: rhoh_out
    logical :: dens_out
    logical :: engi_out

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

    call prof_rapstart('MP_filter', 3)

    rhoh_out = present( rhoh )
    dens_out = present( dens_diff )
    engi_out = present( engi_diff )

    diffq_min = 0.0_rp

    !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
    !$omp reduction(min:diffq_min) &
    !$omp private(i,j,k,iq, &
    !$omp         LHV,LHS,eng,eng0,dens0,diffq,work) &
    !$omp shared(KA,KS,KE,IS,IE,JS,JE,QLA,QIA, &
    !$omp        CVdry,CPdry,CV_VAPOR,CP_VAPOR,CV_WATER,CP_WATER,CV_ICE,CP_ICE, &
    !$omp        DENS,TEMP,CVtot,CPtot,QV,QTRC, &
    !$omp        RHOH,DENS_diff,ENGI_diff,rhoh_out,dens_out,engi_out,limit_negative)
    do j = js, je
    do i = is, ie

       diffq(:) = 0.0_rp

       do k = ks, ke
          eng(k) = cvtot(k,i,j) * temp(k,i,j)
       end do

       if ( rhoh_out ) then
          do k = ks, ke
             eng0(k) = eng(k)
          end do
       end if

       if ( qla > 0 ) then
          call atmos_hydrometeor_lhv( ka, ks, ke, temp(:,i,j), lhv(:) )

          do iq = 1, qla
          do k = ks, ke
             work = - min( qtrc(k,i,j,iq), 0.0_rp ) ! work is positive (vapor to liq)
             if ( work > 0.0_rp ) then
                eng(k) = eng(k) + work * lhv(k)
                diffq(k) = diffq(k) - work
                cvtot(k,i,j) = cvtot(k,i,j) + work * ( cv_water - cv_vapor )
                cptot(k,i,j) = cptot(k,i,j) + work * ( cp_water - cp_vapor )
                ! remove negative value of hydrometeors (mass)
                qtrc(k,i,j,iq) = 0.0_rp
             end if
          enddo
          enddo
       end if

       if ( qia > 0 ) then
          call atmos_hydrometeor_lhs( ka, ks, ke, temp(:,i,j), lhs(:) )

          do iq = qla+1, qla+qia
          do k = ks, ke
             work = - min( qtrc(k,i,j,iq), 0.0_rp ) ! work is positive (vapor to ice)
             if ( work > 0.0_rp ) then
                eng(k) = eng(k) + work * lhs(k)
                diffq(k) = diffq(k) - work
                cvtot(k,i,j) = cvtot(k,i,j) + work * ( cv_ice - cv_vapor )
                cptot(k,i,j) = cptot(k,i,j) + work * ( cp_ice - cp_vapor )
                ! remove negative value of hydrometeors (mass)
                qtrc(k,i,j,iq) = 0.0_rp
             end if
          enddo
          enddo
       end if


       if ( abs(limit_negative) > 0.0_rp ) then
          do k = ks, ke
             if ( diffq(k) < - abs(limit_negative) ) then
                diffq_min = min( diffq_min, diffq(k) )
                log_error("ATMOS_PHY_MP_negative_fixer",*) 'large negative is found'
                log_error_cont(*) 'value = ', diffq(k), ' at (', k, ',', i, ',', j, ')'
             end if
          end do
       end if

       do k = ks, ke
          if ( diffq(k) < 0.0_rp ) then
             ! Compensate for the lack of hydrometeors by the water vapor
             qv(k,i,j) = qv(k,i,j) + diffq(k)
          end if
       end do

       if ( rhoh_out ) then
          do k = ks, ke
             rhoh(k,i,j) = ( eng(k) - eng0(k) ) * dens(k,i,j)
          end do
       end if


       if ( dens_out ) then
          do k = ks, ke
             dens0(k) = dens(k,i,j)
          end do
       end if
       if ( engi_out ) then
          do k = ks, ke
             eng0(k) = eng(k) * dens(k,i,j)
          end do
       end if

       ! fix negative QV
       do k = ks, ke
          work = - min( qv(k,i,j), 0.0_rp )
          if ( work > 0.0_rp ) then
             qv(k,i,j) = 0.0_rp
             dens(k,i,j) = dens(k,i,j) * ( 1.0_rp + work ) ! not to change mass of dry air
             cvtot(k,i,j) = ( cvtot(k,i,j) + work * cv_vapor ) / ( 1.0_rp + work )
             cptot(k,i,j) = ( cptot(k,i,j) + work * cp_vapor ) / ( 1.0_rp + work )
             eng(k) = ( eng(k) + work * cv_vapor * temp(k,i,j) ) / ( 1.0_rp + work )
          end if
       enddo

       do k = ks, ke
          ! update
          temp(k,i,j) = eng(k) / cvtot(k,i,j)
       end do

       if ( dens_out ) then
          do k = ks, ke
             dens_diff(k,i,j) = dens(k,i,j) - dens0(k)
          end do
       end if
       if ( engi_out ) then
          do k = ks, ke
             engi_diff(k,i,j) = eng(k) * dens(k,i,j) - eng0(k)
          end do
       end if

    enddo
    enddo


    if (       abs(limit_negative) > 0.0_rp         &
         .AND. abs(limit_negative) < abs(diffq_min) ) then
       log_error_cont(*) 'maximum negative hydrometeor ', diffq_min, ' < ', - abs(limit_negative)
       call prc_abort
    endif

    call prof_rapend('MP_filter', 3)

    return
  end subroutine atmos_phy_mp_negative_fixer

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_saturation_adjustment_3d( &
         KA, KS, KE, IA, IS, IE, JA, JS, JE, &
         DENS,         &
         flag_liquid,  &
         TEMP,         &
         QV, QC, QI,   &
         CPtot, CVtot, &
         RHOE_d        )
    use scale_prc, only: &
       prc_abort
    use scale_atmos_hydrometeor, only: &
       cp_vapor, &
       cp_water, &
       cp_ice,   &
       cv_vapor, &
       cv_water, &
       cv_ice,   &
       lhv,   &
       lhf
    use scale_atmos_saturation, only: &
       atmos_saturation_moist_conversion_dens_all, &
       atmos_saturation_moist_conversion_dens_liq
    implicit none

    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in) :: DENS(ka,ia,ja)
    logical , intent(in) :: flag_liquid

    real(RP), intent(inout) :: TEMP (ka,ia,ja)
    real(RP), intent(inout) :: QV   (ka,ia,ja)
    real(RP), intent(inout) :: QC   (ka,ia,ja)
    real(RP), intent(inout) :: QI   (ka,ia,ja)
    real(RP), intent(inout) :: CPtot(ka,ia,ja)
    real(RP), intent(inout) :: CVtot(ka,ia,ja)

    real(RP), intent(out) :: RHOE_d(ka,ia,ja)

    ! working
    real(RP) :: QV1
    real(RP) :: QC1
    real(RP) :: QI1

    real(RP) :: Emoist ! moist internal energy

    logical :: converged, error

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

    call prof_rapstart('MP_Saturation_adjustment', 2)

    error = .false.

    if ( flag_liquid ) then ! warm rain

       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(KS,KE,IS,IE,JS,JE, &
       !$omp        CP_VAPOR,CP_WATER,CV_VAPOR,CV_WATER,LHV,LHF, &
       !$omp        DENS,QV,QC,TEMP,CPtot,CVtot,RHOE_d,error) &
       !$omp private(i,j,k, &
       !$omp         QV1,QC1,Emoist,converged)
       do j = js, je
       do i = is, ie
       do k = ks, ke

          qv1 = qv(k,i,j)
          qc1 = qc(k,i,j)

          emoist = temp(k,i,j) * cvtot(k,i,j) &
                 + qv1 * lhv

          call atmos_saturation_moist_conversion_dens_liq( dens(k,i,j), emoist,        & ! [IN]
                                                           temp(k,i,j), qv1, qc1,      & ! [INOUT]
                                                           cptot(k,i,j), cvtot(k,i,j), & ! [INOUT]
                                                           converged                   ) ! [OUT]

          if ( .NOT. converged ) then
             log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged! ', k,i,j
             error = .true.
             exit
          endif

          rhoe_d(k,i,j) = - lhv * ( qv1 - qv(k,i,j) ) * dens(k,i,j)

          qv(k,i,j) = qv1
          qc(k,i,j) = qc1

       end do
       end do
       end do

       if ( error ) call prc_abort

    else ! cold rain

       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp shared (KS,KE,IS,IE,JS,JE, &
       !$omp         CP_VAPOR,CP_WATER,CP_ICE,CV_VAPOR,CV_WATER,CV_ICE,LHV,LHF, &
       !$omp         DENS,QV,QC,QI,TEMP,CPtot,CVtot,RHOE_d,error) &
       !$omp private(i,j,k, &
       !$omp         QV1,QC1,QI1,Emoist,converged)
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qv1 = qv(k,i,j)
          qc1 = qc(k,i,j)
          qi1 = qi(k,i,j)

          emoist = temp(k,i,j) * cvtot(k,i,j) &
                 + qv1 * lhv &
                 - qi1 * lhf

          call atmos_saturation_moist_conversion_dens_all( dens(k,i,j), emoist,        & ! [IN]
                                                           temp(k,i,j), qv1, qc1, qi1, & ! [INOUT]
                                                           cptot(k,i,j), cvtot(k,i,j), & ! [INOUT]
                                                           converged                   ) ! [OUT]

          if ( .NOT. converged ) then
             log_error("ATMOS_PHY_MP_saturation_adjustment_3D",*) 'moist_conversion not converged! ', k,i,j
             error = .true.
             exit
          endif

          rhoe_d(k,i,j) = ( - lhv * ( qv1 - qv(k,i,j) ) &
                            + lhf * ( qi1 - qi(k,i,j) ) ) * dens(k,i,j)

          qv(k,i,j) = qv1
          qc(k,i,j) = qc1
          qi(k,i,j) = qi1

       end do
       end do
       end do

       if ( error ) call prc_abort
    endif

    call prof_rapend  ('MP_Saturation_adjustment', 2)

    return
  end subroutine atmos_phy_mp_saturation_adjustment_3d

  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_phy_mp_precipitation_upwind( &
       KA, KS, KE, QHA, QLA, QIA, &
       TEMP, vterm, FDZ, RCDZ, dt,     &
       i, j,                           &
       DENS, RHOQ, CPtot, CVtot, RHOE, &
       mflx, sflx                      )
    use scale_const, only: &
       grav  => const_grav
    use scale_atmos_hydrometeor, only: &
       cp_water, &
       cp_ice, &
       cv_water, &
       cv_ice
    implicit none
    integer,  intent(in) :: KA, KS, KE
    integer,  intent(in) :: QHA, QLA, QIA

    real(RP), intent(in) :: TEMP (ka)
    real(RP), intent(in) :: vterm(ka,qha) ! terminal velocity of cloud mass
    real(RP), intent(in) :: FDZ  (ka)
    real(RP), intent(in) :: RCDZ (ka)
    real(DP), intent(in) :: dt
    integer,  intent(in) :: i, j         ! for debug

    real(RP), intent(inout) :: DENS (ka)
    real(RP), intent(inout) :: RHOQ (ka,qha)
    real(RP), intent(inout) :: CPtot(ka)
    real(RP), intent(inout) :: CVtot(ka)
    real(RP), intent(inout) :: RHOE (ka)

    real(RP), intent(out)   :: mflx (ka)
    real(RP), intent(out)   :: sflx (2)

    real(RP) :: vtermh(ka)
    real(RP) :: qflx  (ka)
    real(RP) :: eflx  (ka)
    real(RP) :: RHOCP (ka)
    real(RP) :: RHOCV (ka)
    real(RP) :: dDENS
    real(RP) :: CP, CV

    integer  :: k, iq
    !---------------------------------------------------------------------------

    ! tracer/energy transport by falldown
    ! 1st order upwind, forward euler, velocity is always negative

    mflx(:) = 0.0_rp
    sflx(:) = 0.0_rp
    qflx(ke) = 0.0_rp
    eflx(ke) = 0.0_rp

    do k = ks, ke
       rhocp(k) = cptot(k) * dens(k)
       rhocv(k) = cvtot(k) * dens(k)
    end do

    do iq = 1, qha
       do k = ks, ke-1
          vtermh(k) = 0.5_rp * ( vterm(k+1,iq) + vterm(k,iq) )
       enddo
       vtermh(ks-1) = vterm(ks,iq)

       !--- mass flux for each tracer, upwind with vel < 0
       do k = ks-1, ke-1
          qflx(k) = vtermh(k) * rhoq(k+1,iq)
       enddo

       !--- update falling tracer
       do k  = ks, ke
          rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx(k) - qflx(k-1) ) * rcdz(k)
       enddo ! falling (water mass & number) tracer

       ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
       if ( iq > qla + qia ) cycle

       do k = ks-1, ke-1
          mflx(k) = mflx(k) + qflx(k)
       end do

       if ( iq > qla ) then ! ice water
          cp = cp_ice
          cv = cv_ice
          sflx(2) = sflx(2) + qflx(ks-1)
       else                 ! liquid water
          cp = cp_water
          cv = cv_water
          sflx(1) = sflx(1) + qflx(ks-1)
       end if

       !--- update density
       do k = ks, ke
          ddens = - ( qflx(k) - qflx(k-1) ) * rcdz(k) * dt
          rhocp(k) = rhocp(k) + cp * ddens
          rhocv(k) = rhocv(k) + cv * ddens
          dens(k) = dens(k) + ddens
       end do

       ! internal energy flux
       do k = ks-1, ke-1
          eflx(k) = qflx(k) * temp(k+1) * cv &
                  + qflx(k) * fdz(k) * grav               ! potential energy
       end do
       !--- update internal energy
       do k = ks, ke
          rhoe(k) = rhoe(k) - ( eflx(k) - eflx(k-1) ) * rcdz(k) * dt
       end do

    end do

    do k = ks, ke
       cptot(k) = rhocp(k) / dens(k)
       cvtot(k) = rhocv(k) / dens(k)
    end do

    return
  end subroutine atmos_phy_mp_precipitation_upwind

  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_phy_mp_precipitation_semilag( &
       KA, KS, KE, QHA, QLA, QIA, &
       TEMP, vterm, FZ, FDZ, RCDZ, dt, &
       i, j,                           &
       DENS, RHOQ, CPtot, CVtot, RHOE, &
       mflx, sflx                      )
    use scale_const, only: &
       grav  => const_grav
    use scale_atmos_hydrometeor, only: &
       cp_water, &
       cp_ice, &
       cv_water, &
       cv_ice
    implicit none
    integer,  intent(in) :: KA, KS, KE
    integer,  intent(in) :: QHA, QLA, QIA

    real(RP), intent(in) :: TEMP (ka)
    real(RP), intent(in) :: vterm(ka,qha) ! terminal velocity of cloud mass
    real(RP), intent(in) :: FZ   (ka)
    real(RP), intent(in) :: FDZ  (ka)
    real(RP), intent(in) :: RCDZ (ka)
    real(DP), intent(in) :: dt
    integer,  intent(in) :: i, j         ! for debug

    real(RP), intent(inout) :: DENS (ka)
    real(RP), intent(inout) :: RHOQ (ka,qha)
    real(RP), intent(inout) :: CPtot(ka)
    real(RP), intent(inout) :: CVtot(ka)
    real(RP), intent(inout) :: RHOE (ka)

    real(RP), intent(out)   :: mflx (ka)
    real(RP), intent(out)   :: sflx (2)

    real(RP) :: qflx(ka)
    real(RP) :: eflx(ka)
    real(RP) :: RHOCP(ka)
    real(RP) :: RHOCV(ka)
    real(RP) :: dDENS
    real(RP) :: CP, CV

    real(RP) :: vtermh(ka)
    real(RP) :: dvterm(ka)
    real(RP) :: cdz   (ka)
    real(RP) :: rfdz2 (ka)
    real(RP) :: dist  (ka)
    real(RP) :: Z_src
    real(RP) :: flx
    integer  :: k_src (ka)
    integer  :: k_dst

    integer  :: k, iq
    !---------------------------------------------------------------------------

    ! tracer/energy transport by falldown
    ! velocity is always negative

    mflx(:) = 0.0_rp
    sflx(:) = 0.0_rp
    qflx(:) = 0.0_rp
    eflx(:) = 0.0_rp

    do k = ks, ke
       rhocp(k) = cptot(k) * dens(k)
       rhocv(k) = cvtot(k) * dens(k)
    end do

    do k = ks, ke
       cdz(k) = 1.0_rp / rcdz(k)
    end do
    do k = ks, ke-1
       rfdz2(k) = 1.0_rp / ( cdz(k) + cdz(k+1) )
    end do

    do iq = 1, qha
       do k = ks, ke-1
          vtermh(k) = ( cdz(k) * vterm(k+1,iq) + cdz(k+1) * vterm(k,iq) ) * rfdz2(k)
       enddo
       vtermh(ks-1) = vterm(ks,iq)

       do k = ks, ke
          dvterm(k) = vtermh(k) - vtermh(k-1)
       enddo

       ! Movement distance of the cell wall by the fall
       ! the midpoint method (second-order Runge-Kutta)
       ! dz/dt = v(z + v dt/2) ~ v(z) + v dt/2 dv/dz + 1/2 (v dt/2)^2 d^2v/dz^2
       do k = ks, ke-1
          dist(k) = - vtermh(k)    * dt                                                                        &
                    + vtermh(k)    * dt**2 / 2.0_rp * ( dvterm(k+1)+dvterm(k) ) * rfdz2(k)                     &
                    - vtermh(k)**2 * dt**3 / 4.0_rp * ( dvterm(k+1)*rcdz(k+1) - dvterm(k)*rcdz(k) ) * rfdz2(k)
          dist(k) = max( dist(k), 0.0_rp )
       enddo
       dist(ks-1) = - vtermh(ks-1) * dt &
                    + vtermh(ks-1) * dt**2 / 2.0_rp * dvterm(ks)*rcdz(ks)
       dist(ks-1) = max( dist(ks-1), 0.0_rp )

       ! wall cannot overtake
       do k = ke-2, ks-1, -1
          dist(k) = min( dist(k), dist(k+1) + cdz(k+1) )
       end do

!        LOG_INFO_CONT(*) "distance", iq
!        do k = KA, 1, -1
!           LOG_INFO_CONT('(1x,I5,3F9.3,ES15.5)') k, dist(k), vtermh(k), vterm(k,iq), RHOQ(k,iq)
!        enddo

       ! search number of source cell
       do k_dst = ks-1, ke-1
          z_src = fz(k_dst) + dist(k_dst)

          k_src(k_dst) = k_dst
          do k = k_dst, ke-1
             if (       z_src >  fz(k  ) &
                  .AND. z_src <= fz(k+1) ) then
                k_src(k_dst) = k
             endif
          enddo
          if ( z_src > fz(ke) ) k_src(k_dst) = ke
       enddo

!        LOG_INFO_CONT(*) "seek", iq
!        do k = KA, 1, -1
!           LOG_INFO_CONT('(1x,2I5,2F9.3)') k, k_src(k), FZ(k), FZ(k)+dist(k)
!        enddo

       if ( iq > qla ) then ! ice water
          cp = cp_ice
          cv = cv_ice
       else                 ! liquid water
          cp = cp_water
          cv = cv_water
       end if

       do k_dst = ks-1, ke-1
          do k = k_dst, k_src(k_dst)-1
             flx = rhoq(k+1,iq) * cdz(k+1) / dt               ! sum column mass rhoq*dz
             qflx(k_dst) = qflx(k_dst) - flx
             eflx(k_dst) = eflx(k_dst) - flx * temp(k+1) * cv ! internal energy flux
             dist(k_dst) = dist(k_dst) - cdz(k+1)             ! residual
          enddo
          if ( k_src(k_dst) < ke ) then
             ! residual (simple upwind)
             flx = rhoq(k_src(k_dst)+1,iq) * dist(k_dst) / dt
             qflx(k_dst) = qflx(k_dst) - flx                            ! sum column mass rhoq*dz
             eflx(k_dst) = eflx(k_dst) -flx * temp(k_src(k_dst)+1) * cv ! internal energy flux
          end if

          eflx(k_dst) = eflx(k_dst) + qflx(k) * fdz(k_dst) * grav ! potential energy
       enddo

!        LOG_INFO_CONT(*) "flux", iq
!        do k = KA, 1, -1
!           LOG_INFO_CONT('(1x,2I5,F9.3,2ES15.5)') k, k_src(k), dist(k), qflx(k), vtermh(k)*RHOQ(k+1,iq)
!        enddo

       !--- update falling tracer
       do k  = ks, ke
          rhoq(k,iq) = rhoq(k,iq) - dt * ( qflx(k) - qflx(k-1) ) * rcdz(k)
       enddo ! falling (water mass & number) tracer

!        LOG_INFO_CONT(*) "tendency", iq
!        do k = KA, 1, -1
!           LOG_INFO_CONT('(1x,I5,ES15.5)') k, - dt * ( qflx(k) - qflx(k-1) ) * RCDZ(k)
!        enddo

       ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
       if ( iq > qla + qia ) cycle

       do k = ks-1, ke-1
          mflx(k) = mflx(k) + qflx(k)
       end do

       if ( iq > qla ) then ! ice water
          sflx(2) = sflx(2) + qflx(ks-1)
       else                 ! liquid water
          sflx(1) = sflx(1) + qflx(ks-1)
       end if

       !--- update density
       do k = ks, ke
          ddens = - ( qflx(k) - qflx(k-1) ) * rcdz(k) * dt
          rhocp(k) = rhocp(k) + cp * ddens
          rhocv(k) = rhocv(k) + cv * ddens
          dens(k) = dens(k) + ddens
       end do

       !--- update internal energy
       do k = ks, ke
          rhoe(k) = rhoe(k) - ( eflx(k) - eflx(k-1) ) * rcdz(k) * dt
       end do

    end do

    do k = ks, ke
       cptot(k) = rhocp(k) / dens(k)
       cvtot(k) = rhocv(k) / dens(k)
    end do

    return
  end subroutine atmos_phy_mp_precipitation_semilag

  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine atmos_phy_mp_precipitation_momentum( &
       KA, KS, KE, &
       DENS, MOMZ, U, V, mflx, &
       RCDZ, RFDZ,             &
       MOMZ_t, RHOU_t, RHOV_t  )
    implicit none

    integer,  intent(in)  :: KA, KS, KE
    real(RP), intent(in)  :: DENS  (ka)
    real(RP), intent(in)  :: MOMZ  (ka)
    real(RP), intent(in)  :: U     (ka)
    real(RP), intent(in)  :: V     (ka)
    real(RP), intent(in)  :: mflx  (ka)
    real(RP), intent(in)  :: RCDZ  (ka)
    real(RP), intent(in)  :: RFDZ  (ka)
    real(RP), intent(out) :: MOMZ_t(ka)
    real(RP), intent(out) :: RHOU_t(ka)
    real(RP), intent(out) :: RHOV_t(ka)

    real(RP) :: flx(ka)

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

    flx(ke) = 0.0_rp

    !--- momentum z (half level)
    do k = ks, ke-1
       flx(k) = ( mflx(k) + mflx(k-1) ) * momz(k) / ( dens(k+1) + dens(k) )
    enddo
    do k = ks, ke-1
       momz_t(k) = - ( flx(k+1) - flx(k) ) * rfdz(k)
    enddo
    momz_t(ke) = 0.0_rp

    !--- momentum x
    do k = ks-1, ke-1
       flx(k) = mflx(k) * u(k+1)
    enddo
    do k = ks, ke
       rhou_t(k) = - ( flx(k) - flx(k-1) ) * rcdz(k)
    enddo

    !--- momentum y
    do k = ks-1, ke-1
       flx(k) = mflx(k) * v(k+1)
    enddo
    do k = ks, ke
       rhov_t(k) = - ( flx(k) - flx(k-1) ) * rcdz(k)
    enddo

    return
  end subroutine atmos_phy_mp_precipitation_momentum

end module scale_atmos_phy_mp_common