From 67e8c5071d8bc37df1a862baaa4bbcd53ac63a7d Mon Sep 17 00:00:00 2001 From: PraveenKumar-NOAA <63739712+PraveenKumar-NOAA@users.noreply.github.com> Date: Wed, 7 Aug 2024 15:52:41 -0400 Subject: [PATCH] Add specificHumidity to ADPUPA BUFR DUMP Python API (#1246) This PR adds specificHumidity to the ADPUPA BUFR DUMP Python API. The following file has been changed: ~/ush/ioda/bufr2ioda/bufr2ioda_adpupa.py The specificHumidity is calculated using the dewPoint and pressure. The information is extracted from the prepBUFR codes and also documented. The validation of the output obs is performed using the specificHumidity from the prepBUFR data. For more details, please check: #1193 Co-authored-by: Praveen Singh --- ush/ioda/bufr2ioda/bufr2ioda_adpupa.py | 33 ++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/ush/ioda/bufr2ioda/bufr2ioda_adpupa.py b/ush/ioda/bufr2ioda/bufr2ioda_adpupa.py index 53dce26ad..fa0ce63fc 100755 --- a/ush/ioda/bufr2ioda/bufr2ioda_adpupa.py +++ b/ush/ioda/bufr2ioda/bufr2ioda_adpupa.py @@ -37,6 +37,21 @@ def Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd): return uob, vob +def Compute_SpecificHumidity_from_dewPoint_and_Pressure(tmdp, prlc): + + # est is saturated vapor pressure in Mb + est = 6.1078 * np.exp((17.269 * (tmdp - 273.16))/((tmdp - 273.16)+237.3)) + # Change est from Mb to Pa + est = est*100. + + # specificHumidity (qob) is in kg/kg + qob = (0.622 * est)/(prlc - (0.378 * est)) + qob = ma.array(qob) + qob = ma.masked_values(qob, qob.fill_value) + + return qob + + def bufr_to_ioda(config, logger): subsets = config["subsets"] @@ -191,6 +206,7 @@ def bufr_to_ioda(config, logger): dewpointTemperatureOE = np.float32(np.ma.masked_array(np.full((len(tmdp)), 0.0))) windEastwardOE = np.float32(np.ma.masked_array(np.full((len(wspd)), 0.0))) windNorthwardOE = np.float32(np.ma.masked_array(np.full((len(wspd)), 0.0))) + specificHumidityOE = np.float32(np.ma.masked_array(np.full((len(qob)), 0.0))) end_time = time.time() running_time = end_time - start_time @@ -208,6 +224,9 @@ def bufr_to_ioda(config, logger): logger.debug(f' uob min/max = {uob.min()} {uob.max()}') logger.debug(f' vob min/max = {vob.min()} {vob.max()}') + qob = Compute_SpecificHumidity_from_dewPoint_and_Pressure(tmdp, prlc) + logger.debug(f' qob min/max = {qob.min()} {qob.max()}') + end_time = time.time() running_time = end_time - start_time logger.info(f'Processing time for creating derived variables : {running_time} seconds') @@ -236,6 +255,7 @@ def bufr_to_ioda(config, logger): logger.debug(f' wspd shape = {wspd.shape}') logger.debug(f' uob shape = {uob.shape}') logger.debug(f' vob shape = {vob.shape}') + logger.debug(f' qob shape = {qob.shape}') logger.debug(f' qmpr shape = {qmpr.shape}') logger.debug(f' qmat shape = {qmat.shape}') @@ -264,6 +284,7 @@ def bufr_to_ioda(config, logger): logger.debug(f' wspd type = {wspd.dtype}') logger.debug(f' uob type = {uob.dtype}') logger.debug(f' vob type = {vob.dtype}') + logger.debug(f' qob type = {qob.dtype}') logger.debug(f' qmpr type = {qmpr.dtype}') logger.debug(f' qmat type = {qmat.dtype}') @@ -365,6 +386,12 @@ def bufr_to_ioda(config, logger): .write_attr('long_name', 'Northward Wind') \ .write_data(vob) + # Specific Humidity + obsspace.create_var('ObsValue/specificHumidity', dtype=qob.dtype, fillval=qob.fill_value) \ + .write_attr('units', 'kg kg-1') \ + .write_attr('long_name', 'Specific Humidity') \ + .write_data(qob) + # Pressure Quality Marker obsspace.create_var('QualityMarker/pressure', dtype=qmpr.dtype, fillval=qmpr.fill_value) \ .write_attr('long_name', 'Pressure Quality Marker') \ @@ -420,6 +447,12 @@ def bufr_to_ioda(config, logger): .write_attr('long_name', 'Northward Wind Observation Error') \ .write_data(windNorthwardOE) + # Specific Humidity Observation Error + obsspace.create_var('ObsError/specificHumidity', dtype=qob.dtype, fillval=qob.fill_value) \ + .write_attr('units', 'kg kg-1') \ + .write_attr('long_name', 'Specific Humidity Observation Error') \ + .write_data(specificHumidityOE) + end_time = time.time() running_time = end_time - start_time logger.info(f"Running time for splitting and output IODA: {running_time} seconds")