scale_atm_phy_mp_dgm_common Module Reference

module FElib / Atmosphere / Physics cloud microphysics / common More...

Functions/Subroutines
subroutine, public	atm_phy_mp_dgm_common_gen_intweight (intweight, lcmesh)
subroutine, public	atm_phy_mp_dgm_common_gen_vmap (vmapm, vmapp, lmesh, elem)
subroutine, public	atm_phy_mp_dgm_common_precipitation (dens, rhoq, cptot, cvtot, rhoe, flx_hydro, sflx_rain, sflx_snow, esflx, temp, vterm, dt, rnstep, dz, lift, nz, vmapm, vmapp, intweight, qha, qla, qia, lcmesh, elem)
subroutine, public	atm_phy_mp_dgm_common_precipitation_momentum (momu_t, momv_t, momz_t, dens, momu, momv, momz, mflx, dz, lift, nz, vmapm, vmapp, lcmesh, elem)
subroutine, public	atm_phy_mp_dgm_common_negative_fixer (qtrc, ddens, pres, cvtot, cptot, rtot, dens_hyd, pres_hyd, dt, lmesh, elem, qa, qla, qia, drhot)

Detailed Description

module FElib / Atmosphere / Physics cloud microphysics / common

Description

cloud microphysics process common subroutines

To preserve nonnegativity in precipitation process, a limiter proposed by Light and Durran (2016, MWR) is used

Reference

• Light and Durran 2016: Preserving Nonnegativity in Discontinuous Galerkin Approximations to Scalar Transport via Truncation and Mass Aware Rescaling (TMAR). Monthly Weather Review, 144(12), 4771–4786.

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

atm_phy_mp_dgm_common_gen_intweight()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_gen_intweight	(	real(rp), dimension(lcmesh%refelem3d%nfaces,lcmesh%refelem3d%nfptot), intent(out)	intweight,
		class(localmesh3d), target	lcmesh )

Definition at line 77 of file scale_atm_phy_mp_dgm_common.F90.

    use scale_polynominal, only: polynominal_gengausslobattoptintweight
    implicit none
    class(LocalMesh3D), target :: lcmesh
    real(RP), intent(out) :: IntWeight(lcmesh%refElem3D%Nfaces,lcmesh%refElem3D%NfpTot)
 
    class(ElementBase3D), pointer :: elem
    real(RP), allocatable :: intWeight_lgl1DPts_h(:)
    real(RP), allocatable :: intWeight_lgl1DPts_v(:)   
    real(RP), allocatable :: intWeight_h(:) 
    real(RP), allocatable :: intWeight_v(:)  
    
    integer :: f
    integer :: i, j, k, l
    integer :: is, ie
    !--------------------------------------------
 
    elem => lcmesh%refElem3D
    intweight(:,:) = 0.0_rp
 
    allocate( intweight_lgl1dpts_h(elem%Nnode_h1D) )
    allocate( intweight_lgl1dpts_v(elem%Nnode_v) )
    allocate( intweight_h(elem%Nnode_h1D*elem%Nnode_v) )
    allocate( intweight_v(elem%Nnode_h1D**2) )
 
    intweight_lgl1dpts_h(:) = polynominal_gengausslobattoptintweight(elem%PolyOrder_h)
    intweight_lgl1dpts_v(:) = polynominal_gengausslobattoptintweight(elem%PolyOrder_v)
 
    do f=1, elem%Nfaces_h
      do k=1, elem%Nnode_v
      do i=1, elem%Nnode_h1D
        l = i + (k-1)*elem%Nnode_h1D
        intweight_h(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_v(k)
      end do
      end do
 
      is = (f-1)*elem%Nfp_h + 1
      ie = is + elem%Nfp_h - 1
      intweight(f,is:ie) = intweight_h(:)
    end do
 
    do f=1, elem%Nfaces_v
      do j=1, elem%Nnode_h1D
      do i=1, elem%Nnode_h1D
        l = i + (j-1)*elem%Nnode_h1D
        intweight_v(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j)
      end do
      end do
 
      is = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v + 1
      ie = is + elem%Nfp_v - 1
      intweight(elem%Nfaces_h+f,is:ie) = intweight_v(:)
    end do
 
    return

References scale_polynominal::polynominal_gengausslobattoptintweight().

atm_phy_mp_dgm_common_gen_vmap()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_gen_vmap	(	integer, dimension(elem%nfptot,lmesh%nez), intent(out)	vmapm,
		integer, dimension(elem%nfptot,lmesh%nez), intent(out)	vmapp,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem )

Definition at line 137 of file scale_atm_phy_mp_dgm_common.F90.

    implicit none
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    integer, intent(out) :: vmapM(elem%NfpTot,lmesh%NeZ)
    integer, intent(out) :: vmapP(elem%NfpTot,lmesh%NeZ)    
 
    integer :: ke_z
    integer :: f
    integer :: vs, ve
    !-------------------------------------------------------
 
    !$omp parallel private(f, vs, ve)
    !$omp do
    do ke_z=1, lmesh%NeZ
      do f=1, elem%Nfaces_h
        vs = 1 + (f-1)*elem%Nfp_h
        ve = vs + elem%Nfp_h - 1
        vmapm(vs:ve,ke_z) = elem%Fmask_h(:,f) + (ke_z-1)*elem%Np
      end do
      do f=1, elem%Nfaces_v
        vs = elem%Nfp_h*elem%Nfaces_h + 1 + (f-1)*elem%Nfp_v
        ve = vs + elem%Nfp_v - 1
        vmapm(vs:ve,ke_z) = elem%Fmask_v(:,f) + (ke_z-1)*elem%Np
      end do
      vmapp(:,ke_z) = vmapm(:,ke_z)
    end do
    !$omp do
    do ke_z=1, lmesh%NeZ
      vs = elem%Nfp_h*elem%Nfaces_h + 1
      ve = vs + elem%Nfp_v - 1
      if (ke_z > 1) &
        vmapp(vs:ve,ke_z) = elem%Fmask_v(:,2) + (ke_z-2)*elem%Np
 
      vs = elem%Nfp_h*elem%Nfaces_h + elem%Nfp_v + 1
      ve = vs + elem%Nfp_v - 1
      if (ke_z < lmesh%NeZ) &
        vmapp(vs:ve,ke_z) = elem%Fmask_v(:,1) + ke_z*elem%Np
    end do
    !$omp end parallel
 
    return

atm_phy_mp_dgm_common_precipitation()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_precipitation	(	real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout)	dens,
		real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d,qha), intent(inout)	rhoq,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout)	cptot,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout)	cvtot,
		real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout)	rhoe,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout)	flx_hydro,
		real(rp), dimension(elem%nfp_v,lcmesh%ne2da), intent(inout)	sflx_rain,
		real(rp), dimension(elem%nfp_v,lcmesh%ne2da), intent(inout)	sflx_snow,
		real(rp), dimension (elem%nfp_v,lcmesh%ne2da), intent(inout)	esflx,
		real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d), intent(in)	temp,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d,qha), intent(in)	vterm,
		real(rp), intent(in)	dt,
		real(rp), intent(in)	rnstep,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	lift,
		real(rp), dimension(elem%nfptot,lcmesh%nez,lcmesh%ne2d), intent(in)	nz,
		integer, dimension(elem%nfptot,lcmesh%nez), intent(in)	vmapm,
		integer, dimension(elem%nfptot,lcmesh%nez), intent(in)	vmapp,
		real(rp), dimension(elem%nfaces,elem%nfptot), intent(in)	intweight,
		integer, intent(in)	qha,
		integer, intent(in)	qla,
		integer, intent(in)	qia,
		class(localmesh3d), intent(in)	lcmesh,
		class(elementbase3d), intent(in)	elem )

Parameters

[in]

qha

hydrometeor (water + ice)

Definition at line 184 of file scale_atm_phy_mp_dgm_common.F90.

    
    use scale_atmos_hydrometeor, only: &
       cv_water, &
       cp_water, &
       cv_ice,   &
       cp_ice
    
    implicit none
 
    class(LocalMesh3D), intent(in) :: lcmesh
    class(ElementBase3D), intent(in) :: elem
    integer, intent(in) :: QHA
    real(RP), intent(inout) :: DENS (elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(inout) :: RHOQ (elem%Np,lcmesh%NeZ,lcmesh%Ne2D,QHA)
    real(RP), intent(inout) :: CPtot(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(inout) :: CVtot(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(inout) :: RHOE (elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(inout) :: FLX_hydro(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(inout) :: sflx_rain(elem%Nfp_v,lcmesh%Ne2DA)
    real(RP), intent(inout) :: sflx_snow(elem%Nfp_v,lcmesh%Ne2DA)
    real(RP), intent(inout) :: esflx    (elem%Nfp_v,lcmesh%Ne2DA)
    real(RP), intent(in) :: TEMP (elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(in) :: vterm(elem%Np,lcmesh%NeZ,lcmesh%Ne2D,QHA)
    real(RP), intent(in) :: dt
    real(RP), intent(in) :: rnstep
    type(SparseMat), intent(in) :: Dz
    type(SparseMat), intent(in) :: Lift
    real(RP), intent(in) :: nz(elem%NfpTot,lcmesh%NeZ,lcmesh%Ne2D)
    integer, intent(in) :: vmapM(elem%NfpTot,lcmesh%NeZ)
    integer, intent(in) :: vmapP(elem%NfpTot,lcmesh%NeZ)
    real(RP), intent(in) :: IntWeight(elem%Nfaces,elem%NfpTot)
    integer, intent(in) :: QLA, QIA
 
    real(RP) :: qflx(elem%Np)
    real(RP) :: eflx(elem%Np)
    real(RP) :: DENS0(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: RHOCP(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: RHOCV(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: NDcoefEuler(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: DzRHOQ(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: DzRHOE(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: dDENS(elem%Np)
    real(RP) :: CP(QHA)
    real(RP) :: CV(QHA)
 
    real(RP) :: fct_coef(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: RHOQ0, RHOQ1, RHOQ_tmp(elem%Np)
    real(RP) :: netOutwardFlux(lcmesh%NeZ,lcmesh%Ne2D)
    real(RP) :: del_flux(elem%NfpTot,lcmesh%NeZ,lcmesh%Ne2D,2)
 
    real(RP) :: Fz(elem%Np), LiftDelFlx(elem%Np)
    real(RP) :: RHOQ_save(elem%Np)
 
    integer :: ke2D
    integer :: ke_z
    integer :: ke
    integer :: iq
 
    real(RP) :: Q
    real(RP) :: delz
    !-------------------------------------------------------
 
    do iq = 1, qha
      if ( iq > qla + qia ) then
        cp(iq) = undef 
        cv(iq) = undef
      else if ( iq > qla ) then ! ice water
        cp(iq) = cp_ice
        cv(iq) = cv_ice
      else                      ! liquid water
        cp(iq) = cp_water
        cv(iq) = cv_water
      end if
    end do
 
    !$omp parallel do collapse(2)
    do ke2d = 1, lcmesh%Ne2D
    do ke_z = 1, lcmesh%NeZ
      dens0(:,ke_z,ke2d) = dens(:,ke_z,ke2d)
      rhocp(:,ke_z,ke2d) = cptot(:,ke_z,ke2d) * dens(:,ke_z,ke2d)
      rhocv(:,ke_z,ke2d) = cvtot(:,ke_z,ke2d) * dens(:,ke_z,ke2d)
    end do
    end do
 
    do iq = 1, qha
      call atm_phy_mp_dgm_precipitation_get_delflux_dq( &
        del_flux(:,:,:,:),                                                      & ! (out)
        dens0(:,:,:), rhoq(:,:,:,iq), temp(:,:,:), cv(iq), nz(:,:,:), vmapm(:,:), vmapp(:,:), & ! (in)
        lcmesh, elem                                                            ) ! (in)
      
      !$omp parallel do private( &
      !$omp ke2D, ke_z, ke, delz, Fz, LiftDelFlx )
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
        delz = ( lcmesh%pos_ev(lcmesh%EToV(ke,5),3) - lcmesh%pos_ev(lcmesh%EToV(ke,1),3) ) / dble( elem%Nnode_v )
 
        call sparsemat_matmul( dz, rhoq(:,ke_z,ke2d,iq), fz )
        call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
        dzrhoq(:,ke_z,ke2d) = lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
 
        call sparsemat_matmul( dz, rhoq(:,ke_z,ke2d,iq) * cv(iq) * temp(:,ke_z,ke2d), fz )
        call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
        dzrhoe(:,ke_z,ke2d) = lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
        
        ndcoefeuler(:,ke_z,ke2d) = 0.5_rp * delz * abs(vterm(:,ke_z,ke2d,iq))
      end do
      end do
 
      call atm_phy_mp_dgm_netoutwardflux( &
        netoutwardflux(:,:),                                                  & ! (out)
        rhoq(:,:,:,iq), vterm(:,:,:,iq), dzrhoq(:,:,:), ndcoefeuler(:,:,:),   & ! (in) 
        lcmesh%J(:,:), lcmesh%Fscale(:,:),                                    & ! (in)
        nz(:,:,:), vmapm(:,:), vmapp(:,:), lcmesh%VMapM(:,:), intweight(:,:), & ! (in)
        lcmesh, elem                                                          ) ! (in)
      
      !$omp parallel do collapse(2) private(ke, Q)
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
 
        q = sum( lcmesh%J(:,ke) * elem%IntWeight_lgl(:) * rhoq(:,ke_z,ke2d,iq) ) / dt      
        fct_coef(:,ke_z,ke2d) = max( 0.0_rp, min( 1.0_rp, q / ( netoutwardflux(ke_z,ke2d) + 1.0e-10_rp ) ) )  
      end do ! end loop for ke_z
      end do ! end loop for ke2D
 
      call atm_phy_mp_dgm_precipitation_get_delflux( &
        del_flux(:,:,:,:),                                           & ! (out)
        dens0(:,:,:), rhoq(:,:,:,iq), temp(:,:,:), vterm(:,:,:,iq),  & ! (in)
        dzrhoq(:,:,:), dzrhoe(:,:,:), ndcoefeuler(:,:,:),            & ! (in)
        fct_coef(:,:,:),                                             & ! (in)
        cv(iq), lcmesh%J(:,:), lcmesh%Fscale(:,:), nz(:,:,:),        & ! (in)
        vmapm(:,:), vmapp(:,:), lcmesh%vmapM(:,:), intweight(:,:),   & ! (in)
        lcmesh, elem                                                 ) ! (in)
 
      !$omp parallel do collapse(2) private( &
      !$omp ke2D, ke_z, ke,                            &
      !$omp qflx, eflx, dDENS, RHOQ_tmp, RHOQ0, RHOQ1, &
      !$omp RHOQ_save, &
      !$omp Fz, LiftDelFlx )
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
 
        !--- update falling tracer 
 
        rhoq_save(:) = rhoq(:,ke_z,ke2d,iq)
        qflx(:) = vterm(:,ke_z,ke2d,iq) * rhoq(:,ke_z,ke2d,iq)   &
                - ndcoefeuler(:,ke_z,ke2d) * dzrhoq(:,ke_z,ke2d)
 
        call sparsemat_matmul( dz, qflx(:), fz )
        call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
        
        ddens(:) =  - dt * ( &
          lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
        rhoq_tmp(:) = max( 0.0_rp, rhoq(:,ke_z,ke2d,iq) + ddens(:) )
!        RHOQ_tmp(:) = RHOQ(:,ke_z,ke2D,iq) + dDENS(:)
 
        !
        rhoq0 = sum( lcmesh%Gsqrt(:,ke) * lcmesh%J(:,ke) * elem%IntWeight_lgl(:) * ( rhoq(:,ke_z,ke2d,iq) + ddens(:) ) )
        rhoq1 = sum( lcmesh%Gsqrt(:,ke) * lcmesh%J(:,ke) * elem%IntWeight_lgl(:) * rhoq_tmp(:)          )
 
        ddens(:) = rhoq0 / ( rhoq1 + 1.0e-32_rp ) * rhoq_tmp(:) &
                 - rhoq(:,ke_z,ke2d,iq)
        rhoq(:,ke_z,ke2d,iq) = rhoq(:,ke_z,ke2d,iq) + ddens(:)
 
 
        ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
        if ( iq > qla + qia ) cycle
 
        flx_hydro(:,ke_z,ke2d) = flx_hydro(:,ke_z,ke2d) &
                               + qflx(:) * rnstep
        if ( ke_z == 1 ) then
          if ( iq > qla ) then ! ice water
              sflx_snow(:,ke2d) = sflx_snow(:,ke2d)               &
                                + qflx(elem%Hslice(:,1)) * rnstep
          else                 ! liquid water
              sflx_rain(:,ke2d) = sflx_rain(:,ke2d)               &
                                + qflx(elem%Hslice(:,1)) * rnstep
          end if
        end if
 
        !--- update density
 
        rhocp(:,ke_z,ke2d) = rhocp(:,ke_z,ke2d) + cp(iq) * ddens(:)
        rhocv(:,ke_z,ke2d) = rhocv(:,ke_z,ke2d) + cv(iq) * ddens(:)
        dens(:,ke_z,ke2d) = dens(:,ke_z,ke2d) + ddens(:)
  
        !--- update internal energy   
 
        eflx(:) = vterm(:,ke_z,ke2d,iq) * rhoq_save(:) * temp(:,ke_z,ke2d) * cv(iq)  &
                - ndcoefeuler(:,ke_z,ke2d) * dzrhoe(:,ke_z,ke2d)
        
        call sparsemat_matmul( dz, eflx(:), fz )
        call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
 
        rhoe(:,ke_z,ke2d) = rhoe(:,ke_z,ke2d) - dt * ( &
          + lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) &
          + qflx(:) * grav                                  )
 
        if ( ke_z == 1 ) then
          esflx(:,ke2d) = esflx(:,ke2d) &
                        + eflx(elem%Hslice(:,1)) * rnstep
        end if
 
      end do ! end loop for ke_z
      end do ! end loop for ke2D
    end do ! end loop for iq
 
    !$omp parallel do collapse(2)
    do ke2d = 1, lcmesh%Ne2D
    do ke_z = 1, lcmesh%NeZ  
      cptot(:,ke_z,ke2d) = rhocp(:,ke_z,ke2d) / dens(:,ke_z,ke2d)
      cvtot(:,ke_z,ke2d) = rhocv(:,ke_z,ke2d) / dens(:,ke_z,ke2d)
    end do
    end do
 
    return

atm_phy_mp_dgm_common_precipitation_momentum()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_precipitation_momentum	(	real(rp), dimension(elem%np,lcmesh%nea), intent(out)	momu_t,
		real(rp), dimension(elem%np,lcmesh%nea), intent(out)	momv_t,
		real(rp), dimension(elem%np,lcmesh%nea), intent(out)	momz_t,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in)	dens,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in)	momu,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in)	momv,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in)	momz,
		real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in)	mflx,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	lift,
		real(rp), dimension(elem%nfptot,lcmesh%nez,lcmesh%ne2d), intent(in)	nz,
		integer, dimension(elem%nfptot,lcmesh%nez), intent(in)	vmapm,
		integer, dimension(elem%nfptot,lcmesh%nez), intent(in)	vmapp,
		class(localmesh3d), intent(in)	lcmesh,
		class(elementbase3d), intent(in)	elem )

Definition at line 411 of file scale_atm_phy_mp_dgm_common.F90.

    implicit none
 
    class(LocalMesh3D), intent(in) :: lcmesh
    class(ElementBase3D), intent(in) :: elem
    real(RP), intent(out) :: MOMU_t(elem%Np,lcmesh%NeA)
    real(RP), intent(out) :: MOMV_t(elem%Np,lcmesh%NeA)
    real(RP), intent(out) :: MOMZ_t(elem%Np,lcmesh%NeA)
    real(RP), intent(in) :: DENS(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(in) :: MOMU(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(in) :: MOMV(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(in) :: MOMZ(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    real(RP), intent(in) :: mflx(elem%Np,lcmesh%NeZ,lcmesh%Ne2D)
    type(SparseMat), intent(in) :: Dz
    type(SparseMat), intent(in) :: Lift
    real(RP), intent(in) :: nz(elem%NfpTot,lcmesh%NeZ,lcmesh%Ne2D)
    integer, intent(in) :: vmapM(elem%NfpTot,lcmesh%NeZ)
    integer, intent(in) :: vmapP(elem%NfpTot,lcmesh%NeZ)
    
    integer :: ke2D
    integer :: ke_z
    integer :: ke
 
    real(RP) :: Fz(elem%Np), LiftDelFlx(elem%Np)
    real(RP) :: del_flux(elem%NfpTot,lcmesh%NeZ,lcmesh%Ne2D,3)
 
    real(RP) :: RDENS(elem%Np)
    !-------------------------------------------------------
 
    call atm_phy_mp_dgm_precipitation_momentum_get_delflux( &
      del_flux(:,:,:,:),                                  & ! (out)
      dens(:,:,:), momu(:,:,:), momv(:,:,:), momz(:,:,:), & ! (in)
      mflx(:,:,:),                                        & ! (in)
      nz(:,:,:), vmapm(:,:), vmapp(:,:),                  & ! (in)
      lcmesh, elem                                        ) ! (in)
 
    !$omp parallel do collapse(2) private( &
    !$omp ke2D, ke_z, ke, RDENS, Fz, LiftDelFlx )
    do ke2d = 1, lcmesh%Ne2D
    do ke_z = 1, lcmesh%NeZ
      ke = ke2d + (ke_z-1)*lcmesh%Ne2D
      rdens(:) = 1.0_rp / dens(:,ke_z,ke2d)
 
      call sparsemat_matmul( dz, mflx(:,ke_z,ke2d) * momu(:,ke_z,ke2d) * rdens(:), fz )
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
      momu_t(:,ke) = - ( lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
 
      call sparsemat_matmul( dz, mflx(:,ke_z,ke2d) * momv(:,ke_z,ke2d) * rdens(:), fz )
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
      momv_t(:,ke) = - ( lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
 
      call sparsemat_matmul( dz, mflx(:,ke_z,ke2d) * momz(:,ke_z,ke2d) * rdens(:), fz )
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,3), liftdelflx )
      momz_t(:,ke) = - ( lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
    end do
    end do
 
    return

atm_phy_mp_dgm_common_negative_fixer()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_negative_fixer	(	type(localmeshfieldbaselist), dimension(qa), intent(inout)	qtrc,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	ddens,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	pres,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	cvtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	cptot,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd,
		real(rp), intent(in)	dt,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		integer, intent(in)	qa,
		integer, intent(in)	qla,
		integer, intent(in)	qia,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout), optional	drhot )

Definition at line 477 of file scale_atm_phy_mp_dgm_common.F90.

 
    use scale_const, only: &
      cvdry => const_cvdry,  &
      cpdry => const_cpdry,  &
      rdry => const_rdry
 
    use scale_tracer, only: &
      tracer_mass, tracer_r, tracer_cv, tracer_cp
    use scale_atmos_thermodyn, only: &
      atmos_thermodyn_specific_heat
    use scale_localmeshfield_base, only: localmeshfieldbaselist
    implicit none
 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    integer, intent(in) :: QA
    type(LocalMeshFieldBaseList), intent(inout) :: QTRC(QA)
    real(RP), intent(inout) :: DDENS(elem%Np,lmesh%NeA)
    real(RP), intent(inout) :: PRES(elem%Np,lmesh%NeA)
    real(RP), intent(inout) :: CVtot(elem%Np,lmesh%NeA)
    real(RP), intent(inout) :: CPtot(elem%Np,lmesh%NeA)
    real(RP), intent(inout) :: Rtot(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: dt
    integer, intent(in) :: QLA, QIA
    real(RP), intent(inout), optional :: DRHOT(elem%Np,lmesh%NeA)
 
    integer :: ke
    integer :: iq
 
    real(RP) ::  int_w(elem%Np)
 
    real(RP) :: DENS(elem%Np)
    real(RP) :: DDENS0(elem%Np)
 
    real(RP) :: TRCMASS0(elem%Np), TRCMASS1(elem%Np,QA)
    real(RP) :: MASS0_elem, MASS1_elem
    real(RP) :: IntEn0_elem, IntEn_elem
 
    real(RP) :: QTRC_tmp(elem%Np,QA), Qdry(elem%Np)
    real(RP) :: CVtot_old(elem%Np), CPtot_old(elem%Np), Rtot_old(elem%Np)
    real(RP) :: InternalEn(elem%Np), InternalEn0(elem%Np), TEMP(elem%Np)
    real(RP) :: RHOT_hyd(elem%Np)
 
#ifdef SINGLE
    real(RP), parameter :: TRC_EPS = 1e-32_rp
#else
    real(RP), parameter :: TRC_EPS = 1e-128_rp
#endif
    !------------------------------------------------
 
    !$omp parallel do private( &
    !$omp ke, iq, DENS, DDENS0, InternalEn, InternalEn0, TEMP, QTRC_tmp,       &
    !$omp TRCMASS0, TRCMASS1, MASS0_elem, MASS1_elem, IntEn0_elem, IntEn_elem, &
    !$omp Qdry, CVtot_old, CPtot_old, Rtot_old,                                &
    !$omp int_w, RHOT_hyd )
    do ke = lmesh%NeS, lmesh%NeE
 
      do iq = 1, qa
        qtrc_tmp(:,iq) = qtrc(iq)%ptr%val(:,ke)
      end do
      call atmos_thermodyn_specific_heat( & 
        elem%Np, 1, elem%Np, qa,                                           & ! (in)
        qtrc_tmp, tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:), & ! (in)
        qdry, rtot_old, cvtot_old, cptot_old                               ) ! (out)
 
      dens(:) = dens_hyd(:,ke) + ddens(:,ke)
      ddens0(:) = ddens(:,ke)
      
      ! RHOT_hyd(:) = PRES00 / Rdry * ( PRES_hyd(:,ke) / PRES00 )**( CVdry / CPdry )
      ! ( Internal energy ) = Cvtot * RHO * T = CVtot * RHO * ( PT * EXNER )      
      ! InternalEn0(:) = CVtot_old(:) * ( RHOT_hyd(:) + DRHOT(:,ke) ) &
      !                * ( Rtot_old(:) * ( RHOT_hyd(:) + DRHOT(:,ke) ) / PRES00 )**( Rtot_old(:) / CVtot_old(:) ) 
      !InternalEn0(:) = CVtot_old(:) * PRES(:,ke) / Rtot_old(:)
      internalen0(:) = cvtot(:,ke) * pres(:,ke) / rtot(:,ke)
      !TEMP(:) = InternalEn0(:) / ( DENS(:) * CVtot_old(:) )
      temp(:) = internalen0(:) / ( dens(:) * cvtot(:,ke) )
      internalen(:) = internalen0(:)
 
      int_w(:) = lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:)
      do iq = 1, 1 + qla + qia  
        trcmass0(:) = dens(:) * qtrc_tmp(:,iq)
        trcmass1(:,iq) = max( trc_eps, trcmass0(:) )
 
        mass0_elem = sum( int_w(:) * trcmass0(:)    )
        mass1_elem = sum( int_w(:) * trcmass1(:,iq) )
        trcmass1(:,iq) = max(mass0_elem, 0.0e0_rp) / mass1_elem * trcmass1(:,iq)
!        TRCMASS1(:,iq) = MASS0_elem / MASS1_elem * TRCMASS1(:,iq)
 
        ddens(:,ke) = ddens(:,ke) + ( trcmass1(:,iq) - trcmass0(:) )
        internalen(:) = internalen(:) + ( trcmass1(:,iq) - trcmass0(:) ) * tracer_cv(iq) * temp(:)
      end do
 
      !--
 
      dens(:) = dens_hyd(:,ke) + ddens(:,ke)
      do iq = 1, qa
        qtrc_tmp(:,iq) = trcmass1(:,iq) / dens(:)
        qtrc(iq)%ptr%val(:,ke) = qtrc_tmp(:,iq)
      end do
 
      inten0_elem = sum( int_w(:) * internalen0(:) )
      inten_elem  = sum( int_w(:) * internalen(:) )
      internalen(:) = inten0_elem / inten_elem * internalen(:)
 
      call atmos_thermodyn_specific_heat( &
        elem%Np, 1, elem%Np, qa,                                           & ! (in)
        qtrc_tmp, tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:), & ! (in)
        qdry, rtot(:,ke), cvtot(:,ke), cptot(:,ke)                         ) ! (out)
 
      internalen(:) = internalen(:) - ( ddens(:,ke) - ddens0(:) ) * grav * lmesh%zlev(:,ke)
      pres(:,ke) = internalen(:) * rtot(:,ke) / cvtot(:,ke)
 
      if ( present(drhot) ) then
        drhot(:,ke) = pres00 / rtot(:,ke) * ( pres(:,ke) / pres00 )**( cvtot(:,ke) / cptot(:,ke) ) &
                    - pres00 / rdry * ( pres_hyd(:,ke) / pres00 )**( cvdry / cpdry )
      end if
    end do
 
    return