Unrotate 10m winds

src/opsinputs/opsinputs_cxwriter_mod.F90

@@ -1290,13 +1290,17 @@ subroutine opsinputs_cxwriter_unrotatewinds(self, Ob, Cx)
 character(len=80)                       :: ErrorMessage
 integer(integer64)                      :: CxFields(MaxModelCodes)
 
-integer                   :: Ilev         ! Loop variable
-real(real64), allocatable :: EqLat(:)     ! Latitudes on equatorial grid
-real(real64), allocatable :: EqLon(:)     ! Longitudes on equatorial grid
-real(real64), allocatable :: Coeff1(:)    ! Coefficients for rotation
-real(real64), allocatable :: Coeff2(:)    ! Coefficients for rotation
-real(real64), allocatable :: Uunrot(:)    ! Array for unrotated wind u component
-real(real64), allocatable :: Vunrot(:)    ! Array for unrotated wind v component
+integer                   :: Ilev                    ! Loop variable
+real(real64), allocatable :: EqLat(:)                ! Latitudes on equatorial grid
+real(real64), allocatable :: EqLon(:)                ! Longitudes on equatorial grid
+real(real64), allocatable :: Coeff1(:)               ! Coefficients for rotation
+real(real64), allocatable :: Coeff2(:)               ! Coefficients for rotation
+real(real64), allocatable :: Uunrot(:)               ! Array for unrotated wind u component
+real(real64), allocatable :: Vunrot(:)               ! Array for unrotated wind v component
+real(real64), allocatable :: U10unrot(:)             ! Array for unrotated wind u10 component
+real(real64), allocatable :: V10unrot(:)             ! Array for unrotated wind v10 component
+logical                   :: UpperWinds   = .false.  ! Upper air wind u and v components present
+logical                   :: SurfaceWinds = .false.  ! Surface wind u and v components present
 
 ! Body:
 call Ops_ReadCXControlNL(self % obsgroup, CxFields, BGECall = .false._8, ops_call = .false._8)
@@ -1307,61 +1311,91 @@ subroutine opsinputs_cxwriter_unrotatewinds(self, Ob, Cx)
 end if
 
 ! Require both wind components to be present.
-if (all(CxFields /= StashItem_u) .and. all(CxFields /= StashItem_v)) then
-   return
+if (any(CxFields == StashItem_u) .and. any(CxFields == StashItem_v)) then
+  UpperWinds = .true.
+end if
+if (any(CxFields == StashCode_u10) .and. any(CxFields == StashCode_v10)) then
+  SurfaceWinds = .true.
 end if
 
-! Code initially taken from Ops_RotateWinds.
-! The rotation matrix has been modified to perform an un-rotation of the winds.
-if (Cx % header % NumLocal > 0) then
-  allocate (Coeff1(Cx % header % NumLocal))
-  allocate (Coeff2(Cx % header % NumLocal))
-  allocate (EqLon(Cx % header % NumLocal))
-  allocate (EqLat(Cx % header % NumLocal))
-  allocate (Uunrot(Cx % header % NumLocal))
-  allocate (Vunrot(Cx % header % NumLocal))
-
-  ! Calculate longitudes on equatorial grid
-  call Gen_LatLon_to_Eq (Ob % latitude,      & ! in
-                         Ob % longitude,     & ! in
-                         EqLat(:),           & ! out
-                         EqLon(:),           & ! out
-                         self % RC_PoleLat,  & ! in
-                         self % RC_PoleLong)   ! in
-
-  deallocate (EqLat) ! Don't need latitudes
-
-  ! Get rotation coefficients for winds
-  call Ops_WCoeff (Coeff1(:),              & ! out
-                   Coeff2(:),              & ! out
-                   Ob % longitude,         & ! in
-                   EqLon(:),               & ! in
-                   self % RC_PoleLat,      & ! in
-                   self % RC_PoleLong,     & ! in
-                   Cx % header % NumLocal)   ! in
-
-  deallocate (EqLon)
-
-  ! Unrotate model level wind components.
-  ! Note minus sign in front of Coeff2, which reverses the previous rotation.
-  do Ilev = 1, Cx % Header % u % NumLev
-     call Ops_WEq_to_ll (Coeff1(:), & ! in
-          -Coeff2(:),               & ! in
-          Cx % u(:,Ilev),           & ! in
-          Cx % v(:,Ilev),           & ! in
-          Uunrot(:),                & ! out
-          Vunrot(:),                & ! out
-          Cx % header % NumLocal,   & ! in
-          Cx % header % NumLocal)     ! in
-     ! Put unrotated values into Cx structure
-     Cx % u(:,ILev) = Uunrot(:)
-     Cx % v(:,Ilev) = Vunrot(:)
-  end do
-
-  deallocate (Vunrot)
-  deallocate (Uunrot)
-  deallocate (Coeff2)
-  deallocate (Coeff1)
+if (UpperWinds .or. SurfaceWinds) then
+  ! Code initially taken from Ops_RotateWinds.
+  ! The rotation matrix has been modified to perform an un-rotation of the winds.
+  if (Cx % header % NumLocal > 0) then
+    allocate (Coeff1(Cx % header % NumLocal))
+    allocate (Coeff2(Cx % header % NumLocal))
+    allocate (EqLon(Cx % header % NumLocal))
+    allocate (EqLat(Cx % header % NumLocal))
+    if (UpperWinds) then
+      allocate (Uunrot(Cx % header % NumLocal))
+      allocate (Vunrot(Cx % header % NumLocal))
+    end if
+    if (SurfaceWinds) then
+      allocate (U10unrot(Cx % header % NumLocal))
+      allocate (V10unrot(Cx % header % NumLocal))
+    end if
+
+    ! Calculate longitudes on equatorial grid
+    call Gen_LatLon_to_Eq (Ob % latitude,      & ! in
+                           Ob % longitude,     & ! in
+                           EqLat(:),           & ! out
+                           EqLon(:),           & ! out
+                           self % RC_PoleLat,  & ! in
+                           self % RC_PoleLong)   ! in
+
+    deallocate (EqLat) ! Don't need latitudes
+
+    ! Get rotation coefficients for winds
+    call Ops_WCoeff (Coeff1(:),              & ! out
+                     Coeff2(:),              & ! out
+                     Ob % longitude,         & ! in
+                     EqLon(:),               & ! in
+                     self % RC_PoleLat,      & ! in
+                     self % RC_PoleLong,     & ! in
+                     Cx % header % NumLocal)   ! in
+
+    deallocate (EqLon)
+
+    ! Unrotate model level wind components.
+    ! Note minus sign in front of Coeff2, which reverses the previous rotation.
+    if (UpperWinds) then
+      do Ilev = 1, Cx % Header % u % NumLev
+         call Ops_WEq_to_ll (Coeff1(:), & ! in
+              -Coeff2(:),               & ! in
+              Cx % u(:,Ilev),           & ! in
+              Cx % v(:,Ilev),           & ! in
+              Uunrot(:),                & ! out
+              Vunrot(:),                & ! out
+              Cx % header % NumLocal,   & ! in
+              Cx % header % NumLocal)     ! in
+         ! Put unrotated values into Cx structure
+         Cx % u(:,ILev) = Uunrot(:)
+         Cx % v(:,Ilev) = Vunrot(:)
+      end do
+
+      deallocate (Vunrot)
+      deallocate (Uunrot)
+    end if
+
+    if (SurfaceWinds) then
+      call Ops_WEq_to_ll (Coeff1(:), & ! in
+           -Coeff2(:),               & ! in
+           Cx % u10(:),              & ! in
+           Cx % v10(:),              & ! in
+           U10unrot(:),              & ! out
+           V10unrot(:),              & ! out
+           Cx % header % NumLocal,   & ! in
+           Cx % header % NumLocal)     ! in
+      ! Put unrotated values into Cx structure
+      Cx % u10(:) = U10unrot(:)
+      Cx % v10(:) = V10unrot(:)
+      deallocate (V10unrot)
+      deallocate (U10unrot)
+    end if
+
+    deallocate (Coeff2)
+    deallocate (Coeff1)
+  end if
 end if
 
 end subroutine opsinputs_cxwriter_unrotatewinds

test/CMakeLists.txt

@@ -480,6 +480,10 @@ ADD_WRITER_TEST(NAME     cxwriter_UnRotateWinds
                 YAML     CxWriter_UnRotateWinds.yaml
                 NAMELIST CxWriterNamelists_UnRotateWinds/Sonde.nl
                 DATA     CxWriter_UnRotateWinds.nc4 dummy.nc4)
+ADD_WRITER_TEST(NAME     cxwriter_UnRotateWinds10M
+                YAML     CxWriter_UnRotateWinds10M.yaml
+                NAMELIST CxWriterNamelists_UnRotateWinds10M/Scatwind.nl
+                DATA     CxWriter_UnRotateWinds10M.nc4 dummy.nc4)
 
 # Header field tests
 ADD_WRITER_TEST(NAME     cxwriter_FixedHeader_VertCoord

test/generate_unittest_netcdfs.py

@@ -1019,5 +1019,8 @@ def copy_var_to_var(Group, invarname, outvarname, filename):
                               'eastward_wind', 'northward_wind'],
                              'testinput/cx_globalnamelist_surface.nc4')
 
+    # Unrotate 10m winds
+    output_full_cx_to_netcdf         (['uwind_at_10m', 'vwind_at_10m'],[], 'testinput/CxWriter_UnRotateWinds10M.nc4')
+
     output_1d_multi_level_simulated_var_to_netcdf('relativeHumidity', 'testinput/relative_humidity_Sonde.nc4')
     output_2d_geoval_for_multi_level_obs_to_netcdf('relative_humidity', 'testinput/002_UpperAirCxFieldForMultiLevelObs_relative_humidity.nc4')

test/testinput/CxWriterNamelists_UnRotateWinds10M/Scatwind.nl

@@ -0,0 +1,3 @@
+&CXControlNL
+CxFields=3209,3210
+/

test/testinput/CxWriter_UnRotateWinds10M.yaml

@@ -0,0 +1,80 @@
+window begin: 2018-01-01T00:00:00Z
+window end: 2018-01-01T01:00:00Z
+
+observations:
+  # Un-rotate winds on the LamNoWrapEq horizontal grid
+  - obs space:
+      name: Scatwind
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/dummy.nc4
+      simulated variables: [dummy]
+    geovals:
+      filename: Data/CxWriter_UnRotateWinds10M.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: latitude@MetaData
+        minvalue: 0.0
+    - filter: Cx Writer
+      namelist_directory: testinput/CxWriterNamelists_UnRotateWinds10M
+      reject_obs_with_any_variable_failing_qc: true
+      general_mode: debug
+      IC_PLevels: 5
+      FH_HorizGrid: LamNoWrapEq
+    - filter: Cx Checker
+      expected_surface_variables: ["5", "6"] # IndexCxu10, IndexCxv10
+      expected_upper_air_variables: []
+      expected_main_table_columns:
+      - # batch 1
+        - ["11.05", "-14.86"]          # column 1
+        - ["**********", "**********"] # column 2 (the asterisks represent a missing float)
+        - ["11.08","-15.32"]           # 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
+
+  # Un-rotate winds on the default horizontal grid. In this case there is no rotation so the original values are recovered.
+  - obs space:
+      name: Scatwind
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/dummy.nc4
+      simulated variables: [dummy]
+    geovals:
+      filename: Data/CxWriter_UnRotateWinds10M.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: latitude@MetaData
+        minvalue: 0.0
+    - filter: Cx Writer
+      namelist_directory: testinput/CxWriterNamelists_UnRotateWinds10M
+      reject_obs_with_any_variable_failing_qc: true
+      general_mode: debug
+      IC_PLevels: 5
+    - filter: Cx Checker
+      expected_surface_variables: ["5", "6"] # IndexCxu10, IndexCxv10
+      expected_upper_air_variables: []
+      expected_main_table_columns:
+      - # batch 1
+        - ["7.10", "17.10"]            # column 1
+        - ["**********", "**********"] # column 2 (the asterisks represent a missing float)
+        - ["7.40","17.40"]             # 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