scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common_2.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVI / Common
!!
!! @par Description
!!      HEVI DGM scheme for Atmospheric dynamical process
!!      As the thermodynamics equation, a density * potential temperature equation is used.
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common_2
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    grav => const_grav,  &
    rdry => const_rdry,  &
    cpdry => const_cpdry, &
    cvdry => const_cvdry, &
    pres00 => const_pre00
 
  use scale_sparsemat  
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_element_modalfilter, only: modalfilter
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  use scale_element_operation_base, only: elementoperationbase3d
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
    momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id,   &
    momz_vid => prgvar_momz_id,                               &
    prgvar_num
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
  
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd
 
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve_uv
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  logical, public :: vi_use_lapack_flag = .true.
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  private :: vi_cal_del_flux_dyn
  private :: vi_cal_del_flux_dyn_uv
 
contains
  !------------------------------------------------
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
    implicit none
    namelist / param_atmos_dyn_nonhydro3d_rhot_hevi_common / &
      vi_use_lapack_flag
 
    integer :: ierr
    !---------------------------------------------------------
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_dyn_nonhydro3d_rhot_hevi_common,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_DYN_nonhydro3d_rhot_hevi_common_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_DYN_nonhydro3d_rhot_hevi_common_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_NONHYDRO3D_RHOT_HEVI_COMMON. Check!'
      call prc_abort
    endif
    log_nml(param_atmos_dyn_nonhydro3d_rhot_hevi_common)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
    implicit none
    !---------------------------------------------------------
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
 
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax( &
    DENS_t, MOMZ_t, RHOT_t,                                  & ! (out)
    alph,                                                    & ! (in)
    prog_vars, prog_vars0,                                   & ! (in)
    ddens00, momx00, momy00, momz00, drhot00,                & ! (in)
    dens_hyd, pres_hyd,                                      & ! (in)
    rtot, cptot, cvtot, gsqrtv,                              & ! (in)
    impl_fac, dt,                                            & ! (in)
    lmesh, elem,                                             & ! (in)
    vmapm, vmapp,                                            & ! (in)
    element3d_operation,                                     & ! (in)
    im, jm,                                                  & ! (in)
    b,                                                       & ! (out)
    dens, w, wt, pot, dpdrhot                                ) ! (out)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: im, jm
    real(rp), intent(out) :: dens_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: momz_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: rhot_t(elem%np,lmesh%nea)
    real(rp), intent(in) :: alph(elem%nfptot,lmesh%nex*lmesh%ney,lmesh%nez)
    real(rp), intent(in)  :: prog_vars  (elem%np,lmesh%nex*lmesh%ney*lmesh%nez,prgvar_num)
    real(rp), intent(in)  :: prog_vars0 (elem%np,lmesh%nex*lmesh%ney*lmesh%nez,prgvar_num)
    real(rp), intent(in)  :: ddens00(elem%np,lmesh%nea)
    real(rp), intent(in)  :: momx00 (elem%np,lmesh%nea)
    real(rp), intent(in)  :: momy00 (elem%np,lmesh%nea)
    real(rp), intent(in)  :: momz00 (elem%np,lmesh%nea)
    real(rp), intent(in)  :: drhot00(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)  :: rtot(elem%np,lmesh%nea)
    real(rp), intent(in)  :: cptot(elem%np,lmesh%nea)
    real(rp), intent(in)  :: cvtot(elem%np,lmesh%nea)
    real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
    real(rp), intent(in) :: impl_fac
    real(rp), intent(in) :: dt
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)    
    class(elementoperationbase3d), intent(in) :: element3d_operation
    real(rp), intent(out), optional :: b(im,3,elem%nnode_v,jm,lmesh%ne)
    real(rp), intent(out), optional :: dens(im,elem%nnode_v,jm,lmesh%ne)
    real(rp), intent(out), optional :: w(im,elem%nnode_v,jm,lmesh%ne)
    real(rp), intent(out), optional :: wt(im,elem%nnode_v,jm,lmesh%ne)
    real(rp), intent(out), optional :: pot(im,elem%nnode_v,jm,lmesh%ne)
    real(rp), intent(out), optional :: dpdrhot(im,elem%nnode_v,jm,lmesh%ne)
 
    real(rp) :: flux(elem%np,3), dflux(elem%np,2,3)
    real(rp) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
    real(rp) :: dpres(elem%np), rhot(elem%np)
 
    real(rp) :: rdens_, u_, v_, w_, pt_
    real(rp) :: rgsqrtv(elem%np), rgsqrt(elem%np)
    real(rp) :: gsqrt_, gsqrtdpres_, e33
 
    integer :: ke_xy, ke_z
    integer :: ke, ke2d
    real(rp) :: rdt
    real(rp) :: drho(elem%np)
    
    real(rp) :: gamm, rgamm    
    real(rp) :: rp0
    real(rp) :: rovp0, p0ovr
 
    integer :: i, j, p, pp, pv
 
    real(rp) :: cptot_ov_cvtot
 
    logical :: flag_cal_b
    !--------------------------------------------------------
 
    gamm  = cpdry / cvdry
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    rovp0 = rdry * rp0
    p0ovr = pres00 / rdry
 
    rdt = 1.0_rp / dt
 
    if ( present(b) ) then
      flag_cal_b = .true.
    else
      flag_cal_b = .false.
    end if
 
    call vi_cal_del_flux_dyn( del_flux, & ! (out)
      alph, prog_vars, prog_vars0,                                              & ! (in)
      dens_hyd, pres_hyd, rtot, cptot, cvtot,                                   & ! (in)
      lmesh%Gsqrt, lmesh%GsqrtH, lmesh%gam, gsqrtv, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)    
      lmesh%normal_fn(:,:,3), vmapm, vmapp, elem%IndexH2Dto3D_bnd,              & ! (in)
      lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D                     ) ! (in)
 
    !$omp parallel private( ke_xy, ke_z, ke, ke2d, p, i,j,pv, &
    !$omp DPRES, RHOT, drho, pt_, Flux, DFlux,                &
    !$omp E33, RGsqrtV, CPtot_ov_CVtot                        )
 
    !$omp do collapse(2)
    do ke_z=1, lmesh%NeZ
    do ke_xy=1, lmesh%NeX*lmesh%NeY
      ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
      ke2d = lmesh%EMap3Dto2D(ke)
 
      do p=1, elem%Np
        rgsqrtv(p) = 1.0_rp / gsqrtv(p,ke)
      end do
 
      !-
      do p=1, elem%Np
        flux(p,dens_vid) = prog_vars(p,ke,momz_vid) &
                         + gsqrtv(p,ke) * lmesh%GI3(p,ke,1) * prog_vars(p,ke,momx_vid) &
                         + gsqrtv(p,ke) * lmesh%GI3(p,ke,2) * prog_vars(p,ke,momy_vid)
      end do
 
      rhot(:) = p0ovr * (pres_hyd(:,ke)/pres00)**rgamm &
              + prog_vars(:,ke,rhot_vid)
      do p=1, elem%Np
        pt_ = rhot(p) / ( prog_vars(p,ke,dens_vid) + dens_hyd(p,ke) )
        flux(p,rhot_vid) = pt_ * flux(p,dens_vid)
      end do
 
      do p=1, elem%Np
        flux(p,momz_vid) = pres00 * ( rtot(p,ke) * rp0 * rhot(p) )**(cptot(p,ke) / cvtot(p,ke)) &
                         - pres_hyd(p,ke)
      end do
 
      !-
      call element3d_operation%Dz( flux(:,dens_vid), dflux(:,1,dens_vid) )
      call element3d_operation%Lift( del_flux(:,dens_vid,ke), dflux(:,2,dens_vid)  )
 
      call element3d_operation%Dz( flux(:,rhot_vid), dflux(:,1,rhot_vid) )
      call element3d_operation%Lift( del_flux(:,rhot_vid,ke), dflux(:,2,rhot_vid)  )
 
      call element3d_operation%Dz( flux(:,momz_vid), dflux(:,1,momz_vid) )
      call element3d_operation%Lift( del_flux(:,momz_vid,ke), dflux(:,2,momz_vid)  )
 
      !-
      do p=1, elem%Np
        e33 = lmesh%Escale(p,ke,3,3)
 
        dens_t(p,ke) = - ( &
              e33 * dflux(p,1,dens_vid) &
            + dflux(p,2,dens_vid)       ) * rgsqrtv(p)
          
        rhot_t(p,ke) = - ( &
              e33 * dflux(p,1,rhot_vid) &
            + dflux(p,2,rhot_vid)       ) * rgsqrtv(p)
      end do
 
      call element3d_operation%VFilterPM1( prog_vars(:,ke,dens_vid), &
        drho )
      do p=1, elem%Np
        e33 = lmesh%Escale(p,ke,3,3)
 
        momz_t(p,ke) = - ( &
              e33 * dflux(p,1,momz_vid) &
            + dflux(p,2,momz_vid)       ) * rgsqrtv(p) &
            - grav * drho(p) 
      end do
 
      if ( flag_cal_b ) then
        do j=1, jm
        do pv=1, elem%Nnode_v
        do i=1, im
          p = i + (j-1)*im + (pv-1)*im*jm
          dens(i,pv,j,ke) = dens_hyd(p,ke) + prog_vars(p,ke,dens_vid)
          pot(i,pv,j,ke) = rhot(p) / dens(i,pv,j,ke)
          w(i,pv,j,ke) = prog_vars(p,ke,momz_vid) / dens(i,pv,j,ke)
          wt(i,pv,j,ke) = flux(p,dens_vid) / dens(i,pv,j,ke)
        end do
        end do
        end do
        do j=1, jm
        do pv=1, elem%Nnode_v
        do i=1, im
          p = i + (j-1)*im + (pv-1)*im*jm
          cptot_ov_cvtot = cptot(p,ke) / cvtot(p,ke)
          dpdrhot(i,pv,j,ke) = cptot_ov_cvtot &
                            * pres00 * ( rtot(p,ke) / pres00 * rhot(p) )**cptot_ov_cvtot &
                            / rhot(p)
        end do
        end do
        end do
      end if
 
    end do
    end do
    !$omp end do
 
    if ( flag_cal_b ) then
      !$omp do collapse(2)
      do ke=1, lmesh%Ne
      do pv=1, elem%Nnode_v
      do j=1, jm
      do i=1, im
        p  = i + (j-1)*im + (pv-1)*im*jm
        b(i,1,pv,j,ke) = impl_fac * dens_t(p,ke)         &
                       - prog_vars(p,ke,dens_vid)      &
                       + ddens00(p,ke)
        b(i,2,pv,j,ke) = impl_fac * momz_t(p,ke)         &
                       - prog_vars(p,ke,momz_vid)      &
                       + momz00(p,ke)
        b(i,3,pv,j,ke) = impl_fac * rhot_t(p,ke)         &
                       - prog_vars(p,ke,rhot_vid)      &
                       + drhot00(p,ke)
      end do
      end do
      end do
      end do
      !$omp end do
    end if
    !$omp end parallel
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve( PROG_VARS, &
    b, DENS, W, WT, POT, DPDRHOT, alph,                          &
    GsqrtV, IntrpMat_VPOrdM1,                                    &
    impl_fac, dt,                                                &
    vmapM, vmapP, lmesh, elem, im, jm                            )
    use scale_atm_dyn_dgm_hevi_common_linalgebra, only: &
      linkernel_solve_var3 => atm_dyn_dgm_hevi_common_linalgebra_solve_var3
 
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: im
    integer, intent(in) :: jm
    real(rp), intent(inout) :: prog_vars (elem%nnode_h1d**2,elem%nnode_v,lmesh%nex*lmesh%ney,lmesh%nez,prgvar_num)
    real(rp), intent(inout) :: b(im,3*elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: dens(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: w(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: wt(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: pot(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: dpdrhot(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
    real(rp), intent(in) :: intrpmat_vpordm1(elem%np,elem%np)
    real(rp), intent(in) :: impl_fac
    real(rp), intent(in) :: dt
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
 
    real(rp) :: tmp_v3(3)
    real(rp) :: bndmatd(im,3*elem%nnode_v,3*elem%nnode_v,jm,lmesh%ne2d)
    real(rp) :: bndmatl(im,3*elem%nnode_v,3,jm,lmesh%ne2d)
    real(rp) :: g(im,3*elem%nnode_v,3,jm,lmesh%ne2d,lmesh%nez)
    integer :: ipiv(im,3*elem%nnode_v)
 
    integer :: colmask(elem%nnode_v)
    real(rp) :: intrpmat_vpordm1_(elem%nnode_v,elem%nnode_v)
    real(rp) :: dx3(elem%nnode_v,elem%nnode_v)
    real(rp) :: id(elem%nnode_v,elem%nnode_v)
 
    integer :: ke_xy, ke_z
    integer :: i, j, ij, pv, pv2
    integer :: v
    integer :: pp
    integer :: p, p1, p2, p3
    !-----------------------------------------------------------------
 
    ij = 1
    colmask(:) = elem%Colmask(:,ij)
    do pv2=1, elem%Nnode_v
    do pv=1, elem%Nnode_v
      p = ij + (pv-1)*elem%Nnode_h1D**2
      dx3(pv,pv2) = impl_fac * elem%Dx3(p,colmask(pv2))
      intrpmat_vpordm1_(pv,pv2) = impl_fac * grav * intrpmat_vpordm1(p,colmask(pv2))
    end do
    end do
 
    id(:,:) = 0.0_rp
    do pv=1, elem%Nnode_v
      id(pv,pv) = 1.0_rp
    end do
 
    do ke_z=1,lmesh%NeZ
      call atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd( &
        bndmatl, bndmatd, g(:,:,:,:,:,ke_z),                         & ! (out)
        dens, w, wt, pot, dpdrhot, alph,                             &
        gsqrtv, lmesh%normal_fn(:,:,3), id, dx3, intrpmat_vpordm1_,  &
        impl_fac, dt, lmesh, elem, im, jm, ke_z )
 
      call prof_rapstart( 'hevi_cal_vi_lin', 3)        
!      call PROF_rapstart( 'hevi_cal_vi_lin_pre', 3)          
      if ( ke_z > 1 ) then        
        !$omp parallel do private( ke_xy,j,pv,i,v, p1,p2,p3,tmp_v3 ) collapse(2)
        do ke_xy=1, lmesh%NeX * lmesh%NeY
        do j=1, jm
          do pv=1, 3*elem%Nnode_v
          do i=1, im
            p1 = 1 + (elem%Nnode_v-1)*3
            p2 = p1 + 1
            p3 = p1 + 2
 
            tmp_v3(:) = bndmatl(i,pv,:,j,ke_xy)
            bndmatd(i,pv,1,j,ke_xy) = bndmatd(i,pv,1,j,ke_xy) &
                              - tmp_v3(1) * g(i,p1,1,j,ke_xy,ke_z-1) &
                              - tmp_v3(2) * g(i,p2,1,j,ke_xy,ke_z-1) &
                              - tmp_v3(3) * g(i,p3,1,j,ke_xy,ke_z-1)
            bndmatd(i,pv,2,j,ke_xy) = bndmatd(i,pv,2,j,ke_xy) &
                              - tmp_v3(1) * g(i,p1,2,j,ke_xy,ke_z-1) &
                              - tmp_v3(2) * g(i,p2,2,j,ke_xy,ke_z-1) &
                              - tmp_v3(3) * g(i,p3,2,j,ke_xy,ke_z-1)
            bndmatd(i,pv,3,j,ke_xy) = bndmatd(i,pv,3,j,ke_xy) &
                              - tmp_v3(1) * g(i,p1,3,j,ke_xy,ke_z-1) &
                              - tmp_v3(2) * g(i,p2,3,j,ke_xy,ke_z-1) &
                              - tmp_v3(3) * g(i,p3,3,j,ke_xy,ke_z-1)
            b(i,pv,j,ke_xy,ke_z) = b(i,pv,j,ke_xy,ke_z) &
              - tmp_v3(1) * b(i,p1,j,ke_xy,ke_z-1) &
              - tmp_v3(2) * b(i,p2,j,ke_xy,ke_z-1) &
              - tmp_v3(3) * b(i,p3,j,ke_xy,ke_z-1)              
          end do
          end do
        end do
        end do
      end if
!      call PROF_rapend( 'hevi_cal_vi_lin_pre', 3)
 
      call prof_rapstart( 'hevi_cal_vi_lin_core', 3)
      call linkernel_solve_var3( bndmatd, b(:,:,:,:,ke_z), g(:,:,:,:,:,ke_z), &
        elem%Nnode_v, im, jm, lmesh%Ne2D, ke_z == lmesh%NeZ )
 
      call prof_rapend( 'hevi_cal_vi_lin_core', 3)
      call prof_rapend( 'hevi_cal_vi_lin', 3)
    end do
 
    call prof_rapstart( 'hevi_cal_vi_lin', 3)
!    call PROF_rapstart( 'hevi_cal_vi_lin_post', 3)
    do ke_z=lmesh%NeZ-1, 1, -1
      !$omp parallel do private( ke_xy,j,pv,i, tmp_v3 ) collapse(2)          
      do ke_xy=1, lmesh%NeX * lmesh%NeY
      do j=1, jm
        do pv=1, 3*elem%Nnode_v
        do i=1, im
          tmp_v3(:) = g(i,pv,:,j,ke_xy,ke_z)
          b(i,pv,j,ke_xy,ke_z) = b(i,pv,j,ke_xy,ke_z) &
            - tmp_v3(1) * b(i,1,j,ke_xy,ke_z+1) &
            - tmp_v3(2) * b(i,2,j,ke_xy,ke_z+1) &
            - tmp_v3(3) * b(i,3,j,ke_xy,ke_z+1)
        end do
        end do
      end do
      end do
    end do
    !$omp parallel do collapse(3) private(pv,i,j,p1,pp)
    do ke_z=1, lmesh%NeZ
    do ke_xy=1, lmesh%NeX * lmesh%NeY
    do j=1, jm
      do pv=1, elem%Nnode_v
      do i=1, im
        p1 = 1 + (pv-1)*3
        pp = i + (j-1)*im
        prog_vars(pp,pv,ke_xy,ke_z,dens_vid) = prog_vars(pp,pv,ke_xy,ke_z,dens_vid) + b(i,p1  ,j,ke_xy,ke_z)
        prog_vars(pp,pv,ke_xy,ke_z,momz_vid) = prog_vars(pp,pv,ke_xy,ke_z,momz_vid) + b(i,p1+1,j,ke_xy,ke_z)
        prog_vars(pp,pv,ke_xy,ke_z,rhot_vid) = prog_vars(pp,pv,ke_xy,ke_z,rhot_vid) + b(i,p1+2,j,ke_xy,ke_z)
      end do
      end do
    end do
    end do
    end do
!    call PROF_rapend( 'hevi_cal_vi_lin_post', 3)
    call prof_rapend( 'hevi_cal_vi_lin', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve
 
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv( &
    MOMX_t, MOMY_t,                                          & ! (out)
    alph,                                                    & ! (out)
    prog_vars, prog_vars0,                                   & ! (in)
    ddens00, momx00, momy00, momz00, drhot00,                & ! (in)
    dens_hyd, pres_hyd,                                      & ! (in)
    rtot, cptot, cvtot, gsqrtv,                              & ! (in)
    impl_fac, dt,                                            & ! (in)
    lmesh, elem,                                             & ! (in)
    vmapm, vmapp,                                            & ! (in)
    element3d_operation,                                     & ! (in)
    im, jm,                                                  & ! (in)
    b1d_ij_uv                                                ) ! (out)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: momx_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: momy_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in)  :: prog_vars  (elem%np,lmesh%ne,prgvar_num)
    real(rp), intent(in)  :: prog_vars0 (elem%np,lmesh%ne,prgvar_num)
    real(rp), intent(in)  :: ddens00(elem%np,lmesh%nea)
    real(rp), intent(in)  :: momx00 (elem%np,lmesh%nea)
    real(rp), intent(in)  :: momy00 (elem%np,lmesh%nea)
    real(rp), intent(in)  :: momz00 (elem%np,lmesh%nea)
    real(rp), intent(in)  :: drhot00(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)  :: rtot(elem%np,lmesh%nea)
    real(rp), intent(in)  :: cptot(elem%np,lmesh%nea)
    real(rp), intent(in)  :: cvtot(elem%np,lmesh%nea)
    real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
    real(rp), intent(in) :: impl_fac
    real(rp), intent(in) :: dt
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)   
    class(elementoperationbase3d), intent(in) :: element3d_operation
    integer, intent(in) :: im, jm
    real(rp), intent(out), optional :: b1d_ij_uv(im*elem%nnode_v,2,jm,lmesh%ne)
 
    real(rp) :: rgsqrtv
    real(rp) :: liftdelflx(elem%np,2)
    real(rp) :: del_flux(elem%nfptot,2,lmesh%ne)
    integer :: ke_xy, ke_z
    integer :: ke, ke2d
    integer :: fp
 
    integer :: i, j, p, pp, pv
    !--------------------------------------------------------
 
    call vi_cal_del_flux_dyn_uv( del_flux, alph, & ! (out)
      prog_vars, prog_vars0,                                                   & ! (in)
      dens_hyd, pres_hyd, rtot, cptot, cvtot,                                  & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%GsqrtH, lmesh%gam, gsqrtv, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),     & ! (in)    
      lmesh%normal_fn(:,:,3), vmapm, vmapp, elem%IndexH2Dto3D_bnd,             & ! (in)
      lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D                    ) ! (in)
 
    !$omp parallel private( ke_xy, ke_z, ke, ke2D, fp, pv, i, j, p, pp, &
    !$omp LiftDelFlx, RGsqrtV                         )
 
    !$omp do collapse(2)
    do ke_z=1, lmesh%NeZ
    do ke_xy=1, lmesh%NeX*lmesh%NeY
      ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
      ke2d = lmesh%EMap3Dto2D(ke)
 
      call element3d_operation%Lift( del_flux(:,1,ke), liftdelflx(:,1)  )
      call element3d_operation%Lift( del_flux(:,2,ke), liftdelflx(:,2)  )
 
      do fp=1, elem%Np
        rgsqrtv = 1.0_rp / gsqrtv(fp,ke)
        momx_t(fp,ke) = - liftdelflx(fp,1) * rgsqrtv
        momy_t(fp,ke) = - liftdelflx(fp,2) * rgsqrtv
      end do
    end do
    end do
    !$omp end do
 
    if ( present( b1d_ij_uv ) ) then
      !$omp do collapse(2)
      do ke=1, lmesh%Ne
      do pv=1, elem%Nnode_v
      do j=1, jm
      do i=1, im
        pp = i + (pv-1)*im
        p  = i + (j-1)*im + (pv-1)*im*jm
        b1d_ij_uv(pp,1,j,ke) = impl_fac * momx_t(p,ke)    &
                             - prog_vars(p,ke,momx_vid) &
                             + momx00(p,ke)
        b1d_ij_uv(pp,2,j,ke) = impl_fac * momy_t(p,ke)    &
                             - prog_vars(p,ke,momy_vid) &
                             + momy00(p,ke)
      end do
      end do
      end do
      end do
      !$omp end do
    end if
    !$omp end parallel
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve_uv( PROG_VARS, &
    b_uv, alph,                                                  &
    GsqrtV,                                                      &
    impl_fac, dt,                                                &
    vmapM, vmapP, lmesh, elem, im, jm                            )
    use scale_atm_dyn_dgm_hevi_common_linalgebra, only: &
      linkernel_solve_uv => atm_dyn_dgm_hevi_common_linalgebra_solve_uv
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: im
    integer, intent(in) :: jm
    real(rp), intent(inout) :: prog_vars(elem%nnode_h1d**2,elem%nnode_v,lmesh%nex*lmesh%ney,lmesh%nez,prgvar_num)
    real(rp), intent(inout) :: b_uv(im*elem%nnode_v,2,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
    real(rp), intent(in) :: impl_fac
    real(rp), intent(in) :: dt
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
 
    real(rp) :: tmp
    real(rp) :: bndmatd_uv(im*elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d)
    real(rp) :: bndmatl_uv(im*elem%nnode_v,jm,lmesh%ne2d)
    real(rp) :: g_uv(im*elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    integer :: ipiv_uv(im*elem%nnode_v)
 
    integer :: ke_xy, ke_z
    integer :: i, pv, j
    integer :: v
    integer :: p, pp
    integer :: p1, p2, p3
    !-----------------------------------------------------------------
 
    do ke_z=1, lmesh%NeZ
      call atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv( &
        bndmatl_uv, bndmatd_uv,  g_uv(:,:,:,ke_z),                     & ! (out)
        alph, gsqrtv, impl_fac, dt, lmesh, elem, im, jm, vmapm, vmapp, & ! (in)
        ke_z ) ! (in)
 
      call prof_rapstart( 'hevi_cal_vi_lin_uv', 3)        
      if ( ke_z > 1 ) then
        !$omp parallel do private( ke_xy,j,pv,i, pp,p2,tmp ) collapse(2)
        do ke_xy=1, lmesh%NeX * lmesh%NeY
        do j=1, jm
            do pv=1, elem%Nnode_v
            do i=1, im
              pp = i + (pv-1)*im
              p2 = i + (elem%Nnode_v-1)*im
              tmp = bndmatl_uv(pp,j,ke_xy)
              bndmatd_uv(pp,1,j,ke_xy) = bndmatd_uv(pp,1,j,ke_xy) - tmp * g_uv(p2,j,ke_xy,ke_z-1)
              b_uv(pp,1,j,ke_xy,ke_z) = b_uv(pp,1,j,ke_xy,ke_z) - tmp * b_uv(p2,1,j,ke_xy,ke_z-1)
              b_uv(pp,2,j,ke_xy,ke_z) = b_uv(pp,2,j,ke_xy,ke_z) - tmp * b_uv(p2,2,j,ke_xy,ke_z-1)
            end do
            end do
        end do
        end do
      end if
      
      call linkernel_solve_uv( bndmatd_uv, b_uv(:,:,:,:,ke_z), g_uv(:,:,:,ke_z), &
        elem%Nnode_v, im, jm, lmesh%Ne2D, ke_z==lmesh%NeZ )
 
      call prof_rapend( 'hevi_cal_vi_lin_uv', 3) 
    end do
    call prof_rapstart( 'hevi_cal_vi_lin_uv', 3)
    do ke_z=lmesh%NeZ-1, 1, -1
      !$omp parallel do private( ke_xy,j,pv,i, pp,p2,tmp ) collapse(2)
      do ke_xy=1, lmesh%NeX * lmesh%NeY
      do j=1, jm
        do pv=1, elem%Nnode_v
        do i=1, im
          pp = i + (pv-1)*im
          p2 = i ! + (1-1)*im        
          tmp = g_uv(pp,j,ke_xy,ke_z)
          b_uv(pp,1,j,ke_xy,ke_z) = b_uv(pp,1,j,ke_xy,ke_z) - tmp * b_uv(p2,1,j,ke_xy,ke_z+1)
          b_uv(pp,2,j,ke_xy,ke_z) = b_uv(pp,2,j,ke_xy,ke_z) - tmp * b_uv(p2,2,j,ke_xy,ke_z+1)
        end do
        end do
      end do
      end do
    end do
    !$omp parallel do private(ke_z,ke_xy,j,pv,i,pp,p) collapse(2)
    do ke_z=1, lmesh%NeZ
    do ke_xy=1, lmesh%NeX * lmesh%NeY
      do j=1, jm
        do pv=1, elem%Nnode_v
        do i=1, im    
          p = i + (pv-1)*im
          pp = i + (j-1)*im
          prog_vars(pp,pv,ke_xy,ke_z,momx_vid) = prog_vars(pp,pv,ke_xy,ke_z,momx_vid) + b_uv(p,1,j,ke_xy,ke_z)
          prog_vars(pp,pv,ke_xy,ke_z,momy_vid) = prog_vars(pp,pv,ke_xy,ke_z,momy_vid) + b_uv(p,2,j,ke_xy,ke_z)
        end do
        end do
      end do
    end do
    end do
 
    call prof_rapend( 'hevi_cal_vi_lin_uv', 3) 
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve_uv
 
!--
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd( &
    BndMatL, BndMatD, BndMatU,              & ! (out)
    dens, w, wt, pot, dpdrhot,              & ! (in)
    alph, gsqrtv, nz,                       & ! (in)
    id, fac_dx3, fac_intrpmat_vpordm1,      & ! (in)
    impl_fac, dt,                           & ! (in)
    lmesh, elem, im, jm,                    & ! (in)
    ke_z                                    ) ! (in)
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: im, jm
    real(rp), intent(out) :: bndmatl(im,3,elem%nnode_v,3,jm,lmesh%ne2d)
    real(rp), intent(out) :: bndmatd(im,3,elem%nnode_v,3,elem%nnode_v,jm,lmesh%ne2d)
    real(rp), intent(out) :: bndmatu(im,3,elem%nnode_v,3,jm,lmesh%ne2d)
    real(rp), intent(in) :: dens(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: w(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: wt(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: pot(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: dpdrhot(im,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez)
    real(rp), intent(in) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
    real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: id(elem%nnode_v,elem%nnode_v)
    real(rp), intent(in) :: fac_dx3(elem%nnode_v,elem%nnode_v)
    real(rp), intent(in) :: fac_intrpmat_vpordm1(elem%nnode_v,elem%nnode_v)
    real(rp), intent(in) :: impl_fac
    real(rp), intent(in) :: dt
    integer, intent(in) :: ke_z
 
    real(rp) :: fac_dz(im,elem%nnode_v,elem%nnode_v)
    real(rp) :: tmp1(im,elem%nnode_v)
    real(rp) :: tmp2(im,elem%nnode_v)
 
    integer :: p, p2, pv, pv1, pv2
    integer :: fp,fp2
    integer :: ke, ke2d
    integer :: ke_z2, ke_2
    integer :: ij, i, j
    integer :: f1, f2
    integer :: v
 
    real(rp) :: fac
    logical :: bc_flag
 
    integer, parameter :: dens_vid_lc = 1
    integer, parameter :: momz_vid_lc = 2
    integer, parameter :: rhot_vid_lc = 3
 
    !---------------------------------------------------------------------------
 
    call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
 
    !$omp parallel private(ke2D,j,ke, pv,pv1,pv2, p,p2, fp,fp2, ke_z2,ke_2, i,ij, f1,f2, v, fac, bc_flag, &
    !$omp fac_dz, tmp1, tmp2 )
    !$omp do collapse(2)
    do ke2d=1, lmesh%Ne2D
    do j=1, jm
      ke = ke2d + (ke_z-1)*lmesh%Ne2D
 
      do pv2=1, elem%Nnode_v
      do pv=1, elem%Nnode_v
      do i=1, im
        ij = i + (j-1)*im
        p = ij + (pv-1)*im*jm
        fac_dz(i,pv,pv2) = lmesh%Escale(p,ke,3,3) / gsqrtv(p,ke) * fac_dx3(pv,pv2)
      end do
      end do
      end do
 
      do pv2=1, elem%Nnode_v
      do pv=1, elem%Nnode_v
      do i=1, im
        ! DDENS
        bndmatd(i,dens_vid_lc,pv,dens_vid_lc,pv2,j,ke2d) = id(pv,pv2)
        bndmatd(i,dens_vid_lc,pv,momz_vid_lc,pv2,j,ke2d) = fac_dz(i,pv,pv2)
        bndmatd(i,dens_vid_lc,pv,rhot_vid_lc,pv2,j,ke2d) = 0.0_rp
 
        ! MOMZ
        bndmatd(i,momz_vid_lc,pv,momz_vid_lc,pv2,j,ke2d) = id(pv,pv2) ! &
                                          !+ fac_dz_p(i,pv,pv2,j) * 2.0_RP * W0(i,pv2,j)  [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
        bndmatd(i,momz_vid_lc,pv,dens_vid_lc,pv2,j,ke2d) = fac_intrpmat_vpordm1(pv,pv2) ! &
                                          !- fac_dz_p(i,pv,pv2,j) * W0(i,pv2,j)**2        [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
        bndmatd(i,momz_vid_lc,pv,rhot_vid_lc,pv2,j,ke2d) = fac_dz(i,pv,pv2) * dpdrhot(i,pv2,j,ke2d,ke_z)
 
        ! DRHOT
        bndmatd(i,rhot_vid_lc,pv,dens_vid_lc,pv2,j,ke2d) = - fac_dz(i,pv,pv2) * pot(i,pv2,j,ke2d,ke_z) * wt(i,pv2,j,ke2d,ke_z)
        bndmatd(i,rhot_vid_lc,pv,momz_vid_lc,pv2,j,ke2d) =   fac_dz(i,pv,pv2) * pot(i,pv2,j,ke2d,ke_z)
        bndmatd(i,rhot_vid_lc,pv,rhot_vid_lc,pv2,j,ke2d) = id(pv,pv2) + fac_dz(i,pv,pv2) * wt(i,pv2,j,ke2d,ke_z)
      end do
      end do
      end do
 
      do f1=1, 2
        if (f1==1) then
          ke_z2 = max(ke_z-1,1)
          pv1 = 1; pv2 = elem%Nnode_v; f2 = 2
        else
          ke_z2 = min(ke_z+1,lmesh%NeZ)
          pv1 = elem%Nnode_v; pv2 = 1; f2 = 1
        end if
        ke_2 = ke2d + (ke_z2-1)*lmesh%Ne2D
 
        if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
          bc_flag = .true.
          pv2 = pv1; f2 = f1
        else 
          bc_flag = .false.    
        end if
 
        do pv=1, elem%Nnode_v
        do i=1, im
          ij = i + (j-1)*im
          fp  = elem%Nfp_h * elem%Nfaces_h + (f1-1)*elem%Nfp_v + ij
          fp2 = elem%Nfp_h * elem%Nfaces_h + (f2-1)*elem%Nfp_v + ij
          p = ij + (pv-1)*im*jm
 
          !--
          fac  = 0.5_rp * impl_fac / gsqrtv(p,ke) * elem%Lift(p,fp) * lmesh%Fscale(fp,ke)
          tmp1(i,pv) = fac * max( alph(fp,ke), alph(fp2,ke_2) )
          tmp2(i,pv) = fac * nz(fp,ke)
        end do
        end do
 
        !--
 
        if ( bc_flag ) then
          do pv=1, elem%Nnode_v
          do i=1, im
           ! d ( Eqs. ) / d DENS
            bndmatd(i,rhot_vid_lc,pv,dens_vid_lc,pv1,j,ke2d) = bndmatd(i,rhot_vid_lc,pv,dens_vid_lc,pv1,j,ke2d) + 2.0_rp * tmp2(i,pv) * pot(i,pv1,j,ke2d,ke_z) * wt(i,pv1,j,ke2d,ke_z)
 
           ! d ( Eqs. ) / d MOMZ            
            bndmatd(i,dens_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) = bndmatd(i,dens_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) - 2.0_rp * tmp2(i,pv)
            bndmatd(i,momz_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) = bndmatd(i,momz_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) + 2.0_rp * tmp1(i,pv)
            bndmatd(i,rhot_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) = bndmatd(i,rhot_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) - 2.0_rp * tmp2(i,pv) * pot(i,pv1,j,ke2d,ke_z)
 
           ! d ( Eqs. ) / d RHOT
            bndmatd(i,rhot_vid_lc,pv,rhot_vid_lc,pv1,j,ke2d) = bndmatd(i,rhot_vid_lc,pv,rhot_vid_lc,pv1,j,ke2d) - 2.0_rp * tmp2(i,pv) * wt(i,pv1,j,ke2d,ke_z)
 
          end do
          end do
        else
          do pv=1, elem%Nnode_v
          do i=1, im
            ! d ( Eqs. ) / d DENS
            bndmatd(i,dens_vid_lc,pv,dens_vid_lc,pv1,j,ke2d) = bndmatd(i,dens_vid_lc,pv,dens_vid_lc,pv1,j,ke2d) + tmp1(i,pv)
            bndmatd(i,rhot_vid_lc,pv,dens_vid_lc,pv1,j,ke2d) = bndmatd(i,rhot_vid_lc,pv,dens_vid_lc,pv1,j,ke2d) + tmp2(i,pv) * pot(i,pv1,j,ke2d,ke_z) * wt(i,pv1,j,ke2d,ke_z)
 
            ! d ( Eqs. ) / d MOMZ
            bndmatd(i,dens_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) = bndmatd(i,dens_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) - tmp2(i,pv)
            bndmatd(i,momz_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) = bndmatd(i,momz_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) + tmp1(i,pv)
            bndmatd(i,rhot_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) = bndmatd(i,rhot_vid_lc,pv,momz_vid_lc,pv1,j,ke2d) - tmp2(i,pv) * pot(i,pv1,j,ke2d,ke_z)
 
            ! d ( Eqs. ) / d RHOT
            bndmatd(i,momz_vid_lc,pv,rhot_vid_lc,pv1,j,ke2d) = bndmatd(i,momz_vid_lc,pv,rhot_vid_lc,pv1,j,ke2d) - tmp2(i,pv) * dpdrhot(i,pv1,j,ke2d,ke_z)
            bndmatd(i,rhot_vid_lc,pv,rhot_vid_lc,pv1,j,ke2d) = bndmatd(i,rhot_vid_lc,pv,rhot_vid_lc,pv1,j,ke2d) + tmp1(i,pv) &
                                                                                                                - tmp2(i,pv) * wt(i,pv1,j,ke2d,ke_z)
          end do
          end do
          if ( f1 == 1 ) then
            do pv=1, elem%Nnode_v
            do i=1, im
              ! d ( Eqs. ) / d DENS
              bndmatl(i,dens_vid_lc,pv,dens_vid_lc,j,ke2d) = - tmp1(i,pv)
              bndmatl(i,momz_vid_lc,pv,dens_vid_lc,j,ke2d) = 0.0_rp
              bndmatl(i,rhot_vid_lc,pv,dens_vid_lc,j,ke2d) = - tmp2(i,pv) * pot(i,pv2,j,ke2d,ke_z2) * wt(i,pv2,j,ke2d,ke_z2)
              ! d ( Eqs. ) / d MOMZ
              bndmatl(i,dens_vid_lc,pv,momz_vid_lc,j,ke2d) = + tmp2(i,pv)
              bndmatl(i,momz_vid_lc,pv,momz_vid_lc,j,ke2d) = - tmp1(i,pv)            
              bndmatl(i,rhot_vid_lc,pv,momz_vid_lc,j,ke2d) = + tmp2(i,pv) * pot(i,pv2,j,ke2d,ke_z2)
              ! d ( Eqs. ) / d RHOT
              bndmatl(i,dens_vid_lc,pv,rhot_vid_lc,j,ke2d) = 0.0_rp  
              bndmatl(i,momz_vid_lc,pv,rhot_vid_lc,j,ke2d) = + tmp2(i,pv) * dpdrhot(i,pv2,j,ke2d,ke_z2)
              bndmatl(i,rhot_vid_lc,pv,rhot_vid_lc,j,ke2d) = - tmp1(i,pv) &
                                                             + tmp2(i,pv) * wt(i,pv2,j,ke2d,ke_z2)
            end do
            end do            
          else
            do pv=1, elem%Nnode_v
            do i=1, im
              ! d ( Eqs. ) / d DENS
              bndmatu(i,dens_vid_lc,pv,dens_vid_lc,j,ke2d) = - tmp1(i,pv)
              bndmatu(i,momz_vid_lc,pv,dens_vid_lc,j,ke2d) = 0.0_rp
              bndmatu(i,rhot_vid_lc,pv,dens_vid_lc,j,ke2d) = - tmp2(i,pv) * pot(i,pv2,j,ke2d,ke_z2) * wt(i,pv2,j,ke2d,ke_z2)
              ! d ( Eqs. ) / d MOMZ
              bndmatu(i,dens_vid_lc,pv,momz_vid_lc,j,ke2d) = + tmp2(i,pv)
              bndmatu(i,momz_vid_lc,pv,momz_vid_lc,j,ke2d) = - tmp1(i,pv)            
              bndmatu(i,rhot_vid_lc,pv,momz_vid_lc,j,ke2d) = + tmp2(i,pv) * pot(i,pv2,j,ke2d,ke_z2)
              ! d ( Eqs. ) / d RHOT
              bndmatu(i,dens_vid_lc,pv,rhot_vid_lc,j,ke2d) = 0.0_rp  
              bndmatu(i,momz_vid_lc,pv,rhot_vid_lc,j,ke2d) = + tmp2(i,pv) * dpdrhot(i,pv2,j,ke2d,ke_z2)
              bndmatu(i,rhot_vid_lc,pv,rhot_vid_lc,j,ke2d) = - tmp1(i,pv) &
                                                             + tmp2(i,pv) * wt(i,pv2,j,ke2d,ke_z2)
            end do
            end do            
          end if   
        end if
 
      end do ! end do for f
    end do ! end for j
    end do ! end do for ke2D
    !$omp end do
    !$omp end parallel
 
    call prof_rapend( 'hevi_cal_vi_matbnd', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd
 
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv( &
    BndMatL_uv, BndMatD_uv, BndMatU_uv,     & ! (out)
    alph, gsqrtv,                           & ! (in)
    impl_fac, dt,                           & ! (in)
    lmesh, elem, im, jm,                    & ! (in)
    vmapm, vmapp, ke_z                      ) ! (in)
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: im, jm
    real(rp), intent(out) :: bndmatl_uv(im,elem%nnode_v,jm,lmesh%ne2d)
    real(rp), intent(out) :: bndmatd_uv(im,elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d)
    real(rp), intent(out) :: bndmatu_uv(im,elem%nnode_v,jm,lmesh%ne2d)
    real(rp), intent(in) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
    real(rp), intent(in) :: impl_fac
    real(rp), intent(in) :: dt
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)    
    integer, intent(in) :: ke_z
 
    integer :: colmask(elem%nnode_v)
    real(rp) :: id(elem%nnode_v,elem%nnode_v)
    real(rp) :: dd(elem%nnode_v)
    real(rp) :: tmp1(im,elem%nnode_v)
 
    integer :: p, pv, pv1, pv2
    integer :: fp,fp2
    integer :: ke, ke2d
    integer :: ke_z2, ke_2
    integer :: ij, i, j
    integer :: f1, f2
 
    real(rp) :: fac
    logical :: bc_flag
    !---------------------------------------------------------------------------
 
    call prof_rapstart( 'hevi_cal_vi_matbnd_uv', 3)
 
    id(:,:) = 0.0_rp
    do p=1, elem%Nnode_v
      id(p,p) = 1.0_rp
    end do
 
    !$omp parallel private(ke2D,j,ke, pv2,pv,i, f1, ke_z2, pv1,f2, bc_flag, ij,fp,fp2,p, fac,tmp1) 
    !$omp do collapse(2)
    do ke2d=1, lmesh%Ne2D
    do j=1, jm
      ke = ke2d + (ke_z-1)*lmesh%Ne2D
 
      do pv2=1, elem%Nnode_v
      do pv=1, elem%Nnode_v
      do i=1, im
        bndmatd_uv(i,pv,pv2,j,ke2d) = id(pv,pv2)
      end do
      end do
      end do
 
      do f1=1, 2
        if (f1==1) then
          ke_z2 = max(ke_z-1,1)
          pv1 = 1; f2 = 2
        else
          ke_z2 = min(ke_z+1,lmesh%NeZ)
          pv1 = elem%Nnode_v; f2 = 1
        end if
        ke_2 = ke2d + (ke_z2-1)*lmesh%Ne2D
 
        if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
          bc_flag = .true.
          f2 = f1
        else 
          bc_flag = .false.    
        end if
 
        do pv=1, elem%Nnode_v
        do i=1, im
          ij = i + (j-1)*im
          fp  = elem%Nfp_h * elem%Nfaces_h + (f1-1)*elem%Nfp_v + ij
          fp2 = elem%Nfp_h * elem%Nfaces_h + (f2-1)*elem%Nfp_v + ij
          p = ij + (pv-1)*im*jm
 
          !--
          fac  = 0.5_rp * impl_fac / gsqrtv(p,ke)
          tmp1(i,pv) = fac * elem%Lift(p,fp) * lmesh%Fscale(fp,ke) &
              * max( alph(fp,ke), alph(fp2,ke_2) )
        end do
        end do
 
        if (bc_flag) then
        else
          bndmatd_uv(:,:,pv1,j,ke2d) = bndmatd_uv(:,:,pv1,j,ke2d) + tmp1(:,:)
          if (f1 == 1) then
            bndmatl_uv(:,:,j,ke2d) = - tmp1(:,:)
          else
            bndmatu_uv(:,:,j,ke2d) = - tmp1(:,:)
          end if
        end if
      end do ! end do for f
 
    end do ! end do for j
    end do ! end do for ke2D
    !$omp end do
    !$omp end parallel
    call prof_rapend( 'hevi_cal_vi_matbnd_uv', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv
 
!-- private ----------------
 
!OCL SERIAL  
  subroutine vi_cal_del_flux_dyn_uv( del_flux, alph,          & ! (out)
    pvars_, pvars0_,                                          & ! (in)
    dens_hyd, pres_hyd, rtot, cptot, cvtot,                   & ! (in)
    gsqrt, g11, g12, g22, gsqrth, gam, gsqrtv, g13, g23,      & ! (in)
    nz, vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d  ) ! (in)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(out) ::  del_flux(elem%nfptot,2,lmesh%ne)
    real(rp), intent(out) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in) ::  pvars_ (elem%np*lmesh%ne,prgvar_num)
    real(rp), intent(in) ::  pvars0_(elem%np*lmesh%ne,prgvar_num)
    real(rp), intent(in) ::  dens_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  pres_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  rtot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  cptot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  cvtot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  gsqrt(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g11(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g12(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g22(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  gsqrth(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) :: gam(elem%np*lmesh%nea)
    real(rp), intent(in) ::  gsqrtv(elem%np*lmesh%ne)
    real(rp), intent(in) ::  g13(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g23(elem%np*lmesh%nea)
    real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
    integer, intent(in) :: im2dto3d(elem%nfptot)
    
    integer :: fp, ke, ke2d
    integer :: ip(elem%nfptot), im(elem%nfptot)
    real(rp) :: momx0(elem%nfptot,2)
    real(rp) :: momy0(elem%nfptot,2)
    real(rp) :: drhot0(elem%nfptot,2)
    real(rp) :: rhot_hyd(elem%nfptot,2)
    real(rp) :: cptot_ov_cvtot(elem%nfptot,2)
    real(rp) :: rtot_(elem%nfptot,2)
    real(rp) :: wt0(elem%nfptot,2)
    real(rp) :: pres0(elem%nfptot,2)
    real(rp) :: rdens0(elem%nfptot,2)
 
    real(rp) :: gsqrt_(elem%nfptot,2)
    real(rp) :: gsqrtv_(elem%nfptot,2)
    real(rp) :: rgsqrtv(elem%nfptot,2)
    real(rp) :: rgam2(elem%nfptot,2)
    real(rp) :: g13_(elem%nfptot,2)
    real(rp) :: g23_(elem%nfptot,2)
    real(rp) :: gxz_(elem%nfptot,2)
    real(rp) :: gyz_(elem%nfptot,2)
    real(rp) :: g11_, g12_, g22_
    real(rp) :: gnn_m, gnn_p
 
    real(rp) :: tmp1
 
    real(rp) :: gamm, rgamm, pres0ovrdry, rp0
 
    integer, parameter :: in = 1
    integer, parameter :: ex = 2
    !------------------------------------------------------------------------
 
    gamm  = cpdry / cvdry
    rgamm = cvdry / cpdry
    pres0ovrdry = pres00 / rdry
    rp0 = 1.0_rp / pres00
    
    !$omp parallel private( ke, ke2D, fp, iM, iP,                        &
    !$omp MOMX0, MOMY0, DRHOT0, wt0, rdens0, pres0,                      &
    !$omp CPtot_ov_CVtot, Rtot_, rhot_hyd,                               &
    !$omp Gsqrt_, GsqrtV_, RGsqrtV, G11_, G12_, G22_, G13_, G23_, rgam2, &
    !$omp Gxz_, Gyz_, Gnn_M, Gnn_P,                                      &
    !$omp tmp1 )  
 
    !$omp do
    do ke=lmesh%NeS, lmesh%NeE
      ke2d = lmesh%EMap3Dto2D(ke)
      im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
 
      gsqrt_(:,in) = gsqrt(im)
      gsqrt_(:,ex) = gsqrt(ip)
 
      rgam2(:,in) = 1.0_rp / gam(im)**2
      rgam2(:,ex) = 1.0_rp / gam(ip)**2
      gsqrtv_(:,in) = gsqrtv(im)
      gsqrtv_(:,ex) = gsqrtv(ip)
      rgsqrtv(:,:) = 1.0_rp / gsqrtv_(:,:)      
 
      g13_(:,in) = g13(im)
      g13_(:,ex) = g13(ip)
      g23_(:,in) = g23(im)
      g23_(:,ex) = g23(ip)
 
      momx0(:,in) = pvars0_(im(:),momx_vid)
      momx0(:,ex) = pvars0_(ip(:),momx_vid)
      momy0(:,in) = pvars0_(im(:),momy_vid)
      momy0(:,ex) = pvars0_(ip(:),momy_vid)
      drhot0(:,in) = pvars0_(im(:),rhot_vid)
      drhot0(:,ex) = pvars0_(ip(:),rhot_vid)
      
      rtot_(:,in) = rtot(im(:))
      rtot_(:,ex) = rtot(ip(:))
      cptot_ov_cvtot(:,in) = cptot(im(:)) / cvtot(im(:))
      cptot_ov_cvtot(:,ex) = cptot(ip(:)) / cvtot(ip(:))
 
      !-
      rdens0(:,in) = 1.0_rp / ( dens_hyd(im) + pvars0_(im,dens_vid) )
      rdens0(:,ex) = 1.0_rp / ( dens_hyd(ip) + pvars0_(ip,dens_vid) )
 
      wt0(:,in) = ( pvars0_(im(:),momz_vid) * rgsqrtv(:,in) + g13_(:,in) * momx0(:,in)  &
                                                            + g23_(:,in) * momy0(:,in)  ) * rdens0(:,in)     
      wt0(:,ex) = ( pvars0_(ip(:),momz_vid) * rgsqrtv(:,ex) + g13_(:,ex) * momx0(:,ex)  &
                                                            + g23_(:,ex) * momy0(:,ex)  ) * rdens0(:,ex)
      
      rhot_hyd(:,in) = pres0ovrdry * (pres_hyd(im(:))/pres00)**rgamm
      rhot_hyd(:,ex) = pres0ovrdry * (pres_hyd(ip(:))/pres00)**rgamm
      pres0(:,:) = pres00 * ( rtot_(:,:) * rp0 * ( rhot_hyd(:,:) + drhot0(:,:) ) )**cptot_ov_cvtot(:,:)
 
      !-
      do fp=1, elem%NfpTot
        g11_ = g11(im2dto3d(fp),ke2d);  g12_ = g12(im2dto3d(fp),ke2d);  g22_ = g22(im2dto3d(fp),ke2d)
 
        gxz_(fp,in) = rgam2(fp,in) * ( g11_ * g13_(fp,in) + g12_ * g23_(fp,in) )
        gxz_(fp,ex) = rgam2(fp,ex) * ( g11_ * g13_(fp,ex) + g12_ * g23_(fp,ex) )
 
        gyz_(fp,in) = rgam2(fp,in) * ( g12_ * g13_(fp,in) + g22_ * g23_(fp,in) )
        gyz_(fp,ex) = rgam2(fp,ex) * ( g12_ * g13_(fp,ex) + g22_ * g23_(fp,ex) )
      end do
 
      do fp=1, elem%NfpTot
        gnn_m = ( 1.0_rp * rgsqrtv(fp,in)**2 + g13_(fp,in) * gxz_(fp,in) + g23_(fp,in) * gyz_(fp,in) ) * abs( nz(fp,ke) )
        gnn_p = ( 1.0_rp * rgsqrtv(fp,ex)**2 + g13_(fp,ex) * gxz_(fp,ex) + g23_(fp,ex) * gyz_(fp,ex) ) * abs( nz(fp,ke) )
 
        alph(fp,ke) = nz(fp,ke)**2 * max( abs( wt0(fp,in) ) + sqrt( gnn_m * gamm * pres0(fp,in) * rdens0(fp,in) ), &
                                          abs( wt0(fp,ex) ) + sqrt( gnn_p * gamm * pres0(fp,ex) * rdens0(fp,ex) )  )
      end do
 
      !---- 
      do fp=1, elem%NfpTot
        tmp1 = - 0.5_rp * lmesh%Fscale(fp,ke) * alph(fp,ke)
        del_flux(fp,1,ke) = tmp1 * ( pvars_(ip(fp),momx_vid) - pvars_(im(fp),momx_vid) )
        del_flux(fp,2,ke) = tmp1 * ( pvars_(ip(fp),momy_vid) - pvars_(im(fp),momy_vid) )                       
      end do
    end do
    !$omp end do
    !$omp end parallel
 
    return
  end subroutine vi_cal_del_flux_dyn_uv
 
!OCL SERIAL  
  subroutine vi_cal_del_flux_dyn( del_flux, alph,            & ! (out)
    pvars_, pvars0_,                                         & ! (in)
    dens_hyd, pres_hyd, rtot, cptot, cvtot,                  & ! (in)
    gsqrt, gsqrth, gam, gsqrtv, g13, g23,                    & ! (in)
    nz, vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d ) ! (in)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d 
    real(rp), intent(out) ::  del_flux(elem%nfptot,prgvar_num,lmesh%ne)
    real(rp), intent(in) :: alph(elem%nfptot,lmesh%ne)
    real(rp), intent(in) ::  pvars_ (elem%np*lmesh%ne,prgvar_num)
    real(rp), intent(in) ::  pvars0_(elem%np*lmesh%ne,prgvar_num)
    real(rp), intent(in) ::  dens_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  pres_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  rtot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  cptot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  cvtot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  gsqrt(elem%np*lmesh%nea)
    real(rp), intent(in) ::  gsqrth(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) :: gam(elem%np*lmesh%nea)
    real(rp), intent(in) ::  gsqrtv(elem%np*lmesh%ne)
    real(rp), intent(in) ::  g13(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g23(elem%np*lmesh%nea)
    real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
    integer, intent(in) :: im2dto3d(elem%nfptot)
    
    integer :: fp, ke, ke2d, ke_z, ke_xy
    integer :: ip(elem%nfptot), im(elem%nfptot)
    real(rp) :: ddens(elem%nfptot,2)
    real(rp) :: momz(elem%nfptot,2)
    real(rp) :: momw(elem%nfptot,2)
    real(rp) :: drhot(elem%nfptot,2)
    real(rp) :: rhot_hyd(elem%nfptot,2)
    real(rp) :: cptot_ov_cvtot(elem%nfptot,2)
    real(rp) :: rtot_(elem%nfptot,2)
    real(rp) :: pott(elem%nfptot,2)
    real(rp) :: dpres(elem%nfptot,2)
    real(rp) :: dens(elem%nfptot,2)
 
    real(rp) :: gsqrt_(elem%nfptot,2)
    real(rp) :: gsqrtv_(elem%nfptot,2)
    real(rp) :: rgam2(elem%nfptot,2)
    real(rp) :: g13_(elem%nfptot,2)
    real(rp) :: g23_(elem%nfptot,2)
 
    real(rp) :: tmp1
 
    real(rp) :: rgamm, pres0ovrdry, rp0
 
    integer, parameter :: in = 1
    integer, parameter :: ex = 2
    !------------------------------------------------------------------------
 
    rgamm = cvdry / cpdry
    pres0ovrdry = pres00 / rdry
    rp0 = 1.0_rp / pres00
    
    !$omp parallel private( ke_z, ke_xy, ke, ke2D, fp, iM, iP, &
    !$omp DDENS, MOMZ, MOMW, DRHOT,                           &
    !$omp CPtot_ov_CVtot, Rtot_, rhot_hyd, pott, dpres, dens, &
    !$omp Gsqrt_, GsqrtV_, G13_, G23_, rgam2                  )  
 
    !$omp do collapse(2)
    do ke_z=1, lmesh%NeZ
    do ke_xy=1, lmesh%Ne2D
      ke = ke_xy + (ke_z-1)*lmesh%Ne2D
      ke2d = lmesh%EMap3Dto2D(ke)
 
      im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
 
      gsqrt_(:,in) = gsqrt(im)
      gsqrt_(:,ex) = gsqrt(ip)
 
      rgam2(:,in) = 1.0_rp / gam(im)**2
      rgam2(:,ex) = 1.0_rp / gam(ip)**2
      gsqrtv_(:,in) = gsqrtv(im)
      gsqrtv_(:,ex) = gsqrtv(ip)
 
      g13_(:,in) = g13(im)
      g13_(:,ex) = g13(ip)
      g23_(:,in) = g23(im)
      g23_(:,ex) = g23(ip)
 
      ddens(:,in) = pvars_(im(:),dens_vid)
      ddens(:,ex) = pvars_(ip(:),dens_vid)
      momz(:,in) = pvars_(im(:),momz_vid)
      momz(:,ex) = pvars_(ip(:),momz_vid)
      drhot(:,in) = pvars_(im(:),rhot_vid)
      drhot(:,ex) = pvars_(ip(:),rhot_vid)
 
      rtot_(:,in) = rtot(im(:))
      rtot_(:,ex) = rtot(ip(:))
      cptot_ov_cvtot(:,in) = cptot(im(:)) / cvtot(im(:))
      cptot_ov_cvtot(:,ex) = cptot(ip(:)) / cvtot(ip(:))
      !-
 
      dens(:,in) = dens_hyd(im(:)) + ddens(:,in)
      dens(:,ex) = dens_hyd(ip(:)) + ddens(:,ex)
 
      rhot_hyd(:,in) = pres0ovrdry * (pres_hyd(im(:))/pres00)**rgamm
      rhot_hyd(:,ex) = pres0ovrdry * (pres_hyd(ip(:))/pres00)**rgamm      
      pott(:,:) = ( rhot_hyd(:,:) + drhot(:,:) ) / dens(:,:)
 
      dpres(:,:) = pres00 * ( rtot_(:,:) * rp0 * dens(:,:) * pott(:,:) )**cptot_ov_cvtot(:,:)
      dpres(:,in) = dpres(:,in) - pres_hyd(im(:)) 
      dpres(:,ex) = dpres(:,ex) - pres_hyd(ip(:)) 
      
      !-
      momw(:,in) = momz(:,in) &
              + gsqrtv_(:,in) * g13_(:,in) * pvars_(im,momx_vid) &
              + gsqrtv_(:,in) * g23_(:,in) * pvars_(im,momy_vid)
      momw(:,ex) = momz(:,ex) &
              + gsqrtv_(:,ex) * g13_(:,ex) * pvars_(ip,momx_vid) &
              + gsqrtv_(:,ex) * g23_(:,ex) * pvars_(ip,momy_vid)
      
      !----                  
                    
      if ( ke_z == 1 .or. ke_z == lmesh%NeZ) then
        do fp=1, elem%NfpTot
          if ( im(fp) == ip(fp) ) then
            momz(fp,ex) = - momz(fp,in) &
                          - 2.0_rp * gsqrtv_(fp,in) * (  g13_(fp,in) * pvars_(im(fp),momx_vid) &
                                                       + g23_(fp,in) * pvars_(im(fp),momy_vid) )          
            momw(fp,ex) = - momw(fp,in)
          end if
        end do
      else
      end if
 
      do fp=1, elem%NfpTot
        tmp1 = 0.5_rp * lmesh%Fscale(fp,ke)
 
        del_flux(fp,dens_vid,ke) = tmp1 * ( &
                  + ( momw(fp,ex) - momw(fp,in) ) * nz(fp,ke)    &
                  - alph(fp,ke) * ( ddens(fp,ex) - ddens(fp,in) ) )
        
        
        del_flux(fp,momz_vid,ke) = tmp1 * ( &
                  + ( dpres(fp,ex) - dpres(fp,in) ) * nz(fp,ke) &                    
                  - alph(fp,ke) * ( momz(fp,ex) - momz(fp,in) ) )
        
        del_flux(fp,rhot_vid,ke) = tmp1 * ( &
                  + ( pott(fp,ex) * momw(fp,ex)  - pott(fp,in) * momw(fp,in) ) * nz(fp,ke) &
                  - alph(fp,ke) * ( drhot(fp,ex) - drhot(fp,in) ) )
      end do
    end do
    end do
    !$omp end do
    !$omp end parallel
 
    return
  end subroutine vi_cal_del_flux_dyn
end module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common_2