What is Jantzen's formula for the determinant of the Shapovalov form in the case of generalized Verma modules? The best reference I found is
[Kac, Kazhdan '79]
which extends the results of Shapovalov and Jantzen to the case of infinite dimensional Lie algebras. 
Theorem 1 of this paper gives the Shapovalov formula for the determinant of the bilinear Shapovalov form on the modules of a contragredient Lie algebra. 
In the last equation of this paper, the authors give a generalization of the determinant formula for generalized Verma modules based on parabolic subalgebras. This is exactly the result I am looking for. However, the formula given does not seem to be correct because it does not reduce to the standard Shapovalov determinant when we choose the minimal (Borel) parabolic subalgebra. I tried reading the original paper of [Jantzen '77]
but it is written in German. I also looked into the book [Representations of Semisimple Lie Algebras in the BGG Category O, Humphreys, 2008]. Jantzen's determinant formula is mentioned in the beginning of section 9.13 but not given explicitly.
It would be great if you could point me to an english reference that contains the correct version of this determinant formula.
Alternatively, I would also be satisfied with a refinement of Jantzen's simplicity criterion, as given in section 9.13 of the book mentioned above. I would like to know what are the submodules of $M_I(\lambda)$ when $M_{I}(\lambda)$ is not simple. I expect    $M_I(s_\beta \cdot \lambda)$, at least for one $\beta \in \Psi_\lambda^+$, to be a submodule of $M_I(\lambda)$ when $M_{I}(\lambda)$ is not simple. However, I don't know what is the precise statement.
I am interested in the case of the  conformal group in $d>2$   dimensions which is isomorphic to $SO(d+2)$.
The interesting representations for conformal field theory are generalized Verma modules, where the highest weight belongs to an irreducible finite dimensional representation of $SO(d)$ and is an eigenstate of the generator of $SO(2)$ (here I am considering the subgroups $SO(2)\times SO(d)\subset SO(d+2)$) . Then, the other $2d$ generators of $SO(d+2)$ are associated with $d$ positive roots and $d$ negative roots that are used to create the generalized Verma module.
 A: There isn't a lot of follow-up literature (especially in English) concerning the material covered in Jantzen's 1977 Math. Ann. paper, but maybe these remarks will clarify some of the issues you raise.  (Chapter 9 of my 2008 book was a half-way compromise between giving a broad exposition and writing down all the fine details, which would have required even more notation than I already introduced.)
(1) First I'll try to focus your questions in the setting of (reductive) Lie algebras over $\mathbb{C}$.    Consider for any $d \geq 2$ the compact Lie group $SO(d)$, which is somewhat degenerate unless $d \geq 5$.  When $d \geq 5$, the group is (almost) simple, say of Lie rank $\ell$ (the dimension of a maximal torus).  Now pass to the complexification of its Lie algebra, also of rank $\ell$ and simple.   Here there are two slightly different cases:  if $d=2\ell$ is even, the Lie algebra has type $D_\ell$, while if $d = 2\ell+1$ is odd,  the Lie algebra has type $B_\ell$.    (Their root systems and Weyl groups differ.)  
In your set-up, the starting point is actually the group $SO(d+2)$ for some $d$ (say $d \geq 5$ as before) with complexified Lie algebra now called $\mathfrak{g}$.    It has a subgroup: $SO(d) \times SO(2)$, where $SO(2)$ is just a compact torus of dimension 1.  (Maybe the centers of the two factors can overlap?)    If we again write $\ell$ for the rank of $SO(d)$, then $SO(d+2)$ has rank $\ell+1$ in either of the two cases:   If  $d=2\ell$ is even, then $d+2 = 2(\ell+1)$, while if $d=2\ell+1$ is odd, then $d+2 = 2(\ell+1)+1$.   
For the Lie algebras, this translates into having a Levi subalgebra in a maximal standard parabolic subalgebra $\mathfrak{p}_I= \mathfrak{p} = \mathfrak{l} \oplus \mathfrak{u}$.  Here $\mathfrak{g}$ has rank $\ell+1$, while the (simple) derived algebra of $\mathfrak{l}$ has rank $\ell$ and the center of $\mathfrak{l}$ has dimension 1.       The nilradical $\mathfrak{u}$ of $\mathfrak{p}$ then involves the "leftover" positive roots; in each case, a simple computation shows that the number of these (= $\dim \mathfrak{u}$) is just $d$, as stated in the question.   (The number of positive roots in type $B_\ell$ is $\ell^2$ and in type $D_\ell$ is $\ell^2-\ell$.)    For all of this one has fixed a basis of $\ell + 1$ simple roots for the root system and resulting positive roots relative to a Cartan subalgebra of $\mathfrak{g}$ shared with $\mathfrak{p}$.  In each case, the set $I$ of $\ell$ simple roots for the derived algebra of $\mathfrak{l}$ can be obtained by deleting $\alpha_1$ in Bourbaki numbering (the left node of the Dynkin diagram).     
(2) To construct a parabolic Verma module (often called a "generalized" Verma module), start with an integral weight $\lambda$ of $\mathfrak{g}$  whose coordinates relative to the $\ell$ fundamental weights of the derived algebra of $\mathfrak{l}$ are non-negative.    Thus $\lambda$ yields a finite dimensional simple module for $\mathfrak{l}$, which extends to $\mathfrak{p}$ by making $\mathfrak{u}$ act trivially; induce to $\mathfrak{g}$ to get the desired infinite dimensional parabolic Verma module $M_I(\lambda)$.     
(3) A natural problem is to decide when $M_I(\lambda)$ is actually a simple module.     For arbitrary weights the criteria developed by Jantzen (and summarized in $\S9.13$ of my book) are rather intricate, but for a regular weight (one lying inside a Weyl chamber) the precise condition is more straightforward (my Cor. 9.13(d)).   But for the given set-up this still requires a lot of careful bookkeeping with the root systems.     It's useful to start with the much simpler analogous example of $A_1 \hookrightarrow A_2$ discussed in my $\S9.5$, where again one has a maximal parabolic subalgebra.    Here there are only 3 relevant Weyl chambers (out of 6) and 3 simple modules with linked regular weights $\lambda$.     In this example, the corresponding parabolic Verma modules have respectively 2, 2, 1 composition factors.    In particular, note that a weight in the "middle" chamber (below the usual dominant one) defines a parabolic Verma module which is not simple but doesn't include a parabolic Verma module.   
(4) In his 1977 paper Jantzen also adapted to the parabolic case the earlier ideas which he and Shapovalov had independently developed.    Here in particular one can (in principle but seldom in practice) compute explicitly the determinant of the resulting contravariant form on each weight space.    Besides the inductively computed determinant for the Levi part of the parabolic, one has to compute additional factors for the "leftover" positive roots.    Even after unpacking Jantzen's heavy notation for all of this, I'm not confident about being able to carry out the computation in the question here even for a regular weight in each relevant Weyl chamber.   The first "generic" case  seems to arise for $B_3 \hookrightarrow B_4$ (when $d=7$) or $D_4 \hookrightarrow D_5$ (when $d=8$).     (But a couple of smaller "degenerate" cases might also be considered.)    Direct computation in any of these cases looks cumbersome.
A: Let me add some comments and clarifications for the Jantzen-Kac-Kazhdan formula for the determinant of the Shapovalov form of a generalized Verma module. The formula given in [Kac-Kazhdan '79] is:
$$det(F_I^0)_\eta (\lambda) = \prod_{\alpha\in\Delta^+\backslash \Delta_I} \prod_{n>0} \left( (\lambda + \rho) (h_\alpha) - n \frac{(\alpha,\alpha)}{2} \right)^{\chi_I (\lambda-n\alpha)(\eta)}.$$
This gives the determinant for the component of weight $\lambda - \eta$ in the generalized Verma module $M_I (\lambda)$ of highest weight $\lambda$. An important point to keep in mind is that the exponents in the right hand side may be negative. While the determinant itself is a polynomial, the product in the right hand side is a rational expression, where all denominators cancel.
Kac-Kazhdan introduce $\chi_I (\lambda)$ as the character of the generalized Verma module $M_I (\lambda)$ and then give the formula for it. In fact, $\chi_I$ should be defined by the formula. For some reason, the formula for this character in [Kac-Kazhdan '79] is for a lowest weight module, instead of the highest weight. The proper version should look like this: 
$$\chi_I (\mu)= \left(\prod_{\alpha\in\Delta^+} (1 - e^{-\alpha}) \right)^{-1} \sum_{w \in W_I} \epsilon(w) e^{(\mu-\rho)-w(\mu-\rho)},$$ 
with $\chi_I (\mu) (\eta)$ denoting the coefficient at $e^{\mu - \eta}$ (and not at $e^\eta$ as stated in Kac-Kazhdan).
It is emphasized in [Jantzen '77] that $\chi_I (\mu)$ may be negative or zero. It satisfies $\chi_I (w(\mu+\rho) - \rho) = \epsilon(w) \chi_I (\mu)$ for $w \in W_I$, and in particular vanishes if the dot action of some simple reflection in $W_I$ fixes $\mu$.
Because of this behaviour of $\chi_I (\mu)$, when computing $det(F_I^0)_\eta (\lambda)$ one has to consider not only pairs $\alpha$, $n$ with $n \alpha \leq \eta$, but rather those with $\eta - n \alpha$ in ${\mathbb Z}_+ \Delta^+ + {\mathbb Z} \Delta_I$.
