First be aware that, so far, you have been dealing with projectile velocity as a 2-dimensional phenomenon. This is going to have to change.
In order to conform to standard notation, you'll need to refer to altitude (vertical motion) as the z-coordinate, while horizontal will be handled by x and y coordinates. For ease of calculation, you can assume that you initial velocity has a zero component in one axis, let's say the y axis. Now the equations you're familiar with look the same, except that they are in terms of x and z rather than x and y. Let's start by looking at how to do the simulation in the absence of crosswind.
At any point, you can calculate the x and z velocities $V_x$ and $V_z$, and combine them to get the speed of the projectile, V: $$V = \sqrt{(V_x)^2+(V_y)^2}$$
Use this to calculate the drag on the projectile from the answer you linked to, then decompose the drag force into its x and z components. Apply the components to modify the x and z velocities and calculate the next step.
Now for crosswind. The simplest approach is to assume that the crosswind is much less than the projectile velocities. In this case, you can assume that velocities in the y axis don't affect the other 2 axes.
Start by calculating the relative crosswind - that is the speed of the wind minus the actual y-axis speed. In order to do this, you need to calculate the component of the crosswind which is actually perpendicular to the projectile's motion. For instance, if the wind is at 45 degrees to the projectile's direction, the crosswind sideways to the projectile will be .707 times the wind speed. Recognize that the coefficient of drag is different for crosswind than for the previous calculation, since the "regular" calculation looked at the projectile head-on, while crosswind applies sideways. Knowing the relative crosswind and the crosswind drag coefficient, use the linked answer to calculate the sideways force on the projectile, then apply this to calculate the new velocity and then the new y-axis position.
From the correction for wind direction, it should be clear that, for best accuracy, you ought to factor in the "head-on" component of the wind into the overall drag calculation, but you can always point out that the crosswind is assumed to be much lower than the projectile speed.
The air resistance should be in the opposite direction compared to the acceleration due to the gravity. I really hope you can help me with this!
No. The drag force simply points in the opposite direction of the velocity vector.
Now consider the following simplified model:
Assume the dropping plane was flying horizontally and parallel to the $x$-axis, at speed $v_0$, then at the drop point ($t=0$) the parachute has two velocity vectors with scalars:
$$v_x=v_0$$
$$v_y=0$$
As the chute doesn't deploy immediately, in the $y$ direction two forces act: gravity and air drag, so with Newton we can write:
$$ma=mg-\frac12 \rho C_{y,1}A_{y,1}v_{y}^2$$
Set: $\frac12 \rho C_{y,1}A_{y,1}=\alpha_1$
Then after integration between $t=0, v_y=0$ and $t, v_y$:
$$\large{v_y(t)=\sqrt{\frac{mg}{\alpha_1}\big(1-e^{-\frac{2\alpha_1t}{m}}\big)}}$$
Assume the chute opens at $t=\tau$ then for $t>\tau$ we can derive also:
$$\large{v_y(t)=\sqrt{\frac{1}{\alpha2}\big(mg-\big(mg-\alpha_2v_{y,\tau}^2)e^{-\frac{2\alpha_2t}{m}}\big)}}$$
With:
$$\large{v_{y,\tau}=\sqrt{\frac{mg}{\alpha_1}\big(1-e^{-\frac{2\alpha_1\tau}{m}}\big)}}$$
The parachute also experiences drag in the $x$-direction. Prior to deployment of the chute (and assuming no side wind):
$$ma=-\frac12 \rho C_{x,1}A_{x,1}v_{x}^2$$
Or:
$$a=-\alpha_3v_x^2$$
$\frac12 \rho C_{x,1}A_{x,1}=\alpha_3$
On integrating between $t=0, v_x=v_0$ and $t, v_x$
$$v_x(t)=\frac{v_0}{1+v_0\alpha_3t}$$
And for $t>\tau$:
$$v_x(t)=\frac{v_{x,\tau}}{1+v_{x,\tau}\alpha_4t}$$
where:
$$v_{x,\tau}=\frac{v_0}{1+v_0\alpha_3\tau}$$
Best Answer
Yes, but that equation is not suitable for a case with air resistance. Air resistance will make the acceleration depend on the velocity rather than be constant.
Since the velocity of the ball is different on the way up and the way down, the acceleration curves are different on the way up and the way down. In particular on the way up, the downward acceleration will be greater than $g$, while on the way down the downward acceleration will be less than $g$. This produces an asymmetry that is not shown in your equation.