Skip to content

Commit

Permalink
Improved scaling options for climatology
Browse files Browse the repository at this point in the history 
Added options (through input.nml) for:
[1] user-defined scaling factor >> hyb_scaling > 0,
[2] kept previous functionality to compute the scaling >> hyb_scaling < 0,
[3] turning it off >> hyb_scaling = 0
  • Loading branch information
@mgharamti
mgharamti committed Mar 11, 2022
1 parent 8cf5b2b commit 4386b4e
Show file tree
Hide file tree
Showing 5 changed files with 60 additions and 53 deletions.
101 changes: 54 additions & 47 deletions assimilation_code/modules/assimilation/assim_tools_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -400,12 +400,12 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size,
integer, allocatable :: n_close_state_items(:), n_close_obs_items(:)

! Hybrid configuration variables
logical :: do_hybrid, hybrid_scaling, local_single_ss_hybrid
logical :: do_hybrid, local_single_ss_hybrid
integer :: stat_mean_copy, stat_sd_copy
integer :: stat_obs_mean_copy, stat_obs_var_copy
integer :: hybrid_ens_size
real(r8) :: my_hybrid_weight, my_hybrid_weight_sd
real(r8) :: orig_hybrid_weight
real(r8) :: hybrid_scaling, orig_hybrid_weight
real(r8) :: stat_obs_prior_mean, stat_obs_prior_var
real(r8) :: hyb_obs_var, corr_rho
real(r8), allocatable :: stat_obs_prior(:)
Expand Down Expand Up @@ -706,63 +706,70 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size,
! Rescaling factor so that the total variance of the
! rescaled static covariance more appropriately estimated
! the total forecast-error variance
if (do_hybrid .and. hybrid_scaling) then
i_qc = 0

allocate(dtrd(obs_ens_handle%num_vars), tr_R(obs_ens_handle%num_vars), &
tr_B(obs_ens_handle%num_vars))

! loop over all observations
do i = 1, obs_ens_handle%num_vars

call get_obs_from_key(obs_seq, keys(i), observation)
call get_obs_def(observation, obs_def)
if (do_hybrid) then

! value and variance of obs
obs_err_var = get_obs_def_error_variance(obs_def)
call get_obs_values(observation, obs, obs_val_index)

!print *, 'obs: ', obs(1), 'var: ', obs_err_var

call get_var_owner_index(ens_handle, int(i,i8), owner, owners_index)
if(hybrid_scaling < 0.0_r8) then

if(ens_handle%my_pe == owner) then
obs_qc = obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, owners_index)

if(nint(obs_qc)==0) then
i_qc = i_qc + 1

ens_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START:OBS_PRIOR_MEAN_END, owners_index)

dtrd(i_qc) = (obs(1)-ens_mean(1))**2
tr_R(i_qc) = obs_err_var
tr_B(i_qc) = stat_obs_ens_handle%copies(stat_obs_var_copy, owners_index)
i_qc = 0

allocate(dtrd(obs_ens_handle%num_vars), tr_R(obs_ens_handle%num_vars), &
tr_B(obs_ens_handle%num_vars))

! loop over all observations
do i = 1, obs_ens_handle%num_vars

call get_obs_from_key(obs_seq, keys(i), observation)
call get_obs_def(observation, obs_def)

! value and variance of obs
obs_err_var = get_obs_def_error_variance(obs_def)
call get_obs_values(observation, obs, obs_val_index)

call get_var_owner_index(ens_handle, int(i,i8), owner, owners_index)

if(ens_handle%my_pe == owner) then
obs_qc = obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, owners_index)

if(nint(obs_qc)==0) then
i_qc = i_qc + 1

ens_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START:OBS_PRIOR_MEAN_END, owners_index)

dtrd(i_qc) = (obs(1)-ens_mean(1))**2
tr_R(i_qc) = obs_err_var
tr_B(i_qc) = stat_obs_ens_handle%copies(stat_obs_var_copy, owners_index)

endif
endif
endif
enddo
sum_d = sum(dtrd(1:i_qc))
sum_R = sum(tr_R(1:i_qc))
sum_B = sum(tr_B(1:i_qc))

call sum_across_tasks(sum_d, sum_d_all)
call sum_across_tasks(sum_R, sum_R_all)
call sum_across_tasks(sum_B, sum_B_all)
enddo
sum_d = sum(dtrd(1:i_qc))
sum_R = sum(tr_R(1:i_qc))
sum_B = sum(tr_B(1:i_qc))

call sum_across_tasks(sum_d, sum_d_all)
call sum_across_tasks(sum_R, sum_R_all)
call sum_across_tasks(sum_B, sum_B_all)

!fs = (sum(dtrd(1:i_qc)) - sum(tr_R(1:i_qc)) ) / sum(tr_B(1:i_qc))
fs = (sum_d_all - sum_R_all) / sum_B_all
fs = max(0.0_r8, fs)
fs = (sum_d_all - sum_R_all) / sum_B_all
!print *, 'fs: ', fs
if (fs <= 0.0_r8 ) fs = 1.0_r8

deallocate(dtrd, tr_R, tr_B)
elseif(hybrid_scaling > 0.0_r8) then
fs = hybrid_scaling
else
fs = 1.0_r8
endif

!print *, 'scale: ', fs

!print *, 'input scale: ', hybrid_scaling

stat_ens_handle%copies(1:hybrid_ens_size, :) = sqrt(fs) * &
stat_ens_handle%copies(1:hybrid_ens_size, :)
stat_obs_ens_handle%copies(1:hybrid_ens_size, :) = sqrt(fs) * &
stat_obs_ens_handle%copies(1:hybrid_ens_size, :)
stat_obs_ens_handle%copies(stat_obs_var_copy, :) = fs * &
stat_obs_ens_handle%copies(stat_obs_var_copy, :)

deallocate(dtrd, tr_R, tr_B)
endif

! use MLOOP for the overall outer loop times; LG_GRN is for
Expand Down Expand Up @@ -2287,7 +2294,7 @@ subroutine update_from_hybobs_inc(obs, obs_prior_mean, obs_prior_var, obs_inc,
hyb_obs_prior_var = weight_factor * obs_prior_var + (1.0_r8 - weight_factor) * clim_obs_prior_var
hyb_obs_state_cov = weight_factor * obs_state_cov + (1.0_r8 - weight_factor) * clim_obs_state_cov

if (obs_prior_var > 0.0_r8) then
if (hyb_obs_prior_var > 0.0_r8) then
reg_coef = hyb_obs_state_cov / hyb_obs_prior_var
else
reg_coef = 0.0_r8
Expand Down
2 changes: 1 addition & 1 deletion assimilation_code/modules/assimilation/filter_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ module filter_mod
integer :: hyb_ens_size = 100
character(len=256) :: hyb_state_file_list(MAX_NUM_DOMS) = '' ! Name of files containing a list of static (climatology) state files
character(len=256) :: hyb_state_files(MAX_FILES) = ''
logical :: hyb_scaling = .true.
real(r8) :: hyb_scaling = 0.0_r8
logical :: hyb_initial_from_restart = .false.
logical :: hyb_sd_initial_from_restart = .false.
real(r8) :: hyb_weight_initial = 0.5_r8
Expand Down
6 changes: 3 additions & 3 deletions assimilation_code/modules/assimilation/hybrid_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ module hybrid_mod
integer :: flavor
integer :: sample_size
logical :: output_restart = .false.
logical :: scaling = .true.
real(r8) :: scaling = 0.0_r8
real(r8) :: weight, sd
real(r8) :: minmax_mean(2), minmax_sd(2)
logical :: mean_from_restart
Expand Down Expand Up @@ -138,7 +138,7 @@ subroutine adaptive_hybrid_init(hybrid_handle, hyb_flavor, ens_size, mean_from_r
logical, intent(in) :: mean_from_restart
logical, intent(in) :: sd_from_restart
logical, intent(in) :: output_hybrid
logical, intent(in) :: hyb_scale
real(r8), intent(in) :: hyb_scale
real(r8), intent(in) :: hyb_weight_initial, hyb_weight_sd_initial

! Record the module version if this is first initialize call
Expand Down Expand Up @@ -192,7 +192,7 @@ end function hyb_sd_from_restart
function scale_static_background(hybrid_handle)

type(hybrid_type) :: hybrid_handle
logical :: scale_static_background
real(r8) :: scale_static_background

scale_static_background = hybrid_handle%scaling

Expand Down
2 changes: 1 addition & 1 deletion assimilation_code/modules/assimilation/smoother_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ subroutine init_smoother(ens_handle, POST_INF_COPY, POST_INF_SD_COPY)
1.0_r8, 1.0_r8, 0.0_r8, 1.0_r8, ens_handle, allow_missing, "Lag")

! Define a dummy hybrid handle not supported for smoothers
call adaptive_hybrid_init(lag_hybrid, 0, 1, .false., .false., .false., .false., 1.0_r8, 0.0_r8)
call adaptive_hybrid_init(lag_hybrid, 0, 1, .false., .false., .false., 0.0_r8, 1.0_r8, 0.0_r8)

end subroutine init_smoother

Expand Down
2 changes: 1 addition & 1 deletion models/wrf_hydro/work/input.nml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@
hyb_flavor = 2
hyb_ens_size = 1000
hyb_state_file_list = 'static_file_list.txt'
hyb_scaling = .true.
hyb_scaling = 0
hyb_initial_from_restart = .false.
hyb_sd_initial_from_restart = .false.
hyb_weight_initial = 0.5
Expand Down

0 comments on commit 4386b4e

Please sign in to comment.