hbc2d_NoNormalFlow_ESAtmo2D Subroutine

  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