001 /**
002 *
003 */
004 package org.wdssii.nse;
005
006 import org.apache.commons.logging.Log;
007 import org.apache.commons.logging.LogFactory;
008 import org.wdssii.core.DataType;
009
010 /**
011 *
012 * Creates and maintains a NSE profile
013 *
014 * @author lakshman
015 *
016 */
017 public class NSEProfile {
018 private static Log log = LogFactory.getLog(NSEProfile.class);
019 float[] pressureLevels;
020 float[][] sfcElevation;
021 float[][] mslPressure;
022 float[][] sfcUWind;
023 float[][] sfcVWind;
024 float[][] sfcPressure;
025 float[][] sfcTemperature;
026 float[][] sfcDewPoint;
027
028 float[][][] temperatureC, dewpoint, height, u, v, w;
029
030 /**
031 * Invoke this constructor, set the fields sfcElevation, mslPressure, ...
032 * sfcDewPoint and then call finishInit()
033 *
034 * @param pressureLevels
035 */
036 NSEProfile(float[] pressureLevels) {
037 this.pressureLevels = pressureLevels;
038 temperatureC = new float[pressureLevels.length][][];
039 dewpoint = new float[pressureLevels.length][][];
040 height = new float[pressureLevels.length][][];
041 u = new float[pressureLevels.length][][];
042 v = new float[pressureLevels.length][][];
043 w = new float[pressureLevels.length][][];
044 }
045
046 private int numLat, numLon;
047 private float[][][] satVaporPressure, vaporPressure, rh, lclTemp, lclPres,
048 theta, mixingRatio, virtualTempK, thetaE, wetBulbTempK,
049 moistParcelTemp;
050 private int[][] nearSurfaceGridPoint; // index of first pressureLevel above ground
051
052 /** Call after setting all fields. */
053 void finishInit() {
054 long start = System.nanoTime();
055 log.info("Starting to compute 3D NSE points ...");
056
057 numLat = temperatureC[0].length;
058 numLon = temperatureC[0][0].length;
059
060 satVaporPressure = new float[pressureLevels.length][numLat][numLon];
061 vaporPressure = new float[pressureLevels.length][numLat][numLon];
062 rh = new float[pressureLevels.length][numLat][numLon];
063 lclTemp = new float[pressureLevels.length][numLat][numLon];
064 lclPres = new float[pressureLevels.length][numLat][numLon];
065 theta = new float[pressureLevels.length][numLat][numLon];
066 mixingRatio = new float[pressureLevels.length][numLat][numLon];
067 virtualTempK = new float[pressureLevels.length][numLat][numLon];
068 thetaE = new float[pressureLevels.length][numLat][numLon];
069 wetBulbTempK = new float[pressureLevels.length][numLat][numLon];
070 moistParcelTemp = new float[pressureLevels.length][numLat][numLon];
071 nearSurfaceGridPoint = new int[numLat][numLon];
072
073 for (int i = 0; i < numLat; ++i) {
074 for (int j = 0; j < numLon; ++j) {
075 float ground = sfcElevation[i][j];
076 nearSurfaceGridPoint[i][j] = 0;
077 for (int k = 0; k < pressureLevels.length; ++k) {
078 if (height[k][i][j] < ground) {
079 temperatureC[k][i][j] = dewpoint[k][i][j] = u[k][i][j] = v[k][i][j] = DataType.MissingData;
080 satVaporPressure[k][i][j] = vaporPressure[k][i][j] = rh[k][i][j] = lclTemp[k][i][j] = theta[k][i][j] = mixingRatio[k][i][j] = virtualTempK[k][i][j] = thetaE[k][i][j] = wetBulbTempK[k][i][j] = moistParcelTemp[k][i][j] = DataType.MissingData;
081 } else {
082 if (nearSurfaceGridPoint[i][j] == 0 ){
083 nearSurfaceGridPoint[i][j] = k; // index of first point above ground
084 }
085 initNsePoint(k, i, j);
086 }
087 }
088 }
089 }
090 long end = System.nanoTime();
091 log.info("Took " + (end - start) * 0.001 * 0.001
092 + " milliseconds to compute 3D NSE points");
093 initSfcNsePoints();
094 }
095
096 private void initNsePoint(int k, int i, int j) {
097 // order is important as later variables depends on former ones
098 satVaporPressure[k][i][j] = calcSaturationVaporPressure(temperatureC[k][i][j]);
099 vaporPressure[k][i][j] = calcVaporPressure(dewpoint[k][i][j]);
100 rh[k][i][j] = calcRHFromVaporPressures(vaporPressure[k][i][j], satVaporPressure[k][i][j]);
101 lclTemp[k][i][j] = calcLclTemp(temperatureC[k][i][j], dewpoint[k][i][j]);
102 lclPres[k][i][j] = calcLclPres(lclTemp[k][i][j], pressureLevels[k],
103 temperatureC[k][i][j] + 273.15f);
104 theta[k][i][j] = calcPotentialTemperature(pressureLevels[k],
105 temperatureC[k][i][j] + 273.15f);
106 mixingRatio[k][i][j] = calcMixingRatio(vaporPressure[k][i][j],
107 pressureLevels[k]);
108 virtualTempK[k][i][j] = calcVirtualTemperature(temperatureC[k][i][j],
109 mixingRatio[k][i][j]);
110 thetaE[k][i][j] = calcEquivalentPotentialTemperature(
111 temperatureC[k][i][j], mixingRatio[k][i][j], lclTemp[k][i][j],
112 pressureLevels[k]);
113 wetBulbTempK[k][i][j] = calcWetBulb(vaporPressure[k][i][j],
114 satVaporPressure[k][i][j], dewpoint[k][i][j],
115 temperatureC[k][i][j], pressureLevels[k]);
116 moistParcelTemp[k][i][j] = calcMoistParcelTemp(thetaE[k][i][j],
117 pressureLevels[k]);
118 }
119
120 float[][] sfcRelativeHumidity, sfcVirtualTemperature, sfcTheta, sfcThetaE, sfcMixingRatio, sfcLclTemperature, sfcLclPressure, sfcWetBulbTemperatureK;
121 private void initSfcNsePoints() {
122 long start = System.nanoTime();
123 log.info("Starting to compute Surface NSE points ...");
124 sfcRelativeHumidity = new float[numLat][numLon];
125 sfcMixingRatio = new float[numLat][numLon];
126 sfcVirtualTemperature = new float[numLat][numLon];
127 sfcTheta = new float[numLat][numLon];
128 sfcThetaE = new float[numLat][numLon];
129 sfcLclTemperature = new float[numLat][numLon];
130 sfcLclPressure = new float[numLat][numLon];
131 sfcWetBulbTemperatureK = new float[numLat][numLon];
132
133 for (int i = 0; i < numLat; ++i) {
134 for (int j = 0; j < numLon; ++j) {
135 float vaporPressure = calcVaporPressure(sfcDewPoint[i][j]);
136 float saturationVaporPressure = calcSaturationVaporPressure(sfcTemperature[i][j]);
137 sfcRelativeHumidity[i][j] = calcRHFromVaporPressures(vaporPressure, saturationVaporPressure);
138 sfcMixingRatio[i][j] = calcMixingRatio(vaporPressure, sfcPressure[i][j]);
139 sfcVirtualTemperature[i][j] = calcVirtualTemperature(sfcTemperature[i][j], sfcMixingRatio[i][j]);
140 sfcTheta[i][j] = calcPotentialTemperature(sfcPressure[i][j], sfcTemperature[i][j] + 273.15f);
141 sfcThetaE[i][j] = calcEquivalentPotentialTemperature(sfcTemperature[i][j], sfcRelativeHumidity[i][j], sfcPressure[i][j]);
142 sfcLclTemperature[i][j] = calcLclTemp(sfcTemperature[i][j], sfcDewPoint[i][j]);
143 sfcLclPressure[i][j] = calcLclPres(sfcLclTemperature[i][j], sfcPressure[i][j], sfcTemperature[i][j] + 273.15f);
144 }
145 }
146 long end = System.nanoTime();
147 log.info("Took " + (end - start) * 0.001 * 0.001
148 + " milliseconds to compute Surface NSE points");
149 }
150
151 private static float calcDewPoint(float temp, float rh) {
152 if ((temp != DataType.MissingData) && (rh != DataType.MissingData)) {
153 double RH = rh * 0.01;
154 // from Bolton (MWR 1980):
155 double es = 6.112 * Math.exp((17.67 * temp) / (temp + 243.5));
156 double e = es * RH;
157 double dewp = (243.5 * Math.log(e) - 440.8) / (19.48 - Math.log(e));
158 return (float) dewp;
159 }
160 return DataType.MissingData;
161 }
162
163 private static float calcMoistParcelTemp(float thetae, float pres) {
164
165 // calculate the parcel temperature at the given pressure for a
166 // given parcel thetaE, using Newton-Rhapson interation.
167 //
168 if (pres == DataType.MissingData || thetae == DataType.MissingData) {
169 return DataType.MissingData;
170 }
171 double maxt = 0;
172 if ((thetae - 270) > 0)
173 maxt = thetae - 270;
174
175 // generate a first-guess the GEMPAK way (and in degreeC)
176 float tguess = (float) ((thetae - 0.5 * Math.pow(maxt, 1.05))
177 * Math.pow(pres / 1000, 0.2) - 273.15);
178 // set convergence
179 double conv = 0.01;
180 // 297.621 / 900
181 //
182 double oldguess1 = -999999999999999.0;
183 double oldguess2 = -999999999999999.0;
184 // iterate by degrees
185 for (int i = 1; i <= 100; ++i) {
186 float rh = i * 0.01f;
187 float thetae1 = calcEquivalentPotentialTemperature(tguess, rh, pres);
188 float thetae2 = calcEquivalentPotentialTemperature(tguess + 1, rh,
189 pres);
190 if (thetae1 == DataType.MissingData
191 || thetae2 == DataType.MissingData) {
192 return DataType.MissingData;
193 }
194
195 float correction = (thetae - thetae1) / (thetae2 - thetae1);
196 tguess += correction;
197 // ADDED: are we stuck in a loop on either side of the solution?
198 if ((oldguess2 - tguess) < 0.1 && Math.abs(correction) < 0.3) {
199 // not good enough, but need it to help find other bugs
200 tguess += 273.15f;
201 return (float) tguess;
202 }
203 if (Math.abs(correction) < conv) {
204 // return on solution
205 tguess += 273.15f;
206 return (float) tguess;
207 }
208
209 oldguess1 = tguess;
210 oldguess2 = oldguess1;
211 }
212 // failed to converge
213 return DataType.MissingData;
214
215 }
216
217 /*private static float calcRHFromDewpointAndTemperature(float dewpointC, float tempC){
218 float saturationVaporPressure = calcSaturationVaporPressure(tempC);
219 float vaporPressure = calcVaporPressure(dewpointC);
220 return calcRHFromVaporPressures(vaporPressure, saturationVaporPressure);
221 }*/
222
223 private static float calcRHFromVaporPressures(float vaporPressure,
224 float saturationVaporPressure) {
225 if (vaporPressure != DataType.MissingData
226 && saturationVaporPressure != DataType.MissingData) {
227 return (vaporPressure / saturationVaporPressure);
228 }
229 return DataType.MissingData;
230 }
231
232 private static float calcVaporPressure(float dewpointC) {
233 // from Bolton (MWR 1980):
234 if (dewpointC != DataType.MissingData) {
235 return (float) (6.112 * Math.exp((17.67 * dewpointC)
236 / (dewpointC + 243.5)));
237 }
238 return DataType.MissingData;
239 }
240
241 private static float calcSaturationVaporPressure(float _TempC) {
242 // from Bolton (MWR 1980):
243 if (_TempC != DataType.MissingData) {
244 return (float) (6.112f * Math.exp((17.67f * _TempC)
245 / (_TempC + 243.5f)));
246 }
247 return DataType.MissingData;
248 }
249
250 private static float calcWetBulb(float vaporPressure,
251 float satVaporPressure, float DewPointC, float _TempC, float _Pres) {
252 // uses the AWIPS method for Tw calculation at a point
253 double ed = vaporPressure;
254 double es = satVaporPressure;
255 double DewPointK = (DewPointC + 273.16);
256 float _TempK = _TempC + 273.15f;
257
258 double ew, der;
259 double de = .0001;
260 int iter = 0;
261
262 double s = ((es - ed) / (_TempK - DewPointK));
263 double tw = (((.0006355 * _TempK * _Pres) + (DewPointK * s)) / ((.0006355 * _Pres) + s));
264
265 if (Double.isInfinite(tw) || Double.isNaN(tw))
266 return DataType.MissingData;
267
268 while (de >= .0001 && iter <= 10) {
269 ew = Math.exp(26.66082 - (.0091370924 * tw) - (6106.396 / tw));
270 de = ((.0006355 * _Pres) * (_TempK - tw)) - (ew - ed);
271 der = (ew * (.0091379024 - (6106.396 / (tw * tw))) - (.0006355 * _Pres));
272 ++iter;
273 tw = (tw - (de / der));
274 if (Double.isInfinite(tw) || Double.isNaN(tw))
275 return DataType.MissingData;
276 }
277
278 if (Double.isInfinite(tw) || Double.isNaN(tw))
279 return DataType.MissingData;
280
281 return (float) tw;
282 }
283
284 private static float calcVirtualTemperature(float _TempC, float _MixingRatio) {
285 float _TempK = _TempC + 273.15f;
286 if (_TempK != DataType.MissingData
287 && _MixingRatio != DataType.MissingData) {
288 // because moist air is less dense than dry air (Newton 1717)
289 // eq is originally Hess (1959)
290 double Wtemp = _MixingRatio < 0.1 ? 0 : _MixingRatio / 1000.0; // convert
291 // to
292 // kg/kg
293 // (dimensionless)
294
295 double ret_val = _TempK * (1.0 + 1.609 * Wtemp) / (1.0 + Wtemp);
296 return (float) ret_val;
297 }
298 return DataType.MissingData;
299
300 }
301
302 private static float calcMixingRatio(float _VaporPressure, float _Pres) {
303 if (_VaporPressure != DataType.MissingData
304 && _Pres != DataType.MissingData) {
305 // Hess (1959)
306 double epsillon = .62197;
307 double ret_val = (epsillon * _VaporPressure)
308 / (_Pres - _VaporPressure);
309 // convert to g/kg
310 ret_val *= 1000.0;
311 return (float) ret_val;
312 }
313 return DataType.MissingData;
314
315 }
316
317 private static float calcPotentialTemperature(float _Pres, float _TempK) {
318 if (_Pres != DataType.MissingData && _TempK != DataType.MissingData) {
319 return (float) (_TempK * Math.pow((1000.0 / _Pres), 0.286));
320 }
321 return DataType.MissingData;
322 }
323
324 private static float calcLclPres(float _LCLTemp, float _Pres, float _TempK) {
325 if (_LCLTemp != DataType.MissingData && _Pres != DataType.MissingData
326 && _TempK != DataType.MissingData) {
327 return (float) (_Pres * Math.pow(((_LCLTemp + 273.15) / _TempK),
328 (1 / 0.286)));
329 }
330 return DataType.MissingData;
331 }
332
333 private static float calcLclTemp(float tempC, float dewpointC) {
334 if (tempC != DataType.MissingData && dewpointC != DataType.MissingData) {
335 // From Barnes (JAM 1968, p511)
336 return dewpointC - (0.001296f * dewpointC + 0.1963f)
337 * (tempC - dewpointC);
338 } else {
339 return DataType.MissingData;
340 }
341 }
342
343 private static float calcEquivalentPotentialTemperature(float tempC,
344 float rh, float pres) {
345 float dewpoint = calcDewPoint(tempC, rh);
346 float vaporPressure = calcVaporPressure(dewpoint);
347 float mixingRatio = calcMixingRatio(vaporPressure, pres);
348 float lclTemp = calcLclTemp(tempC, dewpoint);
349 return calcEquivalentPotentialTemperature(tempC, mixingRatio, lclTemp,
350 pres);
351 }
352
353 private static float calcEquivalentPotentialTemperature(float tempC,
354 float mixingRatio, float lclTemp, float pres) {
355 float tempK = tempC + 273.16f;
356 // from Bolton (1980), equation #43
357 if (mixingRatio == DataType.MissingData
358 || tempK == DataType.MissingData
359 || lclTemp == DataType.MissingData) {
360 return DataType.MissingData;
361 }
362 float Wtemp = (mixingRatio < 0.1f) ? 0 : mixingRatio;
363 float exponent = 0.2854f * (1 - 0.00028f * Wtemp);
364 double term1 = tempK * Math.pow((1000 / pres), exponent);
365 double term2 = (3.376 / (273.15 + lclTemp)) - 0.00254;
366 double term3 = Wtemp * (1 + 0.00081 * Wtemp);
367 double thetae = term1 * Math.exp(term2 * term3);
368 if (Double.isInfinite(thetae) || Double.isNaN(thetae)) {
369 return DataType.MissingData;
370 }
371 return (float) thetae;
372 }
373
374 /** @return value of y corresponding to x */
375 private static float interpolate(float y1, float y2, float x1, float x2,
376 float x) {
377 if (y1 == DataType.MissingData || y2 == DataType.MissingData
378 || x1 == DataType.MissingData || x2 == DataType.MissingData) {
379 return DataType.MissingData;
380 }
381 return y1 + (x - x1) * (y2 - y1) / (x2 - x1);
382 }
383
384 void fillHeightOfIsotherm(float thresh, float[][] result) {
385 for (int i = 0; i < numLat; ++i) {
386 for (int j = 0; j < numLon; ++j) {
387 result[i][j] = sfcElevation[i][j];
388 for (int k = pressureLevels.length - 2; k >= 0; --k) {
389 if (DataType.MissingData != temperatureC[k + 1][i][j]
390 && temperatureC[k + 1][i][j] < thresh
391 && temperatureC[k][i][j] >= thresh) {
392 result[i][j] = interpolate(height[k][i][j],
393 height[k + 1][i][j], temperatureC[k][i][j],
394 temperatureC[k + 1][i][j], thresh);
395 break;
396 }
397 }
398 }
399 }
400 }
401
402 void fillHeightOfWetBulb(float threshInCelsius, float[][] result) {
403 final float threshInKelvin = 273.16f + threshInCelsius;
404
405 for (int i = 0; i < numLat; ++i) {
406 for (int j = 0; j < numLon; ++j) {
407 float lastWbTemp = DataType.MissingData;
408 for (int k = pressureLevels.length - 1; k >= 0; --k) {
409 float wbTemp = moistParcelTemp[k][i][j];
410 // moistParcelTemperature[k][i][j];
411 if (k != (pressureLevels.length - 1)) {
412 if (wbTemp >= threshInKelvin
413 && lastWbTemp < threshInKelvin) {
414 result[i][j] = interpolate(height[k][i][j],
415 height[k + 1][i][j], wbTemp, lastWbTemp,
416 threshInKelvin);
417 continue;
418 }
419 }
420 }
421 }
422 }
423 }
424
425 public void fillMeanWindSpeed(float botPressure, float topPressure,
426 float[][] result) {
427 for (int i = 0; i < numLat; ++i) {
428 for (int j = 0; j < numLon; ++j) {
429 float meanu = getMeanValuePressureLayer(i, j, botPressure, topPressure, this.u);
430 float meanv = getMeanValuePressureLayer(i, j, botPressure, topPressure, this.v);
431 result[i][j] = getSpeedFromUV(meanu, meanv);
432 }
433 }
434 }
435
436 static float getSpeedFromUV(float u, float v) {
437 if (u == DataType.MissingData || v == DataType.MissingData) {
438 return DataType.MissingData;
439 } else {
440 return (float) Math.sqrt(u * u + v * v);
441 }
442 }
443
444 float getMeanValueLayer(int x, int y, float botHeight, float topHeight,
445 float[][][] param) {
446 assert (param.length == pressureLevels.length);
447
448 boolean found_top = false;
449 boolean found_bot = false;
450
451 float sum_param = 0;
452
453 // FIXME: made this simple -- as a first cut, we won't bother
454 // interpolating to get the endpoints just right. Unless the
455 // model's vertical resolution is really crummy, then this
456 // probably doesn't matter, anyway!
457 //
458 for (int k = 0; k < pressureLevels.length - 1; k++) {
459 if (this.height[k][x][y] >= botHeight && !found_bot) {
460 found_bot = true;
461 }
462 if (this.height[k][x][y] >= topHeight && !found_top) {
463 if (k > 0) {
464 found_top = true;
465 } else {
466 return DataType.MissingData;
467 }
468 }
469
470 // assume that the parameter value represents the value from
471 // the bottom level to just below the next level (we aren't
472 // interpolating to get the exact values at the endpoints)
473
474 if (found_bot && !found_top) {
475 // FIXME: (a) check MISSING? (b) break loop at this point? why remain inside loop?
476 sum_param += param[k][x][y] * (height[k + 1][x][y] - height[k][x][y]);
477 }
478 }
479 if (found_bot && found_top) {
480 return (sum_param / (topHeight - botHeight));
481 } else {
482 return DataType.MissingData;
483 }
484 }
485
486 private float getMeanValuePressureLayer(int x, int y, float botPres,
487 float topPres, float[][][] values) {
488 float botHeight = DataType.MissingData;
489 float topHeight = DataType.MissingData;
490 boolean botFound = false;
491 boolean topFound = false;
492 for (int k = 0; k < pressureLevels.length - 1; ++k) {
493 if (pressureLevels[k] >= botPres && pressureLevels[k + 1] < botPres) {
494 botHeight = interpolate(this.height[k][x][y],
495 this.height[k + 1][x][y], pressureLevels[k],
496 pressureLevels[k + 1], botPres);
497 botFound = true;
498 } else if (pressureLevels[k] >= topPres
499 && pressureLevels[k + 1] < topPres) {
500 topHeight = interpolate(this.height[k][x][y],
501 this.height[k + 1][x][y], pressureLevels[k],
502 pressureLevels[k + 1], topPres);
503 topFound = true;
504 }
505 }
506 if (botFound && topFound) {
507 float val = getMeanValueLayer(x, y, botHeight, topHeight, values);
508 return val;
509 } else {
510 return DataType.MissingData;
511 }
512 }
513
514 public void fillMeanValuePressureLayerVerticalVelocity(int botPressure, int topPressure, float[][] result) {
515 for (int i = 0; i < numLat; ++i) {
516 for (int j = 0; j < numLon; ++j) {
517 result[i][j] = getMeanValuePressureLayer(i, j, botPressure, topPressure, this.w);
518 }
519 }
520 }
521
522 public void fillMeanValuePressureLayerTemperatureC(int botPressure, int topPressure, float[][] result) {
523 for (int i = 0; i < numLat; ++i) {
524 for (int j = 0; j < numLon; ++j) {
525 result[i][j] = getMeanValuePressureLayer(i, j, botPressure, topPressure, this.temperatureC);
526 }
527 }
528 }
529
530 public void fillNumberOfUndergroundGridpoints(float[][] result) {
531 for (int i = 0; i < numLat; ++i) {
532 for (int j = 0; j < numLon; ++j) {
533 result[i][j] = nearSurfaceGridPoint[i][j];
534 }
535 }
536 }
537
538 public void fillSurfaceRelativeHumidityPercent(float[][] result) {
539 for (int i = 0; i < numLat; ++i) {
540 for (int j = 0; j < numLon; ++j) {
541 float v = sfcRelativeHumidity[i][j];
542 result[i][j] = (v == DataType.MissingData)? v : (v*100);
543 }
544 }
545 }
546
547 }