diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index c064e7cbe1..eaf8e189dc 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -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(:) @@ -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 @@ -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 diff --git a/assimilation_code/modules/assimilation/filter_mod.f90 b/assimilation_code/modules/assimilation/filter_mod.f90 index 029bc813f0..e3542d4430 100644 --- a/assimilation_code/modules/assimilation/filter_mod.f90 +++ b/assimilation_code/modules/assimilation/filter_mod.f90 @@ -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 diff --git a/assimilation_code/modules/assimilation/hybrid_mod.f90 b/assimilation_code/modules/assimilation/hybrid_mod.f90 index 017f3e151a..57ebae8b10 100644 --- a/assimilation_code/modules/assimilation/hybrid_mod.f90 +++ b/assimilation_code/modules/assimilation/hybrid_mod.f90 @@ -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 @@ -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 @@ -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 diff --git a/assimilation_code/modules/assimilation/smoother_mod.f90 b/assimilation_code/modules/assimilation/smoother_mod.f90 index 1854e9fdf6..f32419acfc 100644 --- a/assimilation_code/modules/assimilation/smoother_mod.f90 +++ b/assimilation_code/modules/assimilation/smoother_mod.f90 @@ -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 diff --git a/models/wrf_hydro/work/input.nml b/models/wrf_hydro/work/input.nml index 1b7e2e5670..d711fd37f0 100644 --- a/models/wrf_hydro/work/input.nml +++ b/models/wrf_hydro/work/input.nml @@ -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