Bubble code is not in `m_rhs`? about mfc

[sbryngelson]

@anandrdbz reminder of this issue

MFC/src/simulation/m_rhs.f90

Lines 1160 to 1299 in a80152e

	!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp)
	do l = 0, p
	do k = 0, n
	do j = 0, m
	bub_adv_src(j, k, l) = 0d0
	!$acc loop seq
	do q = 1, nb
	bub_r_src(j, k, l, q) = 0d0
	bub_v_src(j, k, l, q) = 0d0
	bub_p_src(j, k, l, q) = 0d0
	bub_m_src(j, k, l, q) = 0d0
	end do
	end do
	end do
	end do
	ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3
	if (id == ndirs) then
	!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp)
	do l = 0, p
	do k = 0, n
	do j = 0, m
	!$acc loop seq
	do q = 1, nb
	Rtmp(q) = q_prim_qp%vf(rs(q))%sf(j, k, l)
	Vtmp(q) = q_prim_qp%vf(vs(q))%sf(j, k, l)
	end do
	call s_comp_n_from_prim(q_prim_qp%vf(alf_idx)%sf(j, k, l), &
	Rtmp, nbub(j, k, l))
	call s_quad((Rtmp**2.d0)*Vtmp, R2Vav)
	bub_adv_src(j, k, l) = 4.d0*pi*nbub(j, k, l)*R2Vav
	end do
	end do
	end do
	!$acc parallel loop collapse(3) gang vector default(present) private(myalpha_rho, myalpha)
	do l = 0, p
	do k = 0, n
	do j = 0, m
	!$acc loop seq
	do q = 1, nb
	bub_r_src(j, k, l, q) = q_cons_qp%vf(vs(q))%sf(j, k, l)
	!$acc loop seq
	do ii = 1, num_fluids
	myalpha_rho(ii) = q_cons_qp%vf(ii)%sf(j, k, l)
	myalpha(ii) = q_cons_qp%vf(advxb + ii - 1)%sf(j, k, l)
	end do
	myRho = 0d0
	n_tait = 0d0
	B_tait = 0d0
	if (mpp_lim .and. (num_fluids > 2)) then
	!$acc loop seq
	do ii = 1, num_fluids
	myRho = myRho + myalpha_rho(ii)
	n_tait = n_tait + myalpha(ii)*gammas(ii)
	B_tait = B_tait + myalpha(ii)*pi_infs(ii)
	end do
	else if (num_fluids > 2) then
	!$acc loop seq
	do ii = 1, num_fluids - 1
	myRho = myRho + myalpha_rho(ii)
	n_tait = n_tait + myalpha(ii)*gammas(ii)
	B_tait = B_tait + myalpha(ii)*pi_infs(ii)
	end do
	else
	myRho = myalpha_rho(1)
	n_tait = gammas(1)
	B_tait = pi_infs(1)
	end if
	n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
	myRho = q_prim_qp%vf(1)%sf(j, k, l)
	myP = q_prim_qp%vf(E_idx)%sf(j, k, l)
	alf = q_prim_qp%vf(alf_idx)%sf(j, k, l)
	myR = q_prim_qp%vf(rs(q))%sf(j, k, l)
	myV = q_prim_qp%vf(vs(q))%sf(j, k, l)
	if (.not. polytropic) then
	pb = q_prim_qp%vf(ps(q))%sf(j, k, l)
	mv = q_prim_qp%vf(ms(q))%sf(j, k, l)
	call s_bwproperty(pb, q)
	vflux = f_vflux(myR, myV, mv, q)
	pbdot = f_bpres_dot(vflux, myR, myV, pb, mv, q)
	bub_p_src(j, k, l, q) = nbub(j, k, l)*pbdot
	bub_m_src(j, k, l, q) = nbub(j, k, l)*vflux*4.d0*pi*(myR**2.d0)
	else
	pb = 0d0; mv = 0d0; vflux = 0d0; pbdot = 0d0
	end if
	if (bubble_model == 1) then
	! Gilmore bubbles
	Cpinf = myP - pref
	Cpbw = f_cpbw(R0(q), myR, myV, pb)
	myH = f_H(Cpbw, Cpinf, n_tait, B_tait)
	c_gas = f_cgas(Cpinf, n_tait, B_tait, myH)
	Cpinf_dot = f_cpinfdot(myRho, myP, alf, n_tait, B_tait, bub_adv_src(j, k, l), divu%sf(j, k, l))
	myHdot = f_Hdot(Cpbw, Cpinf, Cpinf_dot, n_tait, B_tait, myR, myV, R0(q), pbdot)
	rddot = f_rddot(Cpbw, myR, myV, myH, myHdot, c_gas, n_tait, B_tait)
	else if (bubble_model == 2) then
	! Keller-Miksis bubbles
	Cpinf = myP
	Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
	! c_gas = dsqrt( n_tait*(Cpbw+B_tait) / myRho)
	c_liquid = DSQRT(n_tait*(myP + B_tait)/(myRho*(1.d0 - alf)))
	rddot = f_rddot_KM(pbdot, Cpinf, Cpbw, myRho, myR, myV, R0(q), c_liquid)
	else if (bubble_model == 3) then
	! Rayleigh-Plesset bubbles
	Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
	rddot = f_rddot_RP(myP, myRho, myR, myV, R0(q), Cpbw)
	end if
	bub_v_src(j, k, l, q) = nbub(j, k, l)*rddot
	if (alf < 1.d-11) then
	bub_adv_src(j, k, l) = 0d0
	bub_r_src(j, k, l, q) = 0d0
	bub_v_src(j, k, l, q) = 0d0
	if (.not. polytropic) then
	bub_p_src(j, k, l, q) = 0d0
	bub_m_src(j, k, l, q) = 0d0
	end if
	end if
	end do
	end do
	end do
	end do
	end if