Radar reflectivity (#223)

* reflectivity

* exner_levels

Varfields.md

@@ -56,19 +56,19 @@
 | 60 | `VarField_z` | | | | Implement |
 | 61 | `VarField_SBUVozone` | | | | Implement |
 | 62 | `VarField_GeoBriTemp` | | | | Implement |
-| 63 | `VarField_RadialVelocity` | | | | Implement |
-| 64 | `VarField_RadarBeamElev` | | | | Implement |
-| 65 | `VarField_RadarObRange` | | | | Implement |
-| 66 | `VarField_RadarObAzim` | `RadarObAzim(:,:)` | Azimut of radar ob | `MetaData/radarAzimuth` | |
+| 63 | `VarField_RadialVelocity` | `RadialVelocSO(:,:)` | Radar radial velocity | `ObsValue/radialVelocity` | |
+| 64 | `VarField_RadarBeamElev` | `RadarBeamElev(:,:)` | Radar beam tilt | `MetaData/beamTiltAngle` | |
+| 65 | `VarField_RadarObRange` | `RadarObRange(:,:)` | Radar gate range | `MetaData/gateRange` | |
+| 66 | `VarField_RadarObAzim` | `RadarObAzim(:,:)` | Radar beam azimuth | `MetaData/beamAzimuthAngle` | |
 | 67 | `VarField_GPS_Station_Height` | `GPS_Station_Height(:)`| Height of the Ground GNSS stations | `MetaData/stationElevation` | |
 | 68 | `VarField_clw` | CLW(:,:) | retrieved clw from 1Dvar | OneDVar/liquidWaterContent | |
-| 69 | `VarField_RadIdent` | | | | Implement Implement |
-| 70 | `VarField_Reflectivity` | | | | Implement |
+| 69 | `VarField_RadIdent` | `RadIdent(:)` | Radar numerical identifier | `MetaData/stationIdentification` | |
+| 70 | `VarField_Reflectivity` | `ReflectivitySO(:,:)` | Radar reflectivity | `ObsValue/reflectivity` | |
 | 71 | `VarField_BendingAngle` | `BendingAngle(:,:)` or `BendingAngleAll(:,:)` | GPSRO bending angle (L1, L2, LC etc or just LC) | `ObsValue/bendingAngle` | Decide if a separate variable corresponding to BendingAngleAll is needed |
 | 72 | `VarField_ImpactParam` | `ImpactParam(:,:)` or `ImpactParamAll(:,:)` | GPSRO impact parameter (L1, L2, LC etc or just LC) | `ObsValue/impactParameterRO` | Decide if a separate variable corresponding to ImpactParamAll is needed |
 | 73 | `VarField_RO_Rad_Curv` | `RO_Rad_Curv(:)` | GPSRO Earth\'s local radius of curvature (m) | `MetaData/earthRadiusCurvature` | |
 | 74 | `VarField_RO_geoid_und` | `RO_geoid_und(:)` | GPSRO Geoid undulation (above WGS-84 ellipsoid, m) | `MetaData/geoidUndulation` | |
-| 75 | `VarField_RadAltAboveMSL` | | | | Implement |
+| 75 | `VarField_RadAltAboveMSL` | `RadAltAboveMSL(:,:)` | Radar station elevation above mean sea level | `MetaData/stationElevation` | |
 | 76 | `VarField_BriTempVarError` | `BriTempVarError(:,:)` | scene dependent obs errors for BTs (oblocation, channel) | `ObsErrorData/brightnessTemperature` | |
 | 77 | `VarField_AOD` | `AOD(:,:)` | aerosol optical depth | `ObsValue/aerosolOpticalDepth` | |
 | 78 | `VarField_Theta` | `theta(:,:)` | potential temperature | `ObsValue/airTemperature` | |
@@ -97,4 +97,4 @@
 | 101 | `VarField_dWinddP` | | | | Implement |
 | 102 | `VarField_AzimuthCOG` | | | | Implement |
 | 103 | `VarField_HeightCOG` | | | | Implement |
-| 104 | `VarField_RadFlag` | | | | | | Implement |
+| 104 | `VarField_RadFlag` | `RadFlag(:,:)` | Radar preprocessing flag | `QualityInformation/reflectivity` | |

etc/ukv/cx/RadarZ.nl

@@ -0,0 +1,12 @@
+&CXControlNL
+! These STASH codes correspond to the following variables:
+! 4: theta
+! 10: q
+! 12: qcf
+! 33: orography
+! 254: qcl
+! 255: exnerA
+! 272: qrain
+! 407: p
+CxFields=4,10,12,33,254,255,272,407
+/

etc/ukv/varobs/RadarZ.nl

@@ -0,0 +1,11 @@
+&VarobsControlNL
+! These indices correspond to the following variables:
+! 64: Beam tilt (elevation in OPS terminology)
+! 65: Gate range
+! 66: Gate azimuth
+! 69: Station identifier
+! 70: Reflectivity
+! 75: Station elevation (altitude above MSL in OPS terminology)
+! 104: Radar QC flag
+Varfields=64,65,66,69,70,75,104
+/

src/opsinputs/opsinputs_cxfields_mod.F90

@@ -86,8 +86,8 @@ module opsinputs_cxfields_mod
 character(len=*), parameter, public :: opsinputs_cxfields_RH_AMC = opsinputs_cxfields_unknown
 character(len=*), parameter, public :: opsinputs_cxfields_Cl = "liquid_cloud_volume_fraction_in_atmosphere_layer"
 character(len=*), parameter, public :: opsinputs_cxfields_Cf = "ice_cloud_volume_fraction_in_atmosphere_layer"
-character(len=*), parameter, public :: opsinputs_cxfields_qrain = opsinputs_cxfields_unknown
-character(len=*), parameter, public :: opsinputs_cxfields_ExnerA = opsinputs_cxfields_unknown
+character(len=*), parameter, public :: opsinputs_cxfields_qrain = "qrain"
+character(len=*), parameter, public :: opsinputs_cxfields_ExnerA = "exner_levels"
 character(len=*), parameter, public :: opsinputs_cxfields_RichNumber = opsinputs_cxfields_unknown
 character(len=*), parameter, public :: opsinputs_cxfields_SoilMoisture = opsinputs_cxfields_unknown
 character(len=*), parameter, public :: opsinputs_cxfields_SoilTemp = opsinputs_cxfields_unknown

src/opsinputs/opsinputs_varobswriter_mod.F90

@@ -105,6 +105,7 @@ module opsinputs_varobswriter_mod
     ObsGroupAircraft,          &
     ObsGroupGPSRO,             &
     ObsGroupOceanWinds,        &
+    ObsGroupRadarZ,            &
     ObsGroupSatwind,           &
     ObsGroupScatwind,          &
     ObsGroupScatwindChosen,    &
@@ -707,7 +708,7 @@ subroutine opsinputs_varobswriter_populateobservations( &
 integer                                               :: nVarFields
 integer                                               :: iVarField
 integer                                               :: iobs
-
+character(len=200)                                    :: varname
 logical                                               :: FillChanNum = .false.
 logical                                               :: FillNumChans = .false.
 
@@ -1034,8 +1035,14 @@ subroutine opsinputs_varobswriter_populateobservations( &
          Ob % Header % RadialVelocity, "RadialVelocity", JediToOpsLayoutMapping, Ob % RadialVelocity, &
          ObsSpace, self % channels, Flags, ObsErrors, self % VarobsLength, "radialVelocity", "ObsValue")
     case (VarField_Reflectivity)
-      ! TODO(someone): handle this varfield
-      ! call Ops_Alloc(Ob % Header % ReflectivitySO, "ReflectivitySO", Ob % Header % NumObsLocal, Ob % ReflectivitySO)
+       ! Write DerivedObsValue/reflectivity to both Ob % ReflectivitySO and Ob % Reflectivity.
+       ! See the explanation given for RadialVelocity above.
+      call opsinputs_fill_fillelementtype2dfromsimulatedvariable( &
+         Ob % Header % ReflectivitySO, "ReflectivitySO", JediToOpsLayoutMapping, Ob % ReflectivitySO, &
+         ObsSpace, self % channels, Flags, ObsErrors, self % VarobsLength, "reflectivity", "ObsValue")
+      call opsinputs_fill_fillelementtype2dfromsimulatedvariable( &
+         Ob % Header % Reflectivity, "Reflectivity", JediToOpsLayoutMapping, Ob % Reflectivity, &
+         ObsSpace, self % channels, Flags, ObsErrors, self % VarobsLength, "reflectivity", "ObsValue")
     case (VarField_ReflectivityR)
       ! TODO(someone): handle this varfield
       ! call Ops_Alloc(Ob % Header % ReflectivityR, "ReflectivityR", Ob % Header % NumObsLocal, Ob % ReflectivityR)
@@ -1066,8 +1073,14 @@ subroutine opsinputs_varobswriter_populateobservations( &
       ! TODO(someone): handle this varfield
       ! call Ops_Alloc(Ob % Header % RadNoiseLvl, "RadNoiseLvl", Ob % Header % NumObsLocal, Ob % RadNoiseLvl)
     case (VarField_RadFlag)
-      ! TODO(someone): handle this varfield
-      ! call Ops_Alloc(Ob % Header % RadFlag, "RadFlag", Ob % Header % NumObsLocal, Ob % RadFlag)
+      if (Ob % header % ObsGroup == ObsGroupRadarZ) then
+         varname = "reflectivity"
+      else
+         call gen_fail(RoutineName, "Invalid observation name for RadFlag")
+      end if
+      call opsinputs_fill_fillinteger2d( &
+         Ob % Header % RadFlag, "RadFlag", JediToOpsLayoutMapping, Ob % RadFlag, &
+         ObsSpace, self % channels, self % VarobsLength, trim(varname), "QualityInformation")
     case (VarField_clw)
       call opsinputs_fill_fillelementtype2dfromnormalvariablewithlevels( &
         Ob % Header % CLW , "CLW" , Ob % Header % NumObsLocal, ob % CLW, &

test/CMakeLists.txt

@@ -234,6 +234,10 @@ ADD_WRITER_TEST(NAME     varobswriter_068_VarField_clw
                 YAML     068_VarField_clw.yaml
                 NAMELIST VarObsWriterNamelists_068_VarField_clw/AMSUB.nl
                 DATA     068_VarField_clw.nc4)
+ADD_WRITER_TEST(NAME     varobswriter_070_VarField_reflectivity
+                YAML     070_VarField_reflectivity.yaml
+                NAMELIST VarObsWriterNamelists_070_VarField_reflectivity/RadarZ.nl
+                DATA     070_VarField_reflectivity.nc4)
 ADD_WRITER_TEST(NAME     varobswriter_071_VarField_bendingangle
                 YAML     071_VarField_bendingangle.yaml
                 DATA     071_VarField_bendingangle.nc4)
@@ -382,6 +386,10 @@ ADD_WRITER_TEST(NAME     varobswriter_ukvnamelist_radar_doppler_wind
                 YAML     varobswriter_ukvnamelist_radar_doppler_wind.yaml
                 NAMELIST ../../etc/ukv/varobs/Radar.nl
                 DATA     varobs_ukvnamelist_radar_doppler_wind.nc4)
+ADD_WRITER_TEST(NAME     varobswriter_ukvnamelist_radar_reflectivity
+                YAML     varobswriter_ukvnamelist_radar_reflectivity.yaml
+                NAMELIST ../../etc/ukv/varobs/RadarZ.nl
+                DATA     varobs_ukvnamelist_radar_reflectivity.nc4)
 
 ### CxWriter tests
 
@@ -504,7 +512,14 @@ ADD_WRITER_TEST(NAME     cxwriter_035_UpperAirCxField_Cl
                 YAML     035_UpperAirCxField_Cl.yaml
                 NAMELIST CxWriterNamelists_035_UpperAirCxField_Cl/AMSUB.nl
                 DATA     035_UpperAirCxField_Cl.nc4 dummy.nc4)
-
+ADD_WRITER_TEST(NAME     cxwriter_039_UpperAirCxField_Exner
+                YAML     039_UpperAirCxField_Exner.yaml
+                NAMELIST CxWriterNamelists_039_UpperAirCxField_Exner/RadarZ.nl
+                DATA     039_UpperAirCxField_Exner.nc4 dummy.nc4)
+ADD_WRITER_TEST(NAME     cxwriter_040_UpperAirCxField_Qrain
+                YAML     040_UpperAirCxField_Qrain.yaml
+                NAMELIST CxWriterNamelists_040_UpperAirCxField_Qrain/RadarZ.nl
+                DATA     040_UpperAirCxField_Qrain.nc4 dummy.nc4)
 ADD_WRITER_TEST(NAME     cxwriter_041-042_UpperAirCxField_dust1-dust2
                 YAML     041-042_UpperAirCxField_dust1-dust2.yaml
                 NAMELIST CxWriterNamelists_041-042_UpperAirCxField_dust1-dust2/AOD.nl
@@ -704,6 +719,10 @@ ADD_WRITER_TEST(NAME     cxwriter_ukvnamelist_radar_doppler_wind
                 YAML     cxwriter_ukvnamelist_radar_doppler_wind.yaml
                 NAMELIST ../../etc/ukv/cx/Radar.nl
                 DATA     cx_ukvnamelist_radar_doppler_wind.nc4 dummy.nc4)
+ADD_WRITER_TEST(NAME     cxwriter_ukvnamelist_radar_reflectivity
+                YAML     cxwriter_ukvnamelist_radar_reflectivity.yaml
+                NAMELIST ../../etc/ukv/cx/RadarZ.nl
+                DATA     cx_ukvnamelist_radar_reflectivity.nc4 dummy.nc4)
 
 ### Tests of auxiliary classes
 

test/generate_unittest_netcdfs.py

@@ -16,7 +16,8 @@
 missing_float_nc = 9.969209968386869e+36
 
 
-def output_1d_simulated_var_to_netcdf(var_name, file_name, with_bias=False, radar_doppler_wind=False):
+def output_1d_simulated_var_to_netcdf(var_name, file_name, with_bias=False,
+                                      radar_doppler_wind=False, radar_reflectivity=False):
     f = nc4.Dataset(file_name, 'w', format="NETCDF4")
 
     nlocs = 4
@@ -35,10 +36,10 @@ def output_1d_simulated_var_to_netcdf(var_name, file_name, with_bias=False, rada
     var[:] = [1 * minute, 2 * minute, 3 * minute, 4 * minute]
 
     # Station ID
-    if radar_doppler_wind:
-        # The radar Doppler wind processing uses integers for station identification because
-        # `MetaData/stationIdentification` is mapped to the ODB variable `rad_ident`,
-        # which is an integer.
+    if radar_doppler_wind or radar_reflectivity:
+        # The radar Doppler wind and reflectivity processing both use integers for
+        # station identification because `MetaData/stationIdentification` is mapped
+        # to the ODB variable `rad_ident`, which is an integer.
         var = f.createVariable('MetaData/stationIdentification', 'i', ('Location'))
         var[:] = [1, 2, 3, 4]
     else:
@@ -48,7 +49,7 @@ def output_1d_simulated_var_to_netcdf(var_name, file_name, with_bias=False, rada
             var[i] = s
 
     # Extra variables for radar
-    if radar_doppler_wind:
+    if radar_doppler_wind or radar_reflectivity:
         var = f.createVariable('MetaData/radar_family', 'i', ('Location',))
         var[:] = [11, 12, 13, 14]
         var = f.createVariable('MetaData/beamTiltAngle', 'i', ('Location',))
@@ -59,6 +60,9 @@ def output_1d_simulated_var_to_netcdf(var_name, file_name, with_bias=False, rada
         var[:] = [11, 12, 13, 14]
         var = f.createVariable('MetaData/stationElevation', 'i', ('Location',))
         var[:] = [11, 12, 13, 14]
+    if radar_reflectivity:
+        var = f.createVariable('QualityInformation/reflectivity', 'i', ('Location',))
+        var[:] = [11, 12, 13, 14]
 
     var = f.createVariable('ObsValue/' + var_name, 'f', ('Location',))
     obsVal = [1.1, missing_float, 1.3, 1.4]
@@ -726,6 +730,7 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
                                                                  'thickness_200_50hPa_satid_5Predictor', 'thickness_200_50hPa_satid_8Predictor',
                                                                  'legendre_order_1_satid_5Predictor',    'legendre_order_1_satid_8Predictor'],
                                       'testinput/080_VarField_biaspredictors.nc4', predictors=True)
+    output_1d_simulated_var_to_netcdf   ('reflectivity',  'testinput/070_VarField_reflectivity.nc4', radar_reflectivity=True)
     output_2d_simulated_var_to_netcdf('bendingAngle',              'testinput/071_VarField_bendingangle.nc4', add_occulting_satid=True)
     output_2d_normal_var_to_netcdf('impactParameterRO', 'MetaData', 'testinput/072_VarField_impactparam.nc4')
     output_1d_normal_var_to_netcdf('earthRadiusCurvature', 'MetaData',  'testinput/073_VarField_ro_rad_curv.nc4')
@@ -901,6 +906,18 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
                                  ['MetaData/stationIdentification'],
                                  'testinput/varobs_ukvnamelist_radar_doppler_wind.nc4')
 
+    # Radar reflectivity - UKV
+    output_full_varobs_to_netcdf(['MetaData/latitude',
+                                  'MetaData/longitude',
+                                  'MetaData/beamTiltAngle',
+                                  'MetaData/gateRange',
+                                  'MetaData/beamAzimuthAngle',
+                                  'MetaData/stationElevation',
+                                  'ObsValue/reflectivity', 'ObsError/reflectivity'],
+                                 [],
+                                 ['MetaData/stationIdentification', 'QualityInformation/reflectivity',],
+                                 'testinput/varobs_ukvnamelist_radar_reflectivity.nc4')
+
 
 
 
@@ -936,6 +953,8 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
     output_2d_geoval_to_netcdf       ('liquid_cloud_volume_fraction_in_atmosphere_layer',      'testinput/035_UpperAirCxField_Cl.nc4')
     output_2d_geovals_to_netcdf      (['mass_fraction_of_dust00%s_in_air' % i for i in range(1, 7)], 'testinput/041-046_UpperAirCxField_dust1-dust6.nc4')
     output_2d_geovals_to_netcdf      (['eastward_wind', 'northward_wind'], 'testinput/CxWriter_UnRotateWinds.nc4', shift_by_varindex=False)
+    output_2d_geoval_to_netcdf       ('exner_levels', 'testinput/039_UpperAirCxField_Exner.nc4')
+    output_2d_geoval_to_netcdf       ('qrain',      'testinput/040_UpperAirCxField_Qrain.nc4')
 
     # Cx full output for an obsgroup testing
     # list of 1d-variables; list of 2d-variables; filename for output
@@ -1167,3 +1186,11 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
                              ['eastward_wind', 'northward_wind', 'upward_air_velocity'],
                              'testinput/cx_ukvnamelist_radar_doppler_wind.nc4')
 
+    # Radar reflectivity - UKV
+    output_full_cx_to_netcdf(['surface_altitude'],
+                             ['potential_temperature', 'specific_humidity', 'air_pressure_levels',
+                              'cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water',
+                              'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water',
+                              'exner_levels', 'qrain'],
+                             'testinput/cx_ukvnamelist_radar_reflectivity.nc4')
+

test/testinput/039_UpperAirCxField_Exner.yaml

@@ -0,0 +1,41 @@
+time window:
+  begin: 2018-01-01T00:00:00Z
+  end: 2018-01-01T01:00:00Z
+
+observations:
+  - obs space:
+      name: RadarZ
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/dummy.nc4
+      simulated variables: [dummy]
+    geovals:
+      filename: Data/039_UpperAirCxField_Exner.nc4
+    obs filters:
+    # Set the flag of observations with missing values to "pass": we want to check if these
+    # values are encoded correctly in the Cx file.
+    - filter: Reset Flags to Pass
+      flags_to_reset: [10, 15]  # missing, Hfailed
+    # Reject observation 3: we want to check if it is omitted from the Cx file, as expected.
+    - filter: Domain Check
+      where:
+      - variable:
+          name: MetaData/latitude
+        minvalue: 0.0
+    - filter: Cx Writer
+      namelist_directory: testinput/CxWriterNamelists_039_UpperAirCxField_Exner
+      reject_obs_with_any_variable_failing_qc: true
+      general_mode: debug
+      IC_PLevels: 5
+    - filter: Cx Checker
+      expected_surface_variables: []
+      expected_upper_air_variables: ["39"] # IndexCxExner
+      expected_main_table_columns:
+      - # batch 1
+        - ["1.10", "1.20", "1.30"]       # column 1
+        - ["2.10", "**********", "2.30"] # column 2 (the asterisks represent a missing float)
+        - ["4.10", "4.20", "4.30"]       # column 3
+    HofX: ObsValue # just a placeholder -- not used, but needed to force calls to postFilter.
+    benchmarkFlag: 1000 # just to keep the ObsFilters test happy
+    flaggedBenchmark: 0

test/testinput/040_UpperAirCxField_Qrain.yaml

@@ -0,0 +1,41 @@
+time window:
+  begin: 2018-01-01T00:00:00Z
+  end: 2018-01-01T01:00:00Z
+
+observations:
+  - obs space:
+      name: RadarZ
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/dummy.nc4
+      simulated variables: [dummy]
+    geovals:
+      filename: Data/040_UpperAirCxField_Qrain.nc4
+    obs filters:
+    # Set the flag of observations with missing values to "pass": we want to check if these
+    # values are encoded correctly in the Cx file.
+    - filter: Reset Flags to Pass
+      flags_to_reset: [10, 15]  # missing, Hfailed
+    # Reject observation 3: we want to check if it is omitted from the Cx file, as expected.
+    - filter: Domain Check
+      where:
+      - variable:
+          name: MetaData/latitude
+        minvalue: 0.0
+    - filter: Cx Writer
+      namelist_directory: testinput/CxWriterNamelists_040_UpperAirCxField_Qrain
+      reject_obs_with_any_variable_failing_qc: true
+      general_mode: debug
+      IC_PLevels: 5
+    - filter: Cx Checker
+      expected_surface_variables: []
+      expected_upper_air_variables: ["40"] # IndexCxQrain
+      expected_main_table_columns:
+      - # batch 1
+        - ["1.10", "1.20", "1.30"]       # column 1
+        - ["2.10", "**********", "2.30"] # column 2 (the asterisks represent a missing float)
+        - ["4.10", "4.20", "4.30"]       # column 3
+    HofX: ObsValue # just a placeholder -- not used, but needed to force calls to postFilter.
+    benchmarkFlag: 1000 # just to keep the ObsFilters test happy
+    flaggedBenchmark: 0

test/testinput/063_VarField_radialVelocity.yaml

@@ -23,7 +23,7 @@ observations:
         minvalue: 0.0
     - filter: VarObs Writer
       reject_obs_with_any_variable_failing_qc: true
-      namelist_directory: testinput/VarObsWriterNamelists_021_VarField_surface
+      namelist_directory: testinput/VarObsWriterNamelists_063_VarField_radialVelocity
       general_mode: debug
       station_ID_int_to_string: true
     - filter: VarObs Checker