pbc3d_NoStress_ESAtmo3D Subroutine

  subroutine pbc3d_NoStress_ESAtmo3D(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. Mirroring the gradient
    !! directly (grad_ext = grad_int) gives a non-zero flux through the
    !! wall and is NOT a no-stress condition.
    class(BoundaryCondition),intent(in) :: bc
    class(Model),intent(inout) :: mymodel
    ! Local
    integer :: n,i,j,iEl,k,iVar
    real(prec) :: nhat(1:3),gn,g(1:3)

    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)
            do iVar = 1,m%nvar
              g(1:3) = m%solutionGradient%boundary(i,j,k,iEl,iVar,1:3)
              gn = g(1)*nhat(1)+g(2)*nhat(2)+g(3)*nhat(3)
              m%solutionGradient%extBoundary(i,j,k,iEl,iVar,1) = g(1)-2.0_prec*gn*nhat(1)
              m%solutionGradient%extBoundary(i,j,k,iEl,iVar,2) = g(2)-2.0_prec*gn*nhat(2)
              m%solutionGradient%extBoundary(i,j,k,iEl,iVar,3) = g(3)-2.0_prec*gn*nhat(3)
            enddo
          enddo
        enddo
      enddo
    endselect

  endsubroutine pbc3d_NoStress_ESAtmo3D