diff --git a/src/opsinputs/opsinputs_cxwriter_mod.F90 b/src/opsinputs/opsinputs_cxwriter_mod.F90 index 79a1d9be..6767bbbe 100644 --- a/src/opsinputs/opsinputs_cxwriter_mod.F90 +++ b/src/opsinputs/opsinputs_cxwriter_mod.F90 @@ -1291,13 +1291,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) @@ -1308,61 +1312,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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9aa99f0f..5253a8a0 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -486,6 +486,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 diff --git a/test/generate_unittest_netcdfs.py b/test/generate_unittest_netcdfs.py index 2997649b..43df714b 100644 --- a/test/generate_unittest_netcdfs.py +++ b/test/generate_unittest_netcdfs.py @@ -1029,6 +1029,8 @@ def copy_var_to_var(Group, invarname, outvarname, filename): 'liquid_cloud_volume_fraction_in_atmosphere_layer', 'mass_content_of_cloud_ice_in_atmosphere_layer', 'mass_content_of_cloud_liquid_water_in_atmosphere_layer', 'specific_humidity', 'potential_temperature', 'eastward_wind', 'northward_wind'], 'testinput/cx_globalnamelist_scatwind.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') diff --git a/test/testinput/CxWriterNamelists_UnRotateWinds10M/Scatwind.nl b/test/testinput/CxWriterNamelists_UnRotateWinds10M/Scatwind.nl new file mode 100644 index 00000000..a244dc2e --- /dev/null +++ b/test/testinput/CxWriterNamelists_UnRotateWinds10M/Scatwind.nl @@ -0,0 +1,3 @@ +&CXControlNL +CxFields=3209,3210 +/ diff --git a/test/testinput/CxWriter_UnRotateWinds10M.nc4 b/test/testinput/CxWriter_UnRotateWinds10M.nc4 new file mode 100644 index 00000000..e3828340 Binary files /dev/null and b/test/testinput/CxWriter_UnRotateWinds10M.nc4 differ diff --git a/test/testinput/CxWriter_UnRotateWinds10M.yaml b/test/testinput/CxWriter_UnRotateWinds10M.yaml new file mode 100644 index 00000000..6669fb92 --- /dev/null +++ b/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