scale_atm_dyn_dgm_globalnonhydro3d_etot_heve.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVE
!!
!! @par Description
!!      HEVE DGM scheme for Global Atmospheric Dynamical process. 
!!      The governing equations is a fully compressible nonhydrostatic equations, 
!!      which consist of mass, momentum, and thermodynamics (total energy conservation) equations. 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_globalnonhydro3d_etot_heve
  !-----------------------------------------------------------------------------
  !
  !++ 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_hexahedral, only: hexahedralelement
  use scale_element_operation_base, only: elementoperationbase3d
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    atm_dyn_dgm_nonhydro3d_common_init,                       &
    atm_dyn_dgm_nonhydro3d_common_final,                      &
    dens_vid => prgvar_ddens_id, etot_vid => prgvar_etot_id,  &
    momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id,   &
    momz_vid => prgvar_momz_id,                               &
    prgvar_num, intrpmat_vpordm1
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_globalnonhydro3d_etot_heve_init
  public :: atm_dyn_dgm_globalnonhydro3d_etot_heve_final
  public :: atm_dyn_dgm_globalnonhydro3d_etot_heve_cal_tend
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
contains
!OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_etot_heve_init( mesh )
 
    implicit none
    class(meshbase3d), intent(in) :: mesh
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_etot_heve_init
 
!OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_etot_heve_final()
    implicit none
    !--------------------------------------------
    
    call atm_dyn_dgm_nonhydro3d_common_final()    
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_etot_heve_final  
 
  !-------------------------------
 
!OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_etot_heve_cal_tend( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, EnTot_dt,             & ! (out)
    ddens_, momx_, momy_, momz_, etot_, dpres_,               & ! (in)
    dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd,              & ! (in)
    coriolis, rtot, cvtot, cptot, dphyddx, dphyddy,           & ! (in)
    element3d_operation, dx, dy, dz, sx, sy, sz, lift,        & ! (in)
    lmesh, elem, lmesh2d, elem2d )                              ! (in)
 
    use scale_atm_dyn_dgm_nonhydro3d_etot_heve_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_etot_heve_numflux_get_generalhvc
    use scale_const, only: &
      ohm => const_ohm
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    class(elementoperationbase3d), intent(in) :: element3d_operation
    type(sparsemat), intent(in) :: dx, dy, dz, sx, sy, sz, lift
    real(rp), intent(out) :: dens_dt(elem%np,lmesh%nea)
    real(rp), intent(out) :: momx_dt(elem%np,lmesh%nea)
    real(rp), intent(out) :: momy_dt(elem%np,lmesh%nea)
    real(rp), intent(out) :: momz_dt(elem%np,lmesh%nea)
    real(rp), intent(out) :: entot_dt(elem%np,lmesh%nea)
    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)  :: etot_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: dpres_(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)  :: pres_hyd_ref(elem%np,lmesh%nea)
    real(rp), intent(in)  :: therm_hyd(elem%np,lmesh%nea)
    real(rp), intent(in)  :: coriolis(elem2d%np,lmesh2d%nea)
    real(rp), intent(in)  :: rtot (elem%np,lmesh%nea)    
    real(rp), intent(in)  :: cvtot(elem%np,lmesh%nea)
    real(rp), intent(in)  :: cptot(elem%np,lmesh%nea)
    real(rp), intent(in) :: dphyddx(elem%np,lmesh%nea)
    real(rp), intent(in) :: dphyddy(elem%np,lmesh%nea)
 
    real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
    real(rp) :: dpres_hyd(elem%np), gradphyd_x(elem%np), gradphyd_y(elem%np)
    real(rp) :: del_flux(elem%nfptot,lmesh%ne,prgvar_num)
    real(rp) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
    real(rp) :: rdens_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), drho(elem%np)
 
    real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
    real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np)
    real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
    real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
    real(rp) :: cori(elem%np,2)
    logical :: is_panel1to4
    real(rp) :: s
 
    integer :: ke, ke2d
    integer :: p, p12, p3
 
    real(rp) :: gamm, rgamm    
    real(rp) :: rp0
    real(rp) :: rovp0, p0ovr
 
    real(rp) :: enthalpy_(elem%np)
    real(rp) :: u1_(elem%np), u2_(elem%np)
    !------------------------------------------------------------------------
 
    call prof_rapstart('cal_dyn_tend_bndflux', 3)
    call get_ebnd_flux( &
      del_flux, del_flux_hyd,                                                  & ! (out)
      ddens_, momx_, momy_, momz_, etot_, dpres_, dens_hyd, pres_hyd,          & ! (in)
      rtot, cvtot, cptot,                                                      & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%G_ij(:,:,1,1), lmesh%G_ij(:,:,1,2), lmesh%G_ij(:,:,2,2),           & ! (in)
      lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), lmesh%zlev(:,:),       & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),  & ! (in)
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd,                         & ! (in)
      lmesh, elem, lmesh2d, elem2d                                             ) ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
    gamm  = cpdry / cvdry
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    rovp0 = rdry * rp0
    p0ovr = pres00 / rdry
 
    s = 1.0_rp 
    is_panel1to4 = .true.
    if ( lmesh%panelID == 5 ) then
      is_panel1to4 = .false.
    else if ( lmesh%panelID == 6 ) then
      is_panel1to4 = .false.
      s = - 1.0_rp
    end if
 
    !$omp parallel private(                        &
    !$omp rdens_, u_, v_, w_, wt_,                 &
    !$omp enthalpy_, u1_, u2_,                     &
    !$omp Fx, Fy, Fz, LiftDelFlx,                  &
    !$omp drho, DPRES_hyd, GradPhyd_x, GradPhyd_y, &
    !$omp G11, G12, G22, GsqrtV, RGsqrtV,          &
    !$omp X, Y, twoOVdel2,                         &
    !$omp CORI, ke, ke2D                           )
 
    !$omp do
    do ke2d = lmesh2d%NeS, lmesh2d%NeE
      x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
      y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
    end do
 
    !$omp do
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
      g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1)
      g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2)
      g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2)
      gsqrtv(:)  = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
      rgsqrtv(:) = 1.0_rp / gsqrtv(:)
 
      !--
      rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
      u_(:) = momx_(:,ke) * rdens_(:)
      v_(:) = momy_(:,ke) * rdens_(:)
      w_(:) = momz_(:,ke) * rdens_(:)
      wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:) 
      u1_(:) = lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,1,1) * u_(:) + lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,2,1) * v_(:)
      u2_(:) = lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,2,1) * u_(:) + lmesh%G_ij(elem%IndexH2Dto3D(:),ke2d,2,2) * v_(:)
 
      ! DPRES_(:) = ( CPtot(:,ke) / CVtot(:,ke) - 1.0_RP ) &
      !           * (   ETOT_(:,ke) - Grav * ( DDENS_(:,ke) + DENS_hyd(:,ke) ) * lmesh%zlev(:,ke)        &
      !               - 0.5_RP * ( MOMX_(:,ke) * u1_(:) + MOMY_(:,ke) * u2_(:) + MOMZ_(:,ke) * w_(:) ) ) &
      !           - PRES_hyd(:,ke)
      
      enthalpy_(:) = etot_(:,ke) + pres_hyd(:,ke) + dpres_(:,ke)
 
      x(:) = x2d(elem%IndexH2Dto3D,ke2d)
      y(:) = y2d(elem%IndexH2Dto3D,ke2d)
      twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
 
      cori(:,1) = s * ohm * twoovdel2(:) * ( - x(:) * y(:)        * momx_(:,ke) + (1.0_rp + y(:)**2) * momy_(:,ke) )
      cori(:,2) = s * ohm * twoovdel2(:) * ( - (1.0_rp + x(:)**2) * momx_(:,ke) +  x(:) * y(:)       * momy_(:,ke) )
      
      if ( is_panel1to4 ) then
        cori(:,1) = s * y(:) * cori(:,1)
        cori(:,2) = s * y(:) * cori(:,2)
      end if
 
      drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
 
      !-- Gradient hydrostatic pressure
      
      dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
 
      call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
      gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
 
      call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
      gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
      
      !-- DENS
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
 
      dens_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx(:)     &
          + lmesh%Escale(:,ke,2,2) * fy(:)     &
          + lmesh%Escale(:,ke,3,3) * fz(:)     &
          + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
      
      !-- MOMX
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
 
      momx_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)                                         &
            + lmesh%Escale(:,ke,2,2) * fy(:)                                         &
            + lmesh%Escale(:,ke,3,3) * fz(:)                                         &
            + liftdelflx(:)                   ) / lmesh%Gsqrt(:,ke)                  &
          + twoovdel2(:) * y(:) *                                                    &
            ( - x(:) * y(:) * u_(:) + (1.0_rp + y(:)**2) * v_(:) ) * momx_(:,ke)     &
          - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:)         &
          + cori(:,1)              
 
      !-- MOMY
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
 
      momy_dt(:,ke) = &
            - ( lmesh%Escale(:,ke,1,1) * fx(:)                                         &
              + lmesh%Escale(:,ke,2,2) * fy(:)                                         &
              + lmesh%Escale(:,ke,3,3) * fz(:)                                         &
              + liftdelflx(:)                  ) / lmesh%Gsqrt(:,ke)                   &
            + twoovdel2(:) * x(:) *                                                    &
              ( (1.0_rp + x(:)**2) * u_(:) - x(:) * y(:) * v_(:) ) * momy_(:,ke)       &
            - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:)         &
            + cori(:,2) 
 
      !-- MOMZ
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) *   u_(:) * momz_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) *   v_(:) * momz_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
      
      momz_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)       &
            + lmesh%Escale(:,ke,2,2) * fy(:)       &
            + lmesh%Escale(:,ke,3,3) * fz(:)       &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)  &
          - grav * drho(:)
 
      !-- EnTot
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * enthalpy_(:), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * enthalpy_(:), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * enthalpy_(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,etot_vid), liftdelflx)
      
      entot_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)      &
            + lmesh%Escale(:,ke,2,2) * fy(:)      &
            + lmesh%Escale(:,ke,3,3) * fz(:)      &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) 
    end do
    !$omp end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_etot_heve_cal_tend
 
end module scale_atm_dyn_dgm_globalnonhydro3d_etot_heve