I've never heard of anyone using the real gas equation (i.e. Van der waals equation) to derive potential temperature for atmospheric applications. However, if I were to attempt this i would approach it in a similar way as is done for the ideal case. Assuming this form:

(P+a/V2)(V-b) = RT

I would solve for V.

(PV - Pb +a/V + ab/V2) = RT

Now isolating V...

RT - P(V-b) = a/V(1+b/V) => VRT - PV(V-b) = a + b/V => V2RT - PV2(V-b) = aV +ab => -PV3 + (RT-Pb)V2 - aV - ab = 0

or

V3 - (RT/P - b)V2 + a/PV + ab/P = 0

Using the cubic equation (i.e. Ax3+Bx2+Cx+D = 0) you can obtain the roots for V and potentially use that result to solve for V - I will leave that up to you.

As we all know it is now just a matter of using the 1st law of thermodynamics with zero change in entropy:

cp*dT = VdP

Whatever result you get from V can be used to substitute in the equation above, which then is integrated. I've never solved for the roots myself, but this is how I would approach it in the real gas case.

Again, I'm not sure it is necessary to use the real gas EOS since the thermodynamics of the atmosphere itself is modeled well using an ideal gas. It may be different in oceans where salinity is a factor and the fluid is more dense (molecules are in closer proximity). I'm uncertain about rarefied atmospheres. I know there are different considerations there. It is important to note, however, that there is no positive consensus on how to handle real gases that undergo phase transitions between liquid and gas. Apparently the observed behavior does not match theoretical attempts to model this so far. Therefore, I would hesitate to use a real gas here at all given how phase changes between liquid and gas (or gas to liquid) is quite ubiquitous and due to the fact that the atmosphere can be treated as an ideal gas.

Oh, and just to add - since I didn't see that you were interested in the Peng-Robinsion formulation - the approach would be the same, but with a slightly more complicated cubic equation. Without going through all the intermediate details, the cubic form - if I did all the algebra correct - would look something like this:

V3 + (b-RT/P)V2 + (a*alpha - 3b2 -2bRT/P)V + (Pb3 + RTb2 -a*alpha*b) = 0

You might be able to run this through a cubic solver to get the roots. I'm not sure what that will end up looking like but hopefully the integration won't be terribly complicated.

Edward J. Strobach , thanks for the answers! Actually, I had done all that; Matlab's solver gives me the compressibility Z from which I extract the density as function of (p,T), heat capacities etc. My question is rather on the best way to integrate the EOS since dT/T = f(p,v,T) dp or f(Z,p)dp i.e. the choice of isobars/isotherms/isochors.

I see. Then you must get three different roots. I assume that the corresponding valid roots depend on what the temperature is relative to the critical point? What are the conditions you are interested in? Which root(s) have you used in your integrand?