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