From 15459232e7b53dc8898246c60e816b875d3e7b89 Mon Sep 17 00:00:00 2001 From: james-cotton Date: Mon, 27 Nov 2023 13:55:40 +0000 Subject: [PATCH] Unrotate 10m winds --- src/opsinputs/opsinputs_cxwriter_mod.F90 | 154 +++++++++++------- test/CMakeLists.txt | 4 + test/generate_unittest_netcdfs.py | 3 + .../Scatwind.nl | 3 + test/testinput/CxWriter_UnRotateWinds10M.nc4 | Bin 0 -> 6176 bytes test/testinput/CxWriter_UnRotateWinds10M.yaml | 80 +++++++++ 6 files changed, 184 insertions(+), 60 deletions(-) create mode 100644 test/testinput/CxWriterNamelists_UnRotateWinds10M/Scatwind.nl create mode 100644 test/testinput/CxWriter_UnRotateWinds10M.nc4 create mode 100644 test/testinput/CxWriter_UnRotateWinds10M.yaml diff --git a/src/opsinputs/opsinputs_cxwriter_mod.F90 b/src/opsinputs/opsinputs_cxwriter_mod.F90 index e9495ac9..83c06e6a 100644 --- a/src/opsinputs/opsinputs_cxwriter_mod.F90 +++ b/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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 880e0020..7ca8fcae 100644 --- a/test/CMakeLists.txt +++ b/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 diff --git a/test/generate_unittest_netcdfs.py b/test/generate_unittest_netcdfs.py index 66f4a1fe..e941dce2 100644 --- a/test/generate_unittest_netcdfs.py +++ b/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') 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 0000000000000000000000000000000000000000..e38283401ed3400d4d0e59356919f1624ea84078 GIT binary patch literal 6176 zcmeHIUr1A76hGT_Ym947GgcvZ^)bhAZiNuS#@%e<+`{=3#F*|bujppmb@pHs5tQ~A zNLqo^Qx8Ez&})x9R3JpX*FywBQNG(b_q%86QX~*6xF5Sa-=FiH?|i@e`<=U(j1M+@ zPk5ak4{+PGnsQ9c^jZtcQ(7RCjE_3y=gw|(PHM82cuS-8I>A9MJ|F!DT=WlcNf-dh z>WJmXZZRC@WvA%4JIOkltRcvA^LnX}%Nw~+@65av<7fv^GLjP-bKvppkDDLN^qcz~ zj&=&^bT}M*1GZ8-W!&1;X@&z(!0!S7RrGe5cagl5tOV&p6anC&Fu>=#_B>ziAy*09 z;Gx)RDwfT{5u74y=x*_$kwiK>lu75ZYAlri`8?hlH@L)xbZjJH`bRDQu}QtE;4iO~ zG((LKDg}K;D^>M!Ng1Crlv3H)TB_wMdVYLb3o6_qKnXc9{$Y?QpN^KjKdn-bY5#|!KBLPqlOaQf9t`ziA-q5N$X@GK} z3s7TAeK{p~+Ja2RFk3WE4Q0pRIC@Y@etf9Lkmz{#4r>)4np=N$uPdV_hu!9t>KX-C zX}Ua7{e~XYej6xkOXDNgD_4W4@ndQi(jpZmWP(hWT1$s#(qlv-9 zXdT*19*74;nekcV<)Ci&uX0KI4@hcXImQq8vz>u8vz>u8vz>u8vz>u8vz>u8vz>u8-afffpGX! d^wI2<&DE7J(KoNZ;&UrT<6?0s77pKz{Q|m=h=%|G literal 0 HcmV?d00001 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