From 82e510cfb308dee6d6a5caa52707cd33f2e4aae5 Mon Sep 17 00:00:00 2001 From: Moha Gharamti Date: Sun, 19 Feb 2023 17:47:39 -0700 Subject: [PATCH] Hybrid Flavor 3: Adaptive in space and time This is a relatively large commit. Changes to the code took place while running wrf_hydro and DART experiments for hurricane Ian. The major new capability is the support for spatially and temporally varying adaptive hybrid scheme (i.e., flavor #3). Other changes: - Found a bug in filter_mod where extra copies for the static ensemble handle were not being recoreded properly. - Various updates to the hydroDART python scripts for better file movement. - Streamflow obs converter now checks for NaN discharge. This is useful especially for provisional (i.e., near-real-time) data - Add another scalar for mpi communications - Footprint of flavor #3 is mainly in assim_tools_mod and hybrid_mod --- .../modules/assimilation/assim_tools_mod.f90 | 237 ++++++++++-------- .../modules/assimilation/filter_mod.f90 | 8 +- .../modules/assimilation/hybrid_mod.f90 | 60 ++--- .../modules/utilities/mpi_utilities_mod.f90 | 137 ++++++---- .../utilities/null_mpi_utilities_mod.f90 | 40 +-- .../create_identity_streamflow_obs.f90 | 2 +- .../hydrodartpy/core/run_filter_experiment.py | 15 ++ .../hydrodartpy/core/setup_experiment.py | 5 + .../perturb_channel_bucket_only_forcing.py | 1 + .../wrf_hydro/shell_scripts/DART_cleanup.csh | 20 +- 10 files changed, 303 insertions(+), 222 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index eaf8e189dc..4145ce926b 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -60,8 +60,9 @@ module assim_tools_mod inflate_ens, adaptive_inflate_type, & deterministic_inflate, solve_quadratic -use hybrid_mod, only : hybrid_type, do_single_ss_hybrid, get_hybrid_sample_size, & - scale_static_background, update_hybrid +use hybrid_mod, only : hybrid_type, do_single_ss_hybrid, get_hybrid_sample_size, & + scale_static_background, do_constant_hybrid, update_hybrid, & + do_varying_ss_hybrid use time_manager_mod, only : time_type, get_time @@ -400,18 +401,24 @@ 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, local_single_ss_hybrid +logical :: do_hybrid +logical :: local_constant_hybrid +logical :: local_single_ss_hybrid +logical :: local_varying_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) :: hybrid_scaling, orig_hybrid_weight +real(r8) :: hybrid_scaling real(r8) :: stat_obs_prior_mean, stat_obs_prior_var real(r8) :: hyb_obs_var, corr_rho +real(r8) :: vary_ss_hybrid_mean, vary_ss_hybrid_sd +real(r8), allocatable :: orig_hybrid_weight(:) real(r8), allocatable :: stat_obs_prior(:) real(r8), allocatable :: dtrd(:), tr_R(:), tr_B(:) real(r8) :: sum_d, sum_R, sum_B, sum_d_all, sum_R_all, sum_B_all -real(r8) :: fs, ens_mean(1) +real(r8) :: fs, ens_mean(1), orig_hyb_mean integer :: i_qc ! Are we hybridizing the increments? @@ -428,26 +435,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, do_hybrid = .false. endif - -!print *, 'ens_size: ' , ens_size -!print *, 'hybrid_ens_size: ' , hybrid_ens_size -!print *, 'do_hybrid: ' , do_hybrid -!print *, 'weight val: ' , ens_handle%copies(HYB_MEAN_COPY, 1) -!print *, 'weight sd: ' , ens_handle%copies(HYB_SD_COPY, 1) -!print *, 'ENS_MEAN_COPY: ' , ENS_MEAN_COPY -!print *, 'ENS_SD_COPY: ' , ENS_SD_COPY -!print *, 'ENS_INF_COPY: ' , ENS_INF_COPY -!print *, 'ENS_INF_SD_COPY: ' , ENS_INF_SD_COPY -!print *, 'HYB_MEAN_COPY: ' , HYB_MEAN_COPY -!print *, 'HYB_SD_COPY: ' , HYB_SD_COPY -!print *, 'OBS_KEY_COPY: ' , OBS_KEY_COPY -!print *, 'OBS_GLOBAL_QC_COPY: ' , OBS_GLOBAL_QC_COPY -!print *, 'OBS_PRIOR_MEAN_START: ', OBS_PRIOR_MEAN_START -!print *, 'OBS_PRIOR_MEAN_END: ' , OBS_PRIOR_MEAN_END -!print *, 'OBS_PRIOR_VAR_START: ' , OBS_PRIOR_VAR_START -!print *, 'OBS_PRIOR_VAR_END: ' , OBS_PRIOR_VAR_END - - ! timing disabled by default timing(:) = .false. t_base(:) = 0.0_r8 @@ -458,24 +445,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, allocate(n_close_state_items(obs_ens_handle%num_vars), & n_close_obs_items( obs_ens_handle%num_vars)) -! turn these on carefully - they can generate a lot of output! -! also, to be readable - at least with ifort: -! setenv FORT_FMT_RECL 1024 -! so output lines don't wrap. - -!timing(MLOOP) = .true. -!timing(LG_GRN) = .true. - if (timing(MLOOP)) allocate(elapse_array(obs_ens_handle%num_vars)) -! use maxitems limit here or drown in output. -!timing(SM_GRN) = .false. -!t_limit(SM_GRN) = 4_i8 - -!timing(GC) = .true. -!t_limit(GC) = 4_i8 - - ! allocate rather than dump all this on the stack allocate(close_obs_dist( obs_ens_handle%my_num_vars), & last_close_obs_dist(obs_ens_handle%my_num_vars), & @@ -540,10 +511,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ! Hybrid configuration if (do_hybrid) then - hybrid_scaling = scale_static_background(hybrid) - local_single_ss_hybrid = do_single_ss_hybrid(hybrid) + hybrid_scaling = scale_static_background(hybrid) + local_constant_hybrid = do_constant_hybrid(hybrid) + local_single_ss_hybrid = do_single_ss_hybrid(hybrid) + local_varying_ss_hybrid = do_varying_ss_hybrid(hybrid) endif +print *, '' +print *, 'constant_hybrid : ', local_constant_hybrid +print *, 'single_ss_hybrid : ', local_single_ss_hybrid +print *, 'varying_ss_hybrid: ', local_varying_ss_hybrid + ! Default to printing nothing nth_obs = -1 @@ -583,18 +561,18 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, end if ! Hybrid: if adaptive in time, but spatially constant -if (do_hybrid .and. local_single_ss_hybrid) then - ! Any element should work, here we go with the first one - my_hybrid_weight = ens_handle%copies(HYB_MEAN_COPY, 1) - my_hybrid_weight_sd = ens_handle%copies(HYB_SD_COPY , 1) +if (do_hybrid) then + + my_hybrid_weight = ens_handle%copies(HYB_MEAN_COPY, 1) + my_hybrid_weight_sd = ens_handle%copies(HYB_SD_COPY , 1) ! Use the original weight for hybridizing - orig_hybrid_weight = my_hybrid_weight + ! Need to change for spatially-varying form + allocate(orig_hybrid_weight(ens_handle%num_vars)) + orig_hybrid_weight = ens_handle%copies(HYB_MEAN_COPY, :) + !ens_handle%copies(ENS_MEAN_COPY, :) = ens_handle%copies(HYB_MEAN_COPY, :) endif -!print *, 'my_hybrid_weight: ' , my_hybrid_weight -!print *, 'my_hybrid_weight_sd: ', my_hybrid_weight_sd - ! Get info on my number and indices for obs my_num_obs = get_my_num_vars(obs_ens_handle) call get_my_vars(obs_ens_handle, my_obs_indx) @@ -608,9 +586,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ! do the forward operator calculation call get_my_obs_loc(obs_ens_handle, obs_seq, keys, my_obs_loc, my_obs_kind, my_obs_type, obs_time) -!print *, 'my_obs_kind: ', my_obs_kind -!print *, 'my_obs_type: ', my_obs_type - if (convert_all_obs_verticals_first .and. is_doing_vertical_conversion) then ! convert the vertical of all my observations to the localization coordinate if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) @@ -751,7 +726,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, call sum_across_tasks(sum_B, sum_B_all) 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) @@ -761,8 +736,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, fs = 1.0_r8 endif - !print *, 'scale: ', fs - !print *, 'input scale: ', hybrid_scaling + 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, :) @@ -781,8 +756,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ! Loop through all the (global) observations sequentially SEQUENTIAL_OBS: do i = 1, obs_ens_handle%num_vars - !print *, 'obs i: ', i - if (timing(MLOOP)) call start_timer(t_base(MLOOP)) if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) @@ -814,16 +787,21 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, else call get_state_meta_data(-1 * int(base_obs_type,i8), dummyloc, base_obs_kind) ! identity obs endif - - !print *, 'base_obs_type: ', base_obs_type - !print *, 'base_obs_kind: ', base_obs_kind + + ! Weight coefficient at this obs location + if (do_hybrid) then + ! I don't know the exact obs location, take average in space + if (base_obs_type > 0 ) then + orig_hyb_mean = sum(orig_hybrid_weight)/size(orig_hybrid_weight) + else + ! Identity case: + orig_hyb_mean = orig_hybrid_weight(-1 * int(base_obs_type,i8)) + endif + endif ! Get the value of the observation call get_obs_values(observation, obs, obs_val_index) - !print *, 'obs_val: ', obs(1) - !print *, '' - ! Find out who has this observation and where it is call get_var_owner_index(ens_handle, int(i,i8), owner, owners_index) @@ -851,6 +829,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, OBS_PRIOR_MEAN_END, owners_index) orig_obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START: & OBS_PRIOR_VAR_END, owners_index) + print *, 'ensemble mean: ', orig_obs_prior_mean(1) ! Compute observation space increments for each group do group = 1, num_groups @@ -869,18 +848,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ! Find the static variance for this observation hyb_obs_var = stat_obs_ens_handle%copies(stat_obs_var_copy, owners_index) - - call obs_increment(obs_prior(grp_bot:grp_top), grp_size, obs(1), & - obs_err_var, obs_inc(grp_bot:grp_top), inflate, my_inflate, & - my_inflate_sd, net_a(group), orig_hybrid_weight, stat_obs_prior) + + call obs_increment(obs_prior, ens_size, obs(1), & + obs_err_var, obs_inc, inflate, my_inflate, & + my_inflate_sd, net_a(1), orig_hyb_mean, stat_obs_prior) endif end do - !print *, 'obs_inc: ', obs_inc - + ! Update the hybrid weight value - if (do_hybrid .and. local_single_ss_hybrid) then + SINGLE_SS_HYBRID: if (do_hybrid .and. local_single_ss_hybrid .and. .not. inflate_only ) then if(my_hybrid_weight > 0.0_r8 .and. my_hybrid_weight_sd > 0.0_r8) then - !print *, 'Performing the update for alpha' + print *, 'Performing the update for alpha: Spatially-constant' ens_obs_mean = orig_obs_prior_mean(1) ens_obs_var = orig_obs_prior_var(1) @@ -888,22 +866,23 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ! Spatially-constant (correlation is uniform = 1) corr_rho = 1.0_r8 - !print *, 'a = ', my_hybrid_weight + print *, 'm = ', my_hybrid_weight + print *, 'v = ', my_hybrid_weight_sd**2 call update_hybrid(my_hybrid_weight, my_hybrid_weight_sd, & ens_obs_var, hyb_obs_var, obs_err_var, ens_obs_mean, & obs(1), corr_rho) - !print *, 'se2 = ', ens_obs_var - !print *, 'ss2 = ', hyb_obs_var - !print *, 'so2 = ', obs_err_var - !print *, 'd2 = ', (ens_obs_mean - obs(1))**2 - !print *, '' - !print *, 'post my_hybrid_weight: ', my_hybrid_weight - !print *, 'post my_hybrid_weight_sd: ', my_hybrid_weight_sd - !print *, '' + print *, 'se2 = ', ens_obs_var + print *, 'ss2 = ', hyb_obs_var + print *, 'so2 = ', obs_err_var + print *, 'd2 = ', (ens_obs_mean - obs(1))**2 + print *, '' + print *, 'post my_hybrid_weight_mean: ', my_hybrid_weight + print *, 'post my_hybrid_weight_sd : ', my_hybrid_weight_sd + print *, '' endif - endif + endif SINGLE_SS_HYBRID ! Compute updated values for single state space inflation SINGLE_SS_INFLATE: if(local_single_ss_inflate) then @@ -961,17 +940,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, whichvert_real = real(whichvert_obs_in_localization_coord, r8) if(local_varying_ss_inflate) then call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, obs_inc, & - orig_obs_prior_mean, orig_obs_prior_var, net_a, stat_obs_prior, scalar1=obs_qc, & + orig_obs_prior_mean, orig_obs_prior_var, net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real, & - scalar4=my_hybrid_weight, scalar5=my_hybrid_weight_sd) + scalar4=my_hybrid_weight, scalar5=my_hybrid_weight_sd, scalar6=hyb_obs_var) else if(local_single_ss_inflate .or. local_obs_inflate) then call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, obs_inc, & - net_a, stat_obs_prior, scalar1=my_inflate, scalar2=my_inflate_sd, scalar3=obs_qc, & + net_a, scalar1=my_inflate, scalar2=my_inflate_sd, scalar3=obs_qc, & scalar4=vertvalue_obs_in_localization_coord, scalar5=whichvert_real) else call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, obs_inc, & - net_a, stat_obs_prior, scalar1=obs_qc, & + net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real, & scalar4=my_hybrid_weight, scalar5=my_hybrid_weight_sd) endif @@ -986,17 +965,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, !>the cost of sending unneeded values if(local_varying_ss_inflate) then call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, obs_inc, & - orig_obs_prior_mean, orig_obs_prior_var, net_a, stat_obs_prior, scalar1=obs_qc, & + orig_obs_prior_mean, orig_obs_prior_var, net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real, & - scalar4=my_hybrid_weight, scalar5=my_hybrid_weight_sd) + scalar4=my_hybrid_weight, scalar5=my_hybrid_weight_sd, scalar6=hyb_obs_var) else if(local_single_ss_inflate .or. local_obs_inflate) then call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, obs_inc, & - net_a, stat_obs_prior, scalar1=my_inflate, scalar2=my_inflate_sd, scalar3=obs_qc, & + net_a, scalar1=my_inflate, scalar2=my_inflate_sd, scalar3=obs_qc, & scalar4=vertvalue_obs_in_localization_coord, scalar5=whichvert_real) else call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, obs_inc, & - net_a, stat_obs_prior, scalar1=obs_qc, & + net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real, & scalar4=my_hybrid_weight, scalar5=my_hybrid_weight_sd) endif @@ -1217,11 +1196,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ! Loop through to update each of my state variables that is potentially close if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) + STATE_UPDATE: do j = 1, num_close_states state_index = close_state_ind(j) - - !print *, 'j: ', j - !print *, 'state_index: ', state_index ! the "any" is an expensive test when you do it for every ob. don't test ! if we know there aren't going to be missing values in the state. @@ -1231,6 +1208,12 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE endif + ! Current hybrid weight values in space: + if (local_varying_ss_hybrid) then + vary_ss_hybrid_mean = ens_handle%copies(HYB_MEAN_COPY, state_index) + vary_ss_hybrid_sd = ens_handle%copies(HYB_SD_COPY, state_index) + endif + ! Get the initial values of inflation for this variable if state varying inflation if(local_varying_ss_inflate) then varying_ss_inflate = ens_handle%copies(ENS_INF_COPY, state_index) @@ -1239,7 +1222,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, varying_ss_inflate = 0.0_r8 varying_ss_inflate_sd = 0.0_r8 endif - + ! Compute the distance and covariance factor cov_factor = comp_cov_factor(close_state_dist(j), cutoff_rev, & base_obs_loc, base_obs_type, my_state_loc(state_index), my_state_kind(state_index)) @@ -1257,20 +1240,21 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, if(cov_factor <= 0.0_r8) cycle STATE_UPDATE if (timing(SM_GRN)) call start_timer(t_base(SM_GRN), t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) - if (do_hybrid) then ! Regression for the hybrid case: Only support 1 group for now; i.e., num_groups = 1 if(local_varying_ss_inflate .and. varying_ss_inflate > 0.0_r8 .and. varying_ss_inflate_sd > 0.0_r8) then call update_from_hybobs_inc(obs_prior, obs_prior_mean(1), obs_prior_var(1), & obs_inc, ens_handle%copies(1:ens_size, state_index), ens_size, & stat_obs_prior, stat_obs_prior_mean, stat_obs_prior_var, & stat_ens_handle%copies(1:hybrid_ens_size, state_index), hybrid_ens_size, & - orig_hybrid_weight, increment, reg_coef(1), net_a(1), correl(1)) + !ens_handle%copies(ENS_MEAN_COPY, state_index), + orig_hybrid_weight(state_index), increment, reg_coef(1), net_a(1), correl(1)) else call update_from_hybobs_inc(obs_prior, obs_prior_mean(1), obs_prior_var(1), & obs_inc, ens_handle%copies(1:ens_size, state_index), ens_size, & stat_obs_prior, stat_obs_prior_mean, stat_obs_prior_var, & stat_ens_handle%copies(1:hybrid_ens_size, state_index), hybrid_ens_size, & - orig_hybrid_weight, increment, reg_coef(1), net_a(1)) + !ens_handle%copies(ENS_MEAN_COPY, state_index), + orig_hybrid_weight(state_index), increment, reg_coef(1), net_a(1)) endif else ! Loop through groups to update the state variable ensemble members @@ -1314,6 +1298,42 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, ens_handle%copies(1:ens_size, state_index) + reg_factor * increment endif + ! Compute spatially-varying hybrid weighting coefficient + VARYING_SS_HYBRID: if (do_hybrid .and. local_varying_ss_hybrid .and. .not. inflate_only) then + if (vary_ss_hybrid_mean > 0.0_r8 .and. vary_ss_hybrid_sd > 0.0_r8) then + print *, 'Performing the update for alpha: spatially_varying' + + ens_obs_mean = orig_obs_prior_mean(1) + ens_obs_var = orig_obs_prior_var(1) + + corr_rho = reg_factor * abs(correl(1)) + + print *, 'reg_factor: ', reg_factor + print *, 'correl: ', abs(correl(1)) + print *, 'x = ', j + print *, 'm = ', vary_ss_hybrid_mean + print *, 'v = ', vary_ss_hybrid_sd**2 + + call update_hybrid(vary_ss_hybrid_mean, vary_ss_hybrid_sd, & + ens_obs_var, hyb_obs_var, obs_err_var, ens_obs_mean, & + obs(1), corr_rho) + + print *, 'rho = ', corr_rho + print *, 'se2 = ', ens_obs_var + print *, 'ss2 = ', hyb_obs_var + print *, 'so2 = ', obs_err_var + print *, 'd2 = ', (ens_obs_mean - obs(1))**2 + print *, '' + print *, 'post my_hybrid_weight_mean: ', vary_ss_hybrid_mean + print *, 'post my_hybrid_weight_sd : ', vary_ss_hybrid_sd + print *, '' + + !Record the updated weights + ens_handle%copies(HYB_MEAN_COPY, state_index) = vary_ss_hybrid_mean + ens_handle%copies(HYB_SD_COPY, state_index) = vary_ss_hybrid_sd + endif + endif VARYING_SS_HYBRID + ! Compute spatially-varying state space inflation if(local_varying_ss_inflate) then ! base is the initial inflate value for this state variable @@ -1386,14 +1406,12 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, endif !call test_state_copies(ens_handle, 'after_state_updates') - !------------------------------------------------------ - ! Now everybody updates their obs priors (only ones after this one) if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) OBS_UPDATE: do j = 1, num_close_obs obs_index = close_obs_ind(j) - + ! Only have to update obs that have not yet been used if(my_obs_indx(obs_index) > i) then @@ -1417,15 +1435,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, if(cov_factor <= 0.0_r8) cycle OBS_UPDATE if (timing(SM_GRN)) call start_timer(t_base(SM_GRN), t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) - + ! Loop through and update ensemble members in each group if (do_hybrid) then - call update_from_hybobs_inc(obs_prior, obs_prior_mean(1), obs_prior_var(1), & - obs_inc, obs_ens_handle%copies(1:ens_size, obs_index), ens_size, & - stat_obs_prior, stat_obs_prior_mean, stat_obs_prior_var, & + call update_from_hybobs_inc(obs_prior, obs_prior_mean(1), obs_prior_var(1), obs_inc, & + obs_ens_handle%copies(1:ens_size, obs_index), ens_size, & + stat_obs_prior, stat_obs_prior_mean, stat_obs_prior_var, & stat_obs_ens_handle%copies(1:hybrid_ens_size, obs_index), hybrid_ens_size, & - orig_hybrid_weight, increment, reg_coef(1), net_a(1)) + orig_hyb_mean, increment, reg_coef(1), net_a(1)) else - ! Loop through and update ensemble members in each group do group = 1, num_groups grp_bot = grp_beg(group) grp_top = grp_end(group) @@ -1433,7 +1450,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, obs_prior_var(group), obs_inc(grp_bot:grp_top), & obs_ens_handle%copies(grp_bot:grp_top, obs_index), grp_size, & increment(grp_bot:grp_top), reg_coef(group), net_a(group)) - end do + enddo endif if (timing(SM_GRN)) call read_timer(t_base(SM_GRN), 'update_from_obs_inc_O', & @@ -1481,7 +1498,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, end if ! Every pe needs to get the current my_hybrid and my_hybrid_sd back -if(local_single_ss_hybrid) then +if(local_single_ss_hybrid .and. .not. inflate_only) then ens_handle%copies(HYB_MEAN_COPY, :) = my_hybrid_weight ens_handle%copies(HYB_SD_COPY, :) = my_hybrid_weight_sd end if @@ -1582,7 +1599,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, ens_size, deallocate(n_close_state_items, & n_close_obs_items) -if (do_hybrid) deallocate(stat_obs_prior) +if (do_hybrid) deallocate(stat_obs_prior, orig_hybrid_weight) ! end dealloc end subroutine filter_assim diff --git a/assimilation_code/modules/assimilation/filter_mod.f90 b/assimilation_code/modules/assimilation/filter_mod.f90 index e3542d4430..f93b84ea5f 100644 --- a/assimilation_code/modules/assimilation/filter_mod.f90 +++ b/assimilation_code/modules/assimilation/filter_mod.f90 @@ -632,6 +632,9 @@ subroutine filter_main() call set_num_extra_copies(state_ens_handle, num_extras) +! For the static members, we only have 2 extra copies: meand and sd +if (do_hybrid) call set_num_extra_copies(static_state_ens_handle, 2) + call trace_message('After setting up space for ensembles') ! Don't currently support number of processes > model_size @@ -1006,7 +1009,7 @@ subroutine filter_main() if (do_hybrid) then call allocate_single_copy(static_obs_ens_handle, static_prior_qc_copy) - + call get_obs_ens_distrib_state(static_state_ens_handle, static_obs_ens_handle, & static_qc_ens_handle, seq, keys, obs_val_index, input_qc_index, & STATIC_OBS_ERR_VAR_COPY, STATIC_OBS_VAL_COPY, STATIC_OBS_KEY_COPY, STATIC_OBS_GLOBAL_QC_COPY, & @@ -1019,7 +1022,7 @@ subroutine filter_main() !print *, 'obs_hyb: ', static_obs_ens_handle%copies(1:hyb_ens_size, 1) !print *, 'rest: ' , static_obs_ens_handle%copies(hyb_ens_size+1:TOTAL_STATIC_OBS_COPIES, 1) - + call timestamp_message('After computing prior observation values') call trace_message('After computing prior observation values') @@ -1067,7 +1070,6 @@ subroutine filter_main() call timestamp_message('Before observation assimilation') - if (.not. do_hybrid) then call filter_assim(state_ens_handle, obs_fwd_op_ens_handle, seq, keys, & ens_size, num_groups, obs_val_index, prior_inflate, hybridization, & diff --git a/assimilation_code/modules/assimilation/hybrid_mod.f90 b/assimilation_code/modules/assimilation/hybrid_mod.f90 index 57ebae8b10..0d5b4156cf 100644 --- a/assimilation_code/modules/assimilation/hybrid_mod.f90 +++ b/assimilation_code/modules/assimilation/hybrid_mod.f90 @@ -26,14 +26,15 @@ module hybrid_mod set_hybrid_sd_copy, do_single_ss_hybrid, get_hybrid_mean_copy, & get_hybrid_sd_copy, do_ss_hybrid, get_hybrid_mean, get_hybrid_sd, & do_varying_ss_hybrid, hyb_minmax_task_zero, log_hybrid_info, & - update_hybrid, get_hybrid_sample_size, scale_static_background + update_hybrid, get_hybrid_sample_size, scale_static_background, & + do_constant_hybrid character(len=*), parameter :: source = 'hybrid_mod.f90' -integer, parameter :: NO_HYBRID = 0 -integer, parameter :: CONSTANT_HYBRID = 1 -integer, parameter :: SINGLE_SS_HYBRID = 2 -integer, parameter :: VARYING_SS_HYBRID = 3 +integer, parameter :: NO_HYBRID = 0 !only dynamic ensemble, not climatology +integer, parameter :: CONSTANT_HYBRID = 1 !with climatology but weight not adaptive in time +integer, parameter :: SINGLE_SS_HYBRID = 2 !adaptive in time, constant in space +integer, parameter :: VARYING_SS_HYBRID = 3 !adaptive in time, varying in space ! Type to keep track of information for hybrid type hybrid_type @@ -98,13 +99,13 @@ subroutine validate_hybrid_options(inf_flavor, hyb_flavor, hyb_ens_size, & ! Adaptive hybrid scheme is only spatially-constant for now, ! Inflation is relied on to spread the information properly in space ! Only support hybrid with inf-flavors 0, 2, 4, 5 (1 is deprecated) -if (do_hybrid .and. hyb_flavor >= SINGLE_SS_HYBRID) then +if (do_hybrid .and. hyb_flavor <= SINGLE_SS_HYBRID) then if (inf_flavor(1) == SINGLE_SS_INFLATION .or. inf_flavor(2) == SINGLE_SS_INFLATION) then - write(string1, '(A, i2, A)') 'Adaptive hybrid scheme does not support inf_flavor:', SINGLE_SS_HYBRID, ' (i.e., spatially-uniform inflation)' + write(string1, '(A, i2, A, i2, A)') 'Adaptive hybrid scheme flavor:', hyb_flavor, & + ' does not support inf_flavor:', SINGLE_SS_INFLATION, ' (i.e., spatially-uniform inflation)' call error_handler(E_ERR, 'validate_hybrid_options:', string1, source) endif -endif - +endif ! Check the size of the climatological ensemble ! We need this to be fairly large ~O(1000) @@ -116,12 +117,12 @@ subroutine validate_hybrid_options(inf_flavor, hyb_flavor, hyb_ens_size, & if (do_hybrid) output_hybrid = .true. -! Observation space inflation not currently supported -if(hyb_flavor == VARYING_SS_HYBRID) then - call error_handler(E_ERR, 'validate_hybrid_options', & - 'Varying state-space hybridization (type 3) not currently supported', source, & - text2 = 'Look out for this type in future DART releases.') -endif +! Spatially-varying hybrid not currently supported +!if(hyb_flavor == VARYING_SS_HYBRID) then +! call error_handler(E_ERR, 'validate_hybrid_options', & +! 'Varying state-space hybridization (type 3) not currently supported', source, & +! text2 = 'Look out for this type in future DART releases.') +!endif end subroutine validate_hybrid_options @@ -239,6 +240,18 @@ function do_single_ss_hybrid(hybrid_handle) end function do_single_ss_hybrid +!------------------------------------------------------------------------------- +!> Returns true if time-invarian and fixed in space and time + +function do_constant_hybrid(hybrid_handle) + +logical :: do_constant_hybrid +type(hybrid_type), intent(in) :: hybrid_handle + +do_constant_hybrid = (hybrid_handle%flavor == CONSTANT_HYBRID) + +end function do_constant_hybrid + !------------------------------------------------------------------------------- !> Given prior mean and sd values for the hybrid weighting coefficient, @@ -265,7 +278,6 @@ subroutine update_hybrid(weight, weight_sd, se2, ss2, so2, x, obs, rho) real(r8), intent(in) :: obs, so2 real(r8), intent(in) :: rho -!real(r8), parameter :: small_diff = 1.0e-8 real(r8), parameter :: small_diff = epsilon(1.0_r8) integer :: k, find_index(1) @@ -279,15 +291,15 @@ subroutine update_hybrid(weight, weight_sd, se2, ss2, so2, x, obs, rho) ! If the weight_sd not positive, keep everything the same if(weight_sd <= 0.0_r8) return +! Check for bad correlation +if(rho /= rho .or. rho <= 0.0_r8) return !1.0e-5_r8 + m = weight v = weight_sd**2 d2 = (obs - x)**2 ss = se2 - ss2 -!print *, 'd2: ', d2 -!print *, 'ss: ', ss - if (ss == 0.0_r8 .or. abs(ss) <= small_diff) then ! If the ensemble and static variances are very close ! Then just keep the current weighting. @@ -325,15 +337,11 @@ subroutine update_hybrid(weight, weight_sd, se2, ss2, so2, x, obs, rho) disc = R**2 - Q**3 -!print *, 'disc: ', disc -!print *, 'Q: ', Q, ' R: ', R - ! Find cubic roots if (disc < 0.0_r8) then ! 3 distict real roots theta = acos( R / sqrt(Q**3) ) / 3.0_r8 - !print *, 'theta: ', theta mulfac = - 2.0_r8 * sqrt(Q) addfac = - a / 3.0_r8 @@ -345,8 +353,6 @@ subroutine update_hybrid(weight, weight_sd, se2, ss2, so2, x, obs, rho) abssep = abs(sol - m) - !print *, 'sol: ', sol - ! Closest root find_index = minloc(abssep) new_weight = sol(find_index(1)) @@ -361,14 +367,10 @@ subroutine update_hybrid(weight, weight_sd, se2, ss2, so2, x, obs, rho) H = Q / G endif - !print *, 'G: ', G, 'H: ', H - new_weight = G + H - a / 3.0_r8 endif -!print *, 'new_weight = ', new_weight, ' d2 = ', d2, ' se2 = ', se2, ' ss2 = ', ss2, ' so2 = ', so2 - ! Make sure weight satisfies constraints weight = new_weight if(weight < weight_lower_bound) weight = weight_lower_bound diff --git a/assimilation_code/modules/utilities/mpi_utilities_mod.f90 b/assimilation_code/modules/utilities/mpi_utilities_mod.f90 index b604c82be4..57cc89f702 100644 --- a/assimilation_code/modules/utilities/mpi_utilities_mod.f90 +++ b/assimilation_code/modules/utilities/mpi_utilities_mod.f90 @@ -86,18 +86,18 @@ module mpi_utilities_mod ! this directory. It is a sed script that comments in and out the interface ! block below. Please leave the BLOCK comment lines unchanged. -! !!SYSTEM_BLOCK_EDIT START COMMENTED_OUT -! !#if .not. defined (__GFORTRAN__) .and. .not. defined(__NAG__) -! ! interface block for getting return code back from system() routine -! interface -! function system(string) -! character(len=*) :: string -! integer :: system -! end function system -! end interface -! ! end block -! !#endif -! !!SYSTEM_BLOCK_EDIT END COMMENTED_OUT + !!SYSTEM_BLOCK_EDIT START COMMENTED_IN + !#if .not. defined (__GFORTRAN__) .and. .not. defined(__NAG__) + ! interface block for getting return code back from system() routine + interface + function system(string) + character(len=*) :: string + integer :: system + end function system + end interface + ! end block + !#endif + !!SYSTEM_BLOCK_EDIT END COMMENTED_IN ! allow global sum to be computed for integers, r4, and r8s @@ -912,16 +912,16 @@ end function iam_task0 !> intent(in) here, but they call a routine which is intent(inout) so they !> must be the same here. -subroutine broadcast_send(from, array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5) +subroutine broadcast_send(from, array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) integer, intent(in) :: from ! arrays are really only intent(in) here, but must match array_broadcast() call. real(r8), intent(inout) :: array1(:) - real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) - real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) + real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 real(r8) :: packbuf(PACKLIMIT) -real(r8) :: local(5) +real(r8) :: local(6) logical :: doscalar, morethanone integer :: itemcount @@ -938,14 +938,14 @@ subroutine broadcast_send(from, array1, array2, array3, array4, array5, array6, endif ! for relatively small array sizes, pack them into a single send/recv pair. -call countup(array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5, & +call countup(array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6, & itemcount, morethanone, doscalar) if (itemcount <= PACKLIMIT .and. morethanone) then - call packit(packbuf, array1, array2, array3, array4, array5, array6, doscalar, & - scalar1, scalar2, scalar3, scalar4, scalar5) + call packit(packbuf, array1, array2, array3, array4, array5, array6, array7, doscalar, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) call array_broadcast(packbuf, from, itemcount) @@ -959,9 +959,10 @@ subroutine broadcast_send(from, array1, array2, array3, array4, array5, array6, if (present(array4)) call array_broadcast(array4, from) if (present(array5)) call array_broadcast(array5, from) if (present(array6)) call array_broadcast(array6, from) + if (present(array7)) call array_broadcast(array7, from) if (doscalar) then - call packscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5) + call packscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) call array_broadcast(local, from) endif @@ -984,16 +985,16 @@ end subroutine broadcast_send !> intent(out) here, but they call a routine which is intent(inout) so they !> must be the same here. -subroutine broadcast_recv(from, array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5) +subroutine broadcast_recv(from, array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) integer, intent(in) :: from ! arrays are really only intent(out) here, but must match array_broadcast() call. real(r8), intent(inout) :: array1(:) - real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) - real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) + real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 real(r8) :: packbuf(PACKLIMIT) -real(r8) :: local(5) +real(r8) :: local(6) logical :: doscalar, morethanone integer :: itemcount @@ -1010,16 +1011,16 @@ subroutine broadcast_recv(from, array1, array2, array3, array4, array5, array6, endif ! for relatively small array sizes, pack them into a single send/recv pair. -call countup(array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5, & +call countup(array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6, & itemcount, morethanone, doscalar) if (itemcount <= PACKLIMIT .and. morethanone) then call array_broadcast(packbuf, from, itemcount) - call unpackit(packbuf, array1, array2, array3, array4, array5, array6, doscalar, & - scalar1, scalar2, scalar3, scalar4, scalar5) + call unpackit(packbuf, array1, array2, array3, array4, array5, array6, array7, doscalar, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) else @@ -1030,12 +1031,13 @@ subroutine broadcast_recv(from, array1, array2, array3, array4, array5, array6, if (present(array3)) call array_broadcast(array3, from) if (present(array4)) call array_broadcast(array4, from) if (present(array5)) call array_broadcast(array5, from) - if (present(array6)) call array_broadcast(array6, from) + if (present(array6)) call array_broadcast(array6, from) + if (present(array7)) call array_broadcast(array7, from) if (doscalar) then call array_broadcast(local, from) call unpackscalar(local, scalar1, scalar2, scalar3, & - scalar4, scalar5) + scalar4, scalar5, scalar6) endif endif @@ -1050,12 +1052,12 @@ end subroutine broadcast_recv !> also note if there's more than a single array (array1) to send, !> and if there are any scalars specified. -subroutine countup(array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5, & +subroutine countup(array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6, & numitems, morethanone, doscalar) real(r8), intent(in) :: array1(:) - real(r8), intent(in), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) - real(r8), intent(in), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(in), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) + real(r8), intent(in), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 integer, intent(out) :: numitems logical, intent(out) :: morethanone, doscalar @@ -1082,6 +1084,10 @@ subroutine countup(array1, array2, array3, array4, array5, array6, & numitems = numitems + size(array6) morethanone = .true. endif +if (present(array7)) then + numitems = numitems + size(array7) + morethanone = .true. +endif if (present(scalar1)) then numitems = numitems + 1 morethanone = .true. @@ -1107,6 +1113,11 @@ subroutine countup(array1, array2, array3, array4, array5, array6, & morethanone = .true. doscalar = .true. endif +if (present(scalar6)) then + numitems = numitems + 1 + morethanone = .true. + doscalar = .true. +endif end subroutine countup @@ -1114,13 +1125,13 @@ end subroutine countup !> pack multiple small arrays into a single buffer before sending. -subroutine packit(buf, array1, array2, array3, array4, array5, array6, doscalar, & - scalar1, scalar2, scalar3, scalar4, scalar5) +subroutine packit(buf, array1, array2, array3, array4, array5, array6, array7, doscalar, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) real(r8), intent(out) :: buf(:) real(r8), intent(in) :: array1(:) - real(r8), intent(in), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) + real(r8), intent(in), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) logical, intent(in) :: doscalar - real(r8), intent(in), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(in), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 integer :: sindex, eindex @@ -1159,6 +1170,12 @@ subroutine packit(buf, array1, array2, array3, array4, array5, array6, doscalar, sindex = eindex+1 endif +if (present(array7)) then + eindex = sindex + size(array7) - 1 + buf(sindex:eindex) = array7(:) + sindex = eindex+1 +endif + if (doscalar) then if (present(scalar1)) then buf(sindex) = scalar1 @@ -1184,6 +1201,11 @@ subroutine packit(buf, array1, array2, array3, array4, array5, array6, doscalar, buf(sindex) = scalar5 sindex = sindex+1 endif + + if (present(scalar6)) then + buf(sindex) = scalar6 + sindex = sindex+1 + endif endif end subroutine packit @@ -1192,13 +1214,13 @@ end subroutine packit !> unpack multiple small arrays from a single buffer after receiving. -subroutine unpackit(buf, array1, array2, array3, array4, array5, array6, doscalar, & - scalar1, scalar2, scalar3, scalar4, scalar5) +subroutine unpackit(buf, array1, array2, array3, array4, array5, array6, array7, doscalar, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) real(r8), intent(in) :: buf(:) real(r8), intent(out) :: array1(:) - real(r8), intent(out), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) + real(r8), intent(out), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) logical, intent(in) :: doscalar - real(r8), intent(out), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(out), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 integer :: sindex, eindex @@ -1237,6 +1259,12 @@ subroutine unpackit(buf, array1, array2, array3, array4, array5, array6, doscala sindex = eindex+1 endif +if (present(array7)) then + eindex = sindex + size(array7) - 1 + array7(:) = buf(sindex:eindex) + sindex = eindex+1 +endif + if (doscalar) then if (present(scalar1)) then scalar1 = buf(sindex) @@ -1262,6 +1290,11 @@ subroutine unpackit(buf, array1, array2, array3, array4, array5, array6, doscala scalar5 = buf(sindex) sindex = sindex+1 endif + + if (present(scalar6)) then + scalar6 = buf(sindex) + sindex = sindex+1 + endif endif end subroutine unpackit @@ -1270,9 +1303,9 @@ end subroutine unpackit !> for any values specified, pack into a single array -subroutine packscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5) - real(r8), intent(out) :: local(5) - real(r8), intent(in), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 +subroutine packscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) + real(r8), intent(out) :: local(6) + real(r8), intent(in), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 local = 0.0_r8 @@ -1281,6 +1314,7 @@ subroutine packscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5) if (present(scalar3)) local(3) = scalar3 if (present(scalar4)) local(4) = scalar4 if (present(scalar5)) local(5) = scalar5 +if (present(scalar6)) local(6) = scalar6 end subroutine packscalar @@ -1288,15 +1322,16 @@ end subroutine packscalar !> for any values specified, unpack from a single array -subroutine unpackscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5) - real(r8), intent(in) :: local(5) - real(r8), intent(out), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 +subroutine unpackscalar(local, scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) + real(r8), intent(in) :: local(6) + real(r8), intent(out), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 if (present(scalar1)) scalar1 = local(1) if (present(scalar2)) scalar2 = local(2) if (present(scalar3)) scalar3 = local(3) if (present(scalar4)) scalar4 = local(4) if (present(scalar5)) scalar5 = local(5) +if (present(scalar6)) scalar6 = local(6) end subroutine unpackscalar diff --git a/assimilation_code/modules/utilities/null_mpi_utilities_mod.f90 b/assimilation_code/modules/utilities/null_mpi_utilities_mod.f90 index c64dd5cf17..647784e5ce 100644 --- a/assimilation_code/modules/utilities/null_mpi_utilities_mod.f90 +++ b/assimilation_code/modules/utilities/null_mpi_utilities_mod.f90 @@ -59,18 +59,18 @@ module mpi_utilities_mod ! this directory. It is a sed script that comments in and out the interface ! block below. Please leave the BLOCK comment lines unchanged. -! !!SYSTEM_BLOCK_EDIT START COMMENTED_OUT -! !#if .not. defined (__GFORTRAN__) .and. .not. defined(__NAG__) -! ! interface block for getting return code back from system() routine -! interface -! function system(string) -! character(len=*) :: string -! integer :: system -! end function system -! end interface -! ! end block -! !#endif -! !!SYSTEM_BLOCK_EDIT END COMMENTED_OUT + !!SYSTEM_BLOCK_EDIT START COMMENTED_IN + !#if .not. defined (__GFORTRAN__) .and. .not. defined(__NAG__) + ! interface block for getting return code back from system() routine + interface + function system(string) + character(len=*) :: string + integer :: system + end function system + end interface + ! end block + !#endif + !!SYSTEM_BLOCK_EDIT END COMMENTED_IN ! allow global sum to be computed for integers, r4, and r8s @@ -304,13 +304,13 @@ end function iam_task0 !> Returns with nothing to do. Does validate the 'from' task id. !> Not an error to call. -subroutine broadcast_send(from, array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5) +subroutine broadcast_send(from, array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) integer, intent(in) :: from ! arrays are really only intent(in) here, but must match array_broadcast() call. real(r8), intent(inout) :: array1(:) - real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) - real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) + real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 if ( .not. module_initialized ) call initialize_mpi_utilities() @@ -332,13 +332,13 @@ end subroutine broadcast_send !> Returns with nothing to do. Does validate the 'from' task id. !> Not an error to call. -subroutine broadcast_recv(from, array1, array2, array3, array4, array5, array6, & - scalar1, scalar2, scalar3, scalar4, scalar5) +subroutine broadcast_recv(from, array1, array2, array3, array4, array5, array6, array7, & + scalar1, scalar2, scalar3, scalar4, scalar5, scalar6) integer, intent(in) :: from ! arrays are really only intent(out) here, but must match array_broadcast() call. real(r8), intent(inout) :: array1(:) - real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:) - real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5 + real(r8), intent(inout), optional :: array2(:), array3(:), array4(:), array5(:), array6(:), array7(:) + real(r8), intent(inout), optional :: scalar1, scalar2, scalar3, scalar4, scalar5, scalar6 if ( .not. module_initialized ) call initialize_mpi_utilities() diff --git a/models/wrf_hydro/create_identity_streamflow_obs.f90 b/models/wrf_hydro/create_identity_streamflow_obs.f90 index edee0c5735..aeaa661996 100644 --- a/models/wrf_hydro/create_identity_streamflow_obs.f90 +++ b/models/wrf_hydro/create_identity_streamflow_obs.f90 @@ -308,7 +308,7 @@ program create_identity_streamflow_obs OBSLOOP: do n = 1, nobs - if ( discharge(n) < 0.0_r8 ) cycle OBSLOOP + if ( discharge(n) < 0.0_r8 .or. discharge(n) /= discharge(n) ) cycle OBSLOOP ! relate the TimeSlice:station to the RouteLink:gage so we can ! determine the location diff --git a/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/run_filter_experiment.py b/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/run_filter_experiment.py index c784a05d4e..8fde9e21a5 100644 --- a/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/run_filter_experiment.py +++ b/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/run_filter_experiment.py @@ -364,6 +364,21 @@ def manage_filter_output( post_mean_file.rename('input_postinf_mean{0}.nc'.format(domain_tag)) post_sd_file.rename('input_postinf_sd{0}.nc'.format(domain_tag)) + # MEG: Hybridization coefficient section + hybrd_mean_file = run_dir.joinpath('output_hybridweight_mean{0}.nc'.format(domain_tag)) + hybrd_sd_file = run_dir.joinpath('output_hybridweight_sd{0}.nc'.format(domain_tag)) + if hybrd_mean_file.exists(): + shutil.copy2( + str(hybrd_mean_file), + str(output_dir_date / ('output_hybridweight_mean{0}.{1}.nc'.format(domain_tag, datestr))) + ) + shutil.copy2( + str(hybrd_sd_file), + str(output_dir_date / ('output_hybridweight_sd{0}.{1}.nc'.format(domain_tag, datestr))) + ) + hybrd_mean_file.rename('input_hybridweight_mean{0}.nc'.format(domain_tag)) + hybrd_sd_file.rename('input_hybridweight_sd{0}.nc'.format(domain_tag)) + # Other output files. potential_output_file_globs = [ 'forecast*mean*nc', 'forecast*sd*nc', 'forecast_member_*.nc', diff --git a/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/setup_experiment.py b/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/setup_experiment.py index a2b63a2dd2..a6f8b56c47 100644 --- a/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/setup_experiment.py +++ b/models/wrf_hydro/hydro_dart_py/hydrodartpy/core/setup_experiment.py @@ -215,6 +215,11 @@ def setup_experiment(config_file): new_name = po.name.replace('output', 'input') _ = shutil.copy(po, config['experiment']['run_dir'] / new_name) + # Adding the stage of the hybrid weighting coefficient files + hybrid_wght_files = from_filter_path.glob('output_*hybridweight*') + for hyb in hybrid_wght_files: + new_name = hyb.name.replace('output', 'input') + _ = shutil.copy(hyb, config['experiment']['run_dir'] / new_name) # ################################### # Place scripts into the run dir. diff --git a/models/wrf_hydro/python/perturb/perturb_channel_bucket_only_forcing.py b/models/wrf_hydro/python/perturb/perturb_channel_bucket_only_forcing.py index 6681835cee..e5a4386792 100644 --- a/models/wrf_hydro/python/perturb/perturb_channel_bucket_only_forcing.py +++ b/models/wrf_hydro/python/perturb/perturb_channel_bucket_only_forcing.py @@ -49,6 +49,7 @@ def perturb_channel_only_forcing( paths = out_dir / chrtout_basename chrtout.to_netcdf(paths, mode='w') + chrtout.close() return paths diff --git a/models/wrf_hydro/shell_scripts/DART_cleanup.csh b/models/wrf_hydro/shell_scripts/DART_cleanup.csh index dc89ac7632..3c95de7c0a 100755 --- a/models/wrf_hydro/shell_scripts/DART_cleanup.csh +++ b/models/wrf_hydro/shell_scripts/DART_cleanup.csh @@ -41,10 +41,10 @@ # qpeek see output from a running job # #PBS -N DART_cleanup -#PBS -l walltime=00:10:00 +#PBS -l walltime=00:30:00 #PBS -q regular #PBS -l select=1:ncpus=36:mpiprocs=36 -#PBS -A NRAL0017 +#PBS -A P86850054 # #=============================================================================== # LSF directives bsub < script.csh @@ -58,15 +58,17 @@ #BSUB -o DART_cleanup.%J.log #BSUB -q regular #BSUB -n 16 -#BSUB -W 0:30:00 +#BSUB -W 0:50:00 #BSUB -P project # #=============================================================================== # USER SETTINGS HERE -set DARTdir = ${HOME}/WRF_Hydro/wrf_hydro_dart/models/wrf_hydro -set HydroDARTdir = /glade/scratch/${USER}/wrfhydro_dart/flo_cut/runs/da_ln1_n80_op2_lp4_sp4_noloc_ga_ana_sample_var +module swap openmpi mpt + +set DARTdir = /glade/work/${USER}/DART/DART_development/models/wrf_hydro +set HydroDARTdir = /glade/scratch/${USER}/wrfhydro_dart/Ian/runs/h1_N20_a3_R08_infprpo_m0.5_s0.01 set MARKER = all # END USER SETTINGS @@ -176,9 +178,9 @@ echo "DART time support ended at "`date` #=============================================================================== foreach STAGE ( input preassim analysis output ) - foreach TYPE (mean sd priorinf_mean priorinf_sd postinf_mean postinf_sd \ - mean_d01 sd_d01 priorinf_mean_d01 priorinf_sd_d01 postinf_mean_d01 postinf_sd_d01 \ - mean_d02 sd_d02 priorinf_mean_d02 priorinf_sd_d02 postinf_mean_d02 postinf_sd_d02 ) + foreach TYPE (mean sd priorinf_mean priorinf_sd postinf_mean postinf_sd hybridweight_mean hybridweight_sd \ + mean_d01 sd_d01 priorinf_mean_d01 priorinf_sd_d01 postinf_mean_d01 postinf_sd_d01 hybridweight_mean_d01 hybridweight_sd_d01 \ + mean_d02 sd_d02 priorinf_mean_d02 priorinf_sd_d02 postinf_mean_d02 postinf_sd_d02 hybridweight_mean_d02 hybridweight_sd_d02 ) set OUTPUT = ${MARKER}_${STAGE}_${TYPE}.nc ls -1 output/*/${STAGE}_${TYPE}.*.out.nc | sort >! concatlist.txt @@ -195,6 +197,8 @@ end echo "DART diagnostics summarized at "`date` +#exit 0 + #=============================================================================== # Concatenate all ensemble members together #===============================================================================