Update varobs and cx writers for visibility

Varfields.md

@@ -10,7 +10,7 @@
 | 3 | `VarField_rh` | `rh2(:)` or `rh(:,:)` | relative humidity | `ObsValue/relativeHumidity` (`ObsValue/relativeHumidityAt2M` for rh2) | |
 | 4 | `VarField_u` | `u10(:)` or `u(:,:)` | eastward wind | `ObsValue/windEastward` (`ObsValue/windEastwardAt10M` for u10) | |
 | 5 | `VarField_v` | `v10(:)` or `v(:,:)` | northward wind | `ObsValue/windNorthward` (`ObsValue/windNorthwardAt10M` for u10) | |
-| 6 | `VarField_logvis` | | | | Implement |
+| 6 | `VarField_logvis` | `logvis(:)` | Base 10 log of horizontal visibility (m) | `ObsValue/horizonalVisibility` | |
 | 7 | `VarField_tcwv` | `TCWV(:)` | total column water vapour | `ObsValue/precipitableWater` | |
 | 8 | `VarField_windspeed` | `WindSpeed(:)` | 10 metre windspeed | `ObsValue/WindSpeed` | |
 | 9 | `VarField_lwp` | | | | Implement |

etc/ukv/cx/Surface.nl

@@ -11,12 +11,12 @@
 !  4         IndexCxrh2              3245                  Yes
 !  5         IndexCxu10              3209                  Yes
 !  6         IndexCxv10              3210                  Yes
-!  8         IndexCxvis              3247                  No
+!  8         IndexCxvis              3247                  Yes
 ! 13         IndexCxTskinSea         24                    Yes
 ! 16         IndexCxpmsl             16222                 No
 ! 17         IndexCxSeaIce           31                    No
-! 20         IndexCxqt2              3255                  No
-! 21         IndexCxaerosol          90                    No
+! 20         IndexCxqt2              3255                  Yes
+! 21         IndexCxaerosol          90                    Yes
 ! 22         IndexCxPSurfParamA      20000                 No
 ! 23         IndexCxPSurfParamB      20001                 No
 ! 24         IndexCxCloudAmount      9217                  No
@@ -42,5 +42,5 @@
 ! 36         IndexCxRichNumber       3208                  No
 ! 37         IndexCxSoilMoisture     8223                  No
 ! 38         IndexCxSoilTemp         8225                  No
-CxFields=33,1,3236,3245,3209,3210,24,4,2,3,10,407,12,254
+CxFields=33,1,3236,3245,3209,3210,3247,24,3255,90,4,2,3,10,407,12,254
 /

src/opsinputs/opsinputs_cxfields_mod.F90

@@ -29,8 +29,8 @@ module opsinputs_cxfields_mod
 character(len=*), parameter, public :: opsinputs_cxfields_pmsl = var_pmsl
 character(len=*), parameter, public :: opsinputs_cxfields_SeaIce = var_sfc_ifrac
 character(len=*), parameter, public :: opsinputs_cxfields_SnowAmount = opsinputs_cxfields_unknown
-character(len=*), parameter, public :: opsinputs_cxfields_qt2 = opsinputs_cxfields_unknown
-character(len=*), parameter, public :: opsinputs_cxfields_aerosol = opsinputs_cxfields_unknown
+character(len=*), parameter, public :: opsinputs_cxfields_qt2 = "qt_1p5m"
+character(len=*), parameter, public :: opsinputs_cxfields_aerosol = "aerosol"
 character(len=*), parameter, public :: opsinputs_cxfields_PSurfParamA = "surf_param_a"
 character(len=*), parameter, public :: opsinputs_cxfields_PSurfParamB = "surf_param_b"
 character(len=*), parameter, public :: opsinputs_cxfields_LapseRate = opsinputs_cxfields_unknown

src/opsinputs/opsinputs_cxwriter_mod.F90

@@ -643,8 +643,8 @@ subroutine opsinputs_cxwriter_addrequiredgeovars(self, geovars)
       GeoVarName = opsinputs_cxfields_SnowAmount
     case (StashCode_qt2) ! IndexCxqt2
       GeoVarName = opsinputs_cxfields_qt2
-    ! wsmigaj: I haven't been able to identify the stash code associated with this field
-    !  GeoVarName = opsinputs_cxfields_aerosol
+    case (StashItem_aerosol) ! IndexCxaerosol
+     GeoVarName = opsinputs_cxfields_aerosol
     case (StashCode_PsurfParamA) ! IndexCxPsurfParamA
       GeoVarName = opsinputs_cxfields_PSurfParamA
     case (StashCode_PSurfParamB) ! IndexCxPSurfParamB
@@ -775,7 +775,7 @@ subroutine opsinputs_cxwriter_addrequiredgeovars(self, geovars)
       end if
 
   end select
-  
+
   if (GeoVarName /= opsinputs_cxfields_unknown) then
     call geovars % push_back(GeoVarName)
   end if
@@ -973,13 +973,11 @@ subroutine opsinputs_cxwriter_populatecx(self, ReportFlags, Cx)
       call opsinputs_fill_fillrealfromgeoval( &
         Cx % Header % qt2, "qt2", Cx % Header % NumLocal, Cx % qt2, &
         self % GeoVals, opsinputs_cxfields_qt2, self % JediToOpsLayoutMapping)
-    ! wsmigaj: I haven't been able to identify the stash code associated with this field
-    ! case (?) ! IndexCxaerosol
-      !  if (Cx % Header % ObsGroup == ObsGroupSurface) then
-      !    call opsinputs_fill_fillrealfromgeoval( &
-      !      Cx % Header % aerosol, "aerosol", Cx % Header % NumLocal, Cx % aerosol, &
-      !      self % GeoVals, opsinputs_cxfields_aerosol, self % JediToOpsLayoutMapping)
-      !  end if
+    case (StashItem_aerosol) ! IndexCxaerosol
+      ! Note this is treated specially, see opsinputs_fill_fillrealfromgeoval
+      call opsinputs_fill_fillrealfromgeoval( &
+        Cx % Header % aerosol, "aerosol", Cx % Header % NumLocal, Cx % aerosol, &
+        self % GeoVals, opsinputs_cxfields_aerosol, self % JediToOpsLayoutMapping)
     case (StashCode_PsurfParamA) ! IndexCxPsurfParamA
       call opsinputs_fill_fillrealfromgeoval( &
         Cx % Header % PSurfParamA, "PSurfParamA", Cx % Header % NumLocal, Cx % PSurfParamA, &

src/opsinputs/opsinputs_fill_mod.F90

@@ -1088,8 +1088,8 @@ subroutine opsinputs_fill_fillreal2d_norecords( &
 !e.g. HIRS in ATOVS stream
 
 numchans = size(JediVarNamesWithChannels)
-!sizeOfVarobsArray comes from intitial setting of size_of_varobs_array 
-! used to define the size of the channel array to fill. 
+!sizeOfVarobsArray comes from intitial setting of size_of_varobs_array
+! used to define the size of the channel array to fill.
 if (present(sizeOfVarobsArray)) then
   if (sizeOfVarobsArray > 0) then
     numchans = sizeOfVarobsArray
@@ -1121,7 +1121,7 @@ subroutine opsinputs_fill_fillreal2d_norecords( &
 	else
 	  exit
         end if
-      else 
+      else
         if (.not. compressChannels) then
           arrayindex = Channels(iChannel)
         end if
@@ -1362,6 +1362,8 @@ subroutine opsinputs_fill_fillrealfromgeoval( &
 ! Local declarations:
 type(ufo_geoval), pointer                       :: GeoVal
 real(kind_real)                                 :: MissingReal
+integer(integer64), dimension(2)                :: GeoValShape
+integer(integer64)                              :: NumLevs
 
 character(len=*), parameter                     :: &
   RoutineName = "opsinputs_fill_fillrealfromgeoval"
@@ -1380,17 +1382,33 @@ subroutine opsinputs_fill_fillrealfromgeoval( &
     call ufo_geovals_get_var(GeoVals, JediVarName, GeoVal, must_be_found = .false.)
     if (associated(GeoVal)) then
       if (GeoVal % nval /= 1) then
-        write (ErrorMessage, '("GeoVal ",A," contains more than one value per location. &
-          &Only the first of these values will be written to the VarObs file")') JediVarName
-        call gen_warn(RoutineName, ErrorMessage)
+        if (JediVarName /= "aerosol") then ! special case - see below
+          write (ErrorMessage, '("GeoVal ",A," contains ",I0," values per location. &
+            &Only the first of these values will be written to the VarObs file")') JediVarName, GeoVal % nval
+          call gen_warn(RoutineName, ErrorMessage)
+        end if
       end if
     end if
 
     ! Fill the OPS data structures
     call Ops_Alloc(Hdr, OpsVarName, NumObs, Real1)
-    where (GeoVal % vals(1,:) /= MissingReal)
-      Real1 = GeoVal % vals(1,:)
-    end where
+    if (JediVarName == "aerosol") then
+      ! This field "Total Aerosol (for Vis)" exists on 70 levels, but OPS only
+      ! outputs (and VAR expects) a single value at the surface (GeoVal level
+      ! 70 or cx level 1) for the calculation of visibility. See
+      ! Ops_BGEandCXCreate.inc line 166 and Ops_CxComplete.inc line 299. This
+      ! field is only used for surface visibility calculations so no further
+      ! specific handling is required.
+      GeoValShape = shape(GeoVaL % vals)
+      NumLevs = GeoValShape(1)
+      where (GeoVal % vals(NumLevs,:) /= MissingReal)
+        Real1 = GeoVal % vals(NumLevs,:)
+      end where
+    else
+      where (GeoVal % vals(1,:) /= MissingReal)
+        Real1 = GeoVal % vals(1,:)
+      end where
+    end if
 end if
 end subroutine opsinputs_fill_fillrealfromgeoval
 
@@ -2583,13 +2601,13 @@ subroutine opsinputs_fill_setpgefinal(PGE, MissingDouble, PackPGEs, Element)
 end if
 
 if (PackPGEs) then
-  ! For varfields which Ops_VarobPGEs expects PGEs in packed form: 
+  ! For varfields which Ops_VarobPGEs expects PGEs in packed form:
   ! pack consistent with OPS by multiplying by PPF and then truncate.
   Element % PGEFinal = Element % PGEFinal * PPF
   Element % PGEFinal = AINT(Element % PGEFinal)
 else
   ! Varfields which Ops_VarobPGEs expects PGEs in unpacked form,
-  ! To reduce the error (and avoid having to take the truncation error into account when preparing 
+  ! To reduce the error (and avoid having to take the truncation error into account when preparing
   ! known good outputs for tests), round the number after multiplying by PPF and then divide to undo.
   Element % PGEFinal = Element % PGEFinal * PPF
   Element % PGEFinal = NINT(Element % PGEFinal)

src/opsinputs/opsinputs_varobswriter_mod.F90

@@ -856,8 +856,9 @@ subroutine opsinputs_varobswriter_populateobservations( &
           ObsSpace, self % channels, Flags, ObsErrors, self % VarobsLength, "windNorthward", "ObsValue")
       end if
     case (VarField_logvis)
-      ! TODO(someone): handle this varfield
-      ! call Ops_Alloc(Ob % Header % logvis, "logvis", Ob % Header % NumObsLocal, Ob % logvis)
+      call opsinputs_fill_fillelementtypefromsimulatedvariable( &
+        Ob % Header % logvis, "logvis", Ob % Header % NumObsLocal, Ob % logvis, &
+        ObsSpace, Flags, ObsErrors, "horizontalVisibility", "ObsValue")
     case (VarField_tcwv)
       if (Ob % Header % ObsGroup == ObsGroupSatTCWV) then
         call opsinputs_fill_fillelementtypefromsimulatedvariable(Ob % Header % tcwv, "TCWV", Ob % Header % NumObsLocal, Ob % tcwv, &

test/CMakeLists.txt

@@ -145,6 +145,9 @@ ADD_WRITER_TEST(NAME     varobswriter_007_VarField_tcwv
 ADD_WRITER_TEST(NAME     varobswriter_008_VarField_WindSpeed
                 YAML     008_VarField_WindSpeed.yaml
                 DATA     008_VarField_WindSpeed.nc4)
+ADD_WRITER_TEST(NAME     varobswriter_009_VarField_horizontalVisibility
+                YAML     009_VarField_horizontalVisibility.yaml
+                DATA     009_VarField_horizontalVisibility.nc4)
 ADD_WRITER_TEST(NAME     varobswriter_010_VarField_britemp
                 YAML     010_VarField_britemp.yaml
                 NAMELIST VarObsWriterNamelists_010_VarField_britemp/AMSUB.nl

test/generate_unittest_netcdfs.py

@@ -700,6 +700,7 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
     output_simulated_var_profiles_to_netcdf('windNorthward',          'testinput/005_VarField_v_Sonde.nc4')
     output_1d_simulated_var_to_netcdf('precipitableWater',            'testinput/007_VarField_tcwv.nc4')
     output_1d_simulated_var_to_netcdf('windSpeed',                    'testinput/008_VarField_WindSpeed.nc4')
+    output_1d_simulated_var_to_netcdf('horizontalVisibility',         'testinput/009_VarField_horizontalVisibility.nc4')
     output_2d_simulated_var_to_netcdf('brightnessTemperature',      'testinput/010_VarField_britemp.nc4', with_bias=True)
     output_1d_normal_var_to_netcdf   ('skinTemperature', 'OneDVar', 'testinput/011_VarField_tskin.nc4')
     output_2d_normal_var_to_netcdf   ('cloudAmount', 'DerivedObsValue', 'testinput/015_VarField_cloud.nc4', use_chans=True)
@@ -1126,11 +1127,11 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
 
     # Surface - UKV
     output_full_cx_to_netcdf(['skin_temperature', 'surface_altitude', 'surface_pressure', 'uwind_at_10m',
-                              'vwind_at_10m', 'surface_temperature', 'relative_humidity_2m'],
+                              'vwind_at_10m', 'surface_temperature', 'relative_humidity_2m', 'visibility_1p5m','qt_1p5m'],
                              ['potential_temperature', 'specific_humidity', 'cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water',
                               'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water', 'air_pressure_levels',
                               'cloud_volume_fraction_in_atmosphere_layer', 'liquid_cloud_volume_fraction_in_atmosphere_layer', 'ice_cloud_volume_fraction_in_atmosphere_layer',
-                              'eastward_wind', 'northward_wind'],
+                              'eastward_wind', 'northward_wind', 'aerosol'],
                              'testinput/cx_ukvnamelist_surface.nc4')
 
     # Scatwind
@@ -1156,7 +1157,7 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
                              'testinput/cx_globalnamelist_screen.nc4')
 
     # Oceanwinds
-    output_full_cx_to_netcdf(['skin_temperature', 'ice_area_fraction', 'surface_altitude', 'surface_pressure', 
+    output_full_cx_to_netcdf(['skin_temperature', 'ice_area_fraction', 'surface_altitude', 'surface_pressure',
                               'uwind_at_10m', 'vwind_at_10m'],
                              ['air_pressure_levels'],
                              'testinput/cx_globalnamelist_oceanwinds.nc4')

test/testinput/009_VarField_horizontalVisibility.yaml

@@ -0,0 +1,47 @@
+time window:
+  begin: 2018-01-01T00:00:00Z
+  end: 2018-01-01T02:00:00Z
+
+observations:
+  - obs space:
+      name: Surface
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/009_VarField_horizontalVisibility.nc4
+      simulated variables: [horizontalVisibility]
+    obs filters:
+    # Double all observation errors: we want to check if error changes made by filters are
+    # propagated to VarObs files
+    - filter: BlackList
+      action:
+        name: inflate error
+        inflation factor: 2.0
+    # Set the flag of observations with missing values to "pass": we want to check if these
+    # values are encoded correctly in the VarObsFile.
+    - 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 VarObs file, as expected.
+    - filter: Domain Check
+      where:
+      - variable:
+          name: MetaData/latitude
+        minvalue: 0.0
+    - filter: VarObs Writer
+      reject_obs_with_any_variable_failing_qc: true
+      general_mode: debug
+    - filter: VarObs Checker
+      expected_main_table_columns:
+        # Only observations 1, 2 and 4 are passed; observation 3 is rejected by the domain check
+        field:    ["6", "6", "6"]
+        ob value: ["1.10000", "-1073741824.00000", "1.40000"]
+        ob error: ["0.20000", "-1073741824.00000", "0.80000"]
+        # 1.111 is the missing value indicator for PGEs. VarObs files store PGEs multiplied by 10000.
+        pge:      ["4210.00000", "11110.00000", "2380.00000"]
+        lat:      ["21.00000", "22.00000", "24.00000"]
+        lon:      ["31.00000", "32.00000", "34.00000"]
+        time:     ["-3540.00000", "-3480.00000", "-3360.00000"]
+        Callsign: ["station_1", "station_2", "station_4"]
+    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/cxwriter_ukvnamelist_surface.yaml

@@ -29,21 +29,21 @@ observations:
       general_mode: debug
       IC_PLevels: 5
     - filter: Cx Checker
-      expected_surface_variables: ["1","2","3","4","5","6","13"]
+      expected_surface_variables: ["1","2","3","4","5","6","8","13","20","21"]
       expected_upper_air_variables: ["1","3","4","5","11","29","30"]
       expected_main_table_columns:
       - # observation 3 is rejected by the tests above hence only 3 (1,2,4) columns
-        - ["17.10","27.10","57.10","67.10","37.10","47.10",
-        "7.10","1.10","1.20","1.30","81.10","81.20","81.30",
+        - ["17.10","27.10","57.10","67.10","37.10","47.10","77.10",
+        "7.10","87.10","101.10","1.10","1.20","1.30","81.10","81.20","81.30",
         "91.10","91.20","91.30","11.10","11.20","11.30",
         "41.10","41.20","41.30","21.10","21.20","21.30",
         "31.10","31.20","31.30"] # column 1 - 1st observation
-        - ["**********","**********","**********","**********","**********","**********","**********",
-        "2.10","**********","2.30","82.10","**********","82.30","92.10","**********","92.30","12.10",
+        - ["**********","**********","**********","**********","**********","**********","**********","**********",
+        "**********","102.10","2.10","**********","2.30","82.10","**********","82.30","92.10","**********","92.30","12.10",
         "**********","12.30","42.10","**********","42.30","22.10","**********","22.30","32.10","**********",
         "32.30"]  # column 2 - 2nd observation
-        - ["17.40","27.40","57.40","67.40","37.40","47.40",
-        "7.40","4.10","4.20","4.30","84.10","84.20","84.30",
+        - ["17.40","27.40","57.40","67.40","37.40","47.40","77.40",
+        "7.40","87.40","104.10","4.10","4.20","4.30","84.10","84.20","84.30",
         "94.10","94.20","94.30","14.10","14.20","14.30",
         "44.10","44.20","44.30","24.10","24.20","24.30",
         "34.10","34.20","34.30"] # column 3 - observation 4