add tdc uq into HLLC

[Debraheem]

Correctly includes tdc Uq into HLLC as a non-HSE source in the du/dt eqn, and fixes other bugs. Needs further polish before merge.

[warrickball]

I feel like there's some simple factorisation of the constants available here. We've got (√(3/2))² from the pow2 expressions and ½ out front, so is this not the same as

Suggested change

	if (k < s% nz) then
	TDC_eturb_cell_start = 0.5d0*(pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + &
	pow2(s% mlt_vc_old(k+1)/sqrt_2_div_3))
	TDC_eturb_cell = 0.5d0*(pow2(s% mlt_vc_ad(k)/sqrt_2_div_3) + &
	pow2(shift_p1(s% mlt_vc_ad(k+1))/sqrt_2_div_3))
	else ! center cell averaged with 0 for inner face
	TDC_eturb_cell_start = 0.5d0*pow2(s% mlt_vc_old(k)/sqrt_2_div_3)
	TDC_eturb_cell = 0.5d0*pow2(s% mlt_vc(k)/sqrt_2_div_3)
	if (k < s% nz) then
	TDC_eturb_cell_start = 0.75d0*(pow2(s% mlt_vc_old(k)) + &
	pow2(s% mlt_vc_old(k+1)))
	TDC_eturb_cell = 0.75d0*(pow2(s% mlt_vc_ad(k)) + &
	pow2(shift_p1(s% mlt_vc_ad(k+1))))
	else ! center cell averaged with 0 for inner face
	TDC_eturb_cell_start = 0.75d0*pow2(s% mlt_vc_old(k))
	TDC_eturb_cell = 0.75d0*pow2(s% mlt_vc(k))

[Debraheem]

Ahh yes, I didn't even consider reducing this. Good catch. (3/2) * (1/2) -> 3/4

[warrickball]

Do we actually use op_err? It looks like both it and ierr and initialised as zero, then ierr is changed to match op_err whenever the latter is non-zero. But that seems the same as just using ierr instead of op_err.

[Debraheem]

This branch is more or less ready for merge. However, I'm not sure if maybe I should add a subroutine for calculating nonlinear work integrals. Might need to copy some utilites from the RSP functionality over?

something like this:

mesa/star/private/rsp.f90

Lines 860 to 940 in c9e7dca

	subroutine calculate_work_integrals(s)
	type (star_info), pointer :: s
	integer :: i, k
	real(dp) :: dt, dm, dVol, P_tw, Pvsc_tw, Ptrb_tw
	character (len=256) :: fname
	dt = s% dt
	! LAST STEP OF PdV
	if(INSIDE==1.and.IWORK==1) then
	IWORK=0
	do I=1,NZN
	k = NZN+1-i
	dm = s% dm(k)
	dVol = VV0(I) - s% Vol_start(k)
	P_tw = THETA*PPP0(I) + &
	THETA1*(s% Pgas_start(k) + s% Prad_start(k))
	Pvsc_tw = THETA*PPQ0(I) + THETA1*s% Pvsc_start(k)
	Ptrb_tw = THETAT*PPT0(I) + THETAT1*s% Ptrb_start(k)
	WORK(I) = WORK(I) + &
	dVol*s% dm(k)*(P_tw + Pvsc_tw + Ptrb_tw) &
	- dt*dm*s% Eq(k)
	WORKQ(I) = WORKQ(I) + dVol*dm*Pvsc_tw
	WORKT(I) = WORKT(I) + dVol*dm*Ptrb_tw
	WORKC(I) = WORKC(I) - dt*dm*s% Eq(k)
	end do
	if (s% rsp_num_periods == s% RSP_work_period) then
	fname = trim(s% log_directory) // '/' // trim(s% RSP_work_filename)
	write(*,*) 'write work integral data to ' // trim(fname)
	open(78,file=trim(fname),status='unknown')
	do I=1,NZN
	k = NZN+1-i
	write(78,'(I3,tr1,F7.4,4(tr1,d16.10))') &
	I, log10(sum(s% dm(1:k))), &
	WORK(I)/EKDEL, WORKQ(I)/EKDEL, &
	WORKT(I)/EKDEL, WORKC(I)/EKDEL
	end do
	close(78)
	end if
	PDVWORK=0.d0
	do I=1,NZN
	k = NZN+1-i
	PDVWORK=PDVWORK+WORK(I)
	WORK(I)=0.d0
	WORKQ(I)=0.d0
	WORKT(I)=0.d0
	WORKC(I)=0.d0
	end do
	s% rsp_GRPDV=PDVWORK/EKDEL
	end if
	! INITIAL STEP OF PdV:
	if((INSIDE==1.and.IWORK==0).or. &
	(s% rsp_num_periods==0.and.IWORK==0))then
	IWORK=1
	do I=1,NZN
	k = NZN+1-i
	VV0(I) = s% Vol_start(k)
	PPP0(I) = s% Pgas_start(k) + s% Prad_start(k)
	PPQ0(I) = s% Pvsc_start(k)
	PPT0(I) = s% Ptrb_start(k)
	PPC0(I) = s% Chi_start(k)
	end do
	end if
	! FIRST AND NEXT STEPS of PdV:
	if(IWORK==1)then
	do I=1,NZN
	k = NZN+1-i
	dm = s% dm(k)
	dVol = s% Vol(k) - s% Vol_start(k)
	P_tw = THETA*(s% Pgas(k) + s% Prad(k)) &
	+ THETA1*(s% Pgas_start(k) + s% Prad_start(k))
	Pvsc_tw = THETA*s% Pvsc(k) + THETA1*s% Pvsc_start(k)
	Ptrb_tw = THETAT*s% Ptrb(k) + THETAT1*s% Ptrb_start(k)
	WORK(I) = WORK(I) + &
	dm*(dVol*(P_tw + Pvsc_tw + Ptrb_tw) - dt*s% Eq(k))
	WORKQ(I)= WORKQ(I) + dm*dVol*Pvsc_tw
	WORKT(I)= WORKT(I) + dm*dVol*Ptrb_tw
	WORKC(I)= WORKC(I) - dt*dm*s% Eq(k)
	end do
	end if
	end subroutine calculate_work_integrals