DiffusiveBoundaryFlux_ESAtmo2D_t Subroutine

  subroutine DiffusiveBoundaryFlux_ESAtmo2D_t(this)
    !! Fill diffFlux%boundaryNormal with the SIPG-stabilised BR1 flux:
    !!
    !!   f_R^diff(iVar) = -coeff(iVar) * (avg_grad . n) * nmag
    !!                    + tau(iVar) * (uL - uR) * nmag
    !!
    !! tau(iVar) = eta_penalty * coeff(iVar) * (N+1)^2 / length_scale.
    !! avg_grad is solutionGradient%avgBoundary (populated by
    !! AverageSides()); uL, uR are solution%boundary, %extBoundary.
    implicit none
    class(ESAtmo2D_t),intent(inout) :: this
    ! Local
    integer :: i,k,iel
    real(prec) :: nhat(1:2),nmag,gradn,coeff,tau,np2
    real(prec) :: uL,uR

    np2 = real((this%solution%interp%N+1)**2,prec)

    do concurrent(i=1:this%solution%N+1,k=1:4,iel=1:this%mesh%nElem)
      nhat = this%geometry%nHat%boundary(i,k,iel,1,1:2)
      nmag = this%geometry%nScale%boundary(i,k,iel,1)

      ! rho — no mass diffusion, no penalty
      this%diffFlux%boundaryNormal(i,k,iel,1) = 0.0_prec

      ! Momentum equations use nu
      coeff = this%nu
      tau = this%eta_penalty*coeff*np2/this%length_scale

      gradn = this%solutionGradient%avgBoundary(i,k,iel,2,1)*nhat(1)+ &
              this%solutionGradient%avgBoundary(i,k,iel,2,2)*nhat(2)
      uL = this%solution%boundary(i,k,iel,2)
      uR = this%solution%extBoundary(i,k,iel,2)
      this%diffFlux%boundaryNormal(i,k,iel,2) = (-coeff*gradn+tau*(uL-uR))*nmag

      gradn = this%solutionGradient%avgBoundary(i,k,iel,3,1)*nhat(1)+ &
              this%solutionGradient%avgBoundary(i,k,iel,3,2)*nhat(2)
      uL = this%solution%boundary(i,k,iel,3)
      uR = this%solution%extBoundary(i,k,iel,3)
      this%diffFlux%boundaryNormal(i,k,iel,3) = (-coeff*gradn+tau*(uL-uR))*nmag

      ! rho*theta uses kappa
      coeff = this%kappa
      tau = this%eta_penalty*coeff*np2/this%length_scale

      gradn = this%solutionGradient%avgBoundary(i,k,iel,4,1)*nhat(1)+ &
              this%solutionGradient%avgBoundary(i,k,iel,4,2)*nhat(2)
      uL = this%solution%boundary(i,k,iel,4)
      uR = this%solution%extBoundary(i,k,iel,4)
      this%diffFlux%boundaryNormal(i,k,iel,4) = (-coeff*gradn+tau*(uL-uR))*nmag

      ! Geopotential — no diffusion, no penalty.
      this%diffFlux%boundaryNormal(i,k,iel,5) = 0.0_prec
    enddo

  endsubroutine DiffusiveBoundaryFlux_ESAtmo2D_t