TwoPointFluxMethod_ESAtmo2D_t Subroutine

  subroutine TwoPointFluxMethod_ESAtmo2D_t(this)
    !! Pre-projected scalar contravariant two-point Souza et al. (2023) EC
    !! flux. For each node pair (a, b) along reference direction r:
    !!
    !!   Fc^r_v_i = sum_d 0.5*(Ja^r_d(a) + Ja^r_d(b)) * f_d(v_i, sL=s(a), sR=s(b))
    !!
    !! Gravity is NOT split into the pressure flux here; it is handled by
    !! the Souza non-conservative gravity flux differencing in SourceMethod
    !! (which uses the geopotential carried as state variable index 5).
    implicit none
    class(ESAtmo2D_t),intent(inout) :: this
    ! Local
    integer :: nn,i,j,d,iEl,iVar
    real(prec) :: sL(1:this%nvar),sR(1:this%nvar)
    real(prec) :: Fphys(1:this%nvar,1:2)
    real(prec) :: Fc

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

      sL = this%solution%interior(i,j,iEl,1:this%nvar)

      ! ------------- xi^1: pair (i,j)-(nn,j) -------------
      sR = this%solution%interior(nn,j,iEl,1:this%nvar)
      Fphys = this%twopointflux2d(sL,sR)
      do iVar = 1,this%nvar
        Fc = 0.0_prec
        do d = 1,2
          Fc = Fc+0.5_prec*( &
               this%geometry%dsdx%interior(i,j,iEl,1,d,1)+ &
               this%geometry%dsdx%interior(nn,j,iEl,1,d,1))* &
               Fphys(iVar,d)
        enddo
        this%twoPointFlux%interior(nn,i,j,iEl,iVar,1) = Fc
      enddo

      ! ------------- xi^2: pair (i,j)-(i,nn) -------------
      sR = this%solution%interior(i,nn,iEl,1:this%nvar)
      Fphys = this%twopointflux2d(sL,sR)
      do iVar = 1,this%nvar
        Fc = 0.0_prec
        do d = 1,2
          Fc = Fc+0.5_prec*( &
               this%geometry%dsdx%interior(i,j,iEl,1,d,2)+ &
               this%geometry%dsdx%interior(i,nn,iEl,1,d,2))* &
               Fphys(iVar,d)
        enddo
        this%twoPointFlux%interior(nn,i,j,iEl,iVar,2) = Fc
      enddo

    enddo

  endsubroutine TwoPointFluxMethod_ESAtmo2D_t