+ * delta_{id}
!! f_d(rho*th) = _log *
!! f_d(Phi) = 0 (geopotential is carried, not advected)
!!
!! Pressure here is the *total* pressure; gravity is handled by the
!! Souza non-conservative term in SourceMethod, not by a flux split.
class(ESAtmo2D_t),intent(in) :: this
real(prec),intent(in) :: sL(1:this%nvar)
real(prec),intent(in) :: sR(1:this%nvar)
real(prec) :: flux(1:this%nvar,1:2)
! Local
real(prec) :: rhoL,uL,vL,thetaL,pL,rthL
real(prec) :: rhoR,uR,vR,thetaR,pR,rthR
real(prec) :: rho_log,rth_log,u_avg,v_avg,p_avg
real(prec) :: gamma
gamma = this%cp/this%cv
! Left primitive variables
rhoL = sL(1)
rthL = sL(4)
thetaL = rthL/rhoL
uL = sL(2)/rhoL
vL = sL(3)/rhoL
pL = this%p0*(rhoL*this%Rd*thetaL/this%p0)**gamma
! Right primitive variables
rhoR = sR(1)
rthR = sR(4)
thetaR = rthR/rhoR
uR = sR(2)/rhoR
vR = sR(3)/rhoR
pR = this%p0*(rhoR*this%Rd*thetaR/this%p0)**gamma
! Logarithmic and arithmetic means
rho_log = log_mean(rhoL,rhoR)
rth_log = log_mean(rthL,rthR)
u_avg = 0.5_prec*(uL+uR)
v_avg = 0.5_prec*(vL+vR)
p_avg = 0.5_prec*(pL+pR)
! x-direction flux (d=1)
flux(1,1) = rho_log*u_avg
flux(2,1) = rho_log*u_avg*u_avg+p_avg
flux(3,1) = rho_log*v_avg*u_avg
flux(4,1) = rth_log*u_avg
flux(5,1) = 0.0_prec
! y-direction flux (d=2)
flux(1,2) = rho_log*v_avg
flux(2,2) = rho_log*u_avg*v_avg
flux(3,2) = rho_log*v_avg*v_avg+p_avg
flux(4,2) = rth_log*v_avg
flux(5,2) = 0.0_prec
endfunction twopointflux2d_ESAtmo2D_t
subroutine TwoPointFluxMethod_ESAtmo2D_t(this)
!! Pre-projected scalar contravariant two-point Souza et al. (2023) EC
!! flux. For each node pair (a, b) along reference direction r:
!!
!! Fc^r_v_i = sum_d 0.5*(Ja^r_d(a) + Ja^r_d(b)) * f_d(v_i, sL=s(a), sR=s(b))
!!
!! Gravity is NOT split into the pressure flux here; it is handled by
!! the Souza non-conservative gravity flux differencing in SourceMethod
!! (which uses the geopotential carried as state variable index 5).
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: nn,i,j,d,iEl,iVar
real(prec) :: sL(1:this%nvar),sR(1:this%nvar)
real(prec) :: Fphys(1:this%nvar,1:2)
real(prec) :: Fc
do concurrent(nn=1:this%solution%N+1,i=1:this%solution%N+1, &
j=1:this%solution%N+1,iEl=1:this%mesh%nElem)
sL = this%solution%interior(i,j,iEl,1:this%nvar)
! ------------- xi^1: pair (i,j)-(nn,j) -------------
sR = this%solution%interior(nn,j,iEl,1:this%nvar)
Fphys = this%twopointflux2d(sL,sR)
do iVar = 1,this%nvar
Fc = 0.0_prec
do d = 1,2
Fc = Fc+0.5_prec*( &
this%geometry%dsdx%interior(i,j,iEl,1,d,1)+ &
this%geometry%dsdx%interior(nn,j,iEl,1,d,1))* &
Fphys(iVar,d)
enddo
this%twoPointFlux%interior(nn,i,j,iEl,iVar,1) = Fc
enddo
! ------------- xi^2: pair (i,j)-(i,nn) -------------
sR = this%solution%interior(i,nn,iEl,1:this%nvar)
Fphys = this%twopointflux2d(sL,sR)
do iVar = 1,this%nvar
Fc = 0.0_prec
do d = 1,2
Fc = Fc+0.5_prec*( &
this%geometry%dsdx%interior(i,j,iEl,1,d,2)+ &
this%geometry%dsdx%interior(i,nn,iEl,1,d,2))* &
Fphys(iVar,d)
enddo
this%twoPointFlux%interior(nn,i,j,iEl,iVar,2) = Fc
enddo
enddo
endsubroutine TwoPointFluxMethod_ESAtmo2D_t
pure function riemannflux2d_ESAtmo2D_t(this,sL,sR,dsdx,nhat) result(flux)
!! Local Lax-Friedrichs (Rusanov) Riemann flux. Provided as a fallback;
!! the model overrides BoundaryFlux directly with the LMARS solver.
!!
!! F* = 0.5*(fL.n + fR.n) - 0.5*lambda_max*(sR - sL)
!!
!! where lambda_max = max(|vL.n| + cL, |vR.n| + cR)
!! and c = sqrt(gamma * p / rho) is the sound speed.
class(ESAtmo2D_t),intent(in) :: this
real(prec),intent(in) :: sL(1:this%nvar)
real(prec),intent(in) :: sR(1:this%nvar)
real(prec),intent(in) :: dsdx(1:this%nvar,1:2)
real(prec),intent(in) :: nhat(1:2)
real(prec) :: flux(1:this%nvar)
! Local
real(prec) :: rhoL,uL,vL,thetaL,pL,unL,cL
real(prec) :: rhoR,uR,vR,thetaR,pR,unR,cR
real(prec) :: fL(1:4),fR(1:4)
real(prec) :: gamma,lam
gamma = this%cp/this%cv
! Left state
rhoL = sL(1)
uL = sL(2)/rhoL
vL = sL(3)/rhoL
thetaL = sL(4)/rhoL
pL = this%p0*(rhoL*this%Rd*thetaL/this%p0)**gamma
unL = uL*nhat(1)+vL*nhat(2)
cL = sqrt(gamma*pL/rhoL)
! Right state
rhoR = sR(1)
uR = sR(2)/rhoR
vR = sR(3)/rhoR
thetaR = sR(4)/rhoR
pR = this%p0*(rhoR*this%Rd*thetaR/this%p0)**gamma
unR = uR*nhat(1)+vR*nhat(2)
cR = sqrt(gamma*pR/rhoR)
! Normal flux from left state
fL(1) = rhoL*unL
fL(2) = sL(2)*unL+pL*nhat(1)
fL(3) = sL(3)*unL+pL*nhat(2)
fL(4) = sL(4)*unL
! Normal flux from right state
fR(1) = rhoR*unR
fR(2) = sR(2)*unR+pR*nhat(1)
fR(3) = sR(3)*unR+pR*nhat(2)
fR(4) = sR(4)*unR
! Maximum wave speed
lam = max(abs(unL)+cL,abs(unR)+cR)
! LLF flux
flux(1:4) = 0.5_prec*(fL(1:4)+fR(1:4))- &
0.5_prec*lam*(sR(1:4)-sL(1:4))
! Geopotential is carried in the state vector but has no flux.
flux(5) = 0.0_prec
if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning
endfunction riemannflux2d_ESAtmo2D_t
subroutine BoundaryFlux_ESAtmo2D_t(this)
!! LMARS (Low-Mach Approximate Riemann Solver, Chen et al. 2013)
!! interface flux. No hydrostatic pressure split: gravity is folded
!! into SourceMethod via the Souza non-conservative form using the
!! geopotential carried in the state vector (variable index 5).
!!
!! un* = 0.5*(unL + unR) - (pR - pL) / (2 * rho_bar * c_bar)
!! p* = 0.5*(pL + pR) - 0.5 * rho_bar * c_bar * (unR - unL)
!!
!! Upwind on sign of un* selects which side supplies the conserved
!! state for the advective part. Pressure adds to normal momentum.
!! Geopotential (var 5) has zero flux.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: i,k,iel
real(prec) :: nhat(1:2),nmag
real(prec) :: rhoL,uL,vL,thetaL,pL,unL,cL
real(prec) :: rhoR,uR,vR,thetaR,pR,unR,cR
real(prec) :: rho_bar,c_bar,rc,un_star,p_star,gamma
real(prec) :: sL(1:4),sR(1:4),s_up(1:4)
gamma = this%cp/this%cv
do concurrent(i=1:this%solution%N+1,k=1:4,iel=1:this%mesh%nElem)
nhat = this%geometry%nHat%boundary(i,k,iEl,1,1:2)
nmag = this%geometry%nScale%boundary(i,k,iEl,1)
sL = this%solution%boundary(i,k,iel,1:4)
sR = this%solution%extBoundary(i,k,iel,1:4)
rhoL = sL(1)
uL = sL(2)/rhoL
vL = sL(3)/rhoL
thetaL = sL(4)/rhoL
pL = this%p0*(rhoL*this%Rd*thetaL/this%p0)**gamma
unL = uL*nhat(1)+vL*nhat(2)
cL = sqrt(gamma*pL/rhoL)
rhoR = sR(1)
uR = sR(2)/rhoR
vR = sR(3)/rhoR
thetaR = sR(4)/rhoR
pR = this%p0*(rhoR*this%Rd*thetaR/this%p0)**gamma
unR = uR*nhat(1)+vR*nhat(2)
cR = sqrt(gamma*pR/rhoR)
rho_bar = 0.5_prec*(rhoL+rhoR)
c_bar = 0.5_prec*(cL+cR)
rc = rho_bar*c_bar
un_star = 0.5_prec*(unL+unR)-(pR-pL)/(2.0_prec*rc)
p_star = 0.5_prec*(pL+pR)-0.5_prec*rc*(unR-unL)
! Upwind on un_star (advective part of LMARS)
if(un_star >= 0.0_prec) then
s_up = sL
else
s_up = sR
endif
this%flux%boundaryNormal(i,k,iEl,1) = (s_up(1)*un_star)*nmag
this%flux%boundaryNormal(i,k,iEl,2) = (s_up(2)*un_star+p_star*nhat(1))*nmag
this%flux%boundaryNormal(i,k,iEl,3) = (s_up(3)*un_star+p_star*nhat(2))*nmag
this%flux%boundaryNormal(i,k,iEl,4) = (s_up(4)*un_star)*nmag
this%flux%boundaryNormal(i,k,iEl,5) = 0.0_prec
enddo
endsubroutine BoundaryFlux_ESAtmo2D_t
subroutine SourceMethod_ESAtmo2D_t(this)
!! Souza et al. (2023) non-conservative gravity flux differencing.
!!
!! The rho*v equation carries the body force -rho * partial_y Phi where
!! Phi = g*y is the geopotential (state variable index 5; no flux, no
!! source of its own). On a curvilinear mesh, the SBP-EC strong-form
!! flux differencing for the non-conservative product rho * partial_y Phi
!! reads:
!!
!! [rho * d_y Phi]_i = (1/J_i) sum_r sum_j D_split^r[i_r, j]
!! * 0.5*(Ja^r_y(i) + Ja^r_y(j))
!! * _log(s_i, s_j) * (Phi_j - Phi_i)
!!
!! and we set source(rho*v) = - that result. The other variables
!! (rho, rho*u, rho*theta, Phi) all have zero source.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: i,j,iEl,nn
real(prec) :: rho_ijk,phi_ijk,rho_p,phi_p
real(prec) :: rho_log,Ja_avg,acc,jac
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iEl=1:this%mesh%nElem)
rho_ijk = this%solution%interior(i,j,iEl,1)
phi_ijk = this%solution%interior(i,j,iEl,5)
acc = 0.0_prec
do nn = 1,this%solution%N+1
! xi^1 partner (nn, j) — uses Ja^1_y (d=2, r=1)
rho_p = this%solution%interior(nn,j,iEl,1)
phi_p = this%solution%interior(nn,j,iEl,5)
rho_log = log_mean(rho_ijk,rho_p)
Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,iEl,1,2,1)+ &
this%geometry%dsdx%interior(nn,j,iEl,1,2,1))
acc = acc+this%solution%interp%dSplitMatrix(nn,i)* &
Ja_avg*rho_log*(phi_p-phi_ijk)
! xi^2 partner (i, nn) — uses Ja^2_y (d=2, r=2)
rho_p = this%solution%interior(i,nn,iEl,1)
phi_p = this%solution%interior(i,nn,iEl,5)
rho_log = log_mean(rho_ijk,rho_p)
Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,iEl,1,2,2)+ &
this%geometry%dsdx%interior(i,nn,iEl,1,2,2))
acc = acc+this%solution%interp%dSplitMatrix(nn,j)* &
Ja_avg*rho_log*(phi_p-phi_ijk)
enddo
jac = this%geometry%J%interior(i,j,iEl,1)
this%source%interior(i,j,iEl,1) = 0.0_prec
this%source%interior(i,j,iEl,2) = 0.0_prec
this%source%interior(i,j,iEl,3) = -acc/jac
this%source%interior(i,j,iEl,4) = 0.0_prec
this%source%interior(i,j,iEl,5) = 0.0_prec
enddo
endsubroutine SourceMethod_ESAtmo2D_t
subroutine pbc2d_NoStress_ESAtmo2D(bc,mymodel)
!! Parabolic boundary condition: zero diffusive flux normal to the wall
!! (no-stress for momentum, no-heat-flux for rho*theta).
!!
!! Reflects the normal component of the interior solution gradient:
!!
!! grad_ext = grad_int - 2 * (grad_int . n) * n
!!
!! After AverageSides() this gives avgGrad . n = 0 at every wall node,
!! so the BR1 diffusive boundary flux f^diff = -coeff*(avgGrad.n)*nmag
!! vanishes identically on no-normal-flow walls.
class(BoundaryCondition),intent(in) :: bc
class(Model),intent(inout) :: mymodel
! Local
integer :: n,i,iEl,k,iVar
real(prec) :: nhat(1:2),gn,g(1:2)
select type(m => mymodel)
class is(ESAtmo2D_t)
do n = 1,bc%nBoundaries
iEl = bc%elements(n)
k = bc%sides(n)
do i = 1,m%solution%interp%N+1
nhat = m%geometry%nHat%boundary(i,k,iEl,1,1:2)
do iVar = 1,m%nvar
g(1:2) = m%solutionGradient%boundary(i,k,iEl,iVar,1:2)
gn = g(1)*nhat(1)+g(2)*nhat(2)
m%solutionGradient%extBoundary(i,k,iEl,iVar,1) = g(1)-2.0_prec*gn*nhat(1)
m%solutionGradient%extBoundary(i,k,iEl,iVar,2) = g(2)-2.0_prec*gn*nhat(2)
enddo
enddo
enddo
endselect
endsubroutine pbc2d_NoStress_ESAtmo2D
subroutine AdditionalInit_ESAtmo2D_t(this)
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
procedure(SELF_bcMethod),pointer :: bcfunc
bcfunc => hbc2d_NoNormalFlow_ESAtmo2D
call this%hyperbolicBCs%RegisterBoundaryCondition( &
SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc)
! Parabolic BC for the same wall tag: zero diffusive normal flux
! (no-stress / no-heat-flux). hyperbolicBCs and parabolicBCs are
! independent linked lists, so registering the same tag here does
! not clobber the hyperbolic registration above.
bcfunc => pbc2d_NoStress_ESAtmo2D
call this%parabolicBCs%RegisterBoundaryCondition( &
SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc)
! Diffusive-flux scratch buffers. Always allocated (memory cost is
! modest); only used when nu>0 or kappa>0 (and gradient_enabled is
! therefore .true.).
call this%diffFlux%Init(this%solution%interp,this%nvar,this%mesh%nElem)
call this%diffFlux%AssociateGeometry(this%geometry)
call this%diffDiv%Init(this%solution%interp,this%nvar,this%mesh%nElem)
endsubroutine AdditionalInit_ESAtmo2D_t
subroutine AdditionalFree_ESAtmo2D_t(this)
implicit none
class(ESAtmo2D_t),intent(inout) :: this
call this%diffFlux%Free()
call this%diffDiv%Free()
endsubroutine AdditionalFree_ESAtmo2D_t
subroutine SetDiffusion_ESAtmo2D_t(this,nu,kappa,eta_penalty)
!! Set the constant-coefficient Laplacian diffusion coefficients
!! (kinematic momentum diffusivity and thermal diffusivity, both
!! in m^2/s) and the dimensionless SIPG jump penalty. Setting nu
!! or kappa > 0 enables the gradient pipeline so that the diffusive
!! flux methods receive solutionGradient.
!!
!! length_scale is computed from the area Jacobian: for a quad
!! reference cell [-1,1]^2 with area 4, the physical element area
!! is 4* so the characteristic edge length is 2*^(1/2).
!! The smallest length over the mesh is conservative.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
real(prec),intent(in) :: nu
real(prec),intent(in) :: kappa
real(prec),intent(in),optional :: eta_penalty
! Local
real(prec) :: jmin
this%nu = nu
this%kappa = kappa
if(present(eta_penalty)) this%eta_penalty = eta_penalty
if(nu > 0.0_prec .or. kappa > 0.0_prec) then
this%gradient_enabled = .true.
endif
jmin = minval(this%geometry%J%interior)
this%length_scale = 2.0_prec*sqrt(jmin)
endsubroutine SetDiffusion_ESAtmo2D_t
subroutine DiffusiveFluxMethod_ESAtmo2D_t(this)
!! Fill diffFlux%interior with the constant-coefficient Laplacian
!! flux at every interior node:
!!
!! F_d(rho) = 0 (no mass diffusion)
!! F_d(rho*v_i) = -nu * d(rho*v_i)/dx_d
!! F_d(rho*theta)= -kappa * d(rho*theta)/dx_d
!! F_d(Phi) = 0
!!
!! solutionGradient%interior(i,j,iel,iVar,d) holds d(s_iVar)/dx_d.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: i,j,iel,d
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iel=1:this%mesh%nElem,d=1:2)
this%diffFlux%interior(i,j,iel,1,d) = 0.0_prec
this%diffFlux%interior(i,j,iel,2,d) = &
-this%nu*this%solutionGradient%interior(i,j,iel,2,d)
this%diffFlux%interior(i,j,iel,3,d) = &
-this%nu*this%solutionGradient%interior(i,j,iel,3,d)
this%diffFlux%interior(i,j,iel,4,d) = &
-this%kappa*this%solutionGradient%interior(i,j,iel,4,d)
! Geopotential carries no diffusion.
this%diffFlux%interior(i,j,iel,5,d) = 0.0_prec
enddo
endsubroutine DiffusiveFluxMethod_ESAtmo2D_t
subroutine DiffusiveBoundaryFlux_ESAtmo2D_t(this)
!! Fill diffFlux%boundaryNormal with the SIPG-stabilised BR1 flux:
!!
!! f_R^diff(iVar) = -coeff(iVar) * (avg_grad . n) * nmag
!! + tau(iVar) * (uL - uR) * nmag
!!
!! tau(iVar) = eta_penalty * coeff(iVar) * (N+1)^2 / length_scale.
!! avg_grad is solutionGradient%avgBoundary (populated by
!! AverageSides()); uL, uR are solution%boundary, %extBoundary.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: i,k,iel
real(prec) :: nhat(1:2),nmag,gradn,coeff,tau,np2
real(prec) :: uL,uR
np2 = real((this%solution%interp%N+1)**2,prec)
do concurrent(i=1:this%solution%N+1,k=1:4,iel=1:this%mesh%nElem)
nhat = this%geometry%nHat%boundary(i,k,iel,1,1:2)
nmag = this%geometry%nScale%boundary(i,k,iel,1)
! rho — no mass diffusion, no penalty
this%diffFlux%boundaryNormal(i,k,iel,1) = 0.0_prec
! Momentum equations use nu
coeff = this%nu
tau = this%eta_penalty*coeff*np2/this%length_scale
gradn = this%solutionGradient%avgBoundary(i,k,iel,2,1)*nhat(1)+ &
this%solutionGradient%avgBoundary(i,k,iel,2,2)*nhat(2)
uL = this%solution%boundary(i,k,iel,2)
uR = this%solution%extBoundary(i,k,iel,2)
this%diffFlux%boundaryNormal(i,k,iel,2) = (-coeff*gradn+tau*(uL-uR))*nmag
gradn = this%solutionGradient%avgBoundary(i,k,iel,3,1)*nhat(1)+ &
this%solutionGradient%avgBoundary(i,k,iel,3,2)*nhat(2)
uL = this%solution%boundary(i,k,iel,3)
uR = this%solution%extBoundary(i,k,iel,3)
this%diffFlux%boundaryNormal(i,k,iel,3) = (-coeff*gradn+tau*(uL-uR))*nmag
! rho*theta uses kappa
coeff = this%kappa
tau = this%eta_penalty*coeff*np2/this%length_scale
gradn = this%solutionGradient%avgBoundary(i,k,iel,4,1)*nhat(1)+ &
this%solutionGradient%avgBoundary(i,k,iel,4,2)*nhat(2)
uL = this%solution%boundary(i,k,iel,4)
uR = this%solution%extBoundary(i,k,iel,4)
this%diffFlux%boundaryNormal(i,k,iel,4) = (-coeff*gradn+tau*(uL-uR))*nmag
! Geopotential — no diffusion, no penalty.
this%diffFlux%boundaryNormal(i,k,iel,5) = 0.0_prec
enddo
endsubroutine DiffusiveBoundaryFlux_ESAtmo2D_t
subroutine CalculateTendency_ESAtmo2D_t(this)
!! ESAtmo2D tendency = EC inviscid pipeline (parent) + optional
!! constant-coefficient Laplacian diffusion (BR1 weak-form DG).
!!
!! When nu>0 or kappa>0 the parent's CalculateTendency already
!! computes solutionGradient (because gradient_enabled was set by
!! SetDiffusion); we then fill diffFlux from solutionGradient,
!! compute its DG divergence, and accumulate into fluxDivergence
!! before forming dSdt.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
! Local
integer :: i,j,iEl,iVar
real(prec) :: bMi1,bMi2,bMj1,bMj2,qwi,qwj,jac
call this%solution%BoundaryInterp()
call this%solution%SideExchange(this%mesh)
call this%PreTendency()
call this%SetBoundaryCondition()
if(this%gradient_enabled) then
call this%CalculateSolutionGradient()
call this%SetGradientBoundaryCondition()
call this%solutionGradient%AverageSides()
endif
call this%SourceMethod()
call this%BoundaryFlux()
call this%TwoPointFluxMethod()
call this%twoPointFlux%UpdateDevice()
call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior)
! EC surface contribution (1/J) M^{-1} B^T f_Riemann
! Side ordering: 1=South, 2=East, 3=North, 4=West
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)
bMi1 = this%solution%interp%bMatrix(i,2) ! East boundary value at node i
bMi2 = this%solution%interp%bMatrix(i,1) ! West boundary value at node i
bMj1 = this%solution%interp%bMatrix(j,2) ! North boundary value at node j
bMj2 = this%solution%interp%bMatrix(j,1) ! South boundary value at node j
qwi = this%solution%interp%qWeights(i)
qwj = this%solution%interp%qWeights(j)
jac = this%geometry%J%interior(i,j,iEl,1)
this%fluxDivergence%interior(i,j,iEl,iVar) = &
this%fluxDivergence%interior(i,j,iEl,iVar)+ &
(bMi1*this%flux%boundaryNormal(j,2,iEl,iVar)+ &
bMi2*this%flux%boundaryNormal(j,4,iEl,iVar))/(qwi*jac)+ &
(bMj1*this%flux%boundaryNormal(i,3,iEl,iVar)+ &
bMj2*this%flux%boundaryNormal(i,1,iEl,iVar))/(qwj*jac)
enddo
! Add the diffusive contribution if any coefficient is nonzero.
if(this%nu > 0.0_prec .or. this%kappa > 0.0_prec) then
call this%DiffusiveFluxMethod()
call this%DiffusiveBoundaryFlux()
call this%diffFlux%MappedDGDivergence(this%diffDiv%interior)
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)
this%fluxDivergence%interior(i,j,iEl,iVar) = &
this%fluxDivergence%interior(i,j,iEl,iVar)+ &
this%diffDiv%interior(i,j,iEl,iVar)
enddo
endif
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)
this%dSdt%interior(i,j,iEl,iVar) = &
this%source%interior(i,j,iEl,iVar)- &
this%fluxDivergence%interior(i,j,iEl,iVar)
enddo
endsubroutine CalculateTendency_ESAtmo2D_t
subroutine hbc2d_NoNormalFlow_ESAtmo2D(bc,mymodel)
!! No-normal-flow (wall) boundary condition.
!!
!! The normal component of momentum is negated while tangential
!! components are mirrored (preserved). Density and potential
!! temperature (tracers) are mirrored from the interior.
!!
!! Exterior momentum:
!! (rho*v)_ext = (rho*v)_int - 2 * ((rho*v)_int . nhat) * nhat
!!
!! This ensures v_ext . n = -(v_int . n) so the Riemann solver
!! sees zero normal velocity at the wall.
class(BoundaryCondition),intent(in) :: bc
class(Model),intent(inout) :: mymodel
! Local
integer :: n,i,iEl,k
real(prec) :: nhat(1:2),rhovn
select type(m => mymodel)
class is(ESAtmo2D_t)
do n = 1,bc%nBoundaries
iEl = bc%elements(n)
k = bc%sides(n)
do i = 1,m%solution%interp%N+1
nhat = m%geometry%nHat%boundary(i,k,iEl,1,1:2)
! Mirror density (tracer)
m%solution%extBoundary(i,k,iEl,1) = &
m%solution%boundary(i,k,iEl,1)
! Reflect momentum: negate normal component, preserve tangential
! (rho*v)_ext = (rho*v)_int - 2*((rho*v)_int . n)*n
rhovn = m%solution%boundary(i,k,iEl,2)*nhat(1)+ &
m%solution%boundary(i,k,iEl,3)*nhat(2)
m%solution%extBoundary(i,k,iEl,2) = &
m%solution%boundary(i,k,iEl,2)-2.0_prec*rhovn*nhat(1)
m%solution%extBoundary(i,k,iEl,3) = &
m%solution%boundary(i,k,iEl,3)-2.0_prec*rhovn*nhat(2)
! Mirror potential temperature (tracer)
m%solution%extBoundary(i,k,iEl,4) = &
m%solution%boundary(i,k,iEl,4)
! Mirror geopotential — Phi only depends on y, so the
! mirror across a no-normal-flow wall is the same value.
m%solution%extBoundary(i,k,iEl,5) = &
m%solution%boundary(i,k,iEl,5)
enddo
enddo
endselect
endsubroutine hbc2d_NoNormalFlow_ESAtmo2D
subroutine SetHydrostaticBalance_ESAtmo2D_t(this,theta0)
!! Initialise a hydrostatically balanced atmosphere with uniform
!! potential temperature theta0, zero velocity, and the geopotential
!! Phi = g*y carried as state variable index 5.
!!
!! Exner function: pi(y) = 1 - g*y / (cp*theta0)
!! rho(y) = p0/(Rd*theta0) * pi(y)^(cv/Rd)
!! rho*theta(y) = p0/Rd * pi(y)^(cv/Rd)
!!
!! Boundary and extBoundary buffers of the solution are populated
!! analytically from per-face y so that face values are bit-exact
!! with what BoundaryInterp would produce (and extBoundary at walls
!! mirrors). SideExchange will overwrite extBoundary at interior
!! element interfaces with the neighbour's value — for a smooth
!! profile this is a no-op.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
real(prec),intent(in) :: theta0
! Local
integer :: i,j,iEl
integer :: side
real(prec) :: y,exner,rho
print*,__FILE__," : Setting hydrostatic balance with theta0 = ",theta0
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iEl=1:this%mesh%nElem)
y = this%geometry%x%interior(i,j,iEl,1,2)
exner = 1.0_prec-this%g*y/(this%cp*theta0)
rho = this%p0/(this%Rd*theta0)*exner**(this%cv/this%Rd)
this%solution%interior(i,j,iEl,1) = rho
this%solution%interior(i,j,iEl,2) = 0.0_prec
this%solution%interior(i,j,iEl,3) = 0.0_prec
this%solution%interior(i,j,iEl,4) = rho*theta0
this%solution%interior(i,j,iEl,5) = this%g*y
enddo
! Populate boundary + extBoundary buffers analytically from the
! per-face y-coordinates. Geopotential at the boundary mirrors
! by construction (Phi = g*y is single-valued in y).
do concurrent(i=1:this%solution%N+1,side=1:4,iEl=1:this%mesh%nElem)
y = this%geometry%x%boundary(i,side,iEl,1,2)
exner = 1.0_prec-this%g*y/(this%cp*theta0)
rho = this%p0/(this%Rd*theta0)*exner**(this%cv/this%Rd)
this%solution%boundary(i,side,iEl,1) = rho
this%solution%boundary(i,side,iEl,2) = 0.0_prec
this%solution%boundary(i,side,iEl,3) = 0.0_prec
this%solution%boundary(i,side,iEl,4) = rho*theta0
this%solution%boundary(i,side,iEl,5) = this%g*y
this%solution%extBoundary(i,side,iEl,1) = rho
this%solution%extBoundary(i,side,iEl,2) = 0.0_prec
this%solution%extBoundary(i,side,iEl,3) = 0.0_prec
this%solution%extBoundary(i,side,iEl,4) = rho*theta0
this%solution%extBoundary(i,side,iEl,5) = this%g*y
enddo
call this%ReportMetrics()
call this%solution%UpdateDevice()
endsubroutine SetHydrostaticBalance_ESAtmo2D_t
subroutine AddThermalBubble_ESAtmo2D_t(this,dtheta,r0,x0,y0)
!! Adds a pressure-balanced warm bubble perturbation.
!!
!! The potential temperature perturbation has a cos^2 profile:
!!
!! theta'(x,y) = dtheta * cos^2(pi*r / (2*r0)) for r <= r0
!! theta'(x,y) = 0 for r > r0
!!
!! where r = sqrt((x-x0)^2 + (y-y0)^2).
!!
!! To maintain pressure balance (p unchanged), the density is
!! adjusted while rho*theta remains constant:
!!
!! theta_new = theta_old + theta'
!! rho_new = rho_old * theta_old / theta_new
!! (rho*theta)_new = rho_old * theta_old = (rho*theta)_old
!!
!! The buoyancy force arises from the density deficit in the
!! gravitational source term -rho*g.
implicit none
class(ESAtmo2D_t),intent(inout) :: this
real(prec),intent(in) :: dtheta ! Perturbation amplitude [K]
real(prec),intent(in) :: r0 ! Bubble radius [m]
real(prec),intent(in) :: x0,y0 ! Bubble center [m]
! Local
integer :: i,j,iEl
real(prec) :: x,y,r,rho_old,theta_old,theta_new,thetap
print*,__FILE__," : Adding thermal bubble perturbation"
print*,__FILE__," : dtheta = ",dtheta
print*,__FILE__," : r0 = ",r0
print*,__FILE__," : center = ",x0,y0
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
iEl=1:this%mesh%nElem)
x = this%geometry%x%interior(i,j,iEl,1,1)-x0
y = this%geometry%x%interior(i,j,iEl,1,2)-y0
r = sqrt(x*x+y*y)
if(r <= r0) then
rho_old = this%solution%interior(i,j,iEl,1)
theta_old = this%solution%interior(i,j,iEl,4)/rho_old
thetap = dtheta*cos(pi*r/(2.0_prec*r0))**2
theta_new = theta_old+thetap
! Adjust density to maintain pressure balance: rho*theta = const
this%solution%interior(i,j,iEl,1) = rho_old*theta_old/theta_new
! rho*theta is unchanged (pressure-balanced perturbation)
endif
enddo
call this%ReportMetrics()
call this%solution%UpdateDevice()
endsubroutine AddThermalBubble_ESAtmo2D_t
endmodule SELF_ESAtmo2D_t