AIRS Code Collection
variance.c
Go to the documentation of this file.
1/*
2 This file is part of the AIRS Code Collection.
3
4 the AIRS Code Collections is free software: you can redistribute it
5 and/or modify it under the terms of the GNU General Public License
6 as published by the Free Software Foundation, either version 3 of
7 the License, or (at your option) any later version.
8
9 The AIRS Code Collection is distributed in the hope that it will be
10 useful, but WITHOUT ANY WARRANTY; without even the implied warranty
11 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with the AIRS Code Collection. If not, see
16 <http://www.gnu.org/licenses/>.
17
18 Copyright (C) 2019-2025 Forschungszentrum Juelich GmbH
19*/
20
26#include "libairs.h"
27
28/* ------------------------------------------------------------
29 Dimensions...
30 ------------------------------------------------------------ */
31
32/* Number of latitudes for threshold tables. */
33#define NLAT_GW 19
34#define NLAT_SURF 6
35#define NLAT_TROP 73
36
37/* Number of months for threshold tables. */
38#define NMON 12
39
40/* Maximum number of longitudes. */
41#define NX 3600
42
43/* Maximum number of latitudes. */
44#define NY 1800
45
46/* ------------------------------------------------------------
47 Global variables...
48 ------------------------------------------------------------ */
49
50/* Latitudes for gravity wave variance thresholds. */
51static double t_gw_lat[NLAT_GW]
52 = { -90, -80, -70, -60, -50, -40, -30, -20, -10, 0,
53 10, 20, 30, 40, 50, 60, 70, 80, 90
54};
55
56/* Gravity wave variance thresholds (ascending orbits). */
57static double t_gw_asc[NMON][NLAT_GW]
58 = { {0.00387, 0.00422, 0.00633, 0.0124, 0.0216, 0.0324,
59 0.0553, 0.0791, 0.0501, 0.0136, 0.0134, 0.0151,
60 0.0522, 0.321, 0.697, 0.776, 0.696, 0.764, 0.771},
61{0.00913, 0.00942, 0.00867, 0.00897, 0.0112, 0.0168,
62 0.0314, 0.0484, 0.032, 0.0128, 0.0122, 0.0134,
63 0.0382, 0.124, 0.345, 0.404, 0.545, 1.16, 1.18},
64{0.0845, 0.0664, 0.0384, 0.0227, 0.0147, 0.0118,
65 0.0141, 0.0184, 0.0162, 0.0123, 0.0124, 0.0124,
66 0.0159, 0.0509, 0.085, 0.103, 0.188, 0.367, 0.529},
67{0.265, 0.297, 0.216, 0.106, 0.0666, 0.0299,
68 0.0169, 0.0129, 0.0116, 0.012, 0.0135, 0.0141,
69 0.0134, 0.0137, 0.017, 0.0268, 0.0259, 0.0319, 0.0323},
70{0.326, 0.44, 0.628, 0.567, 0.434, 0.235,
71 0.0601, 0.0214, 0.0132, 0.0113, 0.0144, 0.0185,
72 0.0179, 0.0142, 0.0116, 0.00945, 0.00865, 0.00918, 0.00878},
73{0.537, 0.73, 1.39, 1.75, 1.35, 0.528,
74 0.188, 0.0311, 0.0133, 0.0124, 0.0205, 0.0313,
75 0.0297, 0.0216, 0.0166, 0.0131, 0.00983, 0.00606, 0.0049},
76{0.382, 1.15, 1.57, 2.13, 1.66, 0.851,
77 0.126, 0.0204, 0.0133, 0.0135, 0.0281, 0.0385,
78 0.0375, 0.0312, 0.0223, 0.0143, 0.00949, 0.0061, 0.00493},
79{0.226, 0.697, 1.68, 1.56, 1.14, 0.496,
80 0.0616, 0.0143, 0.0126, 0.013, 0.0216, 0.0252,
81 0.0241, 0.0206, 0.0152, 0.0106, 0.00976, 0.0105, 0.00998},
82{0.236, 0.489, 0.648, 0.553, 0.524, 0.21,
83 0.033, 0.0129, 0.0116, 0.0129, 0.0163, 0.0165,
84 0.0153, 0.014, 0.0141, 0.0185, 0.0301, 0.0591, 0.0745},
85{0.046, 0.082, 0.112, 0.0806, 0.0516, 0.0469,
86 0.0225, 0.0139, 0.0127, 0.0121, 0.0125, 0.0138,
87 0.0176, 0.0357, 0.0563, 0.062, 0.133, 0.327, 0.3},
88{0.00669, 0.00867, 0.0117, 0.0117, 0.014, 0.015,
89 0.0203, 0.0213, 0.0144, 0.0116, 0.0124, 0.0179,
90 0.0574, 0.185, 0.346, 0.442, 0.54, 0.669, 0.664},
91{0.00355, 0.00381, 0.00658, 0.0125, 0.0217, 0.0304,
92 0.0424, 0.0515, 0.0315, 0.0139, 0.0137, 0.0161,
93 0.0582, 0.306, 0.999, 1.2, 1.14, 0.621, 0.448}
94};
95
96/* Gravity wave variance thresholds (descending orbits). */
97static double t_gw_dsc[NMON][NLAT_GW]
98 = { {0.00383, 0.00458, 0.00866, 0.019, 0.0348, 0.0598,
99 0.144, 0.234, 0.135, 0.0373, 0.0325, 0.0377,
100 0.0858, 0.497, 1.4, 1.32, 0.808, 0.771, 0.773},
101{0.00999, 0.0123, 0.0141, 0.0148, 0.0177, 0.0286,
102 0.0626, 0.102, 0.0717, 0.0302, 0.0261, 0.03,
103 0.086, 0.268, 0.631, 0.716, 1.17, 1.24, 1.21},
104{0.103, 0.096, 0.0715, 0.0535, 0.0343, 0.0245,
105 0.025, 0.0315, 0.0303, 0.0233, 0.023, 0.0257,
106 0.0353, 0.118, 0.197, 0.359, 0.541, 0.585, 0.586},
107{0.272, 0.293, 0.276, 0.226, 0.146, 0.0689,
108 0.0373, 0.0245, 0.0232, 0.0232, 0.0224, 0.0217,
109 0.0242, 0.031, 0.0441, 0.0664, 0.0623, 0.053, 0.0361},
110{0.331, 0.44, 0.641, 0.868, 0.824, 0.47,
111 0.115, 0.0444, 0.0269, 0.0223, 0.0274, 0.0332,
112 0.0273, 0.023, 0.0191, 0.0172, 0.0138, 0.0107, 0.00894},
113{0.554, 0.716, 1.31, 2.29, 2.43, 1.05,
114 0.41, 0.0651, 0.0269, 0.0257, 0.0447, 0.0622,
115 0.0497, 0.0357, 0.0258, 0.0182, 0.0117, 0.00697, 0.00502},
116{0.427, 0.905, 1.44, 2.78, 2.76, 1.52,
117 0.278, 0.041, 0.0279, 0.0296, 0.0629, 0.0818,
118 0.0758, 0.0534, 0.0356, 0.0227, 0.012, 0.00692, 0.00513},
119{0.245, 0.74, 1.88, 2.32, 1.89, 0.883,
120 0.122, 0.0292, 0.0264, 0.0289, 0.0516, 0.059,
121 0.0495, 0.0373, 0.0268, 0.0185, 0.0163, 0.0131, 0.0103},
122{0.272, 0.551, 0.812, 0.844, 0.852, 0.486,
123 0.0842, 0.0269, 0.0225, 0.0239, 0.0322, 0.0324,
124 0.0307, 0.0304, 0.035, 0.0484, 0.0692, 0.0956, 0.0948},
125{0.0644, 0.125, 0.177, 0.135, 0.0922, 0.0899,
126 0.0524, 0.0249, 0.0214, 0.0218, 0.0251, 0.0293,
127 0.0403, 0.0903, 0.168, 0.246, 0.358, 0.378, 0.288},
128{0.00676, 0.00923, 0.0148, 0.0195, 0.0261, 0.0286,
129 0.0302, 0.0343, 0.0298, 0.024, 0.0252, 0.0403,
130 0.131, 0.448, 0.681, 0.923, 0.839, 0.684, 0.629},
131{0.00347, 0.00412, 0.00995, 0.0221, 0.0363, 0.0531,
132 0.104, 0.168, 0.112, 0.0365, 0.0335, 0.0382,
133 0.128, 0.563, 1.62, 1.87, 1.47, 0.652, 0.408}
134};
135
136/* Latitudes for zonal mean tropopause temperatures. */
137static double t_trop_lat[NLAT_TROP]
138 = { 90, 87.5, 85, 82.5, 80, 77.5, 75, 72.5, 70, 67.5, 65, 62.5, 60,
139 57.5, 55, 52.5, 50, 47.5, 45, 42.5, 40, 37.5, 35, 32.5, 30, 27.5,
140 25, 22.5, 20, 17.5, 15, 12.5, 10, 7.5, 5, 2.5, 0, -2.5, -5, -7.5,
141 -10, -12.5, -15, -17.5, -20, -22.5, -25, -27.5, -30, -32.5, -35,
142 -37.5, -40, -42.5, -45, -47.5, -50, -52.5, -55, -57.5, -60, -62.5,
143 -65, -67.5, -70, -72.5, -75, -77.5, -80, -82.5, -85, -87.5, -90
144};
145
146/* Zonal mean tropopause temperatures. */
147static double t_trop[NMON][NLAT_TROP]
148 = { {211.152, 211.237, 211.434, 211.549, 211.614, 211.776, 211.974,
149 212.234, 212.489, 212.808, 213.251, 213.692, 214.193, 214.591,
150 214.985, 215.327, 215.658, 215.956, 216.236, 216.446, 216.738,
151 216.836, 216.032, 213.607, 209.281, 205, 201.518, 198.969,
152 197.123, 195.869, 195.001, 194.409, 193.985, 193.734, 193.617,
153 193.573, 193.6, 193.642, 193.707, 193.856, 194.131, 194.558,
154 195.121, 195.907, 196.91, 198.192, 199.744, 201.583, 203.672,
155 206.012, 208.542, 211.135, 213.681, 216.085, 218.317, 220.329,
156 222.071, 223.508, 224.612, 225.357, 225.761, 225.863, 225.657,
157 225.287, 224.813, 224.571, 224.385, 224.3, 224.257, 224.173,
158 223.786, 222.713, 222.11},
159{212.593, 212.621, 212.801, 212.888, 212.912, 213.054, 213.245,
160 213.512, 213.726, 213.962, 214.259, 214.508, 214.823, 215.037,
161 215.297, 215.545, 215.808, 216.063, 216.323, 216.539, 216.867,
162 217.051, 216.532, 214.512, 210.371, 205.658, 201.758, 198.937,
163 197.047, 195.817, 194.96, 194.386, 193.993, 193.771, 193.673,
164 193.635, 193.658, 193.691, 193.744, 193.872, 194.126, 194.54,
165 195.085, 195.847, 196.8, 198.013, 199.489, 201.261, 203.298,
166 205.596, 208.082, 210.628, 213.156, 215.563, 217.822, 219.903,
167 221.745, 223.311, 224.566, 225.451, 225.947, 226.079, 225.849,
168 225.406, 224.889, 224.643, 224.431, 224.246, 224.079, 223.884,
169 223.42, 222.402, 221.871},
170{215.529, 215.491, 215.539, 215.621, 215.691, 215.808, 215.847,
171 215.881, 215.878, 215.907, 216.02, 216.113, 216.297, 216.342,
172 216.38, 216.369, 216.342, 216.284, 216.185, 215.989, 215.855,
173 215.626, 215.023, 213.432, 209.979, 205.886, 202.212, 199.414,
174 197.488, 196.216, 195.327, 194.732, 194.347, 194.158, 194.095,
175 194.079, 194.116, 194.154, 194.195, 194.302, 194.534, 194.922,
176 195.461, 196.253, 197.288, 198.644, 200.309, 202.293, 204.553,
177 207.033, 209.538, 211.911, 214.016, 215.862, 217.572, 219.179,
178 220.655, 221.959, 223.052, 223.867, 224.344, 224.451, 224.179,
179 223.706, 223.163, 222.876, 222.613, 222.385, 222.154, 221.842,
180 221.304, 220.402, 220.06},
181{219.921, 219.916, 219.99, 219.989, 219.916, 219.867, 219.73,
182 219.522, 219.16, 218.765, 218.448, 218.144, 217.99, 217.756,
183 217.553, 217.311, 217.025, 216.684, 216.241, 215.649, 215.05,
184 214.302, 213.219, 211.496, 208.729, 205.649, 202.594, 200.066,
185 198.144, 196.733, 195.687, 194.991, 194.586, 194.429, 194.418,
186 194.443, 194.492, 194.534, 194.59, 194.718, 194.997, 195.481,
187 196.165, 197.159, 198.462, 200.142, 202.154, 204.533, 207.208,
188 209.848, 212.088, 213.845, 215.222, 216.348, 217.384, 218.383,
189 219.313, 220.131, 220.799, 221.271, 221.479, 221.405, 221.012,
190 220.4, 219.702, 219.227, 218.827, 218.434, 217.977, 217.477,
191 216.783, 215.974, 215.707},
192{225.363, 225.255, 225.064, 224.745, 224.351, 224, 223.551,
193 222.966, 222.195, 221.435, 220.802, 220.245, 219.871, 219.424,
194 218.99, 218.529, 218.013, 217.445, 216.76, 215.859, 214.723,
195 213.049, 211.032, 208.767, 206.449, 204.302, 202.113, 200.187,
196 198.501, 197.153, 196.117, 195.441, 195.121, 195.073, 195.146,
197 195.212, 195.261, 195.288, 195.343, 195.485, 195.772, 196.284,
198 197.018, 198.125, 199.624, 201.604, 204.073, 207.036, 210.193,
199 212.853, 214.611, 215.635, 216.287, 216.801, 217.284, 217.716,
200 218.057, 218.253, 218.282, 218.115, 217.729, 217.15, 216.376,
201 215.449, 214.428, 213.574, 212.847, 212.281, 211.718, 211.211,
202 210.616, 210.112, 210.056},
203{228.431, 228.261, 227.966, 227.457, 226.812, 226.208, 225.518,
204 224.71, 223.701, 222.762, 222.045, 221.486, 221.142, 220.761,
205 220.361, 219.896, 219.34, 218.646, 217.626, 215.983, 213.624,
206 210.817, 208.017, 205.73, 203.8, 202.363, 200.96, 199.778,
207 198.695, 197.845, 197.166, 196.743, 196.6, 196.66, 196.809,
208 196.925, 196.985, 196.996, 197.033, 197.135, 197.335, 197.754,
209 198.367, 199.335, 200.693, 202.564, 205.001, 208.084, 211.473,
210 214.407, 216.208, 217.018, 217.314, 217.394, 217.371, 217.234,
211 216.961, 216.517, 215.878, 215.027, 213.952, 212.697, 211.274,
212 209.736, 208.172, 206.872, 205.84, 205.093, 204.32, 203.816,
213 203.55, 203.49, 203.606},
214{229.01, 228.807, 228.45, 227.839, 227.084, 226.377, 225.589,
215 224.712, 223.665, 222.724, 222.058, 221.658, 221.519, 221.376,
216 221.136, 220.673, 219.926, 218.742, 216.744, 214.028, 210.994,
217 208.374, 206.131, 204.563, 203.251, 202.328, 201.313, 200.411,
218 199.531, 198.876, 198.356, 198.104, 198.088, 198.21, 198.385,
219 198.502, 198.57, 198.601, 198.652, 198.731, 198.869, 199.207,
220 199.737, 200.595, 201.802, 203.491, 205.771, 208.765, 212.241,
221 215.403, 217.439, 218.251, 218.297, 217.988, 217.533, 216.941,
222 216.161, 215.154, 213.887, 212.35, 210.525, 208.481, 206.287,
223 204.068, 202.033, 200.405, 199.106, 198.225, 197.435, 197.02,
224 197.133, 197.527, 197.808},
225{226.525, 226.354, 225.996, 225.433, 224.842, 224.358, 223.818,
226 223.202, 222.426, 221.723, 221.266, 220.98, 220.893, 220.707,
227 220.392, 219.928, 219.182, 218.015, 216.051, 213.399, 210.617,
228 208.318, 206.311, 204.838, 203.515, 202.527, 201.397, 200.423,
229 199.494, 198.848, 198.385, 198.212, 198.294, 198.49, 198.707,
230 198.853, 198.933, 198.967, 199.01, 199.079, 199.207, 199.537,
231 200.081, 200.968, 202.215, 203.946, 206.254, 209.291, 212.876,
232 216.262, 218.487, 219.387, 219.436, 219.048, 218.405, 217.527,
233 216.372, 214.919, 213.152, 211.096, 208.767, 206.247, 203.609,
234 201.029, 198.763, 196.961, 195.578, 194.635, 193.923, 193.54,
235 193.632, 193.944, 193.912},
236{223.293, 223.158, 222.945, 222.571, 222.126, 221.749, 221.362,
237 220.946, 220.404, 219.946, 219.704, 219.599, 219.611, 219.429,
238 219.124, 218.702, 218.063, 217.157, 215.827, 213.879, 211.352,
239 208.833, 206.504, 204.728, 203.168, 201.992, 200.735, 199.74,
240 198.833, 198.213, 197.801, 197.661, 197.765, 197.963, 198.182,
241 198.336, 198.42, 198.456, 198.505, 198.609, 198.794, 199.19,
242 199.796, 200.758, 202.089, 203.915, 206.262, 209.295, 212.807,
243 216.083, 218.329, 219.47, 219.877, 219.846, 219.507, 218.85,
244 217.84, 216.448, 214.652, 212.509, 210.083, 207.534, 204.982,
245 202.596, 200.463, 198.769, 197.441, 196.546, 195.902, 195.472,
246 195.193, 195.066, 195.006},
247{219.564, 219.492, 219.415, 219.191, 218.926, 218.801, 218.691,
248 218.561, 218.298, 218.06, 217.982, 217.956, 218.038, 217.954,
249 217.81, 217.532, 217.08, 216.439, 215.549, 214.31, 212.725,
250 210.573, 208.019, 205.585, 203.459, 201.779, 200.162, 198.879,
251 197.771, 196.987, 196.459, 196.19, 196.172, 196.274, 196.435,
252 196.544, 196.601, 196.644, 196.727, 196.904, 197.184, 197.696,
253 198.42, 199.497, 200.934, 202.825, 205.151, 208.005, 211.279,
254 214.441, 216.87, 218.493, 219.498, 220.072, 220.353, 220.336,
255 219.991, 219.271, 218.142, 216.636, 214.804, 212.776, 210.636,
256 208.535, 206.516, 204.825, 203.383, 202.281, 201.365, 200.561,
257 199.896, 199.415, 199.382},
258{215.926, 215.884, 215.897, 215.814, 215.689, 215.692, 215.707,
259 215.767, 215.815, 215.92, 216.138, 216.327, 216.588, 216.668,
260 216.664, 216.553, 216.373, 216.112, 215.711, 215.025, 214.106,
261 212.596, 210.346, 207.503, 204.604, 202.251, 200.231, 198.607,
262 197.228, 196.174, 195.382, 194.87, 194.61, 194.54, 194.579,
263 194.615, 194.66, 194.709, 194.82, 195.074, 195.487, 196.103,
264 196.904, 198.01, 199.43, 201.246, 203.431, 206.007, 208.905,
265 211.81, 214.34, 216.36, 217.918, 219.141, 220.159, 220.965,
266 221.514, 221.754, 221.637, 221.135, 220.226, 218.986, 217.475,
267 215.879, 214.251, 212.918, 211.84, 211.026, 210.288, 209.553,
268 208.791, 208.132, 208.053},
269{212.893, 212.911, 213.03, 213.109, 213.224, 213.453, 213.653,
270 213.836, 213.98, 214.166, 214.481, 214.787, 215.179, 215.435,
271 215.688, 215.908, 216.084, 216.217, 216.262, 216.123, 215.819,
272 214.977, 213.173, 210.214, 206.619, 203.437, 200.836, 198.843,
273 197.271, 196.078, 195.164, 194.509, 194.057, 193.82, 193.742,
274 193.723, 193.762, 193.813, 193.903, 194.121, 194.49, 195.016,
275 195.698, 196.627, 197.82, 199.359, 201.204, 203.355, 205.78,
276 208.414, 211.057, 213.521, 215.662, 217.504, 219.133, 220.544,
277 221.723, 222.631, 223.274, 223.649, 223.737, 223.547, 223.053,
278 222.357, 221.52, 220.948, 220.527, 220.247, 220.013, 219.726,
279 219.273, 218.506, 218.144}
280};
281
282/* ------------------------------------------------------------
283 Main...
284 ------------------------------------------------------------ */
285
286int main(
287 int argc,
288 char *argv[]) {
289
290 static pert_t *pert;
291
292 static wave_t *wave;
293
294 static char pertname[LEN], set[LEN];
295
296 const double dc_hlat = 25, dc_tlim = 250;
297
298 static double bt[NX][NY], bt_8mu[NX][NY], bt_8mu_min[NX][NY],
299 bt_8mu_max[NX][NY], dt[NX][NY], mtime[NX][NY], glat[NY], glon[NX],
300 fdc[NX][NY], fwg[NX][NY], fgw[NX][NY], fcw[NX][NY], mean[NX][NY],
301 min[NX][NY], max[NX][NY], var[NX][NY], t_dc, t_gw, help[NX * NY];
302
303 const int nmin = 10, iradius = 30;
304
305 static int n[NX][NY], ndc[NX][NY], ngw[NX][NY], ncw[NX][NY], nwg[NX][NY],
306 det_gw, det_cw, det_dc, det_wg, ilat, ncid, varid, minid, maxid, lonid,
307 latid, npid, dimid[10], help2[NX * NY];
308
309 /* Check arguments... */
310 if (argc < 4)
311 ERRMSG("Give parameters: <ctl> <var.tab> <pert1.nc> [<pert2.nc> ...]");
312
313 /* Get control parameters... */
314 scan_ctl(argc, argv, "SET", -1, "full", set);
315 scan_ctl(argc, argv, "PERTNAME", -1, "4mu", pertname);
316 const int nx = (int) scan_ctl(argc, argv, "NX", -1, "360", NULL);
317 const double lon0 = scan_ctl(argc, argv, "LON0", -1, "-180", NULL);
318 const double lon1 = scan_ctl(argc, argv, "LON1", -1, "180", NULL);
319 const int ny = (int) scan_ctl(argc, argv, "NY", -1, "180", NULL);
320 const double lat0 = scan_ctl(argc, argv, "LAT0", -1, "-90", NULL);
321 const double lat1 = scan_ctl(argc, argv, "LAT1", -1, "90", NULL);
322 const int bg_poly_x =
323 (int) scan_ctl(argc, argv, "BG_POLY_X", -1, "0", NULL);
324 const int bg_poly_y =
325 (int) scan_ctl(argc, argv, "BG_POLY_Y", -1, "0", NULL);
326 const int bg_smooth_x =
327 (int) scan_ctl(argc, argv, "BG_SMOOTH_X", -1, "0", NULL);
328 const int bg_smooth_y =
329 (int) scan_ctl(argc, argv, "BG_SMOOTH_Y", -1, "0", NULL);
330 const double gauss_fwhm = scan_ctl(argc, argv, "GAUSS_FWHM", -1, "0", NULL);
331 const double var_dh = scan_ctl(argc, argv, "VAR_DH", -1, "0", NULL);
332 const double thresh_gw =
333 scan_ctl(argc, argv, "THRESH_GW", -1, "-999", NULL);
334 const double thresh_dc =
335 scan_ctl(argc, argv, "THRESH_DC", -1, "-999", NULL);
336 const double dt_trop = scan_ctl(argc, argv, "DT_TROP", -1, "0", NULL);
337 const double dt230 = scan_ctl(argc, argv, "DT230", -1, "0.16", NULL);
338 const double nu = scan_ctl(argc, argv, "NU", -1, "2345.0", NULL);
339 const int output = (int) scan_ctl(argc, argv, "OUTPUT", -1, "1", NULL);
340
341 /* Allocate... */
342 ALLOC(pert, pert_t, 1);
343
344 /* Check grid dimensions... */
345 if (nx < 1 || nx > NX)
346 ERRMSG("Set 1 <= NX <= MAX!");
347 if (ny < 1 || ny > NY)
348 ERRMSG("Set 1 <= NY <= MAX!");
349
350 /* Loop over perturbation files... */
351 for (int iarg = 3; iarg < argc; iarg++) {
352
353 /* Read perturbation data... */
354 FILE *in;
355 if (!(in = fopen(argv[iarg], "r")))
356 continue;
357 else {
358 fclose(in);
359 read_pert(argv[iarg], pertname, pert);
360 }
361
362 /* Recalculate background and perturbations... */
363 if (bg_poly_x > 0 || bg_poly_y > 0 ||
364 bg_smooth_x > 0 || bg_smooth_y > 0 || gauss_fwhm > 0 || var_dh > 0) {
365
366 /* Allocate... */
367 ALLOC(wave, wave_t, 1);
368
369 /* Convert to wave analysis struct... */
370 pert2wave(pert, wave, 0, pert->ntrack - 1, 0, pert->nxtrack - 1);
371
372 /* Estimate background... */
373 background_poly(wave, bg_poly_x, bg_poly_y);
374 background_smooth(wave, bg_smooth_x, bg_smooth_y);
375
376 /* Gaussian filter... */
377 gauss(wave, gauss_fwhm);
378
379 /* Compute variance... */
380 variance(wave, var_dh);
381
382 /* Copy data... */
383 for (int ix = 0; ix < wave->nx; ix++)
384 for (int iy = 0; iy < wave->ny; iy++) {
385 pert->pt[iy][ix] = wave->pt[ix][iy];
386 pert->var[iy][ix] = wave->var[ix][iy];
387 }
388
389 /* Free... */
390 free(wave);
391 }
392
393 /* Detection... */
394 for (int itrack = 0; itrack < pert->ntrack; itrack++)
395 for (int ixtrack = 0; ixtrack < pert->nxtrack; ixtrack++) {
396
397 /* Check data... */
398 if (pert->time[itrack][ixtrack] < 0
399 || pert->lon[itrack][ixtrack] < -180
400 || pert->lon[itrack][ixtrack] > 180
401 || pert->lat[itrack][ixtrack] < -90
402 || pert->lat[itrack][ixtrack] > 90
403 || pert->pt[itrack][ixtrack] < -100
404 || pert->pt[itrack][ixtrack] > 100
405 || !gsl_finite(pert->bt[itrack][ixtrack])
406 || !gsl_finite(pert->pt[itrack][ixtrack])
407 || !gsl_finite(pert->var[itrack][ixtrack])
408 || !gsl_finite(pert->dc[itrack][ixtrack]))
409 continue;
410
411 /* Get and check ascending/descending flag... */
412 int asc =
413 (pert->lat[itrack > 0 ? itrack : itrack + 1][pert->nxtrack / 2] >
414 pert->lat[itrack > 0 ? itrack - 1 : itrack][pert->nxtrack / 2]);
415 if (((set[0] == 'a' || set[0] == 'A') && !asc)
416 || ((set[0] == 'd' || set[0] == 'D') && asc))
417 continue;
418
419 /* Check am/pm flag... */
420 double lt = fmod(pert->time[itrack][ixtrack], 86400.) / 3600.;
421 if (((set[0] == 'm' || set[0] == 'M') && lt > 12.)
422 || ((set[0] == 'n' || set[0] == 'N') && lt < 12.))
423 continue;
424
425 /* Get grid indices... */
426 int ix =
427 (int) ((pert->lon[itrack][ixtrack] - lon0) / (lon1 -
428 lon0) * (double) nx);
429 int iy =
430 (int) ((pert->lat[itrack][ixtrack] - lat0) / (lat1 -
431 lat0) * (double) ny);
432 if (ix < 0 || ix >= nx || iy < 0 || iy >= ny)
433 continue;
434
435 /* Get month index... */
436 int imon =
437 (int) (fmod(pert->time[0][0] / 60. / 60. / 24. / 365.25, 1.) *
438 NMON);
439 if (imon < 0 || imon >= NMON)
440 continue;
441
442 /* Get gravity wave detection threshold... */
443 if (thresh_gw <= 0.0) {
444 ilat = locate_irr(t_gw_lat, NLAT_GW, pert->lat[itrack][ixtrack]);
445 if (asc)
446 t_gw = LIN(t_gw_lat[ilat], t_gw_asc[imon][ilat],
447 t_gw_lat[ilat + 1], t_gw_asc[imon][ilat + 1],
448 pert->lat[itrack][ixtrack]);
449 else
450 t_gw = LIN(t_gw_lat[ilat], t_gw_dsc[imon][ilat],
451 t_gw_lat[ilat + 1], t_gw_dsc[imon][ilat + 1],
452 pert->lat[itrack][ixtrack]);
453 } else
454 t_gw = thresh_gw;
455
456 /* Get deep convection detection threshold... */
457 if (thresh_dc <= 0.0) {
458 ilat =
459 locate_irr(t_trop_lat, NLAT_TROP, pert->lat[itrack][ixtrack]);
460 t_dc =
461 LIN(t_trop_lat[ilat], t_trop[imon][ilat], t_trop_lat[ilat + 1],
462 t_trop[imon][ilat + 1], pert->lat[itrack][ixtrack]) + dt_trop;
463 } else
464 t_dc = thresh_dc + dt_trop;
465
466 /* Detection of gravity waves... */
467 det_gw = (pert->var[itrack][ixtrack] >= t_gw);
468
469 /* Detection of convective waves... */
470 det_cw = 0;
471 if (det_gw)
472 for (int itrack2 = GSL_MAX(itrack - iradius, 0);
473 itrack2 <= GSL_MIN(itrack + iradius, pert->ntrack - 1);
474 itrack2++)
475 for (int ixtrack2 = GSL_MAX(ixtrack - iradius, 0);
476 ixtrack2 <= GSL_MIN(ixtrack + iradius, pert->nxtrack - 1);
477 ixtrack2++) {
478 if (det_cw)
479 break;
480 det_cw = (pert->dc[itrack2][ixtrack2] <= t_dc);
481 }
482
483 /* Detection of deep convection... */
484 det_dc = (pert->dc[itrack][ixtrack] <= t_dc);
485
486 /* Detection of wave generation... */
487 det_wg = 0;
488 if (det_dc)
489 for (int itrack2 = GSL_MAX(itrack - iradius, 0);
490 itrack2 <= GSL_MIN(itrack + iradius, pert->ntrack - 1);
491 itrack2++)
492 for (int ixtrack2 = GSL_MAX(ixtrack - iradius, 0);
493 ixtrack2 <= GSL_MIN(ixtrack + iradius, pert->nxtrack - 1);
494 ixtrack2++) {
495 if (det_wg)
496 break;
497 det_wg = (pert->var[itrack2][ixtrack2] >= t_gw);
498 }
499
500 /* Count events... */
501 n[ix][iy]++;
502 if (det_dc)
503 ndc[ix][iy]++;
504 if (det_wg)
505 nwg[ix][iy]++;
506 if (det_gw)
507 ngw[ix][iy]++;
508 if (det_cw)
509 ncw[ix][iy]++;
510
511 /* Get statistics of perturbations... */
512 mean[ix][iy] += pert->pt[itrack][ixtrack];
513 var[ix][iy] += gsl_pow_2(pert->pt[itrack][ixtrack]);
514 max[ix][iy] = GSL_MAX(max[ix][iy], pert->pt[itrack][ixtrack]);
515 min[ix][iy] = GSL_MIN(min[ix][iy], pert->pt[itrack][ixtrack]);
516
517 /* Get statistics of brightness temperatures... */
518 bt[ix][iy] += pert->bt[itrack][ixtrack];
519 bt_8mu[ix][iy] += pert->dc[itrack][ixtrack];
520 if (n[ix][iy] > 1) {
521 bt_8mu_min[ix][iy]
522 = GSL_MIN(bt_8mu_min[ix][iy], pert->dc[itrack][ixtrack]);
523 bt_8mu_max[ix][iy]
524 = GSL_MAX(bt_8mu_max[ix][iy], pert->dc[itrack][ixtrack]);
525 } else {
526 bt_8mu_min[ix][iy] = pert->dc[itrack][ixtrack];
527 bt_8mu_max[ix][iy] = pert->dc[itrack][ixtrack];
528 }
529
530 /* Get mean time... */
531 mtime[ix][iy] += pert->time[itrack][ixtrack];
532 }
533 }
534
535 /* Analyze results... */
536 for (int ix = 0; ix < nx; ix++)
537 for (int iy = 0; iy < ny; iy++) {
538
539 /* Get geolocation... */
540 mtime[ix][iy] /= (double) n[ix][iy];
541 glon[ix]
542 = lon0 + (ix + 0.5) / (double) nx *(
543 lon1 - lon0);
544 glat[iy]
545 = lat0 + (iy + 0.5) / (double) ny *(
546 lat1 - lat0);
547
548 /* Normalize brightness temperatures... */
549 bt[ix][iy] /= (double) n[ix][iy];
550 bt_8mu[ix][iy] /= (double) n[ix][iy];
551
552 /* Get fractions... */
553 fdc[ix][iy] = (double) ndc[ix][iy] / (double) n[ix][iy] * 100.;
554 fwg[ix][iy] = (double) nwg[ix][iy] / (double) ndc[ix][iy] * 100.;
555 fgw[ix][iy] = (double) ngw[ix][iy] / (double) n[ix][iy] * 100.;
556 fcw[ix][iy] = (double) ncw[ix][iy] / (double) ngw[ix][iy] * 100.;
557
558 /* Check number of observations... */
559 if (n[ix][iy] < nmin) {
560 fdc[ix][iy] = GSL_NAN;
561 fwg[ix][iy] = GSL_NAN;
562 fgw[ix][iy] = GSL_NAN;
563 fcw[ix][iy] = GSL_NAN;
564 bt_8mu[ix][iy] = GSL_NAN;
565 bt_8mu_min[ix][iy] = GSL_NAN;
566 bt_8mu_max[ix][iy] = GSL_NAN;
567 }
568
569 /* Check detections of deep convection at high latitudes... */
570 if (fabs(glat[iy]) > dc_hlat && bt_8mu[ix][iy] <= dc_tlim) {
571 fdc[ix][iy] = GSL_NAN;
572 fwg[ix][iy] = GSL_NAN;
573 fcw[ix][iy] = GSL_NAN;
574 }
575
576 /* Estimate noise... */
577 if (dt230 > 0) {
578 const double nesr = PLANCK(230.0 + dt230, nu) - PLANCK(230.0, nu);
579 dt[ix][iy] = BRIGHT(PLANCK(bt[ix][iy], nu) + nesr, nu) - bt[ix][iy];
580 }
581
582 /* Get mean perturbation and variance... */
583 mean[ix][iy] /= (double) n[ix][iy];
584 var[ix][iy] =
585 var[ix][iy] / (double) n[ix][iy] - gsl_pow_2(mean[ix][iy]);
586 }
587
588 /* Write ASCII file... */
589 if (output == 1) {
590
591 /* Write info... */
592 printf("Write variance statistics: %s\n", argv[2]);
593
594 /* Create file... */
595 FILE *out;
596 if (!(out = fopen(argv[2], "w")))
597 ERRMSG("Cannot create file!");
598
599 /* Write header... */
600 fprintf(out,
601 "# $1 = time [s]\n"
602 "# $2 = longitude [deg]\n"
603 "# $3 = latitude [deg]\n"
604 "# $4 = number of footprints\n"
605 "# $5 = fraction of convection events [%%]\n"
606 "# $6 = fraction of wave generating events [%%]\n"
607 "# $7 = fraction of gravity wave events [%%]\n"
608 "# $8 = fraction of convective wave events [%%]\n"
609 "# $9 = mean perturbation [K]\n"
610 "# $10 = minimum perturbation [K]\n");
611 fprintf(out,
612 "# $11 = maximum perturbation [K]\n"
613 "# $12 = variance [K^2]\n"
614 "# $13 = mean surface temperature [K]\n"
615 "# $14 = minimum surface temperature [K]\n"
616 "# $15 = maximum surface temperature [K]\n"
617 "# $16 = mean background temperature [K]\n"
618 "# $17 = noise estimate [K]\n");
619
620 /* Write results... */
621 for (int iy = 0; iy < ny; iy++) {
622 if (iy == 0 || nx > 1)
623 fprintf(out, "\n");
624 for (int ix = 0; ix < nx; ix++)
625 fprintf(out, "%.2f %g %g %d %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
626 mtime[ix][iy], glon[ix], glat[iy], n[ix][iy],
627 fdc[ix][iy], fwg[ix][iy], fgw[ix][iy], fcw[ix][iy],
628 mean[ix][iy], min[ix][iy], max[ix][iy], var[ix][iy],
629 bt_8mu[ix][iy], bt_8mu_min[ix][iy], bt_8mu_max[ix][iy],
630 bt[ix][iy], dt[ix][iy]);
631 }
632
633 /* Close file... */
634 fclose(out);
635 }
636
637 /* Write netCDF file... */
638 else if (output == 2) {
639
640 /* Create netCDF file... */
641 printf("Write variance statistics: %s\n", argv[2]);
642 NC(nc_create(argv[2], NC_CLOBBER, &ncid));
643
644 /* Set dimensions... */
645 NC(nc_def_dim(ncid, "lat", (size_t) ny, &dimid[0]));
646 NC(nc_def_dim(ncid, "lon", (size_t) nx, &dimid[1]));
647
648 /* Add variables... */
649 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 1, &dimid[0], &latid));
650 add_att(ncid, latid, "deg", "latitude");
651 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 1, &dimid[1], &lonid));
652 add_att(ncid, lonid, "deg", "longitude");
653 NC(nc_def_var(ncid, "var", NC_FLOAT, 2, dimid, &varid));
654 add_att(ncid, varid, "K^2", "brightness temperature variance");
655 NC(nc_def_var(ncid, "min", NC_FLOAT, 2, dimid, &minid));
656 add_att(ncid, minid, "K", "brightness temperature minimum");
657 NC(nc_def_var(ncid, "max", NC_FLOAT, 2, dimid, &maxid));
658 add_att(ncid, maxid, "K", "brightness temperature maximum");
659 NC(nc_def_var(ncid, "np", NC_INT, 2, dimid, &npid));
660 add_att(ncid, npid, "1", "number of footprints");
661
662 /* Leave define mode... */
663 NC(nc_enddef(ncid));
664
665 /* Write data... */
666 NC(nc_put_var_double(ncid, latid, glat));
667 NC(nc_put_var_double(ncid, lonid, glon));
668 for (int ix = 0; ix < nx; ix++)
669 for (int iy = 0; iy < ny; iy++)
670 help[iy * nx + ix] = var[ix][iy] - POW2(dt[ix][iy]);
671 NC(nc_put_var_double(ncid, varid, help));
672 for (int ix = 0; ix < nx; ix++)
673 for (int iy = 0; iy < ny; iy++)
674 help[iy * nx + ix] = min[ix][iy];
675 NC(nc_put_var_double(ncid, minid, help));
676 for (int ix = 0; ix < nx; ix++)
677 for (int iy = 0; iy < ny; iy++)
678 help[iy * nx + ix] = max[ix][iy];
679 NC(nc_put_var_double(ncid, maxid, help));
680 for (int ix = 0; ix < nx; ix++)
681 for (int iy = 0; iy < ny; iy++)
682 help2[iy * nx + ix] = n[ix][iy];
683 NC(nc_put_var_int(ncid, npid, help2));
684
685 /* Close file... */
686 NC(nc_close(ncid));
687 }
688
689 else
690 ERRMSG("Unknown output format!");
691
692 /* Free... */
693 free(pert);
694
695 return EXIT_SUCCESS;
696}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35
locate_irr
int locate_irr(const double *xx, const int n, const double x)
Find array index for irregular grid.
Definition: jurassic.c:4103
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5114
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:126
POW2
#define POW2(x)
Compute x^2.
Definition: jurassic.h:187
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory.
Definition: jurassic.h:121
PLANCK
#define PLANCK(T, nu)
Compute Planck function.
Definition: jurassic.h:183
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:164
background_smooth
void background_smooth(wave_t *wave, int npts_x, int npts_y)
Smooth background.
Definition: libairs.c:176
add_att
void add_att(int ncid, int varid, const char *unit, const char *long_name)
Add variable attributes to netCDF file.
Definition: libairs.c:30
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libairs.c:126
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libairs.c:1584
pert2wave
void pert2wave(pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
Convert radiance perturbation data to wave analysis struct.
Definition: libairs.c:999
read_pert
void read_pert(char *filename, char *pertname, pert_t *pert)
Read radiance perturbation data.
Definition: libairs.c:1137
gauss
void gauss(wave_t *wave, double fwhm)
Apply Gaussian filter to perturbations...
Definition: libairs.c:579
libairs.h
AIRS Code Collection library declarations.
pert_t
Perturbation data.
Definition: libairs.h:205
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libairs.h:232
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libairs.h:229
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libairs.h:214
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libairs.h:223
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libairs.h:208
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libairs.h:211
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libairs.h:217
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libairs.h:226
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libairs.h:220
wave_t
Wave analysis data.
Definition: libairs.h:287
wave_t::nx
int nx
Number of across-track values.
Definition: libairs.h:290
wave_t::ny
int ny
Number of along-track values.
Definition: libairs.h:293
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libairs.h:323
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libairs.h:320
main
int main(int argc, char *argv[])
Definition: variance.c:286
NLAT_GW
#define NLAT_GW
Definition: variance.c:33
NLAT_TROP
#define NLAT_TROP
Definition: variance.c:35
NX
#define NX
Definition: variance.c:41
NY
#define NY
Definition: variance.c:44
NMON
#define NMON
Definition: variance.c:38
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4