Initialise a hydrostatically balanced atmosphere with uniform potential temperature theta0, zero velocity, and the geopotential Phi = g*z carried as state variable index 6.
Exner function: pi(z) = 1 - gz / (cptheta0) rho(z) = p0/(Rdtheta0) * pi(z)^(cv/Rd) rhotheta(z) = p0/Rd * pi(z)^(cv/Rd)
Boundary and extBoundary buffers of the solution are populated analytically from per-face z 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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ESAtmo3D_t), | intent(inout) | :: | this | |||
| real(kind=prec), | intent(in) | :: | theta0 |
subroutine SetHydrostaticBalance_ESAtmo3D_t(this,theta0)
!! Initialise a hydrostatically balanced atmosphere with uniform
!! potential temperature theta0, zero velocity, and the geopotential
!! Phi = g*z carried as state variable index 6.
!!
!! Exner function: pi(z) = 1 - g*z / (cp*theta0)
!! rho(z) = p0/(Rd*theta0) * pi(z)^(cv/Rd)
!! rho*theta(z) = p0/Rd * pi(z)^(cv/Rd)
!!
!! Boundary and extBoundary buffers of the solution are populated
!! analytically from per-face z 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(ESAtmo3D_t),intent(inout) :: this
real(prec),intent(in) :: theta0
! Local
integer :: i,j,k,iEl
integer :: side
real(prec) :: z,exner,rho
print*,__FILE__," : Setting hydrostatic balance with theta0 = ",theta0
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
k=1:this%solution%N+1,iEl=1:this%mesh%nElem)
z = this%geometry%x%interior(i,j,k,iEl,1,3)
exner = 1.0_prec-this%g*z/(this%cp*theta0)
rho = this%p0/(this%Rd*theta0)*exner**(this%cv/this%Rd)
this%solution%interior(i,j,k,iEl,1) = rho
this%solution%interior(i,j,k,iEl,2) = 0.0_prec
this%solution%interior(i,j,k,iEl,3) = 0.0_prec
this%solution%interior(i,j,k,iEl,4) = 0.0_prec
this%solution%interior(i,j,k,iEl,5) = rho*theta0
this%solution%interior(i,j,k,iEl,6) = this%g*z
enddo
! Populate boundary + extBoundary buffers analytically from the
! per-face z-coordinates. Geopotential at the boundary mirrors
! by construction (Phi = g*z is single-valued in z).
do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
side=1:6,iEl=1:this%mesh%nElem)
z = this%geometry%x%boundary(i,j,side,iEl,1,3)
exner = 1.0_prec-this%g*z/(this%cp*theta0)
rho = this%p0/(this%Rd*theta0)*exner**(this%cv/this%Rd)
this%solution%boundary(i,j,side,iEl,1) = rho
this%solution%boundary(i,j,side,iEl,2) = 0.0_prec
this%solution%boundary(i,j,side,iEl,3) = 0.0_prec
this%solution%boundary(i,j,side,iEl,4) = 0.0_prec
this%solution%boundary(i,j,side,iEl,5) = rho*theta0
this%solution%boundary(i,j,side,iEl,6) = this%g*z
this%solution%extBoundary(i,j,side,iEl,1) = rho
this%solution%extBoundary(i,j,side,iEl,2) = 0.0_prec
this%solution%extBoundary(i,j,side,iEl,3) = 0.0_prec
this%solution%extBoundary(i,j,side,iEl,4) = 0.0_prec
this%solution%extBoundary(i,j,side,iEl,5) = rho*theta0
this%solution%extBoundary(i,j,side,iEl,6) = this%g*z
enddo
call this%ReportMetrics()
call this%solution%UpdateDevice()
endsubroutine SetHydrostaticBalance_ESAtmo3D_t