CalculateTendency_ESAtmo2D_t Subroutine

  subroutine CalculateTendency_ESAtmo2D_t(this)
    !! ESAtmo2D tendency = EC inviscid pipeline (parent) + optional
    !! constant-coefficient Laplacian diffusion (BR1 weak-form DG).
    !!
    !! When nu>0 or kappa>0 the parent's CalculateTendency already
    !! computes solutionGradient (because gradient_enabled was set by
    !! SetDiffusion); we then fill diffFlux from solutionGradient,
    !! compute its DG divergence, and accumulate into fluxDivergence
    !! before forming dSdt.
    implicit none
    class(ESAtmo2D_t),intent(inout) :: this
    ! Local
    integer :: i,j,iEl,iVar
    real(prec) :: bMi1,bMi2,bMj1,bMj2,qwi,qwj,jac

    call this%solution%BoundaryInterp()
    call this%solution%SideExchange(this%mesh)

    call this%PreTendency()
    call this%SetBoundaryCondition()

    if(this%gradient_enabled) then
      call this%CalculateSolutionGradient()
      call this%SetGradientBoundaryCondition()
      call this%solutionGradient%AverageSides()
    endif

    call this%SourceMethod()
    call this%BoundaryFlux()

    call this%TwoPointFluxMethod()
    call this%twoPointFlux%UpdateDevice()
    call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior)

    ! EC surface contribution (1/J) M^{-1} B^T f_Riemann
    ! Side ordering: 1=South, 2=East, 3=North, 4=West
    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)

      bMi1 = this%solution%interp%bMatrix(i,2) ! East boundary value at node i
      bMi2 = this%solution%interp%bMatrix(i,1) ! West boundary value at node i
      bMj1 = this%solution%interp%bMatrix(j,2) ! North boundary value at node j
      bMj2 = this%solution%interp%bMatrix(j,1) ! South boundary value at node j
      qwi = this%solution%interp%qWeights(i)
      qwj = this%solution%interp%qWeights(j)
      jac = this%geometry%J%interior(i,j,iEl,1)

      this%fluxDivergence%interior(i,j,iEl,iVar) = &
        this%fluxDivergence%interior(i,j,iEl,iVar)+ &
        (bMi1*this%flux%boundaryNormal(j,2,iEl,iVar)+ &
         bMi2*this%flux%boundaryNormal(j,4,iEl,iVar))/(qwi*jac)+ &
        (bMj1*this%flux%boundaryNormal(i,3,iEl,iVar)+ &
         bMj2*this%flux%boundaryNormal(i,1,iEl,iVar))/(qwj*jac)
    enddo

    ! Add the diffusive contribution if any coefficient is nonzero.
    if(this%nu > 0.0_prec .or. this%kappa > 0.0_prec) then
      call this%DiffusiveFluxMethod()
      call this%DiffusiveBoundaryFlux()
      call this%diffFlux%MappedDGDivergence(this%diffDiv%interior)
      do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                    iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)
        this%fluxDivergence%interior(i,j,iEl,iVar) = &
          this%fluxDivergence%interior(i,j,iEl,iVar)+ &
          this%diffDiv%interior(i,j,iEl,iVar)
      enddo
    endif

    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar)

      this%dSdt%interior(i,j,iEl,iVar) = &
        this%source%interior(i,j,iEl,iVar)- &
        this%fluxDivergence%interior(i,j,iEl,iVar)
    enddo

  endsubroutine CalculateTendency_ESAtmo2D_t