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.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/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 22f8ecc..7eceaf8 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, GSL }; struct SfcCorrectionTypeParameterTraitsHelper { typedef SfcCorrectionType EnumType; @@ -30,7 +30,7 @@ struct SfcCorrectionTypeParameterTraitsHelper { static constexpr util::NamedEnumerator namedValues[] = { { SfcCorrectionType::UKMO, "UKMO" }, { SfcCorrectionType::WRFDA, "WRFDA" }, - { SfcCorrectionType::RRFS_GSL, "RRFS_GSL" } + { SfcCorrectionType::GSL, "GSL" } }; }; @@ -64,23 +64,66 @@ 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 /// 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, 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 + {"lapse_rate", + "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, + this}; + + /// 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, " + "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}; + + /// 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 " + "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, \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 " + "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 73c6d40..8ae2ee1 100644 --- a/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 +++ b/src/ufo/operators/sfccorrected/ufo_sfccorrected_mod.F90 @@ -9,12 +9,14 @@ module ufo_sfccorrected_mod use oops_variables_mod + use obs_variables_mod use ufo_vars_mod use missing_values_mod use iso_c_binding use kinds use ufo_constants_mod, only : grav, rd, Lclr, t2tv use gnssro_mod_transform, only : geop2geometric, geometric2geop + use mpi implicit none private @@ -23,19 +25,24 @@ 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 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_level + logical :: threshold + real(kind_real) :: min_threshold, max_threshold 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(6) :: geovars_list = (/ var_ps, var_geomz, var_sfc_geomz, var_ts, var_prs, var_sfc_t2m /) contains @@ -46,8 +53,13 @@ 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) :: constant_lapse_rate +integer :: local_lapse_rate_level +logical :: threshold +real(kind_real) :: min_threshold, max_threshold + character(max_string) :: debug_msg !> In case where a user wants to specify geoVaLs variable name of model @@ -72,6 +84,33 @@ 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."GSL") then + 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 * 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("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) + call abor1_ftn(debug_msg) + end select +endif + end subroutine ufo_sfccorrected_setup ! ------------------------------------------------------------------------------ @@ -90,20 +129,23 @@ 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 +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 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 :: 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 +integer :: rank, ierr, unit_lr +character(len=30) :: filename + missing = missing_value(missing) nobs = obsspace_get_nlocs(obss) @@ -130,20 +172,22 @@ 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 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 @@ -156,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. @@ -229,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,:) @@ -296,9 +362,63 @@ subroutine ufo_sfccorrected_simobs(self, geovals, obss, nvars, nlocs, hofx) deallocate(model_t_2000) !---------------- -case ("RRFS_GSL") +case ("GSL") !---------------- -! to be implemented + + model_ts = model_t2m%vals(1,:) + + allocate(lr(nobs)) + + if (self%lapse_rate_option.eq."Local") then + + 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,:)) + + 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 + + else if (self%lapse_rate_option.eq."Constant") then + + ! Local lapse rate in K/m + lr = self%lapse_rate + + else + + write(err_msg, '(a)') & + "No terrain adjustment applied to temperature observations." + call fckit_log%debug(err_msg) + + endif + + 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 @@ -329,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) @@ -448,5 +570,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: 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_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, T_m2o +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_o +end where + + cor_tsfc = T_o2m + +end subroutine da_int_lr + ! ------------------------------------------------------------------------------ end module ufo_sfccorrected_mod