SetHydrostaticBalance_ESAtmo2D_t Subroutine

  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