Calculation of the Vertical Velocity in the Asian Summer Monsoon Anticyclone Region Using the Thermodynamic Method With in situ and Satellite Data

Correctly calculating the vertical velocity of the Asian summer monsoon anticyclone (ASMA) region is helpful for accurately knowing the ozone stratosphere–troposphere exchange, so as to explore the variation of ozone in the ASMA region. Therefore, the vertical velocity over the ASMA in June, July, August, and September 2012 and 2016 was calculated using the thermodynamic method, which may avoid the deviations produced by the kinematics method using the mass continuity equation. In order to improve the accuracy, we used high-resolution heating rate datasets obtained via the radiation model in Canadian Atmospheric Global Climate Model called CanAM4.3_RAD based on in situ observations and revised satellite data from MLS/AIRS. The vertical velocity calculated by the thermodynamic method (VT) is then compared with the data from ERA-Interim (VERA–I). In the daytime, values of VT were similar to VERA–I and were dominated by ascending motion, although VT showed descending motion at the western edge of the ASMA below 100 hPa. The intensity of VT was slightly smaller than that of VERA–I at lower levels (200–100 hPa) over the ASMA region and significantly weaker above 100 hPa. The situation was more complex at night. Both VT and VERA–I showed the convergence of vertical wind at 150 hPa and the divergence at 80 hPa, but VT had a smaller standard deviation. VT showed descending in the western and northern ASMA, but VERA–I only descended in the west. The descending motion in the west, seen in both VT and VERA–I, is produced by the heating difference between the Qinghai–Tibet Plateau and the Iranian Plateau. The difference of the two vertical velocities in the northern ASMA may indicate the different understandings of the local Hadley Circulation and local Brewer-Dobson Circulation.

The Qinghai-Tibet Plateau, with an altitude of about 4,000 m, has an important role in global stratospheric-tropospheric exchange (Liu, 1998(Liu, , 2007Fu et al., 2006;Fan et al., 2008;Bian et al., 2011bBian et al., , 2020Tian et al., 2017). More research is required into upward transport over the Qinghai-Tibet Plateau as this is the most representative component of ozone in stratospherictropospheric exchange (Vaughan et al., 1994;Holton et al., 1995;Lv et al., 2009). The Asian summer monsoon anticyclone (ASMA) region, a stable high-pressure system over the Qinghai-Tibet Plateau in the summer, is the most important region of stratosphere-troposphere exchange (Zhou et al., 1995a;Cong et al., 2001;Chen et al., 2010;Deng et al., 2015;Fan et al., 2017a,b;Li et al., 2017Yu et al., 2017;Shi et al., 2018;Yan et al., 2019) and accounts for almost 75% of this exchange process (Rosenlof, 1995;Stohl et al., 2003;Hu et al., 2009;Tian et al., 2009). Because of its existence, the structure and variation of ozone concentrations in the Qinghai-Tibet Plateau region are significantly different from the global average (Zhou and Luo, 1994;Zhou et al., 1995b;Zou and Gao, 1997;Zheng et al., 1998;Gettelman et al., 2004;Park et al., 2007;Vernier et al., 2009;Tang et al., 2019). Using satellite observations, Randel et al. (2010) confirmed that the ASMA region is an important channel for the mass transport from the troposphere to the stratosphere in summer. Ploeger et al. (2010) found that the timescales, paths, and diffusion of the mass flux are mainly determined by vertical transport. Bergman et al. (2012) proved that this vertical transport mainly relies on the vertical motion caused by the deviation of large-scale horizontal circulation over the ASMA region. Therefore, accurate calculations of the vertical velocity over the ASMA region are required to accurately estimate the transport of ozone.
The reliable calculation of the vertical velocity is a longstanding issue in numerical weather prediction (Krishnamurti and Bounoua, 1996). The vertical velocity is typically inferred from the mass continuity equation, which has been called the kinematics method by previous researchers (Eluszkiewicz et al., 2000;Mahowald et al., 2002;Schoeberl et al., 2003;Chipperfield, 2006). This method has a more simplified form (Cong et al., 2001;Gettelman et al., 2004;Chen et al., 2010;Zhan and Li, 2012). Sherman (1953) estimated the vertical velocity using the vertical equation and Lord and Jacqueline (1988), Uma and Rao (2009), Ploeger et al. (2010), and de Roode et al. (2012 assessed the impact of the vertical velocity scheme on the modeling of transport. Davis et al. (2007) used Doppler radar to estimate the mean vertical velocity of air.
The magnitude of the horizontal divergence is partially canceled in the ASMA region as a result of the opposite signs of the zonal and meridional divergence of the horizontal and vertical winds (Guo et al., 2012) and the divergence is insensitive to small perturbations in the large-scale horizontal circulation (Rex et al., 2008;Ploeger et al., 2010). Therefore, the deviations in calculating the divergence by the kinematic method will be carried over. The thermodynamic method has been used to avoid deviations from the kinematic method in the upper troposphere and lower stratosphere (UTLS) (Kistler et al., 2001;Rex et al., 2008;Ploeger et al., 2010). This method uses the diabatic heating rate to reflect the vertical motion of the atmosphere across the isentropic (θ) surface, and the vertical velocity is determined in a θ coordinate system (Rosenlof, 1995;Eluszkiewicz et al., 1996;Randel et al., 2002;Fueglistaler et al., 2005;Rex et al., 2008). Bergman et al. (2013) showed that numerical solutions are affected by the accuracy and resolution of the heating rates in the thermodynamic method. Some researchers have used satellite data and assimilation datasets to calculate heating rates (Rosenlof, 1995;Eluszkiewicz et al., 1996;Randel et al., 2002) and the results were clearly different (Liu et al., 2003;Yang et al., 2004).
Although some studies have applied the thermodynamic method, there is a lack of such research in the ASMA region. We used the revised satellite detection data from the microwave limb sounder (MLS) and the atmospheric infrared sounder (AIRS) for 2012 and 2016 reported by Shi et al. (2017) and the radiation model in Canadian Atmospheric Global Climate Model called CanAM4.3_RAD (Appendix A) to calculate the heating rate. We then applied the thermodynamic equation (Zhou et al., 2014) in isentropic coordinates to introduce a method to calculate the vertical velocity over the ASMA region.
Section "Materials and Methods" of this paper presents the data and methodology. Section "Results" compares the thermodynamic vertical velocity with the reanalysis data in different ways. We discuss the possible reasons for the deviations and the potential factors influencing the results in section "Discussion and Conclusions."

ERA-Interim Reanalysis Data
The ERA-Interim (ERA-I) reanalysis dataset has a higher credibility in the Qinghai-Tibet Plateau region than other reanalysis datasets, especially in the UTLS region, and realizes a more accurate description of the interannual variability (Monge-Sanz et al., 2007;Hu et al., 2009;Guo et al., 2017;Huang et al., 2018;Shi et al., 2018). The vertical velocity in the ERA-I dataset has previously been used as the kinematic vertical velocity in this region. We downloaded the ERA-I data directly from the European Centre for Medium-Range Weather Forecasts' website 1 , including the horizontal wind (u, v), vertical velocity (ω), and temperature (T) fields. The u, v, T, and p were used to determine the vertical velocity by thermodynamic method, and the ω was compared with the vertical velocity we calculated. The vertical resolution in the isobaric surfaces is 13 layers from 300 to 10 hPa. The horizontal grid resolution is 0.75 • × 0.75 • and the temporal resolution is 12 h. We used the dataset for the time period June, July, August, and September (JJAS) in 2012 and 2016.

Heating Rate Data
We substituted the heating rate with only the radiative heating rate data and ignored the sensible heat and latent heat. This is because the water vapor phase transition over the ASMA region is not significant as a result of the small amount of water vapor present and the low temperatures. The exchange of sensible heat is weak due to the small vertical temperature lapse rate and convection in the UTLS region. We used the CanAM4.3_RAD developed by Salzen et al. (2013) to calculate the radiative heating rate. This model can provide a high-accuracy simulation of the cloud radiation effect and the Earth's radiation system . The concentration of ozone, the temperature, and the water vapor content were required to drive the model and these data were obtained from the in situ observational profile of Ali station and the revised satellite dataset from MLS/AIRS. The Ali dataset was obtained using a balloon-borne electrochemical concentration cell sonde at the Indus River (32.5 • N, 80 • E) in the central region of the ASMA. The electrochemical concentration cell measured the profiles with an accuracy within 10% in the troposphere and 5% in the stratosphere (Shi et al., 2017). The sounding balloons were launched at midnight from June to September 2016 to an altitude of 30 km. A total of 17 valid ozone profiles were obtained. In order to reduce the deviation taken by the time, human operation, and weather process, the JJAS average of Ali data was used to represent the summer mean.
The revised satellite dataset from the MLS/AIRS was corrected using the Ali dataset based on the MLS Level 2 Version 4 1 https://apps.ecmwf.int/datasets/ and the AIRS Version 6 Level 2 (L2) Support by different methods in JJAS 2012 and 2016. The original AIRS data had 55 layers from 515 to 10 hPa, and the original MLS data had 33 layers in total from upper troposphere to 0.1 hPa in the vertical direction (Yan et al., 2015;Gamelin and Carvalho, 2017;Herman et al., 2017;Shi et al., 2017). According to the spatial and temporal characteristics of the two data themselves, MLS was observed only in the longitudinal direction, along the horizontal plane of the flight path, while AIRS was observed in both the latitudinal and longitudinal directions, so we processed the MLS/AIRS into different analytical data. In the center of ASMA (about 25 • N-42 • N, 40 • E-100 • E), MLS and AIRS were revised by subtracting the summer mean deviation profile of Ali, and in the other area of ASMA, data were corrected by subtracting the half of the summer mean deviation profile of Ali (Shi et al., 2017). Their horizontal resolutions were 1.9 • × 2.5 • and 1 • × 1 • , respectively. In this way, we obtained the average heating rate at station Ali in JJAS 2016, the mean semi-monthly MLS heating rate, and the daily AIRS heating rate in JJAS 2012 and 2016.

Data Processing
The heating rate and ERA-I reanalysis data were interpolated in the horizontal and vertical directions to obtain the same temporal and spatial resolution before carrying out our calculations. We used typical linear interpolation for u, v, T, ω in the ERA-Interim dataset in the new horizontal grids of 1.9 • × 2.5 • and 1 • × 1 • , which are consistent with the MLS and AIRS heating rate data. We then interpolated the ERA-Interim data and the heating rate data into the same isentropic layers in the vertical direction. The layers were selected based on a combination of the vertical resolution between the two datasets. In this way, the resolution advantage of each dataset was preserved. The area was limited to

Calculation of Vertical Velocity
Previous studies (Eluszkiewicz et al., 2000;Mahowald et al., 2002;Schoeberl et al., 2003;Chipperfield, 2006) have shown that the vertical velocity can be determined using the kinematics equation in isobaric coordinates: which indicates that ω mainly depends on changes in the divergence term ∂u ∂x + ∂v ∂y . However, the term ∂u ∂x + ∂v ∂y is insensitive to small perturbations in u and v and deviations from the observational data are common. We therefore applied the isentropic (θ) coordinate system: Eq. 2 gives a method to determine the vertical wind by directly solving the pressure change rate in the isentropic coordinate system. This is then expanded in θ-coordinates in Eqs 3 and 4. ω θ is the vertical wind in θ-coordinates: The general form of the thermodynamic energy equation gives: whereQ (J s −1 ) is the diabatic heating rate, which means the change of adiabatic heating over time, c P = 1,000 J/K·kg, is the specific heat at constant pressure, T (K) is the temperature, and p (Pa) is the atmosphere pressure. Logarithmic differentiation of the potential temperature equation θ = T P0 P R/Cp gives c P d ln θ = c P d ln T-Rd ln p where R is the gas constant. Combine Eq. 7 with Eq. 6 and let s = c P ln θ, so that the thermodynamic equation can be written as Transforming this with θ gives dθ dt =q T whereq = dQ dt = θ C PQ . Substituting Eq. 9 into Eq. 4, the diagnostic equation for the vertical velocity in isentropic coordinate system is obtained: The first term on the right-hand side is the local pressure tendency term, which represents the vertical velocity caused by the isentropic surface moving up and down locally. The second is the airflow advection term along the isentropic surface, which can be used to judge the direction of vertical velocity when the air flows from high to low pressure, producing an ascending motion and vice versa. The last diabatic term indicates that because ∂P ∂θ < 0, there is ascending motion when there is diabatic heating (Q > 0) and descending motion when there is diabatic cooling (Q < 0) (Zhou et al., 2014).

Statistical Analysis
We introduced some statistical variables for quantitative analysis (Rodgers, 1990;Von Clarmann, 2006;Shi et al., 2017): Eq. 11 is the standard deviation (σ) of a variable where V i is one of the vertical velocity values, the subscript i indicates the different times, andV is the mean value. This variable can illustrate the degree of discretization of each vertical velocity. The mean bias (b) in Eq. 12 was applied to describe the degree of deviation between the vertical velocity of the AIRS and ERA-I data. The standard error (se) is used in Eq. 13 to reduce the inaccuracy caused by abnormal or accidental conditions in Eq. 12. It is stipulated that when the interval (b -se, b + se) does not contain the value 0, then the difference between the two datasets described by b will be accepted and the result is considered to be statistically significant.

RESULTS
Substituting three different heating rates into Eq. 10, we calculated the average vertical velocity profiles for Ali station in JJAS (V Ali ), the daily vertical velocity data of the AIRS (V AIRS ), and the monthly average vertical velocity of the MLS (V MLS ). We then compared these results with the vertical velocity in the ERA-Interim (V ERA−I ) dataset, which had been interpolated at the same spatiotemporal resolution. We used V Ali , V AIRS , and V ERA−I to analyze the results for Ali station. We compared V MLS with V ERA−I for the ASMA region distribution analysis. Figure 1 shows the profiles of the nighttime vertical velocity at Ali station. It can be seen that V Ali and V ERA−I had a similar tendency, with a decreasing rate with height below 150 hPa, between 70 and 50 hPa, and an increasing rate with height in the other layers. In terms of the absolute values, the values in V Ali and V ERA−I approach each other below 80 hPa, but are different above 80 hPa. To make the difference clearer, we reduced the range of abscissa in (-0.1 to 0.4) in Figure 1A. V ERA shows weaker ascending motion in the range of 70-50 hPa than that below 80 hPa, while V Ali is mainly concentrated near the zero line, indicating insignificant vertical movement in the upper levels. They both indicate the vertical motions in the lower stratosphere (70-50 hPa) is weaker. For two vertical velocities, the wind reversed from ascending to descending around 200 hPa and the intensity of descending motion decreased gradually with increasing height. The vertical velocity reached zero near 80 hPa and the wind direction reversed again to a dominant ascending movement. The characteristics of low convergence and high divergence over the Qinghai-Tibet Plateau are described simultaneously by these two vertical velocities.

Results for Ali Station
To avoid the influence of a single station, we compared V AIRS with V ERA−I in different layers (Figure 2). We chose the regional average over (32.5±2 • N, 80.5±2 • E) around Ali station to represent the single-point data. Figure 2 clearly shows that there are significant differences between day and night.
In the daytime, V ERA−I (blue) describes a more pronounced upward movement in the region around Ali station that decreases with height. By contrast, V AIRS only shows a clear ascending motion in the lower layers (200 and 150 hPa) and the values are much lower than V ERA−I . V AIRS approaches zero at 70 hPa.
At night, V ERA−I shows an alternation from ascending to descending motion at 200 hPa and is dominated by descending motion in the range 150-100 hPa. V ERA−I shows a remarkable upward motion in the higher layers around 70 hPa most of the time. V AIRS shows the same variation as V ERA−I , but with smaller values. This is because the radiation heating rate of the atmosphere in the UTLS in the daytime is mainly caused by shortwave solar radiation, which has a strong positive heating effect, leading to upward motion. Most of the heating at night is caused by longwave radiation from the ground, so the heating rate is negative, representing descending motion. The nighttime characteristic represented by V AIRS is consistent with the results of V Ali .
We examined the mean and standard deviation of V AIRS and V ERA−I in JJAS to conduct a more intuitive quantitative analysis and comparison (Figure 3). For the vertical velocity in 2012, V AIRS and V ERA−I are mostly positive in all the layers in the daytime ( Figure 3A) and V AIRS shows a much weaker updraft than V ERA−I . V AIRS reverses from positive to negative between 150 and 100 hPa, leading to wind convergence with a lower standard deviation; but this is not seen for V ERA−I . The consistency of the distribution of V AIRS and V ERA−I increases during the night-time (Figure 3C), especially in the upper stratosphere. The reversal of winds appears near 150 and 80 hPa, consistent with V Ali (Figure 1). The standard deviations also show that the dispersion of both V AIRS and V ERA−I gradually decreases with height in the daytime.
The overall distribution of the mean and standard deviation in 2016 were similar to those in 2012, but the differences between two vertical velocities were more significant than in 2012, especially in the lower stratosphere. The results in 2016 had a greater deviation, especially for V AIRS at night.
We substituted V AIRS and V ERA−I into Eq. 13 and obtained the range (b -se, b + se) to investigate the mean bias and standard error of the vertical velocity (Figure 4). According to the daytime results for 2012 ( Figure 4A) and 2016 (Figure 4B), the mean bias is concentrated near zero below 250 hPa, but is almost positive above 250 hPa. This shows that V ERA−I differs slightly from V AIRS in the lower layers, but that the difference is more pronounced at higher levels in the daytime. Above 250 hPa, the range (b -se, b + se) with the error bars does not contain 0. These results show that, in the daytime, the deviations all have more statistical significance above 250 hPa. In other words, V ERA−I is generally higher than V AIRS above 250 hPa.
The results for nighttime are different. In 2012, the mean bias was positive at all layers, which means that V AIRS was smaller than V ERA−I . But combined with the values of error bars, the results were statistically significant above 180 hPa, but insignificant at lower levels. Therefore V ERA−I is larger than V AIRS in the lower stratosphere and the deviation between them in the upper troposphere has more uncertainty. The situation was the opposite in 2016 and was statistically significant at lower levels (250-110 hPa), but not at higher levels.

Regional Comparison Analysis
Using (32.5 • N, 80.5 • E) as the center point, V MLS was expanded along 32.5 • N and 80.5 • E to obtain the longitudinal and latitudinal vertical profiles. Figure 5 shows the characteristics of V MLS and V ERA−I during the daytime in 2012. The latitudinal profiles in Figures 5A,B show that both V MLS and V ERA−I have clear ascending motion in most areas of the ASMA and convergence at 80 hPa near 91.25 • E (the center of the Qinghai-Tibet Plateau). But V MLS shows descending motion in the west of 91.25 • E and below 100 hPa, which is not shown in V ERA−I . In terms of the deviation in Figure 5C, V MLS has a weaker latitudinal ascending motion than V ERA−I , especially above 100 hPa, where the differences pass the significance test. According to the longitudinal section in Figures 5D,E, upward motion appears simultaneously in most areas of V ERA−I and V MLS . Below 150 hPa, the ascending motion for V MLS (Figure 5D) is as strong as V ERA−I (Figure 5E), but above 150 hPa, V MLS is clearly weaker than V ERA−I . These differences also have passed the significance test ( Figure 5F). In contrast, V MLS is negative north of 34.4 • N and near 100 hPa, although the intensity is not as strong as for the ascending motion.
The nighttime conditions (Figure 6) are different from the daytime conditions (Figure 5). In the latitudinal profiles (Figures 6A,B), V MLS shows strong descending motion below 150 hPa, whereas V ERA−I shows significant ascending. Above 150 hPa, V MLS is negative west of 91.25 • E and positive to the east. Compared with the daytime, the area in which V MLS is descending has expanded to lower and higher layers in the vertical direction. For V ERA−I , descending motion is seen to the west of 91.25 • E with a center of intensity near 100 hPa, but the vertical depth of the descending air is less than for V MLS. The latitudinal profiles for V MLS at nighttime and daytime both show descending motion in the west over the ASMA region. That probably indicates there is a difference in the descending branch of the Hadley Circulation in this region between V MLS and V ERA−I . The Hadley Circulation originally leads to a descending motion, but the stronger heating of the Qinghai-Tibet Plateau in the east and the weaker heating of the Iranian Plateau in the west (Zhang F. et al., 2017) finally introduced a more clearly upward motion over the Qinghai-Tibet Plateau. Thus, the atmosphere in the west descends while the atmosphere in the east ascends. V MLS is more consistent with this explanation both day and night. Meanwhile, V MLS and V ERA−I both clearly show ascending movement in the east of the ASMA region. In the east of 99.25 • E, between 100 and 70 hPa, V MLS ascends significantly slowly than V ERA−I , as illustrated in the test of deviation ( Figure 6C).
The longitudinal sections in Figures 6D,E show V ERA−I ascending, but V MLS descending in the northern part of the ASMA and in the lower layers below 150 hPa. With the local Hadley Circulation and the local Brewer-Dobson Circulation     existing here, a downdraft occurs in this region. We suggest that V MLS can describe the descending motion more clearly. The ascending movement of V MLS is significantly weaker than V ERA−I in the same region. Figure 6 shows that V MLS is mainly subsiding in the northwest and reaches the upper troposphere, whereas V ERA−I is mainly subsiding in the lower layers in the west with a stronger intensity center. The descending motion of V MLS below 150 hPa shows that there is a vertical divergence at 150 hPa in the west and north of the ASMA region. This result is consistent with that of V Ali and V AIRS at night. Figure 7 shows the overall daily mean vertical velocity of JJAS in 2012. The distributions of both V MLS and V ERA−I are more consistent with the characteristics in daytime (Figure 5), indicating that the vertical motion in the ASMA region depends mainly on the vertical velocity in the daytime. The ascending movement of V MLS is weaker than that of V ERA−I , but still reflects the subsiding circulation in the western and northern ASMA.

Conclusion
The thermodynamic method has theoretical advantages over the kinematic method for calculating the vertical velocity over the ASMA region. We established a method based on the thermodynamic equation and introduced heating rate data obtained from in situ observations and revised satellite data from MLS/AIRS to determine the vertical velocity in JJAS 2012 and 2106. By comparing V AIRS , V Ali , and V MLS with V ERA−I , we found that V Ali and V ERA−I present similar behavior in upper troposphere, with a reversal of the wind direction near 150 and 80 hPa during the night in JJAS 2012. A comparison of V AIRS and V ERA−I showed that V AIRS behaved differently from V ERA−I at daytime. V AIRS ascended weakly below 150 hPa and reversed to descending then reaches zero at 100 hPa, while V ERA−I was clearly shown ascending in all layers. At nighttime, V AIRS has similar behaviors of V ERA−I with a weaker vertical motion intensity, which are also consistent with V Ali. The differences between them have statistical significance with the degree of credibility in one standard error. We consider that V AIRS is more in line with V ERA−I in 2012 and better describes the vertical motion in the ASMA region.
The vertical profiles both day and night of V MLS along (32.5 • N, 80.5 • E) illustrated that upward motions are in the east and south, and descending motions are in the west and north over the ASMA region, while V ERA−I only showed descending motions in the west and north at the nighttime. The local flow of the Hadley Circulation and the Brewer-Dobson Circulation in this region may show a downdraft directly in the north. Besides, as a result of the weaker heating of the Iranian Plateau in the west than the Qinghai-Tibet Plateau in the east, the ascending motion over the Iranian Plateau is weaker than that over the Qinghai-Tibet Plateau (Zhang F. et al., 2017). Therefore, the vertical motions over the Iranian Plateau mainly present descending. We suggest that V MLS can describe the vertical motions over the ASMA more accurately-that is, descending in the west and north and ascending in the east and south. The two vertical velocity profiles were similar in the ASMA region during the night-time. Based on the daily mean vertical velocity, we found that the vertical motion over the ASMA region mainly depends on the vertical velocity in the daytime.

Discussion
The thermodynamic vertical velocity we calculated is slightly weaker than the reanalysis data in some areas over the ASMA, which is inconsistent with the expectation that the velocity obtained by the kinematics method is insufficient. The probable reasons for this deviation are as follows.
(1) The radiative heating rate used in this paper resulted in an inaccurate diabatic calculation.
(2) For consistency, the variables u, v, T, andQ were interpolated vertically to fixed isentropic (θ) levels. However, a variable in adjacent grids has different staggered changes in the vertical direction. After interpolation, the original vertical distribution characteristics would be eliminated, thereby introducing a deviation (Rex et al., 2008). (3) The deviations in the ERA-I dataset cannot be neglected and new reanalysis data should be analyzed in the future (Gao et al., 2014;Gianpaolo et al., 2015).

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
DG initiated and coordinated the work. DG and PS provided the calculation and analysis of vertical velocities. DG, PS, and MW wrote the manuscript. CS, YL, and MW gave valuable suggestions for revisions. CS and CZ revised the MLS and AIRS satellite data. WL calculated the heating rates.