The article is very well written and discusses a problem with a well-established methodology by subdiving the perturbation equations into several layers whene the analytical solution is found, reducing the differential problem to an algebraic linear system.

I have  few comments to improve the analysis:

1) Section 2.1 and 2.2 do not have any graphical schematic to help the reader to understand the layers subdivision. I think that adding those (at least for one section) will facilitate the understanding

2) Equation 14. The total derivative operator is undefined. The authors are also focusing on stationary waves, right?

3) Equation 20. Are the derivatives evaluated at Z=H?

4) The method described by equation 20 is very common in acoustics (see the book of Salomons, Computational Atmospheric Acoustics about the FFP method)

5) around line 245: since the authors give importance to the computational time, it is worth to state what solver was used to solve the banded matrix? Was that the native numpy routines or did they use a home-made algorithm?

6) Line 254. I would replace frequency with wavenumber since frequency is more related to temporal variations, while your method is for stationary waves. This applies to the entire manuscript

7) Line 289. The agreement is qualitatively well but not perfect. How can one improve the agremeent? By adding more layers? or there is a limitation in the original data from Wells and Vosper?

There are some typos here and there. I have found two at page 5 at rows 2 and 6 where coefficients and inversions should be singular.

Comments on

“Extending the applicability of a wind-farm gravity-wave model to vertically non-uniform atmospheres” by  Koen Devesse , Luca Lanzilao , Sebastiaan Jamaer, Nicole van Lipzig , and Johan Meyer

This paper adopts the gravity wave formulation of Smith (2010) and the two layer ABL of Allearts and Meyers (2019) while examining variable temperature and wind structure in the troposphere. This extension is fairly work intensive as they allow for almost any vertical distribution of wind and stability. They do this using a matrix solution to the interface conditions between each thin layer up through the troposphere. At the top, they (rather crudely) define a sudden jump in conditions to a uniform stratosphere. Another assumption is that waves are absorbed when they encounter a critical level where wind turning creates a zero intrinsic frequency. A nice feature of the paper is that they check their solution methodology rather carefully against some earlier mountain wave solutions. It seems to me that everything checks out, although I cannot verify every detail. They apply their model to a large set of actual atmospheric profiles taken from the ERA5 re-analysis. Doing this, they find a significant sensitivity daily variation in the profile.

My general suggestions are:

1. Because of the complexity of the tropospheric structure, they do not seem to have isolated the causal relationship between wind response and particular profile features. This random relationship is illustrated in Figure 8 and 9. This is a little disappointing.
2. On the above point, it might be worth checking the following idea. The impact of profile details on the lower troposphere is probably caused either by the way the waves are launched ( the low level N and U) or by the way that waves are reflected downwards. If the latter is true, then their idealized abrupt tropopause might be important as it is probably the main reflector (Fig 2?). An early paper by Klemp and Lilly (1975, referenced here), tried to explain severe downslope wind based on a tuning related to tropopause reflection. If this is the case, one can define another “Froude number” using the critical speeds for deep wave resonance of this type.  Keep in mind however, that tropopause reflection is probably overdone in this model due the assumed sharp tropopause.
3. One problem is that the authors discuss a large number of different model runs but the prrofiles are imprecisely described. I found it difficult to know the properties of each run. The problem begins with the Fig 1 where these plots do not match the equations just below (33 and 34). The problem gets  worse from there as the reference to different wind and stability profiles are too casual and imprecise. This problem must be fixed for reader to follow the logic of the paper. Perhaps a table of run characteristics would help.
4. I am not sure I see the point of figures 4 and 5. I think they are trying to show the impact of the Froude number based on the inversion strength (g’). It seems from these plots that that Froude number makes little difference. I kind of expected this. When tropospheric stability (N) is very small the inversion Froude number makes of big difference  but for realistic values of N, the effect of the inversion is much less. There results in Fig 4 and 5 just seem to verify this general property. The properties aloft are more important than the inversion (g’).
5. In broad terms I think this paper is valuable and significant as it points out that variable tropospheric wind and stability profiles significantly impact wind farm disturbance patterns.
6. Minor Points
   1. A little more explanation of Fig. 3 would be helpful. Why so many peaks?
   2. Line 122 Why does hydrostatic require the inversion?
   3. Line 130 How does it couple to an actuator disk model? Gaussain filter with L=1km
   4. Line 306 gives N=0.013 while the figure 3 caption gives N=0.0113. Are both correct?
   5. Line 227 to 230 are unclear.
   6. The title could be shorter and more attractive

Recommendation: Publish with revisions.