scale_atm_phy_tb_dgm_dns.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics turbulence
!!
!! @par Description
!!      Turbulence process for DNS
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_tb_dgm_dns
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    eps  => const_eps,     &
    grav => const_grav,    &
    rdry => const_rdry,    &
    cpdry => const_cpdry,  &
    cvdry => const_cvdry,  &
    pres00 => const_pre00, &
    karman  => const_karman, &
    rplanet => const_radius    
 
  use scale_sparsemat  
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  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
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_phy_tb_dgm_dns_init
  public :: atm_phy_tb_dgm_dns_final
  public :: atm_phy_tb_dgm_dns_cal_grad
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: cal_del_flux_grad
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private, parameter   :: OneOverThree  = 1.0_rp / 3.0_rp
  real(RP), private, parameter   :: twoOverThree  = 2.0_rp / 3.0_rp
  real(RP), private, parameter   :: FourOverThree = 4.0_rp / 3.0_rp
 
  real(RP), private              :: DNS_NU        = 1.512e-5_rp !< kinematic viscosity coefficient [m2/s] for air at 20degC
  real(RP), private              :: DNS_MU        = 1.8e-5_rp   !< molecular diffusive coefficient [m2/s] for air at 20degC
 
contains
!OCL SERIAL
  subroutine atm_phy_tb_dgm_dns_init( mesh )
    implicit none    
    class(meshbase3d), intent(in) :: mesh
 
    logical  :: consistent_tke = .true.
 
    namelist / param_atmos_phy_tb_dgm_dns / &
      dns_nu, dns_mu
        
    integer :: ierr
    !--------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_TB_dgm_dns_setup",*) 'Setup'
    log_info("ATMOS_PHY_TB_dgm_dns_setup",*) 'Eddy Viscocity Model for DNS'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_tb_dgm_dns,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_TB_dgm_dns_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_TB_dgm_dns_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_DGM_DNS. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_tb_dgm_dns)
 
    return
  end subroutine atm_phy_tb_dgm_dns_init
 
!OCL SERIAL
  subroutine atm_phy_tb_dgm_dns_final()
    implicit none
    !--------------------------------------------------------------------
 
    return
  end subroutine atm_phy_tb_dgm_dns_final
 
!> Calculate parameterized stress tensor and eddy heat flux with turbulent model
!!
!OCL SERIAL  
  subroutine atm_phy_tb_dgm_dns_cal_grad( &
    T11, T12, T13, T21, T22, T23, T31, T32, T33,                & ! (out)
    df1, df2, df3,                                              & ! (out)
    tke, nu, kh,                                                & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd,    & ! (in)
    pres, pt,                                                   & ! (in)
    dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d, & ! (in)
    is_bound                                                    ) ! (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) :: t11(elem%np,lmesh%nea)      !< (1,1) component of stress tensor
    real(rp), intent(out) :: t12(elem%np,lmesh%nea)      !< (1,2) component of stress tensor
    real(rp), intent(out) :: t13(elem%np,lmesh%nea)      !< (1,3) component of stress tensor
    real(rp), intent(out) :: t21(elem%np,lmesh%nea)      !< (2,1) component of stress tensor
    real(rp), intent(out) :: t22(elem%np,lmesh%nea)      !< (2,2) component of stress tensor
    real(rp), intent(out) :: t23(elem%np,lmesh%nea)      !< (2,3) component of stress tensor
    real(rp), intent(out) :: t31(elem%np,lmesh%nea)      !< (3,1) component of stress tensor
    real(rp), intent(out) :: t32(elem%np,lmesh%nea)      !< (3,2) component of stress tensor
    real(rp), intent(out) :: t33(elem%np,lmesh%nea)      !< (3,3) component of stress tensor
    real(rp), intent(out)  :: df1(elem%np,lmesh%nea)     !< Diffusive heat flux in x1 direction / density
    real(rp), intent(out)  :: df2(elem%np,lmesh%nea)     !< Diffusive heat flux in x2 direction / density
    real(rp), intent(out)  :: df3(elem%np,lmesh%nea)     !< Diffusive heat flux in x3 direction / density
    real(rp), intent(out) :: tke(elem%np,lmesh%nea)      !< Parameterized turbulent kinetic energy
    real(rp), intent(out) :: nu(elem%np,lmesh%nea)       !< Eddy viscosity
    real(rp), intent(out) :: kh(elem%np,lmesh%nea)       !< Eddy diffusivity
    real(rp), intent(in)  :: ddens_(elem%np,lmesh%nea)   !< Density perturbation
    real(rp), intent(in)  :: momx_ (elem%np,lmesh%nea)   !< Momentum in x1 direction
    real(rp), intent(in)  :: momy_ (elem%np,lmesh%nea)   !< Momentum in x2 direction
    real(rp), intent(in)  :: momz_ (elem%np,lmesh%nea)   !< Momentum in x3 direction
    real(rp), intent(in)  :: drhot_(elem%np,lmesh%nea)   !< Density x potential temperature perturbation
    real(rp), intent(in)  :: dens_hyd(elem%np,lmesh%nea) !< Reference pressure in hydrostatic balance
    real(rp), intent(in)  :: pres_hyd(elem%np,lmesh%nea) !< Reference density in hydrostatic balance
    real(rp), intent(in)  :: pres(elem%np,lmesh%nea)     !< Pressure
    real(rp), intent(in)  :: pt(elem%np,lmesh%nea)       !< Potential temperature
    type(sparsemat), intent(in) :: dx, dy, dz             !< Differential matrix managed by sparse matrix type
    type(sparsemat), intent(in) :: sx, sy, sz             !< Stiffness matrix managed by sparse matrix type
    type(sparsemat), intent(in) :: lift                   !< Lifting matrix managed by sparse matrix type
    logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne) !< Flag whether nodes are located at domain boundaries
 
    real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
    real(rp) :: dens(elem%np), rdens(elem%np), rhot(elem%np), q(elem%np)
    real(rp) :: ddensdxi(elem%np,3)
    real(rp) :: dveldxi(elem%np,3,3)
    real(rp) :: del_flux_rho (elem%nfptot,lmesh%ne,3)
    real(rp) :: del_flux_mom (elem%nfptot,lmesh%ne,3,3)
    real(rp) :: del_flux_rhot(elem%nfptot,lmesh%ne,3)
 
    real(rp) :: s11(elem%np), s12(elem%np), s22(elem%np), s23(elem%np), s31(elem%np), s33(elem%np)
    real(rp) :: skkovthree
    real(rp) :: coef
 
    integer :: ke
    integer :: p
    !--------------------------------------------------------------------
 
    call cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot,        & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, pt,            & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (in)
      lmesh, elem, is_bound )                                                   ! (in)
  
    !$omp parallel do private( &
    !$omp Fx, Fy, Fz, LiftDelFlx,                  &
    !$omp DENS, RHOT, RDENS, Q, DdensDxi, DVelDxi, &
    !$omp S11, S12, S22, S23, S31, S33,            &
    !$omp coef, SkkOvThree,                        &
    !$omp p )
    do ke=lmesh%NeS, lmesh%NeE
      !---
      dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
      rdens(:) = 1.0_rp / dens(:)
      rhot(:) = dens(:) * pt(:,ke)
 
      ! gradient of density
      call sparsemat_matmul( dx, dens, fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,1), liftdelflx )
      ddensdxi(:,1) = lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:)
 
      call sparsemat_matmul( dy, dens, fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,2), liftdelflx )
      ddensdxi(:,2) = lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:)
 
      call sparsemat_matmul( dz, dens, fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,3), liftdelflx )
      ddensdxi(:,3) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
      
      ! gradient of u
      q(:) = momx_(:,ke) * rdens(:)
 
      call sparsemat_matmul( dx, momx_(:,ke), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,1), liftdelflx )
      dveldxi(:,1,1) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, momx_(:,ke), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,1), liftdelflx )
      dveldxi(:,2,1) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, momx_(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,1), liftdelflx )
      dveldxi(:,3,1) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
      ! gradient of v
      q(:) = momy_(:,ke) * rdens(:)
 
      call sparsemat_matmul( dx, momy_(:,ke), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,2), liftdelflx )
      dveldxi(:,1,2) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, momy_(:,ke), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,2), liftdelflx )
      dveldxi(:,2,2) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, momy_(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,2), liftdelflx )
      dveldxi(:,3,2) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
      ! gradient of w
      q(:) = momz_(:,ke) * rdens(:)
 
      call sparsemat_matmul( dx, momz_(:,ke), fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,3), liftdelflx )
      dveldxi(:,1,3) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, momz_(:,ke), fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,3), liftdelflx )
      dveldxi(:,2,3) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, momz_(:,ke), fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,3), liftdelflx )
      dveldxi(:,3,3) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
      ! gradient of pt
      q(:) = rhot(:) * rdens(:)
 
      call sparsemat_matmul( dx, rhot, fx )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,1), liftdelflx )
      df1(:,ke) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
 
      call sparsemat_matmul( dy, rhot, fy )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,2), liftdelflx )
      df2(:,ke) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
 
      call sparsemat_matmul( dz, rhot, fz )
      call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,3), liftdelflx )
      df3(:,ke) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
 
 
      ! Calculate the component of strain velocity tensor
      s11(:) = dveldxi(:,1,1)
      s12(:) = 0.5_rp * ( dveldxi(:,1,2) + dveldxi(:,2,1) )
      s22(:) = dveldxi(:,2,2)
      s23(:) = 0.5_rp * ( dveldxi(:,2,3) + dveldxi(:,3,2) )
      s31(:) = 0.5_rp * ( dveldxi(:,1,3) + dveldxi(:,3,1) )
      s33(:) = dveldxi(:,3,3)
 
      ! Caclulate eddy viscosity & eddy diffusivity
 
      nu(:,ke) = dns_nu
      kh(:,ke) = dns_mu
 
      !---
 
      do p=1, elem%Np
        skkovthree = ( s11(p) + s22(p) + s33(p) ) * oneoverthree
        coef = 2.0_rp * nu(p,ke)
 
        t11(p,ke) = dens(p) * ( coef * ( s11(p) - skkovthree )  )
        t12(p,ke) = dens(p) * coef * s12(p)
        t13(p,ke) = dens(p) * coef * s31(p)
 
        t21(p,ke) = dens(p) * coef * s12(p)
        t22(p,ke) = dens(p) * ( coef * ( s22(p) - skkovthree )  )
        t23(p,ke) = dens(p) * coef * s23(p)
 
        t31(p,ke) = dens(p) * coef * s31(p)
        t32(p,ke) = dens(p) * coef * s23(p)
        t33(p,ke) = dens(p) * ( coef * ( s33(p) - skkovthree )  )
      end do
 
      df1(:,ke) = kh(:,ke) * df1(:,ke)
      df2(:,ke) = kh(:,ke) * df2(:,ke)
      df3(:,ke) = kh(:,ke) * df3(:,ke)
    end do
  
    return
  end subroutine atm_phy_tb_dgm_dns_cal_grad
 
!OCL SERIAL  
  subroutine cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot,  & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, pt_,           & ! (in)
    nx, ny, nz, vmapm, vmapp, lmesh, elem, is_bound                         ) ! (in)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    real(rp), intent(out) ::  del_flux_rho(elem%nfptot*lmesh%ne,3)
    real(rp), intent(out) ::  del_flux_mom(elem%nfptot*lmesh%ne,3,3)
    real(rp), intent(out) ::  del_flux_rhot(elem%nfptot*lmesh%ne,3)
    real(rp), intent(in) ::  ddens_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  momx_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  momy_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  momz_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  drhot_(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) ::  pt_(elem%np*lmesh%nea)
    real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
    real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
    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)
    logical, intent(in) :: is_bound(elem%nfptot*lmesh%ne)
    
    integer :: i, ip, im
    real(rp) :: densm, densp
    real(rp) :: del
    real(rp) :: facx, facy, facz
    !------------------------------------------------------------------------
    
    !$omp parallel do private ( iM, iP,   &
    !$omp densM, densP,                   &
    !$omp del, facx, facy, facz           )
    do i=1, elem%NfpTot * lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
 
      densm = ddens_(im) + dens_hyd(im)
      densp = ddens_(ip) + dens_hyd(ip)
 
      if ( is_bound(i) ) then
        facx = 1.0_rp
        facy = 1.0_rp 
        facz = 1.0_rp
      else
        ! facx = 1.0_RP - sign(1.0_RP,nx(i))
        ! facy = 1.0_RP - sign(1.0_RP,ny(i))
        ! facz = 1.0_RP - sign(1.0_RP,nz(i))
        facx = 1.0_rp
        facy = 1.0_rp
        facz = 1.0_rp
      end if
 
      del = 0.5_rp * ( densp - densm )
      del_flux_rho(i,1) = facx * del * nx(i)
      del_flux_rho(i,2) = facy * del * ny(i)
      del_flux_rho(i,3) = facz * del * nz(i)
 
      del = 0.5_rp * ( momx_(ip) - momx_(im) )
      del_flux_mom(i,1,1) = facx * del * nx(i)
      del_flux_mom(i,2,1) = facy * del * ny(i)
      del_flux_mom(i,3,1) = facz * del * nz(i)
 
      del = 0.5_rp * ( momy_(ip) - momy_(im) )
      del_flux_mom(i,1,2) = facx * del * nx(i)
      del_flux_mom(i,2,2) = facy * del * ny(i)
      del_flux_mom(i,3,2) = facz * del * nz(i)
 
      del = 0.5_rp * ( momz_(ip) - momz_(im) )
      del_flux_mom(i,1,3) = facx * del * nx(i)
      del_flux_mom(i,2,3) = facy * del * ny(i)
      del_flux_mom(i,3,3) = facz * del * nz(i)
 
      del = 0.5_rp * ( densp * pt_(ip) - densm * pt_(im) )
      del_flux_rhot(i,1) = facx * del * nx(i)
      del_flux_rhot(i,2) = facy * del * ny(i)
      del_flux_rhot(i,3) = facz * del * nz(i)
    end do
 
    return
  end subroutine cal_del_flux_grad
 
!-- private --------------------------------------------------------
 
end module scale_atm_phy_tb_dgm_dns