SourceMethod_ESAtmo3D_t Subroutine

  subroutine SourceMethod_ESAtmo3D_t(this)
    !! Souza et al. (2023) non-conservative gravity flux differencing.
    !!
    !! The rho*w equation carries the body force -rho * partial_z Phi where
    !! Phi = g*z is the geopotential, which we keep as state variable index
    !! 6 (no flux, no source of its own — its tendency is identically zero).
    !! On a curvilinear mesh, the SBP-EC strong-form flux differencing for
    !! the non-conservative product rho * partial_z Phi reads:
    !!
    !!   [rho * d_z Phi]_i = (1/J_i) sum_r sum_j D_split^r[i_r, j]
    !!                       * 0.5*(Ja^r_z(i) + Ja^r_z(j))
    !!                       * <rho>_log(s_i, s_j) * (Phi_j - Phi_i)
    !!
    !! and we set source(rho*w) = - that result. The other variables
    !! (rho, rho*u, rho*v, rho*theta, Phi) all have zero source.
    !!
    !! Why this form replaces the WB hydrostatic split: gravity now lives
    !! on the same node-pair carrier as the EC volume flux (every node
    !! pair, including those that span an element interface), so the
    !! gravity / pressure-flux topology mismatch that produced
    !! z-element-aligned banding is gone. Well-balancing is recovered
    !! automatically: in the hydrostatic state the EC volume flux gives
    !! d_z p_hyd and the source gives -rho_hyd*g, which cancel pointwise.
    implicit none
    class(ESAtmo3D_t),intent(inout) :: this
    ! Local
    integer :: i,j,k,iEl,nn
    real(prec) :: rho_ijk,phi_ijk,rho_p,phi_p
    real(prec) :: rho_log,Ja_avg,acc,jac

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

      rho_ijk = this%solution%interior(i,j,k,iEl,1)
      phi_ijk = this%solution%interior(i,j,k,iEl,6)
      acc = 0.0_prec

      do nn = 1,this%solution%N+1
        ! xi^1 partner (nn, j, k)
        rho_p = this%solution%interior(nn,j,k,iEl,1)
        phi_p = this%solution%interior(nn,j,k,iEl,6)
        rho_log = log_mean(rho_ijk,rho_p)
        Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,k,iEl,1,3,1)+ &
                           this%geometry%dsdx%interior(nn,j,k,iEl,1,3,1))
        acc = acc+this%solution%interp%dSplitMatrix(nn,i)* &
              Ja_avg*rho_log*(phi_p-phi_ijk)

        ! xi^2 partner (i, nn, k)
        rho_p = this%solution%interior(i,nn,k,iEl,1)
        phi_p = this%solution%interior(i,nn,k,iEl,6)
        rho_log = log_mean(rho_ijk,rho_p)
        Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,k,iEl,1,3,2)+ &
                           this%geometry%dsdx%interior(i,nn,k,iEl,1,3,2))
        acc = acc+this%solution%interp%dSplitMatrix(nn,j)* &
              Ja_avg*rho_log*(phi_p-phi_ijk)

        ! xi^3 partner (i, j, nn)
        rho_p = this%solution%interior(i,j,nn,iEl,1)
        phi_p = this%solution%interior(i,j,nn,iEl,6)
        rho_log = log_mean(rho_ijk,rho_p)
        Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,k,iEl,1,3,3)+ &
                           this%geometry%dsdx%interior(i,j,nn,iEl,1,3,3))
        acc = acc+this%solution%interp%dSplitMatrix(nn,k)* &
              Ja_avg*rho_log*(phi_p-phi_ijk)
      enddo

      jac = this%geometry%J%interior(i,j,k,iEl,1)
      this%source%interior(i,j,k,iEl,1) = 0.0_prec
      this%source%interior(i,j,k,iEl,2) = 0.0_prec
      this%source%interior(i,j,k,iEl,3) = 0.0_prec
      this%source%interior(i,j,k,iEl,4) = -acc/jac
      this%source%interior(i,j,k,iEl,5) = 0.0_prec
      this%source%interior(i,j,k,iEl,6) = 0.0_prec

    enddo

  endsubroutine SourceMethod_ESAtmo3D_t