Add work integrals from RSP into star

[Debraheem]

We should copy some utilites from RSP functionality over to star for calculating work integrals.

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

Originally posted by @Debraheem in #883 (comment)