hbc3d_NoNormalFlow_ESAtmo3D Subroutine

  subroutine hbc3d_NoNormalFlow_ESAtmo3D(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,j,iEl,k
    real(prec) :: nhat(1:3),rhovn

    select type(m => mymodel)
    class is(ESAtmo3D_t)
      do n = 1,bc%nBoundaries
        iEl = bc%elements(n)
        k = bc%sides(n)
        do j = 1,m%solution%interp%N+1
          do i = 1,m%solution%interp%N+1

            nhat = m%geometry%nHat%boundary(i,j,k,iEl,1,1:3)

            ! Mirror density (tracer)
            m%solution%extBoundary(i,j,k,iEl,1) = &
              m%solution%boundary(i,j,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,j,k,iEl,2)*nhat(1)+ &
                    m%solution%boundary(i,j,k,iEl,3)*nhat(2)+ &
                    m%solution%boundary(i,j,k,iEl,4)*nhat(3)

            m%solution%extBoundary(i,j,k,iEl,2) = &
              m%solution%boundary(i,j,k,iEl,2)-2.0_prec*rhovn*nhat(1)
            m%solution%extBoundary(i,j,k,iEl,3) = &
              m%solution%boundary(i,j,k,iEl,3)-2.0_prec*rhovn*nhat(2)
            m%solution%extBoundary(i,j,k,iEl,4) = &
              m%solution%boundary(i,j,k,iEl,4)-2.0_prec*rhovn*nhat(3)

            ! Mirror potential temperature (tracer)
            m%solution%extBoundary(i,j,k,iEl,5) = &
              m%solution%boundary(i,j,k,iEl,5)

            ! Mirror geopotential — Phi only depends on z, so the
            ! mirror across a no-normal-flow wall is the same value.
            m%solution%extBoundary(i,j,k,iEl,6) = &
              m%solution%boundary(i,j,k,iEl,6)

          enddo
        enddo
      enddo
    endselect

  endsubroutine hbc3d_NoNormalFlow_ESAtmo3D