Merge branch 'develop' into feature/bufr2ioda

NOAA-EMC · Aug 8, 2024 · 182086d · 182086d
2 parents 4e43edd + 67e8c50
commit 182086d
Showing 1 changed file with 33 additions and 0 deletions.
diff --git 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")