From c12a53a3fa14597152136a9fd5d9d4d3a2e8581d Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Tue, 20 Aug 2024 22:27:07 +0000 Subject: [PATCH 01/13] =?UTF-8?q?Add=20functionality=20to=20adjust=20obser?= =?UTF-8?q?ved=20surface=20temperature=20to=20model=20surface=20height=20u?= =?UTF-8?q?sing=20a=20lapse=20rate=20(local/constant)=20=E2=80=94=20implem?= =?UTF-8?q?entation=20pending.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../sfccorrected/ObsSfcCorrectedParameters.h | 11 +++++++++-- .../operators/sfccorrected/ufo_sfccorrected_mod.F90 | 7 ++++++- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index 22f8ecc..c161a0e 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -22,7 +22,7 @@ namespace ufo { /// enum type for surface correction type, and ParameterTraitsHelper for it enum class SfcCorrectionType { - UKMO, WRFDA, RRFS_GSL + UKMO, WRFDA, RRFS_GSL, RRFS_GSL_LOCAL }; struct SfcCorrectionTypeParameterTraitsHelper { typedef SfcCorrectionType EnumType; @@ -30,7 +30,8 @@ struct SfcCorrectionTypeParameterTraitsHelper { static constexpr util::NamedEnumerator namedValues[] = { { SfcCorrectionType::UKMO, "UKMO" }, { SfcCorrectionType::WRFDA, "WRFDA" }, - { SfcCorrectionType::RRFS_GSL, "RRFS_GSL" } + { SfcCorrectionType::CONSTANT_LAPSE_RATE, "CONSTANT_LAPSE_RATE" }, + { SfcCorrectionType::GSL, "GSL" } }; }; @@ -81,6 +82,12 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { /// Note: "station_altitude" default value is "stationElevation" oops::Parameter ObsHeightName{"station_altitude", "stationElevation", this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to CONSTANT_LAPSE_RATE, "lapse_rate" default value is adiabatic lapse rate 9.8 K/km + oops::Parameter LapseRateValue{"lapse_rate", 9.8, this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to GSL, the size of vector "local_lapse_rate_levels" can only be 2. + oops::Parameter> LocalLapseRateLevels{"local_lapse_rate_levels", this}; }; // ----------------------------------------------------------------------------- diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 73c6d40..5ac4298 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -296,7 +296,12 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) deallocate(model_t_2000) !---------------- -case ("RRFS_GSL") +case ("CONSTANT_LAPSE_RATE") +!---------------- +! to be implemented + +!---------------- +case ("GSL") !---------------- ! to be implemented From cf9dd7f1e40dc7f60390e78fd1877b862b65e372 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Tue, 20 Aug 2024 22:31:35 +0000 Subject: [PATCH 02/13] Change scheme name --- src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index c161a0e..ba85e5a 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -22,7 +22,7 @@ namespace ufo { /// enum type for surface correction type, and ParameterTraitsHelper for it enum class SfcCorrectionType { - UKMO, WRFDA, RRFS_GSL, RRFS_GSL_LOCAL + UKMO, WRFDA, CONSTANT_LAPSE_RATE, GSL }; struct SfcCorrectionTypeParameterTraitsHelper { typedef SfcCorrectionType EnumType; From 4f1f7b4b9a4cb4e6fe18b66b49cf7aa7a8401427 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Mon, 26 Aug 2024 21:40:22 +0000 Subject: [PATCH 03/13] Fix to use required parameter for LocalLapseRateLevels --- src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index ba85e5a..aec7b8b 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -87,7 +87,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { oops::Parameter LapseRateValue{"lapse_rate", 9.8, this}; /// Note: Only relevant if \c SfcCorrectionType is set to GSL, the size of vector "local_lapse_rate_levels" can only be 2. - oops::Parameter> LocalLapseRateLevels{"local_lapse_rate_levels", this}; + oops::RequiredParameter> LocalLapseRateLevels{"local_lapse_rate_levels", this}; }; // ----------------------------------------------------------------------------- From f0dfa9f8815483c9e399416c244f53a373a5bb8d Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Mon, 26 Aug 2024 22:27:15 +0000 Subject: [PATCH 04/13] Implement constant lapse rate method for surface temperature adjustment --- .../sfccorrected/ObsSfcCorrectedParameters.h | 4 +- .../sfccorrected/ufo_sfccorrected_mod.F90 | 68 +++++++++++++++++-- 2 files changed, 66 insertions(+), 6 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index aec7b8b..f59859c 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -86,8 +86,8 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { /// Note: Only relevant if \c SfcCorrectionType is set to CONSTANT_LAPSE_RATE, "lapse_rate" default value is adiabatic lapse rate 9.8 K/km oops::Parameter LapseRateValue{"lapse_rate", 9.8, this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL, the size of vector "local_lapse_rate_levels" can only be 2. - oops::RequiredParameter> LocalLapseRateLevels{"local_lapse_rate_levels", this}; + /// Note: Only relevant if \c SfcCorrectionType is set to GSL. + oops::Parameter LocalLapseRateLevels{"local_lapse_rate_levels", 5, this}; }; // ----------------------------------------------------------------------------- diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 5ac4298..bfa5298 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -30,12 +30,14 @@ module ufo_sfccorrected_mod type(oops_variables), public :: geovars character(len=MAXVARLEN) :: da_sfc_scheme character(len=MAXVARLEN) :: station_altitude + real(kind_real) :: lapse_rate + integer, allocatable :: local_lapse_rate_levels contains procedure :: setup => ufo_sfccorrected_setup procedure :: simobs => ufo_sfccorrected_simobs end type ufo_sfccorrected - character(len=MAXVARLEN), dimension(5) :: geovars_list = (/ var_ps, var_geomz, var_sfc_geomz, var_ts, var_prs /) + character(len=MAXVARLEN), dimension(5) :: geovars_list = (/ var_ps, var_geomz, var_sfc_geomz, var_ts, var_prs, var_sfc_t2m /) contains @@ -48,6 +50,9 @@ subroutine ufo_sfccorrected_setup(self, f_conf) type(fckit_configuration), intent(in) :: f_conf character(len=:), allocatable :: str_sfc_scheme, str_var_sfc_geomz, str_var_geomz character(len=:), allocatable :: str_obs_height +real(kind_real), allocatable :: constant_lapse_rate +integer, allocatable :: local_lapse_rate_levels + character(max_string) :: debug_msg !> In case where a user wants to specify geoVaLs variable name of model @@ -72,6 +77,16 @@ subroutine ufo_sfccorrected_setup(self, f_conf) call f_conf%get_or_die("station_altitude", str_obs_height) self%station_altitude = str_obs_height +if (self%da_sfc_scheme.eq."CONSTANT_LAPSE_RATE") then + call f_conf%get_or_die("lapse_rate", constant_lapse_rate) + self%lapse_rate = constant_lapse_rate +endif + +if (self%da_sfc_scheme.eq."GSL") then + call f_conf%get_or_die("local_lapse_rate_levels", local_lapse_rate_levels) + self%local_lapse_rate_levels = local_lapse_rate_levels +endif + end subroutine ufo_sfccorrected_setup ! ------------------------------------------------------------------------------ @@ -92,7 +107,7 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) real(kind_real) :: H2000 = 2000.0 integer :: nobs, iobs, ivar, iobsvar, k, kbot, idx_geop real(kind_real), allocatable :: cor_tsfc(:) -type(ufo_geoval), pointer :: model_ps, model_p, model_sfc_geomz, model_t, model_geomz +type(ufo_geoval), pointer :: model_ps, model_p, model_sfc_geomz, model_t, model_geomz, model_t2m character(len=*), parameter :: myname_="ufo_sfccorrected_simobs" character(max_string) :: err_msg real(kind_real) :: wf @@ -130,13 +145,15 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) write(err_msg,'(a)') 'ufo_sfccorrected:'//new_line('a')// & 'retrieving GeoVaLs with names: '//trim(geovars_list(1))// & ', '//trim(geovars_list(2))//', '//trim(geovars_list(3))// & - ', '//trim(geovars_list(4))//', '//trim(geovars_list(5)) + ', '//trim(geovars_list(4))//', '//trim(geovars_list(5))// & + ', '//trim(geovars_list(6)) call fckit_log%debug(err_msg) call ufo_geovals_get_var(geovals, trim(geovars_list(1)), model_ps) call ufo_geovals_get_var(geovals, trim(geovars_list(2)), model_geomz) call ufo_geovals_get_var(geovals, trim(geovars_list(3)), model_sfc_geomz) call ufo_geovals_get_var(geovals, trim(geovars_list(4)), model_t) call ufo_geovals_get_var(geovals, trim(geovars_list(5)), model_p) +call ufo_geovals_get_var(geovals, trim(geovars_list(6)), model_t2m) ! discover if the model vertical profiles are ordered top-bottom or not kbot = 1 @@ -298,7 +315,10 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) !---------------- case ("CONSTANT_LAPSE_RATE") !---------------- -! to be implemented + + model_ts = model_t2m%vals(:) + + call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, self%lapse_rate) !---------------- case ("GSL") @@ -453,5 +473,45 @@ subroutine da_int_ukmo(nobs, missing, cor_tsfc, H_o, T_o, H_m, P_m, T_2000, P_20 end subroutine da_int_ukmo +! ----------------------------------------------------------- +!> \Conduct terrain height correction for sfc temp +!! +!! \Method: Lapse Rate (adiabatic by default) +!! +!! Where: +!! H_m = model sfc height +!! H_o = obs station height +!! T_o = obs temp at station height +!! T_lr = temperature lapse rate +!! T_o2m = obs temp interpolated from station height to model sfc level + +subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) +implicit none +integer, intent (in) :: nobs +real(c_double), intent (in) :: missing +real(kind_real), dimension(nobs), intent (in) :: H_o +real(kind_real), dimension(nobs), intent (in) :: H_m +real(kind_real), dimension(nobs), intent (in) :: T_o +real(kind_real), dimension(nobs), intent (in) :: T_lr +real(kind_real), dimension(nobs), intent (out) :: cor_tsfc +real(kind_real), dimension(nobs) :: T_o2m +integer i + +! T_o2m : obs temp interpolated to model sfc + +! Adjust temp from station height to model sfc based on lapse rate +! ------------------------------------------------- +where ( H_o /= missing .and. T_o /= missing ) + + T_o2m = T_o - T_lr * ( H_m - H_o) + +elsewhere + T_o2m = T_m ! to give zero analysis increment +end where + + cor_tsfc = T_o2m + +end subroutine da_int_lr + ! ------------------------------------------------------------------------------ end module ufo_sfccorrected_mod From b2031be74d0872aa9afcf78b2775bf964cc8a427 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Tue, 27 Aug 2024 05:23:49 +0000 Subject: [PATCH 05/13] Implement GSL local lapse rate method --- .../sfccorrected/ufo_sfccorrected_mod.F90 | 64 +++++++++++++++++-- 1 file changed, 57 insertions(+), 7 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index bfa5298..6fa1880 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -114,7 +114,7 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) integer :: wi logical :: variable_present, variable_present_t, variable_present_q real(kind_real), dimension(:), allocatable :: obs_height, obs_t, obs_q, obs_psfc, obs_lat, obs_tv -real(kind_real), dimension(:), allocatable :: model_ts, model_zs, model_level1, model_p_2000, model_t_2000, model_psfc +real(kind_real), dimension(:), allocatable :: model_ts, model_zs, model_level1, model_p_2000, model_t_2000, model_psfc, llr real(kind_real), dimension(:), allocatable :: H2000_geop real(kind_real), dimension(:), allocatable :: avg_tv real(kind_real) :: model_znew @@ -318,12 +318,20 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) model_ts = model_t2m%vals(:) - call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, self%lapse_rate) + call da_int_clr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, self%lapse_rate) !---------------- case ("GSL") !---------------- -! to be implemented + + allocate(llr(nobs)) + + llr = (model_t%vals(self%local_lapse_rate_levels,:) - model_t%vals(1,:)) / & + (model_geomz%vals(self%local_lapse_rate_levels,:) - model_geomz%vals(1,:)) + + model_ts = model_t2m%vals(:) + + call da_int_llr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, llr) !----------- case default @@ -476,7 +484,7 @@ end subroutine da_int_ukmo ! ----------------------------------------------------------- !> \Conduct terrain height correction for sfc temp !! -!! \Method: Lapse Rate (adiabatic by default) +!! \Method: Constant Lapse Rate (adiabatic by default) !! !! Where: !! H_m = model sfc height @@ -485,7 +493,48 @@ end subroutine da_int_ukmo !! T_lr = temperature lapse rate !! T_o2m = obs temp interpolated from station height to model sfc level -subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) +subroutine da_int_clr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) +implicit none +integer, intent (in) :: nobs +real(c_double), intent (in) :: missing +real(kind_real), dimension(nobs), intent (in) :: H_o +real(kind_real), dimension(nobs), intent (in) :: H_m +real(kind_real), dimension(nobs), intent (in) :: T_o +real(kind_real), intent (in) :: T_lr +real(kind_real), dimension(nobs), intent (out) :: cor_tsfc +real(kind_real), dimension(nobs) :: T_o2m +integer i + +! T_o2m : obs temp interpolated to model sfc + +! Adjust temp from station height to model sfc based on lapse rate +! ------------------------------------------------- +where ( H_o /= missing .and. T_o /= missing ) + + T_o2m = T_o - T_lr * ( H_m - H_o) + +elsewhere + !T_o2m = T_m ! to give zero analysis increment + T_o2m = T_o +end where + + cor_tsfc = T_o2m + +end subroutine da_int_clr + +! ----------------------------------------------------------- +!> \Conduct terrain height correction for sfc temp +!! +!! \Method: Constant Lapse Rate (adiabatic by default) +!! +!! Where: +!! H_m = model sfc height +!! H_o = obs station height +!! T_o = obs temp at station height +!! T_lr = temperature lapse rate +!! T_o2m = obs temp interpolated from station height to model sfc level + +subroutine da_int_llr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) implicit none integer, intent (in) :: nobs real(c_double), intent (in) :: missing @@ -506,12 +555,13 @@ subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) T_o2m = T_o - T_lr * ( H_m - H_o) elsewhere - T_o2m = T_m ! to give zero analysis increment + !T_o2m = T_m ! to give zero analysis increment + T_o2m = T_o end where cor_tsfc = T_o2m -end subroutine da_int_lr +end subroutine da_int_llr ! ------------------------------------------------------------------------------ end module ufo_sfccorrected_mod From a79709957cc931d038dde4c6cd7f602a486c3100 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Tue, 27 Aug 2024 05:29:38 +0000 Subject: [PATCH 06/13] Bug fix for allocation --- src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 6fa1880..dc5b031 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -31,13 +31,13 @@ module ufo_sfccorrected_mod character(len=MAXVARLEN) :: da_sfc_scheme character(len=MAXVARLEN) :: station_altitude real(kind_real) :: lapse_rate - integer, allocatable :: local_lapse_rate_levels + integer :: local_lapse_rate_levels contains procedure :: setup => ufo_sfccorrected_setup procedure :: simobs => ufo_sfccorrected_simobs end type ufo_sfccorrected - character(len=MAXVARLEN), dimension(5) :: geovars_list = (/ var_ps, var_geomz, var_sfc_geomz, var_ts, var_prs, var_sfc_t2m /) + character(len=MAXVARLEN), dimension(6) :: geovars_list = (/ var_ps, var_geomz, var_sfc_geomz, var_ts, var_prs, var_sfc_t2m /) contains From 846e6140312befe0665cfeb769ec0e19a94ff952 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Tue, 27 Aug 2024 21:31:02 +0000 Subject: [PATCH 07/13] Combing Constant and Local Lapse Rate options --- .../sfccorrected/ObsSfcCorrectedParameters.h | 12 +-- .../sfccorrected/ufo_sfccorrected_mod.F90 | 89 ++++++------------- 2 files changed, 35 insertions(+), 66 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index f59859c..3f74cb7 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -22,7 +22,7 @@ namespace ufo { /// enum type for surface correction type, and ParameterTraitsHelper for it enum class SfcCorrectionType { - UKMO, WRFDA, CONSTANT_LAPSE_RATE, GSL + UKMO, WRFDA, GSL }; struct SfcCorrectionTypeParameterTraitsHelper { typedef SfcCorrectionType EnumType; @@ -30,7 +30,6 @@ struct SfcCorrectionTypeParameterTraitsHelper { static constexpr util::NamedEnumerator namedValues[] = { { SfcCorrectionType::UKMO, "UKMO" }, { SfcCorrectionType::WRFDA, "WRFDA" }, - { SfcCorrectionType::CONSTANT_LAPSE_RATE, "CONSTANT_LAPSE_RATE" }, { SfcCorrectionType::GSL, "GSL" } }; }; @@ -65,7 +64,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { this}; oops::Parameter correctionType{"da_sfc_scheme", - "Scheme used for surface temperature correction (UKMO or WRFDA)", + "Scheme used for surface temperature correction (UKMO, WRFDA or GSL)", SfcCorrectionType::UKMO, this}; /// Note: "height" default value has to be consistent with var_geomz defined @@ -82,11 +81,14 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { /// Note: "station_altitude" default value is "stationElevation" oops::Parameter ObsHeightName{"station_altitude", "stationElevation", this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to GSL, "lapse_rate_option" default value is "LOCAL" + oops::Parameter LapseRateOption{"lapse_rate_option", "Lapse rate option for surface temperature correction (CONSTANT or LOCAL)", "LOCAL", this}; - /// Note: Only relevant if \c SfcCorrectionType is set to CONSTANT_LAPSE_RATE, "lapse_rate" default value is adiabatic lapse rate 9.8 K/km + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "CONSTANT", "lapse_rate" default value is adiabatic lapse rate 9.8 K/km oops::Parameter LapseRateValue{"lapse_rate", 9.8, this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL. + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "LOCAL" oops::Parameter LocalLapseRateLevels{"local_lapse_rate_levels", 5, this}; }; diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index dc5b031..8b2f816 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -30,6 +30,7 @@ module ufo_sfccorrected_mod type(oops_variables), public :: geovars character(len=MAXVARLEN) :: da_sfc_scheme character(len=MAXVARLEN) :: station_altitude + character(len=MAXVARLEN) :: lapse_rate_option real(kind_real) :: lapse_rate integer :: local_lapse_rate_levels contains @@ -48,7 +49,7 @@ subroutine ufo_sfccorrected_setup(self, f_conf) implicit none class(ufo_sfccorrected), intent(inout) :: self type(fckit_configuration), intent(in) :: f_conf -character(len=:), allocatable :: str_sfc_scheme, str_var_sfc_geomz, str_var_geomz +character(len=:), allocatable :: str_sfc_scheme, str_var_sfc_geomz, str_var_geomz, str_lapse_rate_option character(len=:), allocatable :: str_obs_height real(kind_real), allocatable :: constant_lapse_rate integer, allocatable :: local_lapse_rate_levels @@ -77,14 +78,21 @@ subroutine ufo_sfccorrected_setup(self, f_conf) call f_conf%get_or_die("station_altitude", str_obs_height) self%station_altitude = str_obs_height -if (self%da_sfc_scheme.eq."CONSTANT_LAPSE_RATE") then - call f_conf%get_or_die("lapse_rate", constant_lapse_rate) - self%lapse_rate = constant_lapse_rate -endif - if (self%da_sfc_scheme.eq."GSL") then - call f_conf%get_or_die("local_lapse_rate_levels", local_lapse_rate_levels) - self%local_lapse_rate_levels = local_lapse_rate_levels + call f_conf%get_or_die("lapse_rate_option", str_lapse_rate_option) + self%lapse_rate_option = str_lapse_rate_option + select case (trim(self%lapse_rate_option)) + case ("CONSTANT") + call f_conf%get_or_die("lapse_rate", constant_lapse_rate) + self%lapse_rate = constant_lapse_rate + case ("LOCAL") + call f_conf%get_or_die("local_lapse_rate_levels", local_lapse_rate_levels) + self%local_lapse_rate_levels = local_lapse_rate_levels + case default + write(err_msg,*) "ufo_sfccorrected: lapse_rate_option not recognized" + call fckit_log%debug(err_msg) + call abor1_ftn(err_msg) + end select endif end subroutine ufo_sfccorrected_setup @@ -313,25 +321,25 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) deallocate(model_t_2000) !---------------- -case ("CONSTANT_LAPSE_RATE") +case ("GSL") !---------------- model_ts = model_t2m%vals(:) - call da_int_clr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, self%lapse_rate) + allocate(lr(nobs)) -!---------------- -case ("GSL") -!---------------- + if (self%lapse_rate_option.eq."LOCAL") then - allocate(llr(nobs)) + lr = (model_t%vals(self%local_lapse_rate_levels,:) - model_t%vals(1,:)) / & + (model_geomz%vals(self%local_lapse_rate_levels,:) - model_geomz%vals(1,:)) - llr = (model_t%vals(self%local_lapse_rate_levels,:) - model_t%vals(1,:)) / & - (model_geomz%vals(self%local_lapse_rate_levels,:) - model_geomz%vals(1,:)) + else: - model_ts = model_t2m%vals(:) + lr = self%lapse_rate + + endif - call da_int_llr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, llr) + call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, lr) !----------- case default @@ -493,48 +501,7 @@ end subroutine da_int_ukmo !! T_lr = temperature lapse rate !! T_o2m = obs temp interpolated from station height to model sfc level -subroutine da_int_clr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) -implicit none -integer, intent (in) :: nobs -real(c_double), intent (in) :: missing -real(kind_real), dimension(nobs), intent (in) :: H_o -real(kind_real), dimension(nobs), intent (in) :: H_m -real(kind_real), dimension(nobs), intent (in) :: T_o -real(kind_real), intent (in) :: T_lr -real(kind_real), dimension(nobs), intent (out) :: cor_tsfc -real(kind_real), dimension(nobs) :: T_o2m -integer i - -! T_o2m : obs temp interpolated to model sfc - -! Adjust temp from station height to model sfc based on lapse rate -! ------------------------------------------------- -where ( H_o /= missing .and. T_o /= missing ) - - T_o2m = T_o - T_lr * ( H_m - H_o) - -elsewhere - !T_o2m = T_m ! to give zero analysis increment - T_o2m = T_o -end where - - cor_tsfc = T_o2m - -end subroutine da_int_clr - -! ----------------------------------------------------------- -!> \Conduct terrain height correction for sfc temp -!! -!! \Method: Constant Lapse Rate (adiabatic by default) -!! -!! Where: -!! H_m = model sfc height -!! H_o = obs station height -!! T_o = obs temp at station height -!! T_lr = temperature lapse rate -!! T_o2m = obs temp interpolated from station height to model sfc level - -subroutine da_int_llr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) +subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) implicit none integer, intent (in) :: nobs real(c_double), intent (in) :: missing @@ -561,7 +528,7 @@ subroutine da_int_llr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) cor_tsfc = T_o2m -end subroutine da_int_llr +end subroutine da_int_lr ! ------------------------------------------------------------------------------ end module ufo_sfccorrected_mod From 5c01b78b77bf3f876eabb11b574134432a0db5cb Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Tue, 27 Aug 2024 22:10:15 +0000 Subject: [PATCH 08/13] Variable name change based on Hui's suggestion --- .../sfccorrected/ObsSfcCorrectedParameters.h | 2 +- .../operators/sfccorrected/ufo_sfccorrected_mod.F90 | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index 3f74cb7..dd240e5 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -89,7 +89,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { oops::Parameter LapseRateValue{"lapse_rate", 9.8, this}; /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "LOCAL" - oops::Parameter LocalLapseRateLevels{"local_lapse_rate_levels", 5, this}; + oops::Parameter LocalLapseRateLevels{"local_lapse_rate_level", 5, this}; }; // ----------------------------------------------------------------------------- diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 8b2f816..9012ff1 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -32,7 +32,7 @@ module ufo_sfccorrected_mod character(len=MAXVARLEN) :: station_altitude character(len=MAXVARLEN) :: lapse_rate_option real(kind_real) :: lapse_rate - integer :: local_lapse_rate_levels + integer :: local_lapse_rate_level contains procedure :: setup => ufo_sfccorrected_setup procedure :: simobs => ufo_sfccorrected_simobs @@ -52,7 +52,7 @@ subroutine ufo_sfccorrected_setup(self, f_conf) character(len=:), allocatable :: str_sfc_scheme, str_var_sfc_geomz, str_var_geomz, str_lapse_rate_option character(len=:), allocatable :: str_obs_height real(kind_real), allocatable :: constant_lapse_rate -integer, allocatable :: local_lapse_rate_levels +integer, allocatable :: local_lapse_rate_level character(max_string) :: debug_msg @@ -86,8 +86,8 @@ subroutine ufo_sfccorrected_setup(self, f_conf) call f_conf%get_or_die("lapse_rate", constant_lapse_rate) self%lapse_rate = constant_lapse_rate case ("LOCAL") - call f_conf%get_or_die("local_lapse_rate_levels", local_lapse_rate_levels) - self%local_lapse_rate_levels = local_lapse_rate_levels + call f_conf%get_or_die("local_lapse_rate_level", local_lapse_rate_level) + self%local_lapse_rate_level = local_lapse_rate_level case default write(err_msg,*) "ufo_sfccorrected: lapse_rate_option not recognized" call fckit_log%debug(err_msg) @@ -330,8 +330,8 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) if (self%lapse_rate_option.eq."LOCAL") then - lr = (model_t%vals(self%local_lapse_rate_levels,:) - model_t%vals(1,:)) / & - (model_geomz%vals(self%local_lapse_rate_levels,:) - model_geomz%vals(1,:)) + lr = (model_t%vals(self%local_lapse_rate_level,:) - model_t%vals(1,:)) / & + (model_geomz%vals(self%local_lapse_rate_level,:) - model_geomz%vals(1,:)) else: From 1dd19f995890ab1c7a1ca54d2982480395da8dfc Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Thu, 12 Dec 2024 01:47:28 -0600 Subject: [PATCH 09/13] Variables rename and bug fixes --- .../operators/sfccorrected/ObsSfcCorrected.h | 4 +-- .../sfccorrected/ObsSfcCorrected.interface.h | 2 +- .../sfccorrected/ObsSfcCorrectedParameters.h | 30 +++++++++++++------ .../sfccorrected/ufo_sfccorrected_mod.F90 | 24 +++++++-------- 4 files changed, 36 insertions(+), 24 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrected.h b/src/ufo/operators/sfccorrected/ObsSfcCorrected.h index 317e1d2..8e90442 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrected.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrected.h @@ -48,7 +48,7 @@ class ObsSfcCorrected : public ObsOperatorBase, // Other const oops::Variables & requiredVars() const override {return varin_;} - oops::Variables simulatedVars() const override {return operatorVars_;} + oops::ObsVariables simulatedVars() const override {return operatorVars_;} int & toFortran() {return keyOper_;} const int & toFortran() const {return keyOper_;} @@ -58,7 +58,7 @@ class ObsSfcCorrected : public ObsOperatorBase, F90hop keyOper_; const ioda::ObsSpace& odb_; oops::Variables varin_; - oops::Variables operatorVars_; + oops::ObsVariables operatorVars_; }; // ----------------------------------------------------------------------------- diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.h b/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.h index e11127c..98552b2 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.h @@ -50,7 +50,7 @@ extern "C" { // ----------------------------------------------------------------------------- void ufo_sfccorrected_setup_f90(F90hop &, const eckit::Configuration &, - const oops::Variables&operatorVars, + const oops::ObsVariables&operatorVars, const int *operatorVarIndices, const int numOperatorVarIndices, oops::Variables &requiredVars); void ufo_sfccorrected_delete_f90(F90hop &); diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index dd240e5..e215df6 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -71,25 +71,37 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { /// in ufo_variables_mod.F90 oops::Parameter geovarGeomZ{"geovar_geomz", "Model variable for height of vertical levels", - "height", this}; + "height_above_mean_sea_level", this}; /// Note: "surface_altitude" default value has to be consistent with var_sfc_geomz /// in ufo_variables_mod.F90 oops::Parameter geovarSfcGeomZ{"geovar_sfc_geomz", "Model variable for surface height", - "surface_altitude", this}; + "height_above_mean_sea_level_at_surface", this}; /// Note: "station_altitude" default value is "stationElevation" oops::Parameter ObsHeightName{"station_altitude", "stationElevation", this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL, "lapse_rate_option" default value is "LOCAL" - oops::Parameter LapseRateOption{"lapse_rate_option", "Lapse rate option for surface temperature correction (CONSTANT or LOCAL)", "LOCAL", this}; - - /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "CONSTANT", "lapse_rate" default value is adiabatic lapse rate 9.8 K/km - oops::Parameter LapseRateValue{"lapse_rate", 9.8, this}; + /// Note: Only relevant if \c SfcCorrectionType is set to GSL, "lapse_rate_option" default value is "Local" + oops::Parameter LapseRateOption{"lapse_rate_option", "Lapse rate option for surface temperature correction (Constant or Local)", "Local", this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "Constant", "lapse_rate" default value is adiabatic lapse rate 9.8 K/km + oops::Parameter LapseRateValue + {"lapse_rate", + "The lapse rate used to adjust the observed surface temperature to " + "the model's surface level. Used if lapse rate option is set to constant, " + "otherwise ignored.", + 9.8, + this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "LOCAL" - oops::Parameter LocalLapseRateLevels{"local_lapse_rate_level", 5, this}; + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "Local" + oops::Parameter LocalLapseRateLevel + {"local_lapse_rate_level", + "The highest model level used to calculate the local lapse rate, " + "which adjusts the observed surface temperature to the model's surface level. " + "Used if lapse rate option is set to local, otherwise ignored.", + 5, + this}; }; // ----------------------------------------------------------------------------- diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 9012ff1..441295b 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -51,8 +51,8 @@ subroutine ufo_sfccorrected_setup(self, f_conf) type(fckit_configuration), intent(in) :: f_conf character(len=:), allocatable :: str_sfc_scheme, str_var_sfc_geomz, str_var_geomz, str_lapse_rate_option character(len=:), allocatable :: str_obs_height -real(kind_real), allocatable :: constant_lapse_rate -integer, allocatable :: local_lapse_rate_level +real(kind_real) :: constant_lapse_rate +integer :: local_lapse_rate_level character(max_string) :: debug_msg @@ -82,16 +82,16 @@ subroutine ufo_sfccorrected_setup(self, f_conf) call f_conf%get_or_die("lapse_rate_option", str_lapse_rate_option) self%lapse_rate_option = str_lapse_rate_option select case (trim(self%lapse_rate_option)) - case ("CONSTANT") + case ("Constant") call f_conf%get_or_die("lapse_rate", constant_lapse_rate) self%lapse_rate = constant_lapse_rate - case ("LOCAL") + case ("Local") call f_conf%get_or_die("local_lapse_rate_level", local_lapse_rate_level) self%local_lapse_rate_level = local_lapse_rate_level case default - write(err_msg,*) "ufo_sfccorrected: lapse_rate_option not recognized" - call fckit_log%debug(err_msg) - call abor1_ftn(err_msg) + write(debug_msg,*) "ufo_sfccorrected: lapse_rate_option not recognized" + call fckit_log%debug(debug_msg) + call abor1_ftn(debug_msg) end select endif @@ -122,7 +122,7 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) integer :: wi logical :: variable_present, variable_present_t, variable_present_q real(kind_real), dimension(:), allocatable :: obs_height, obs_t, obs_q, obs_psfc, obs_lat, obs_tv -real(kind_real), dimension(:), allocatable :: model_ts, model_zs, model_level1, model_p_2000, model_t_2000, model_psfc, llr +real(kind_real), dimension(:), allocatable :: model_ts, model_zs, model_level1, model_p_2000, model_t_2000, model_psfc, lr real(kind_real), dimension(:), allocatable :: H2000_geop real(kind_real), dimension(:), allocatable :: avg_tv real(kind_real) :: model_znew @@ -168,7 +168,7 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) do iobs = 1, nlocs if (model_geomz%vals(1,iobs) .ne. missing) then if (model_geomz%vals(1,iobs) .gt. model_geomz%vals(model_geomz%nval,iobs)) then - write(err_msg,'(a)') ' ufo_sfcpcorrected:'//new_line('a')// & + write(err_msg,'(a)') ' ufo_sfccorrected:'//new_line('a')// & ' Model vertical height profile is from top to bottom' call fckit_log%debug(err_msg) kbot = model_geomz%nval @@ -324,16 +324,16 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) case ("GSL") !---------------- - model_ts = model_t2m%vals(:) + model_ts = model_t2m%vals(1,:) allocate(lr(nobs)) if (self%lapse_rate_option.eq."LOCAL") then - lr = (model_t%vals(self%local_lapse_rate_level,:) - model_t%vals(1,:)) / & + lr = -1. * (model_t%vals(self%local_lapse_rate_level,:) - model_t%vals(1,:)) / & (model_geomz%vals(self%local_lapse_rate_level,:) - model_geomz%vals(1,:)) - else: + else lr = self%lapse_rate From 2ff27ef284ffb1c5b0de85ffcfa131eca9281494 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Fri, 13 Dec 2024 13:05:45 -0600 Subject: [PATCH 10/13] Bug fixes for inconsistant type for obsvars, and adjusted temperature calculation --- .../ObsSfcCorrected.interface.F90 | 3 ++- .../sfccorrected/ufo_sfccorrected_mod.F90 | 26 ++++++++++++------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.F90 b/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.F90 index 239918c..a768cee 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.F90 +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrected.interface.F90 @@ -36,6 +36,7 @@ subroutine ufo_sfccorrected_setup_c(c_key_self, c_conf, c_obsvars, c_obsvarindic c_geovars) bind(c,name='ufo_sfccorrected_setup_f90') use fckit_configuration_module, only: fckit_configuration use oops_variables_mod +use obs_variables_mod implicit none integer(c_int), intent(inout) :: c_key_self type(c_ptr), value, intent(in) :: c_conf @@ -50,7 +51,7 @@ subroutine ufo_sfccorrected_setup_c(c_key_self, c_conf, c_obsvars, c_obsvarindic call ufo_sfccorrected_registry%setup(c_key_self, self) f_conf = fckit_configuration(c_conf) -self%obsvars = oops_variables(c_obsvars) +self%obsvars = obs_variables(c_obsvars) allocate(self%obsvarindices(self%obsvars%nvars())) self%obsvarindices(:) = c_obsvarindices(:) + 1 ! Convert from C to Fortran indexing self%geovars = oops_variables(c_geovars) diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 441295b..46ac61b 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -9,6 +9,7 @@ module ufo_sfccorrected_mod use oops_variables_mod + use obs_variables_mod use ufo_vars_mod use missing_values_mod use iso_c_binding @@ -23,7 +24,7 @@ module ufo_sfccorrected_mod ! Fortran derived type for observation type type, public :: ufo_sfccorrected private - type(oops_variables), public :: obsvars ! Variables to be simulated + type(obs_variables), public :: obsvars ! Variables to be simulated integer, allocatable, public :: obsvarindices(:) ! Indices of obsvars in the list of all ! simulated variables in the ObsSpace. ! allocated/deallocated at interface layer @@ -339,7 +340,10 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) endif - call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, lr) + ! Convert the unit of the lapse rate from K/km to K/m + lr = lr * 0.001 + + call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, model_ts, lr) !----------- case default @@ -359,12 +363,12 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) ! cor_tsfc is corrected obs temp at model sfc ! do iobs = 1, nlocs - if (cor_tsfc(iobs) /= missing) then + !if (cor_tsfc(iobs) /= missing) then ! for T_o2m, adjusted hofx - O = model_ts - T_o2m, OBS not adjusted. - hofx(ivar,iobs) = model_ts(iobs) + obs_t(iobs) - cor_tsfc(iobs) - else + ! hofx(ivar,iobs) = model_ts(iobs) + obs_t(iobs) - cor_tsfc(iobs) + !else hofx(ivar,iobs) = model_ts(iobs) - end if + !end if enddo enddo @@ -501,16 +505,17 @@ end subroutine da_int_ukmo !! T_lr = temperature lapse rate !! T_o2m = obs temp interpolated from station height to model sfc level -subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) +subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_m, T_lr) implicit none integer, intent (in) :: nobs real(c_double), intent (in) :: missing real(kind_real), dimension(nobs), intent (in) :: H_o real(kind_real), dimension(nobs), intent (in) :: H_m real(kind_real), dimension(nobs), intent (in) :: T_o +real(kind_real), dimension(nobs), intent (in) :: T_m real(kind_real), dimension(nobs), intent (in) :: T_lr real(kind_real), dimension(nobs), intent (out) :: cor_tsfc -real(kind_real), dimension(nobs) :: T_o2m +real(kind_real), dimension(nobs) :: T_o2m, T_m2o integer i ! T_o2m : obs temp interpolated to model sfc @@ -520,10 +525,11 @@ subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) where ( H_o /= missing .and. T_o /= missing ) T_o2m = T_o - T_lr * ( H_m - H_o) + T_m2o = T_m + T_lr * ( H_m - H_o) elsewhere - !T_o2m = T_m ! to give zero analysis increment - T_o2m = T_o + T_o2m = T_m + T_m2o = T_o end where cor_tsfc = T_o2m From 2da845a2ccd2cf5b677e5cf9c4a8a007c86ef7e9 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Mon, 16 Dec 2024 16:53:09 -0600 Subject: [PATCH 11/13] Bug fix to properly account for FV3 inverse vertical coordinate indices; added upper and lower limits for local lapse rate; included additional debug information. --- .../sfccorrected/ObsSfcCorrectedParameters.h | 20 +++- .../sfccorrected/ufo_sfccorrected_mod.F90 | 98 +++++++++++++++---- 2 files changed, 97 insertions(+), 21 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index e215df6..1fcd8f4 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -88,7 +88,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "Constant", "lapse_rate" default value is adiabatic lapse rate 9.8 K/km oops::Parameter LapseRateValue {"lapse_rate", - "The lapse rate used to adjust the observed surface temperature to " + "The lapse rate (K/km) used to adjust the observed surface temperature to " "the model's surface level. Used if lapse rate option is set to constant, " "otherwise ignored.", 9.8, @@ -102,6 +102,24 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { "Used if lapse rate option is set to local, otherwise ignored.", 5, this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "Local" + oops::Parameter MinThreshold + {"min_threshold", + "The minimum lapse rate (K/km) can be applied to adjust the " + "observed surface temperature to the model's surface level. " + "Used if lapse rate option is set to local, otherwise ignored.", + 0.5, + this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "Local" + oops::Parameter MaxThreshold + {"max_threshold", + "The maximum lapse rate (K/km) can be applied to adjust the " + "observed surface temperature to the model's surface level. " + "Used if lapse rate option is set to local, otherwise ignored.", + 10.0, + this}; }; // ----------------------------------------------------------------------------- diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 46ac61b..4622e36 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -34,6 +34,7 @@ module ufo_sfccorrected_mod character(len=MAXVARLEN) :: lapse_rate_option real(kind_real) :: lapse_rate integer :: local_lapse_rate_level + real(kind_real) :: min_threshold, max_threshold contains procedure :: setup => ufo_sfccorrected_setup procedure :: simobs => ufo_sfccorrected_simobs @@ -54,6 +55,7 @@ subroutine ufo_sfccorrected_setup(self, f_conf) character(len=:), allocatable :: str_obs_height real(kind_real) :: constant_lapse_rate integer :: local_lapse_rate_level +real(kind_real) :: min_threshold, max_threshold character(max_string) :: debug_msg @@ -85,10 +87,14 @@ subroutine ufo_sfccorrected_setup(self, f_conf) select case (trim(self%lapse_rate_option)) case ("Constant") call f_conf%get_or_die("lapse_rate", constant_lapse_rate) - self%lapse_rate = constant_lapse_rate + self%lapse_rate = constant_lapse_rate * 0.001 case ("Local") call f_conf%get_or_die("local_lapse_rate_level", local_lapse_rate_level) self%local_lapse_rate_level = local_lapse_rate_level + call f_conf%get_or_die("min_threshold", min_threshold) + self%min_threshold = min_threshold * 0.001 + call f_conf%get_or_die("max_threshold", max_threshold) + self%max_threshold = max_threshold * 0.001 case default write(debug_msg,*) "ufo_sfccorrected: lapse_rate_option not recognized" call fckit_log%debug(debug_msg) @@ -114,7 +120,7 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) ! Local variables real(c_double) :: missing real(kind_real) :: H2000 = 2000.0 -integer :: nobs, iobs, ivar, iobsvar, k, kbot, idx_geop +integer :: nobs, iobs, ivar, iobsvar, k, kbot, ktop_lr, idx_geop real(kind_real), allocatable :: cor_tsfc(:) type(ufo_geoval), pointer :: model_ps, model_p, model_sfc_geomz, model_t, model_geomz, model_t2m character(len=*), parameter :: myname_="ufo_sfccorrected_simobs" @@ -128,6 +134,8 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) real(kind_real), dimension(:), allocatable :: avg_tv real(kind_real) :: model_znew +real(kind_real) :: avg_lr, avg_temp_high, avg_temp_low, avg_height_high, avg_height_low, n_valid_obs, n_obs_within_thresh + missing = missing_value(missing) nobs = obsspace_get_nlocs(obss) @@ -329,21 +337,74 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) allocate(lr(nobs)) - if (self%lapse_rate_option.eq."LOCAL") then + if (self%lapse_rate_option.eq."Local") then - lr = -1. * (model_t%vals(self%local_lapse_rate_level,:) - model_t%vals(1,:)) / & - (model_geomz%vals(self%local_lapse_rate_level,:) - model_geomz%vals(1,:)) + if (kbot == 1) then + ktop_lr = self%local_lapse_rate_level + else + ktop_lr = kbot - self%local_lapse_rate_level + 1 + end if + + ! Local lapse rate in K/m + lr = -1. * (model_t%vals(ktop_lr,:) - model_t%vals(kbot,:)) / & + (model_geomz%vals(ktop_lr,:) - model_geomz%vals(kbot,:)) + + lr = min(self%max_threshold, max(self%min_threshold, lr)) + + n_valid_obs = 0. + n_obs_within_thresh = 0. + avg_lr = 0. + avg_temp_high = 0. + avg_temp_low = 0. + avg_height_high = 0. + avg_height_low = 0. + do iobs=1, nobs + if (model_ts(iobs) > 0) then + n_valid_obs = n_valid_obs + 1 + avg_lr = avg_lr + lr(iobs) + avg_temp_high = avg_temp_high + model_t%vals(ktop_lr,iobs) + avg_temp_low = avg_temp_low + model_t%vals(kbot,iobs) + avg_height_high = avg_height_high + model_geomz%vals(ktop_lr,iobs) + avg_height_low = avg_height_low + model_geomz%vals(kbot,iobs) + if (lr(iobs) > 0.0005 .and. lr(iobs) < 0.0100) then + n_obs_within_thresh = n_obs_within_thresh + 1 + end if + end if + end do + avg_lr = avg_lr * 1000. / n_valid_obs + avg_temp_high = avg_temp_high / n_valid_obs + avg_temp_low = avg_temp_low / n_valid_obs + avg_height_high = avg_height_high / n_valid_obs + avg_height_low = avg_height_low / n_valid_obs + + write(err_msg, '(2a, 5(a,f8.3),a)') & + "Lapse rate option: ", trim(self%lapse_rate_option), & + ", average Local Lapse Rate: ", avg_lr, & + " K/km, average temperature (higher level): ", avg_temp_high, & + " K, average temperature (lower level): ", avg_temp_low, & + " K, average AMSL (higher level): ", avg_height_high, & + " m, average AMSL (lower level): ", avg_height_low, " m." + call fckit_log%info(err_msg) + write(err_msg, '(3(a, f5.0))') & + "Number of METAR observations: ", n_valid_obs, & + ", number of observations within the thresholds: ", n_obs_within_thresh, & + ", number of observations outside the thresholds: ", n_valid_obs-n_obs_within_thresh + call fckit_log%info(err_msg) else + ! Local lapse rate in K/m lr = self%lapse_rate + avg_lr = self%lapse_rate * 1000. - endif + write(err_msg, '(3a, f8.3, a)') & + "Lapse rate option: ", trim(self%lapse_rate_option), & + "Constant Lapse Rate: ", avg_lr, " K/km" + call fckit_log%info(err_msg) - ! Convert the unit of the lapse rate from K/km to K/m - lr = lr * 0.001 + endif - call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, model_ts, lr) + call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, lr) !----------- case default @@ -363,12 +424,12 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) ! cor_tsfc is corrected obs temp at model sfc ! do iobs = 1, nlocs - !if (cor_tsfc(iobs) /= missing) then -! for T_o2m, adjusted hofx - O = model_ts - T_o2m, OBS not adjusted. - ! hofx(ivar,iobs) = model_ts(iobs) + obs_t(iobs) - cor_tsfc(iobs) - !else - hofx(ivar,iobs) = model_ts(iobs) - !end if + if (cor_tsfc(iobs) /= missing) then + ! for T_o2m, adjusted hofx - O = model_ts - T_o2m, OBS not adjusted. + hofx(ivar,iobs) = model_ts(iobs) + obs_t(iobs) - cor_tsfc(iobs) + else + hofx(ivar,iobs) = model_ts(iobs) + end if enddo enddo @@ -505,14 +566,13 @@ end subroutine da_int_ukmo !! T_lr = temperature lapse rate !! T_o2m = obs temp interpolated from station height to model sfc level -subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_m, T_lr) +subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_lr) implicit none integer, intent (in) :: nobs real(c_double), intent (in) :: missing real(kind_real), dimension(nobs), intent (in) :: H_o real(kind_real), dimension(nobs), intent (in) :: H_m real(kind_real), dimension(nobs), intent (in) :: T_o -real(kind_real), dimension(nobs), intent (in) :: T_m real(kind_real), dimension(nobs), intent (in) :: T_lr real(kind_real), dimension(nobs), intent (out) :: cor_tsfc real(kind_real), dimension(nobs) :: T_o2m, T_m2o @@ -525,11 +585,9 @@ subroutine da_int_lr(nobs, missing, cor_tsfc, H_o, T_o, H_m, T_m, T_lr) where ( H_o /= missing .and. T_o /= missing ) T_o2m = T_o - T_lr * ( H_m - H_o) - T_m2o = T_m + T_lr * ( H_m - H_o) elsewhere - T_o2m = T_m - T_m2o = T_o + T_o2m = T_o end where cor_tsfc = T_o2m From 7eb1d63d76a56656ffd499f7fecb1bb14d514d41 Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Mon, 6 Jan 2025 11:51:36 -0600 Subject: [PATCH 12/13] Add option to disable terrain adjustment for T, and include settings to limit the maximum/minimum lapse rate. --- .../sfccorrected/ObsSfcCorrectedParameters.h | 12 +++-- .../sfccorrected/ufo_sfccorrected_mod.F90 | 52 +++++++++++++------ 2 files changed, 44 insertions(+), 20 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h index 1fcd8f4..7eceaf8 100644 --- a/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h +++ b/src/ufo/operators/sfccorrected/ObsSfcCorrectedParameters.h @@ -83,7 +83,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { oops::Parameter ObsHeightName{"station_altitude", "stationElevation", this}; /// Note: Only relevant if \c SfcCorrectionType is set to GSL, "lapse_rate_option" default value is "Local" - oops::Parameter LapseRateOption{"lapse_rate_option", "Lapse rate option for surface temperature correction (Constant or Local)", "Local", this}; + oops::Parameter LapseRateOption{"lapse_rate_option", "Lapse rate option for surface temperature correction (Constant, Local or NoAdjustment)", "Local", this}; /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "Constant", "lapse_rate" default value is adiabatic lapse rate 9.8 K/km oops::Parameter LapseRateValue @@ -94,7 +94,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { 9.8, this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "Local" + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "Local" oops::Parameter LocalLapseRateLevel {"local_lapse_rate_level", "The highest model level used to calculate the local lapse rate, " @@ -103,7 +103,11 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { 5, this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "Local" + /// Should the local lapse rate be restricted to a specific range + /// Note: Only relevant if \c SfcCorrectionType is set to GSL and \c LapseRateOption is set to "Local" + oops::Parameter Threshold{"threshold", true, this}; + + /// Note: Only relevant if \c SfcCorrectionType is set to GSL, \c LapseRateOption is set to "Local", and \c Threshold is set to true oops::Parameter MinThreshold {"min_threshold", "The minimum lapse rate (K/km) can be applied to adjust the " @@ -112,7 +116,7 @@ class ObsSfcCorrectedParameters : public ObsOperatorParametersBase { 0.5, this}; - /// Note: Only relevant if \c SfcCorrectionType is set to GSL and and \c LapseRateOption is set to "Local" + /// Note: Only relevant if \c SfcCorrectionType is set to GSL, \c LapseRateOption is set to "Local", and \c Threshold is set to true oops::Parameter MaxThreshold {"max_threshold", "The maximum lapse rate (K/km) can be applied to adjust the " diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index 4622e36..f1cc36d 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -34,6 +34,7 @@ module ufo_sfccorrected_mod character(len=MAXVARLEN) :: lapse_rate_option real(kind_real) :: lapse_rate integer :: local_lapse_rate_level + logical :: threshold real(kind_real) :: min_threshold, max_threshold contains procedure :: setup => ufo_sfccorrected_setup @@ -55,6 +56,7 @@ subroutine ufo_sfccorrected_setup(self, f_conf) character(len=:), allocatable :: str_obs_height real(kind_real) :: constant_lapse_rate integer :: local_lapse_rate_level +logical :: threshold real(kind_real) :: min_threshold, max_threshold character(max_string) :: debug_msg @@ -91,10 +93,16 @@ subroutine ufo_sfccorrected_setup(self, f_conf) case ("Local") call f_conf%get_or_die("local_lapse_rate_level", local_lapse_rate_level) self%local_lapse_rate_level = local_lapse_rate_level - call f_conf%get_or_die("min_threshold", min_threshold) - self%min_threshold = min_threshold * 0.001 - call f_conf%get_or_die("max_threshold", max_threshold) - self%max_threshold = max_threshold * 0.001 + call f_conf%get_or_die("threshold", threshold) + self%threshold = threshold + if (self%threshold) then + call f_conf%get_or_die("min_threshold", min_threshold) + self%min_threshold = min_threshold * 0.001 + call f_conf%get_or_die("max_threshold", max_threshold) + self%max_threshold = max_threshold * 0.001 + end if + case ("NoAdjustment") + ! Do nothing in this case case default write(debug_msg,*) "ufo_sfccorrected: lapse_rate_option not recognized" call fckit_log%debug(debug_msg) @@ -349,7 +357,9 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) lr = -1. * (model_t%vals(ktop_lr,:) - model_t%vals(kbot,:)) / & (model_geomz%vals(ktop_lr,:) - model_geomz%vals(kbot,:)) - lr = min(self%max_threshold, max(self%min_threshold, lr)) + if (self%threshold) then + lr = min(self%max_threshold, max(self%min_threshold, lr)) + end if n_valid_obs = 0. n_obs_within_thresh = 0. @@ -384,14 +394,14 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) " K, average temperature (lower level): ", avg_temp_low, & " K, average AMSL (higher level): ", avg_height_high, & " m, average AMSL (lower level): ", avg_height_low, " m." - call fckit_log%info(err_msg) + call fckit_log%debug(err_msg) write(err_msg, '(3(a, f5.0))') & "Number of METAR observations: ", n_valid_obs, & ", number of observations within the thresholds: ", n_obs_within_thresh, & ", number of observations outside the thresholds: ", n_valid_obs-n_obs_within_thresh - call fckit_log%info(err_msg) + call fckit_log%debug(err_msg) - else + else if (self%lapse_rate_option.eq."Constant") then ! Local lapse rate in K/m lr = self%lapse_rate @@ -400,11 +410,21 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) write(err_msg, '(3a, f8.3, a)') & "Lapse rate option: ", trim(self%lapse_rate_option), & "Constant Lapse Rate: ", avg_lr, " K/km" - call fckit_log%info(err_msg) + call fckit_log%debug(err_msg) + + else + + write(err_msg, '(a)') & + "No terrain adjustment applied to temperature observations." + call fckit_log%debug(err_msg) endif - call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, lr) + if (self%lapse_rate_option.eq."NoAdjustment") then + cor_tsfc = obs_t + else + call da_int_lr(nobs, missing, cor_tsfc, obs_height, obs_t, model_zs, lr) + end if !----------- case default @@ -424,12 +444,12 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) ! cor_tsfc is corrected obs temp at model sfc ! do iobs = 1, nlocs - if (cor_tsfc(iobs) /= missing) then - ! for T_o2m, adjusted hofx - O = model_ts - T_o2m, OBS not adjusted. - hofx(ivar,iobs) = model_ts(iobs) + obs_t(iobs) - cor_tsfc(iobs) - else - hofx(ivar,iobs) = model_ts(iobs) - end if + if (cor_tsfc(iobs) /= missing) then +! for T_o2m, adjusted hofx - O = model_ts - T_o2m, OBS not adjusted. + hofx(ivar,iobs) = model_ts(iobs) + obs_t(iobs) - cor_tsfc(iobs) + else + hofx(ivar,iobs) = model_ts(iobs) + end if enddo enddo From bfc73178d105fd520c75004ed3c1f2596865188b Mon Sep 17 00:00:00 2001 From: Sijie Pan Date: Wed, 5 Mar 2025 14:27:18 -0600 Subject: [PATCH 13/13] write out lapse rate diag file --- .../sfccorrected/ufo_sfccorrected_mod.F90 | 94 +++++++++---------- 1 file changed, 45 insertions(+), 49 deletions(-) diff --git a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 index f1cc36d..8ae2ee1 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -16,6 +16,7 @@ module ufo_sfccorrected_mod use kinds use ufo_constants_mod, only : grav, rd, Lclr, t2tv use gnssro_mod_transform, only : geop2geometric, geometric2geop + use mpi implicit none private @@ -136,13 +137,14 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) real(kind_real) :: wf integer :: wi logical :: variable_present, variable_present_t, variable_present_q -real(kind_real), dimension(:), allocatable :: obs_height, obs_t, obs_q, obs_psfc, obs_lat, obs_tv +real(kind_real), dimension(:), allocatable :: obs_height, obs_t, obs_q, obs_psfc, obs_lat, obs_lon, obs_tv real(kind_real), dimension(:), allocatable :: model_ts, model_zs, model_level1, model_p_2000, model_t_2000, model_psfc, lr real(kind_real), dimension(:), allocatable :: H2000_geop real(kind_real), dimension(:), allocatable :: avg_tv real(kind_real) :: model_znew -real(kind_real) :: avg_lr, avg_temp_high, avg_temp_low, avg_height_high, avg_height_low, n_valid_obs, n_obs_within_thresh +integer :: rank, ierr, unit_lr +character(len=30) :: filename missing = missing_value(missing) nobs = obsspace_get_nlocs(obss) @@ -198,6 +200,28 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) allocate(model_level1(nobs)) allocate(model_psfc(nobs)) +! Sijie Pan, read lat and lon +if (.not. allocated(obs_lat)) then + variable_present = obsspace_has(obss, "MetaData", "latitude") + if (variable_present) then + call fckit_log%debug(' allocating obs_lat array') + allocate(obs_lat(nobs)) + call obsspace_get_db(obss, "MetaData", "latitude", obs_lat) + else + call abor1_ftn('Variable latitude@MetaData does not exist, aborting') + endif +endif +if (.not. allocated(obs_lon)) then + variable_present = obsspace_has(obss, "MetaData", "longitude") + if (variable_present) then + call fckit_log%debug(' allocating obs_lon array') + allocate(obs_lon(nobs)) + call obsspace_get_db(obss, "MetaData", "longitude", obs_lon) + else + call abor1_ftn('Variable longitude@MetaData does not exist, aborting') + endif +endif + ! If needed, we can convert geopotential heights to geometric altitude ! for full model vertical column using gnssro_mod_transform. We need ! to get the latitude of observation to do this. @@ -271,7 +295,7 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) enddo endif -if (allocated(obs_lat)) deallocate(obs_lat) +!if (allocated(obs_lat)) deallocate(obs_lat) model_psfc = model_ps%vals(1,:) @@ -357,60 +381,30 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) lr = -1. * (model_t%vals(ktop_lr,:) - model_t%vals(kbot,:)) / & (model_geomz%vals(ktop_lr,:) - model_geomz%vals(kbot,:)) + call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) + write(filename, '(a, i4.4, a)') "local_lapse_rate_", rank, ".txt" + unit_lr = 100+rank + open(unit=unit_lr, file=filename, status='unknown') + do iobs = 1, nobs + if (model_p%vals(kbot,iobs) > 0 .and. model_p%vals(ktop_lr,iobs) > 0) then + write(unit_lr, '(i5, 9(f10.3))') & + iobs, obs_lat(iobs), obs_lon(iobs), model_p%vals(kbot,iobs)/100.0, & + model_p%vals(ktop_lr,iobs)/100.0, model_geomz%vals(ktop_lr,iobs), & + model_geomz%vals(kbot,iobs), & + (model_p%vals(kbot,iobs) - model_p%vals(ktop_lr,iobs)) / 100.0, & + model_geomz%vals(ktop_lr,iobs) - model_geomz%vals(kbot,iobs), & + lr(iobs) * 1000.0 + end if + end do + if (self%threshold) then lr = min(self%max_threshold, max(self%min_threshold, lr)) end if - n_valid_obs = 0. - n_obs_within_thresh = 0. - avg_lr = 0. - avg_temp_high = 0. - avg_temp_low = 0. - avg_height_high = 0. - avg_height_low = 0. - do iobs=1, nobs - if (model_ts(iobs) > 0) then - n_valid_obs = n_valid_obs + 1 - avg_lr = avg_lr + lr(iobs) - avg_temp_high = avg_temp_high + model_t%vals(ktop_lr,iobs) - avg_temp_low = avg_temp_low + model_t%vals(kbot,iobs) - avg_height_high = avg_height_high + model_geomz%vals(ktop_lr,iobs) - avg_height_low = avg_height_low + model_geomz%vals(kbot,iobs) - if (lr(iobs) > 0.0005 .and. lr(iobs) < 0.0100) then - n_obs_within_thresh = n_obs_within_thresh + 1 - end if - end if - end do - avg_lr = avg_lr * 1000. / n_valid_obs - avg_temp_high = avg_temp_high / n_valid_obs - avg_temp_low = avg_temp_low / n_valid_obs - avg_height_high = avg_height_high / n_valid_obs - avg_height_low = avg_height_low / n_valid_obs - - write(err_msg, '(2a, 5(a,f8.3),a)') & - "Lapse rate option: ", trim(self%lapse_rate_option), & - ", average Local Lapse Rate: ", avg_lr, & - " K/km, average temperature (higher level): ", avg_temp_high, & - " K, average temperature (lower level): ", avg_temp_low, & - " K, average AMSL (higher level): ", avg_height_high, & - " m, average AMSL (lower level): ", avg_height_low, " m." - call fckit_log%debug(err_msg) - write(err_msg, '(3(a, f5.0))') & - "Number of METAR observations: ", n_valid_obs, & - ", number of observations within the thresholds: ", n_obs_within_thresh, & - ", number of observations outside the thresholds: ", n_valid_obs-n_obs_within_thresh - call fckit_log%debug(err_msg) - else if (self%lapse_rate_option.eq."Constant") then ! Local lapse rate in K/m lr = self%lapse_rate - avg_lr = self%lapse_rate * 1000. - - write(err_msg, '(3a, f8.3, a)') & - "Lapse rate option: ", trim(self%lapse_rate_option), & - "Constant Lapse Rate: ", avg_lr, " K/km" - call fckit_log%debug(err_msg) else @@ -455,6 +449,8 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) enddo deallocate(obs_height) +if (allocated(obs_lat)) deallocate(obs_lat) +if (allocated(obs_lon)) deallocate(obs_lon) if (variable_present_t) deallocate(obs_t) if (variable_present_q) deallocate(obs_q)