From 07e74ec51e144d91981c1686f146526f69facb97 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Wed, 27 Oct 2021 15:46:30 -0600 Subject: [PATCH 01/18] Moved computations for single state space inflation from filter_assim to subroutine update_single_state_space_inflation in adaptive_inflate_mod. Results are bitwise for both prior and posterior adaptive state space inflation general tests. --- .../assimilation/adaptive_inflate_mod.f90 | 60 +++++++++++++++++++ .../modules/assimilation/assim_tools_mod.f90 | 42 ++----------- 2 files changed, 64 insertions(+), 38 deletions(-) diff --git a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 index ae05dd597e..be1797566a 100644 --- a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 +++ b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 @@ -20,6 +20,7 @@ module adaptive_inflate_mod private public :: update_inflation, do_obs_inflate, & + update_single_state_space_inflation, & do_varying_ss_inflate, do_single_ss_inflate, inflate_ens, & adaptive_inflate_init, adaptive_inflate_type, & deterministic_inflate, solve_quadratic, & @@ -95,6 +96,10 @@ module adaptive_inflate_mod ! Flag indicating whether module has been initialized logical :: initialized = .false. +! Used for precision tests in inflation update routines +real(r8), parameter :: small = epsilon(1.0_r8) ! threshold for avoiding NaNs/Inf + + !=============================================================================== contains @@ -645,6 +650,61 @@ subroutine update_inflation(inflate_handle, inflate, inflate_sd, prior_mean, pri end subroutine update_inflation +!------------------------------------------------------------------------------- +!> Computes updatea inflation mean and inflation sd for single state space inflation + +subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd, & + ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, & + ens_size, inflate_only) + +! Computes update values for the inflation (inflate_mean) and its standard +! deviation (inflate_sd) + +type(adaptive_inflate_type), intent(in) :: inflate +real(r8), intent(inout) :: inflate_mean +real(r8), intent(inout) :: inflate_sd +real(r8), intent(in) :: ss_inflate_base +real(r8), intent(in) :: orig_obs_prior_mean +real(r8), intent(in) :: orig_obs_prior_var +real(r8), intent(in) :: obs +real(r8), intent(in) :: obs_err_var +integer, intent(in) :: ens_size +logical, intent(in) :: inflate_only + +real(r8) :: gamma, ens_var_deflate, r_var, r_mean + +! If either inflation or sd is not positive, not really doing inflation +if(inflate_mean > 0.0_r8 .and. inflate_sd > 0.0_r8) then + ! For case with single spatial inflation, use gamma = 1.0_r8 + gamma = 1.0_r8 + ! Deflate the inflated variance; required for efficient single pass + ! This is one of many places that assumes linear state/obs relation + ! over range of ensemble; Essentially, we are removing the inflation + ! which has already been applied in filter to see what inflation should + ! have been needed. + ens_var_deflate = orig_obs_prior_var / & + (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 + + ! If this is inflate_only (i.e. posterior) remove impact of this obs. + ! This is simulating independent observation by removing its impact. + if(inflate_only .and. & + ens_var_deflate > small .and. & + obs_err_var > small .and. & + obs_err_var - ens_var_deflate > small ) then + r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) + r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var) + else + r_var = ens_var_deflate + r_mean = orig_obs_prior_mean + endif + + ! Update the inflation mean value and standard deviation + call update_inflation(inflate, inflate_mean, inflate_sd, & + r_mean, r_var, ens_size, obs, obs_err_var, gamma) +endif + +end subroutine update_single_state_space_inflation + !------------------------------------------------------------------------------- !> Uses one of 2 algorithms in references on DART web site to update the diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 7198b40158..df551f983d 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -56,7 +56,7 @@ module assim_tools_mod use adaptive_inflate_mod, only : do_obs_inflate, do_single_ss_inflate, do_ss_inflate, & do_varying_ss_inflate, & - update_inflation, & + update_inflation, update_single_state_space_inflation, & inflate_ens, adaptive_inflate_type, & deterministic_inflate, solve_quadratic @@ -703,45 +703,11 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Compute updated values for single state space inflation SINGLE_SS_INFLATE: if(local_single_ss_inflate) then - ss_inflate_base = ens_handle%copies(ENS_SD_COPY, 1) ! Update for each group separately do group = 1, num_groups - ! If either inflation or sd is not positive, not really doing inflation - if(my_inflate > 0.0_r8 .and. my_inflate_sd > 0.0_r8) then - ! For case with single spatial inflation, use gamma = 1.0_r8 - ! See adaptive inflation module for details - gamma = 1.0_r8 - ! Deflate the inflated variance; required for efficient single pass - ! This is one of many places that assumes linear state/obs relation - ! over range of ensemble; Essentially, we are removing the inflation - ! which has already been applied in filter to see what inflation should - ! have been needed. - ens_obs_mean = orig_obs_prior_mean(group) - ens_obs_var = orig_obs_prior_var(group) - ! gamma is hardcoded as 1.0, so no test is needed here. - ens_var_deflate = ens_obs_var / & - (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 - - ! If this is inflate_only (i.e. posterior) remove impact of this obs. - ! This is simulating independent observation by removing its impact. - if(inflate_only .and. & - ens_var_deflate > small .and. & - obs_err_var > small .and. & - obs_err_var - ens_var_deflate > small ) then - r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) - r_mean = r_var *(ens_obs_mean / ens_var_deflate - obs(1) / obs_err_var) - else - r_var = ens_var_deflate - r_mean = ens_obs_mean - endif - - if (timing(SM_GRN)) call start_timer(t_base(SM_GRN), t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) - ! Update the inflation value - call update_inflation(inflate, my_inflate, my_inflate_sd, & - r_mean, r_var, grp_size, obs(1), obs_err_var, gamma) - if (timing(SM_GRN)) call read_timer(t_base(SM_GRN), 'update_inflation_C', & - t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) - endif + call update_single_state_space_inflation(inflate, my_inflate, my_inflate_sd, & + ens_handle%copies(ENS_SD_COPY, 1), orig_obs_prior_mean(group), & + orig_obs_prior_var(group), obs(1), obs_err_var, grp_size, inflate_only) end do endif SINGLE_SS_INFLATE From 51377346553643c8923cc54845d44705549269cd Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Thu, 28 Oct 2021 13:52:11 -0600 Subject: [PATCH 02/18] Moved the observation increments computation and the single state space inflation computations to the code done by all processes. This leads to redundant computation but should retain the overall run time. The obs_incs no longer are broadcast but there is a need for the prior ensemble mean and variance to be broadcast for single state space inflation. --- .../modules/assimilation/assim_tools_mod.f90 | 76 +++++++++---------- 1 file changed, 37 insertions(+), 39 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index df551f983d..46d3a2f229 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -308,6 +308,10 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & type(obs_sequence_type), intent(in) :: obs_seq integer, intent(in) :: keys(:) integer, intent(in) :: ens_size, num_groups, obs_val_index +! JLA: At present, this only needs to be inout because of the possible use of +! non-determinstic obs_space adaptive inflation that is not currently supported. +! Implementing that would require communication of the info about the inflation +! values as each observation updated them. type(adaptive_inflate_type), intent(inout) :: inflate integer, intent(in) :: ENS_MEAN_COPY, ENS_SD_COPY, ENS_INF_COPY integer, intent(in) :: ENS_INF_SD_COPY @@ -393,6 +397,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & integer, allocatable :: n_close_state_items(:), n_close_obs_items(:) +! Just to make sure multiple tasks are running for tests +write(*, *) 'my_task_id ', my_task_id() + ! timing disabled by default timing(:) = .false. t_base(:) = 0.0_r8 @@ -559,8 +566,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & end do if (timing(LG_GRN)) call read_timer(t_base(LG_GRN), 'get_state_meta_data') -!call test_get_state_meta_data(my_state_loc, ens_handle%my_num_vars) - !> optionally convert all state location verticals if (convert_all_state_verticals_first .and. is_doing_vertical_conversion) then if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) @@ -571,10 +576,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & if (timing(LG_GRN)) call read_timer(t_base(LG_GRN), 'convert_vertical_state') endif -! PAR: MIGHT BE BETTER TO HAVE ONE PE DEDICATED TO COMPUTING -! INCREMENTS. OWNING PE WOULD SHIP IT'S PRIOR TO THIS ONE -! BEFORE EACH INCREMENT. - ! Get mean and variance of each group's observation priors for adaptive inflation ! Important that these be from before any observations have been used if(local_ss_inflate) then @@ -682,56 +683,35 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif obs_qc = obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, owners_index) + ! Only value of 0 for DART QC field should be assimilated IF_QC_IS_OKAY: if(nint(obs_qc) ==0) then obs_prior = obs_ens_handle%copies(1:ens_size, owners_index) ! Compute the prior mean and variance for this observation + ! Note that these are before DA starts, so can be different from current obs_prior orig_obs_prior_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START: & OBS_PRIOR_MEAN_END, owners_index) orig_obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START: & OBS_PRIOR_VAR_END, owners_index) - ! Compute observation space increments for each group - do group = 1, num_groups - grp_bot = grp_beg(group) - grp_top = grp_end(group) - 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)) - end do - - ! Compute updated values for single state space inflation - SINGLE_SS_INFLATE: if(local_single_ss_inflate) then - ! Update for each group separately - do group = 1, num_groups - call update_single_state_space_inflation(inflate, my_inflate, my_inflate_sd, & - ens_handle%copies(ENS_SD_COPY, 1), orig_obs_prior_mean(group), & - orig_obs_prior_var(group), obs(1), obs_err_var, grp_size, inflate_only) - end do - endif SINGLE_SS_INFLATE - endif IF_QC_IS_OKAY !Broadcast the info from this obs to all other processes ! What gets broadcast depends on what kind of inflation is being done - !>@todo it should also depend on if vertical is being converted. the last - !>two values aren't needed unless vertical conversion is happening. - !>@todo FIXME: this is messy, but should we have 6 different broadcasts, - !>the three below and three more which omit the 2 localization values? - !>how much does this cost in time? time this and see. + !@todo it should also depend on if vertical is being converted. 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, & + call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, & orig_obs_prior_mean, orig_obs_prior_var, net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) - 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, & + call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, & + orig_obs_prior_mean, orig_obs_prior_var, & 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, & + call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, & net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) endif @@ -742,18 +722,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! I don't store this obs; receive the obs prior and increment from broadcast ! Also get qc and inflation information if needed ! also a converted vertical coordinate if needed - !>@todo FIXME see the comment in the broadcast_send() section about - !>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, & + call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & orig_obs_prior_mean, orig_obs_prior_var, net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) 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, & + call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & + orig_obs_prior_mean, orig_obs_prior_var, & 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, & + call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & net_a, scalar1=obs_qc, & scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) endif @@ -775,6 +754,25 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & !> because that's what we pass into the get_close_xxx() routines below. if (is_doing_vertical_conversion) & call set_vertical(base_obs_loc, vertvalue_obs_in_localization_coord, whichvert_obs_in_localization_coord) + + ! Compute observation space increments for each group + do group = 1, num_groups + grp_bot = grp_beg(group) + grp_top = grp_end(group) + 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)) + end do + + ! Compute updated values for single state space inflation + if(local_single_ss_inflate) then + ! Update for each group separately + do group = 1, num_groups + call update_single_state_space_inflation(inflate, my_inflate, my_inflate_sd, & + ens_handle%copies(ENS_SD_COPY, 1), orig_obs_prior_mean(group), & + orig_obs_prior_var(group), obs(1), obs_err_var, grp_size, inflate_only) + end do + endif ! Can compute prior mean and variance of obs for each group just once here do group = 1, num_groups From 81fde3becd6018e9cdeec737f7e473bf69cd11e9 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Sat, 30 Oct 2021 10:19:55 -0600 Subject: [PATCH 03/18] Removed all the timing code to try to make code clearer for mods. Combined loops computing group prior statistics. --- .../modules/assimilation/assim_tools_mod.f90 | 151 +----------------- 1 file changed, 7 insertions(+), 144 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 46d3a2f229..144e5ba969 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -382,52 +382,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & logical :: local_ss_inflate logical :: local_obs_inflate -! timing related vars: -! set timing(N) true to collect and print timing info -integer, parameter :: Ntimers = 5 -integer, parameter :: MLOOP = 1 ! main assimilation loop -integer, parameter :: LG_GRN = 2 ! large section timings -integer, parameter :: SM_GRN = 3 ! inner loops - use carefully! -integer, parameter :: GC = 4 ! get_close() related loops -logical :: timing(Ntimers) ! enable or disable w/ this -real(digits12) :: t_base(Ntimers) ! storage for time info -integer(i8) :: t_items(Ntimers) ! count of number of calls -integer(i8) :: t_limit(Ntimers) ! limit on number printed -real(digits12), allocatable :: elapse_array(:) - integer, allocatable :: n_close_state_items(:), n_close_obs_items(:) ! Just to make sure multiple tasks are running for tests write(*, *) 'my_task_id ', my_task_id() -! timing disabled by default -timing(:) = .false. -t_base(:) = 0.0_r8 -t_items(:) = 0_i8 -t_limit(:) = 0_i8 - ! how about this? look for imbalances in the tasks 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), & @@ -541,7 +504,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & 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)) if (obs_ens_handle%my_num_vars > 0) then call convert_vertical_obs(ens_handle, obs_ens_handle%my_num_vars, my_obs_loc, & my_obs_kind, my_obs_type, get_vertical_localization_coord(), vstatus) @@ -552,7 +514,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif enddo endif - if (timing(LG_GRN)) call read_timer(t_base(LG_GRN), 'convert_vertical_obs') endif ! Get info on my number and indices for state @@ -560,20 +521,16 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & call get_my_vars(ens_handle, my_state_indx) ! Get the location and kind of all my state variables -if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) do i = 1, ens_handle%my_num_vars call get_state_meta_data(my_state_indx(i), my_state_loc(i), my_state_kind(i)) end do -if (timing(LG_GRN)) call read_timer(t_base(LG_GRN), 'get_state_meta_data') !> optionally convert all state location verticals if (convert_all_state_verticals_first .and. is_doing_vertical_conversion) then - if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) if (ens_handle%my_num_vars > 0) then call convert_vertical_state(ens_handle, ens_handle%my_num_vars, my_state_loc, my_state_kind, & my_state_indx, get_vertical_localization_coord(), istatus) endif - if (timing(LG_GRN)) call read_timer(t_base(LG_GRN), 'convert_vertical_state') endif ! Get mean and variance of each group's observation priors for adaptive inflation @@ -623,18 +580,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & allow_missing_in_state = get_missing_ok_status() -! use MLOOP for the overall outer loop times; LG_GRN is for -! sections inside the overall loop, including the total time -! for the state_update and obs_update loops. use SM_GRN for -! sections inside those last 2 loops and be careful - they will -! be called nobs * nstate * ntasks. - ! Loop through all the (global) observations sequentially SEQUENTIAL_OBS: do i = 1, obs_ens_handle%num_vars - - if (timing(MLOOP)) call start_timer(t_base(MLOOP)) - if (timing(LG_GRN)) call start_timer(t_base(LG_GRN)) - ! Some compilers do not like mod by 0, so test first. if (print_every_nth_obs > 0) nth_obs = mod(i, print_every_nth_obs) @@ -687,14 +634,11 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Only value of 0 for DART QC field should be assimilated IF_QC_IS_OKAY: if(nint(obs_qc) ==0) then obs_prior = obs_ens_handle%copies(1:ens_size, owners_index) - - ! Compute the prior mean and variance for this observation ! Note that these are before DA starts, so can be different from current obs_prior orig_obs_prior_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START: & OBS_PRIOR_MEAN_END, owners_index) orig_obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START: & OBS_PRIOR_VAR_END, owners_index) - endif IF_QC_IS_OKAY !Broadcast the info from this obs to all other processes @@ -742,13 +686,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & !----------------------------------------------------------------------- ! Everybody is doing this section, cycle if qc is bad - if(nint(obs_qc) /= 0) then - if (timing(MLOOP)) then - write(msgstring, '(A32,I7)') 'sequential obs cycl: obs', keys(i) - call read_timer(t_base(MLOOP), msgstring, elapsed = elapse_array(i)) - endif - cycle SEQUENTIAL_OBS - endif + if(nint(obs_qc) /= 0) cycle SEQUENTIAL_OBS !> all tasks must set the converted vertical values into the 'base' version of this loc !> because that's what we pass into the get_close_xxx() routines below. @@ -762,6 +700,12 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & 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)) + + ! Also compute prior mean and variance of obs for efficiency here + obs_prior_mean(group) = sum(obs_prior(grp_bot:grp_top)) / grp_size + obs_prior_var(group) = sum((obs_prior(grp_bot:grp_top) - obs_prior_mean(group))**2) / & + (grp_size - 1) + if (obs_prior_var(group) < 0.0_r8) obs_prior_var(group) = 0.0_r8 end do ! Compute updated values for single state space inflation @@ -774,16 +718,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & end do endif - ! Can compute prior mean and variance of obs for each group just once here - do group = 1, num_groups - grp_bot = grp_beg(group) - grp_top = grp_end(group) - obs_prior_mean(group) = sum(obs_prior(grp_bot:grp_top)) / grp_size - obs_prior_var(group) = sum((obs_prior(grp_bot:grp_top) - obs_prior_mean(group))**2) / & - (grp_size - 1) - if (obs_prior_var(group) < 0.0_r8) obs_prior_var(group) = 0.0_r8 - end do - ! If we are doing adaptive localization then we need to know the number of ! other observations that are within the localization radius. We may need ! to shrink it, and so we need to know this before doing get_close() for the @@ -794,15 +728,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & if (.not. close_obs_caching) then - if (timing(GC)) call start_timer(t_base(GC), t_items(GC), t_limit(GC), do_sync=.false.) call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & my_obs_loc, my_obs_kind, my_obs_type, & num_close_obs, close_obs_ind, close_obs_dist, ens_handle) - if (timing(GC)) then - write(msgstring, '(A32,3I7)') 'gc_ob_NC:nobs,tot,obs# ', num_close_obs, obs_ens_handle%my_num_vars, keys(i) - call read_timer(t_base(GC), msgstring, t_items(GC), t_limit(GC), do_sync=.false.) - endif - else if (base_obs_loc == last_base_obs_loc) then @@ -811,14 +739,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & close_obs_dist(:) = last_close_obs_dist(:) num_close_obs_cached = num_close_obs_cached + 1 else - if (timing(GC)) call start_timer(t_base(GC), t_items(GC), t_limit(GC), do_sync=.false.) call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & my_obs_loc, my_obs_kind, my_obs_type, & num_close_obs, close_obs_ind, close_obs_dist, ens_handle) - if (timing(GC)) then - write(msgstring, '(A32,3I7)') 'gc_ob_C: nobs,tot,obs# ', num_close_obs, obs_ens_handle%my_num_vars, keys(i) - call read_timer(t_base(GC), msgstring, t_items(GC), t_limit(GC), do_sync=.false.) - endif last_base_obs_loc = base_obs_loc last_num_close_obs = num_close_obs @@ -829,9 +752,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif n_close_obs_items(i) = num_close_obs - !print*, 'base_obs _oc', base_obs_loc, 'rank ', my_task_id() - !call test_close_obs_dist(close_obs_dist, num_close_obs, i) - !print*, 'num close ', num_close_obs ! set the cutoff default, keep a copy of the original value, and avoid ! looking up the cutoff in a list if the incoming obs is an identity ob @@ -848,13 +768,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! For adaptive localization, need number of other obs close to the chosen observation if(adaptive_localization_threshold > 0) then - if (timing(GC)) call start_timer(t_base(GC), t_items(GC), t_limit(GC), do_sync=.false.) - ! this does a cross-task sum, so all tasks must make this call. total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & close_obs_dist, cutoff_rev*2.0_r8) - if (timing(GC)) call read_timer(t_base(GC), 'count_close', t_items(GC), t_limit(GC), do_sync=.false.) - ! Want expected number of close observations to be reduced to some threshold; ! accomplish this by cutting the size of the cutoff distance. @@ -921,14 +837,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Now everybody updates their close states ! Find state variables on my process that are close to observation being assimilated if (.not. close_obs_caching) then - if (timing(GC)) call start_timer(t_base(GC), t_items(GC), t_limit(GC), do_sync=.false.) call get_close_state(gc_state, base_obs_loc, base_obs_type, & my_state_loc, my_state_kind, my_state_indx, & num_close_states, close_state_ind, close_state_dist, ens_handle) - if (timing(GC)) then - write(msgstring, '(A32,3I7)') 'gc_st_NC:nsts,tot,obs# ', num_close_states, ens_handle%my_num_vars, keys(i) - call read_timer(t_base(GC), msgstring, t_items(GC), t_limit(GC), do_sync=.false.) - endif else if (base_obs_loc == last_base_states_loc) then num_close_states = last_num_close_states @@ -936,14 +847,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & close_state_dist(:) = last_close_state_dist(:) num_close_states_cached = num_close_states_cached + 1 else - if (timing(GC)) call start_timer(t_base(GC), t_items(GC), t_limit(GC), do_sync=.false.) call get_close_state(gc_state, base_obs_loc, base_obs_type, & my_state_loc, my_state_kind, my_state_indx, & num_close_states, close_state_ind, close_state_dist, ens_handle) - if (timing(GC)) then - write(msgstring, '(A32,3I7)') 'gc_st_C: nsts,tot,obs# ', num_close_states, ens_handle%my_num_vars, keys(i) - call read_timer(t_base(GC), msgstring, t_items(GC), t_limit(GC), do_sync=.false.) - endif last_base_states_loc = base_obs_loc last_num_close_states = num_close_states @@ -958,13 +864,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & !call test_close_obs_dist(close_state_dist, num_close_states, i) !call test_state_copies(ens_handle, 'beforeupdates') - if (timing(LG_GRN)) then - write(msgstring, '(A32,I7)') 'before_state_update: obs', keys(i) - call read_timer(t_base(LG_GRN), msgstring) - endif - ! 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) @@ -1001,7 +901,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! If no weight is indicated, no more to do with this state variable 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.) ! Loop through groups to update the state variable ensemble members do group = 1, num_groups grp_bot = grp_beg(group) @@ -1020,8 +919,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & increment(grp_bot:grp_top), reg_coef(group), net_a(group)) endif end do - if (timing(SM_GRN)) call read_timer(t_base(SM_GRN), 'update_from_obs_inc_S', & - t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) ! Compute an information factor for impact of this observation on this state if(num_groups == 1) then @@ -1078,11 +975,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! IS A TABLE LOOKUP POSSIBLE TO ACCELERATE THIS? ! Update the inflation values - if (timing(SM_GRN)) call start_timer(t_base(SM_GRN), t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) call update_inflation(inflate, varying_ss_inflate, varying_ss_inflate_sd, & r_mean, r_var, grp_size, obs(1), obs_err_var, gamma) - if (timing(SM_GRN)) call read_timer(t_base(SM_GRN), 'update_inflation_V', & - t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) else ! if we don't go into the previous if block, make sure these ! have good values going out for the block below @@ -1108,17 +1002,12 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif end do STATE_UPDATE - if (timing(LG_GRN)) then - write(msgstring, '(A32,I7)') 'state_update: obs', keys(i) - call read_timer(t_base(LG_GRN), msgstring) - 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) @@ -1144,7 +1033,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & 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 do group = 1, num_groups grp_bot = grp_beg(group) @@ -1154,8 +1042,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & 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 - if (timing(SM_GRN)) call read_timer(t_base(SM_GRN), 'update_from_obs_inc_O', & - t_items(SM_GRN), t_limit(SM_GRN), do_sync=.false.) ! FIXME: could we move the if test for inflate only to here? @@ -1179,17 +1065,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif endif end do OBS_UPDATE - if (timing(LG_GRN)) then - write(msgstring, '(A32,I7)') 'obs_update: obs', keys(i) - call read_timer(t_base(LG_GRN), msgstring) - endif !call test_state_copies(ens_handle, 'after_obs_updates') - if (timing(MLOOP)) then - write(msgstring, '(A32,I7)') 'sequential obs loop: obs', keys(i) - call read_timer(t_base(MLOOP), msgstring, elapsed = elapse_array(i)) - endif end do SEQUENTIAL_OBS ! Every pe needs to get the current my_inflate and my_inflate_sd back @@ -1203,21 +1081,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & call get_close_destroy(gc_state) call get_close_destroy(gc_obs) -! print some stats about the assimilation -! (if interesting, could print exactly which obs # was fastest and slowest) -if (my_task_id() == 0 .and. timing(MLOOP)) then - write(msgstring, *) 'average assim time: ', sum(elapse_array) / size(elapse_array) - call error_handler(E_MSG,'filter_assim:',msgstring) - - write(msgstring, *) 'minimum assim time: ', minval(elapse_array) - call error_handler(E_MSG,'filter_assim:',msgstring) - - write(msgstring, *) 'maximum assim time: ', maxval(elapse_array) - call error_handler(E_MSG,'filter_assim:',msgstring) -endif - -if (timing(MLOOP)) deallocate(elapse_array) - ! do some stats - being aware that unless we do a reduce() operation ! this is going to be per-task. so only print if something interesting ! shows up in the stats? maybe it would be worth a reduce() call here? From 54f8f735470a223d7603ba7412d6aed2e0ed895b Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Sat, 30 Oct 2021 15:30:59 -0600 Subject: [PATCH 04/18] Removed computation of inflation updates for spatially varying state space from assim_tools into adaptive_inflate. Reduced the broadcast options from 3 to 1. A small amount of undefined values may be broadcast and received for certain options but clarity seems more important for now. Removed an option for not computing the correlation in update_state to clarify code. Could turn correlation into an optional argument if this efficiency was really important (which I doubt). --- .../assimilation/adaptive_inflate_mod.f90 | 106 +++++++++++-- .../modules/assimilation/assim_tools_mod.f90 | 146 ++++-------------- 2 files changed, 119 insertions(+), 133 deletions(-) diff --git a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 index be1797566a..e20f3c1c3e 100644 --- a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 +++ b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 @@ -19,18 +19,18 @@ module adaptive_inflate_mod implicit none private -public :: update_inflation, do_obs_inflate, & - update_single_state_space_inflation, & - do_varying_ss_inflate, do_single_ss_inflate, inflate_ens, & - adaptive_inflate_init, adaptive_inflate_type, & - deterministic_inflate, solve_quadratic, & - log_inflation_info, get_minmax_task_zero, mean_from_restart, & - sd_from_restart, & - output_inf_restart, get_inflate_mean, get_inflate_sd, & - get_is_prior, get_is_posterior, do_ss_inflate, & - set_inflation_mean_copy, set_inflation_sd_copy, get_inflation_mean_copy, & +public :: update_inflation, do_obs_inflate, & + update_single_state_space_inflation, update_varying_state_space_inflation, & + do_varying_ss_inflate, do_single_ss_inflate, inflate_ens, & + adaptive_inflate_init, adaptive_inflate_type, & + deterministic_inflate, solve_quadratic, & + log_inflation_info, get_minmax_task_zero, mean_from_restart, & + sd_from_restart, & + output_inf_restart, get_inflate_mean, get_inflate_sd, & + get_is_prior, get_is_posterior, do_ss_inflate, & + set_inflation_mean_copy, set_inflation_sd_copy, get_inflation_mean_copy, & get_inflation_sd_copy, do_rtps_inflate, validate_inflate_options, & - print_inflation_restart_filename, & + print_inflation_restart_filename, & PRIOR_INF, POSTERIOR_INF, NO_INFLATION, OBS_INFLATION, VARYING_SS_INFLATION, & SINGLE_SS_INFLATION, RELAXATION_TO_PRIOR_SPREAD, ENHANCED_SS_INFLATION @@ -657,9 +657,6 @@ subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, & ens_size, inflate_only) -! Computes update values for the inflation (inflate_mean) and its standard -! deviation (inflate_sd) - type(adaptive_inflate_type), intent(in) :: inflate real(r8), intent(inout) :: inflate_mean real(r8), intent(inout) :: inflate_sd @@ -705,6 +702,87 @@ subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd end subroutine update_single_state_space_inflation +!------------------------------------------------------------------------------- +!> Computes updatea inflation mean and inflation sd for varying state space inflation + +subroutine update_varying_state_space_inflation(inflate, inflate_mean_inout, inflate_sd_inout, & + ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, & + ens_size, reg_factor, correl, inflate_only) + +type(adaptive_inflate_type), intent(in) :: inflate +real(r8), intent(inout) :: inflate_mean_inout +real(r8), intent(inout) :: inflate_sd_inout +real(r8), intent(in) :: ss_inflate_base +real(r8), intent(in) :: orig_obs_prior_mean +real(r8), intent(in) :: orig_obs_prior_var +real(r8), intent(in) :: obs +real(r8), intent(in) :: obs_err_var +integer, intent(in) :: ens_size +real(r8), intent(in) :: reg_factor +real(r8), intent(in) :: correl +logical, intent(in) :: inflate_only + +real(r8) :: inflate_mean, inflate_sd, gamma, ens_var_deflate, r_var, r_mean +real(r8) :: diff_sd, outlier_ratio +logical :: do_adapt_inf_update + +! Copy the inflate_mean and inflate_sd which can be changed by update_inflation +inflate_mean = inflate_mean_inout +inflate_sd = inflate_sd_inout + +if(inflate_mean > 0.0_r8 .and. inflate_sd > 0.0_r8) then + ! Gamma is less than 1 for varying ss, see adaptive inflate module + gamma = reg_factor * abs(correl) + + ! Remove the impact of inflation to allow efficient single pass with assim. + if ( abs(gamma) > small ) then + ens_var_deflate = orig_obs_prior_var / & + (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 + else + ens_var_deflate = orig_obs_prior_var + endif + + ! If this is inflate only (i.e. posterior) remove impact of this obs. + if(inflate_only .and. & + ens_var_deflate > small .and. & + obs_err_var > small .and. & + obs_err_var - ens_var_deflate > small ) then + r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) + r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var) + else + r_var = ens_var_deflate + r_mean = orig_obs_prior_mean + endif + + ! IS A TABLE LOOKUP POSSIBLE TO ACCELERATE THIS? + ! Update the inflation values + call update_inflation(inflate, inflate_mean, inflate_sd, & + r_mean, r_var, ens_size, obs, obs_err_var, gamma) +else + ! if we don't go into the previous if block, make sure these + ! have good values going out for the block below + r_mean = orig_obs_prior_mean + r_var = orig_obs_prior_var +endif + +! Update adaptive values if posterior outlier_ratio test doesn't fail. +! Match code in obs_space_diags() in filter.f90 +do_adapt_inf_update = .true. +if (inflate_only) then + diff_sd = sqrt(obs_err_var + r_var) + if (diff_sd > 0.0_r8) then + outlier_ratio = abs(obs - r_mean) / diff_sd +! JLA: ********** Magic number of 3 sd outlier here needs to be addressed ********* + do_adapt_inf_update = (outlier_ratio <= 3.0_r8) + endif +endif +if (do_adapt_inf_update) then + inflate_mean_inout = inflate_mean + inflate_sd_inout = inflate_sd +endif + +end subroutine update_varying_state_space_inflation + !------------------------------------------------------------------------------- !> Uses one of 2 algorithms in references on DART web site to update the diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 144e5ba969..e4abdf834d 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -57,6 +57,7 @@ module assim_tools_mod use adaptive_inflate_mod, only : do_obs_inflate, do_single_ss_inflate, do_ss_inflate, & do_varying_ss_inflate, & update_inflation, update_single_state_space_inflation, & + update_varying_state_space_inflation, & inflate_ens, adaptive_inflate_type, & deterministic_inflate, solve_quadratic @@ -326,13 +327,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & real(r8) :: reg_factor, impact_factor real(r8) :: net_a(num_groups), reg_coef(num_groups), correl(num_groups) real(r8) :: cov_factor, obs(1), obs_err_var, my_inflate, my_inflate_sd -real(r8) :: varying_ss_inflate, varying_ss_inflate_sd -real(r8) :: ss_inflate_base, obs_qc, cutoff_rev, cutoff_orig -real(r8) :: gamma, ens_obs_mean, ens_obs_var, ens_var_deflate -real(r8) :: r_mean, r_var +real(r8) :: obs_qc, cutoff_rev, cutoff_orig real(r8) :: orig_obs_prior_mean(num_groups), orig_obs_prior_var(num_groups) real(r8) :: obs_prior_mean(num_groups), obs_prior_var(num_groups) -real(r8) :: diff_sd, outlier_ratio real(r8) :: vertvalue_obs_in_localization_coord, whichvert_real real(r8), allocatable :: close_obs_dist(:) real(r8), allocatable :: close_state_dist(:) @@ -375,7 +372,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & type(obs_def_type) :: obs_def type(time_type) :: obs_time, this_obs_time -logical :: do_adapt_inf_update logical :: allow_missing_in_state logical :: local_single_ss_inflate logical :: local_varying_ss_inflate @@ -642,44 +638,23 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif IF_QC_IS_OKAY !Broadcast the info from this obs to all other processes - ! What gets broadcast depends on what kind of inflation is being done - !@todo it should also depend on if vertical is being converted. + ! orig_obs_prior_mean and orig_obs_prior_var only used with adaptive inflation + ! my_inflate and my_inflate_sd only used with single state space inflation + ! vertvalue_obs_in_localization_coord and whichvert_real only used for vertical + ! coordinate transformation 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, & - orig_obs_prior_mean, orig_obs_prior_var, net_a, scalar1=obs_qc, & - scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) - else if(local_single_ss_inflate .or. local_obs_inflate) then - call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, & - orig_obs_prior_mean, orig_obs_prior_var, & - 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, & - net_a, scalar1=obs_qc, & - scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) - endif + call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, & + net_a, orig_obs_prior_mean, orig_obs_prior_var, & + scalar1=obs_qc, scalar2=vertvalue_obs_in_localization_coord, & + scalar3=whichvert_real, scalar4=my_inflate, scalar5=my_inflate_sd) ! Next block is done by processes that do NOT own this observation !----------------------------------------------------------------------- else - ! I don't store this obs; receive the obs prior and increment from broadcast - ! Also get qc and inflation information if needed - ! also a converted vertical coordinate if needed - if(local_varying_ss_inflate) then - call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & - orig_obs_prior_mean, orig_obs_prior_var, net_a, scalar1=obs_qc, & - scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) - else if(local_single_ss_inflate .or. local_obs_inflate) then - call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & - orig_obs_prior_mean, orig_obs_prior_var, & - 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, & - net_a, scalar1=obs_qc, & - scalar2=vertvalue_obs_in_localization_coord, scalar3=whichvert_real) - endif + call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & + net_a, orig_obs_prior_mean, orig_obs_prior_var, & + scalar1=obs_qc, scalar2=vertvalue_obs_in_localization_coord, & + scalar3=whichvert_real, scalar4=my_inflate, scalar5=my_inflate_sd) whichvert_obs_in_localization_coord = nint(whichvert_real) endif @@ -876,15 +851,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE 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) - varying_ss_inflate_sd = ens_handle%copies(ENS_INF_SD_COPY, state_index) - else - 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)) @@ -905,19 +871,11 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & do group = 1, num_groups grp_bot = grp_beg(group) grp_top = grp_end(group) - ! Do update of state, correl only needed for varying ss inflate - if(local_varying_ss_inflate .and. varying_ss_inflate > 0.0_r8 .and. & - varying_ss_inflate_sd > 0.0_r8) then - call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & - obs_prior_var(group), obs_inc(grp_bot:grp_top), & - ens_handle%copies(grp_bot:grp_top, state_index), grp_size, & - increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) - else - call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & - obs_prior_var(group), obs_inc(grp_bot:grp_top), & - ens_handle%copies(grp_bot:grp_top, state_index), grp_size, & - increment(grp_bot:grp_top), reg_coef(group), net_a(group)) - endif + ! Do update of state, correl only needed for varying ss inflate but compute for all + call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & + obs_prior_var(group), obs_inc(grp_bot:grp_top), & + ens_handle%copies(grp_bot:grp_top, state_index), grp_size, & + increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) end do ! Compute an information factor for impact of this observation on this state @@ -941,64 +899,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Compute spatially-varying state space inflation if(local_varying_ss_inflate) then - ! base is the initial inflate value for this state variable - ss_inflate_base = ens_handle%copies(ENS_SD_COPY, state_index) - ! Loop through each group to update inflation estimate - GroupInflate: do group = 1, num_groups - if(varying_ss_inflate > 0.0_r8 .and. varying_ss_inflate_sd > 0.0_r8) then - ! Gamma is less than 1 for varying ss, see adaptive inflate module - gamma = reg_factor * abs(correl(group)) - ! Deflate the inflated variance using the INITIAL state inflate - ! value (before these obs started gumming it up). - ens_obs_mean = orig_obs_prior_mean(group) - ens_obs_var = orig_obs_prior_var(group) - - ! Remove the impact of inflation to allow efficient single pass with assim. - if ( abs(gamma) > small ) then - ens_var_deflate = ens_obs_var / & - (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 - else - ens_var_deflate = ens_obs_var - endif - - ! If this is inflate only (i.e. posterior) remove impact of this obs. - if(inflate_only .and. & - ens_var_deflate > small .and. & - obs_err_var > small .and. & - obs_err_var - ens_var_deflate > small ) then - r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) - r_mean = r_var *(ens_obs_mean / ens_var_deflate - obs(1) / obs_err_var) - else - r_var = ens_var_deflate - r_mean = ens_obs_mean - endif - - ! IS A TABLE LOOKUP POSSIBLE TO ACCELERATE THIS? - ! Update the inflation values - call update_inflation(inflate, varying_ss_inflate, varying_ss_inflate_sd, & - r_mean, r_var, grp_size, obs(1), obs_err_var, gamma) - else - ! if we don't go into the previous if block, make sure these - ! have good values going out for the block below - r_mean = orig_obs_prior_mean(group) - r_var = orig_obs_prior_var(group) - endif - - ! Update adaptive values if posterior outlier_ratio test doesn't fail. - ! Match code in obs_space_diags() in filter.f90 - do_adapt_inf_update = .true. - if (inflate_only) then - diff_sd = sqrt(obs_err_var + r_var) - if (diff_sd > 0.0_r8) then - outlier_ratio = abs(obs(1) - r_mean) / diff_sd - do_adapt_inf_update = (outlier_ratio <= 3.0_r8) - endif - endif - if (do_adapt_inf_update) then - ens_handle%copies(ENS_INF_COPY, state_index) = varying_ss_inflate - ens_handle%copies(ENS_INF_SD_COPY, state_index) = varying_ss_inflate_sd - endif - end do GroupInflate + do group = 1, num_groups + call update_varying_state_space_inflation(inflate, & + ens_handle%copies(ENS_INF_COPY, state_index), & + ens_handle%copies(ENS_INF_SD_COPY, state_index), & + ens_handle%copies(ENS_SD_COPY, state_index), & + orig_obs_prior_mean(group), orig_obs_prior_var(group), obs(1), & + obs_err_var, grp_size, reg_factor, correl(group), inflate_only) + end do endif end do STATE_UPDATE From 4cd06cffc28cb3fd7504a0c16fc27c8b9a09a21c Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Mon, 1 Nov 2021 10:29:28 -0600 Subject: [PATCH 05/18] Compacting the adaptive localization code prior to pulling it out in a subroutine. localization_diagnostic files are identical to baseline. --- .../modules/assimilation/assim_tools_mod.f90 | 62 +++++-------------- 1 file changed, 17 insertions(+), 45 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index e4abdf834d..eba2c876a2 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -742,7 +742,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! For adaptive localization, need number of other obs close to the chosen observation if(adaptive_localization_threshold > 0) then - ! this does a cross-task sum, so all tasks must make this call. total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & close_obs_dist, cutoff_rev*2.0_r8) @@ -750,62 +749,35 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Want expected number of close observations to be reduced to some threshold; ! accomplish this by cutting the size of the cutoff distance. if(total_num_close_obs > adaptive_localization_threshold) then - cutoff_rev = revised_distance(cutoff_rev*2.0_r8, adaptive_localization_threshold, & total_num_close_obs, base_obs_loc, & adaptive_cutoff_floor*2.0_r8) / 2.0_r8 - - if ( output_localization_diagnostics ) then - - ! to really know how many obs are left now, you have to - ! loop over all the obs, again, count how many kinds are - ! going to be assim, and explicitly check the distance and - ! see if it's closer than the new cutoff ( times 2 ), and - ! then do a global sum to get the total. since this costs, - ! do it only when diagnostics are requested. - - ! this does a cross-task sum, so all tasks must make this call. - rev_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & - close_obs_dist, cutoff_rev*2.0_r8) - - - ! GSR output the new cutoff - ! Here is what we might want: - ! time, ob index #, ob location, new cutoff, the assimilate obs count, owner (which process has this ob) - ! obs_time, obs_val_index, base_obs_loc, cutoff_rev, total_num_close_obs, owner - ! break up the time into secs and days, and break up the location into lat, lon and height - ! nsc - the min info here that can't be extracted from the obs key is: - ! key (obs#), total_num_close_obs (close w/ original cutoff), revised cutoff & new count - if (my_task_id() == 0) then - call get_obs_def(observation, obs_def) - this_obs_time = get_obs_def_time(obs_def) - call get_time(this_obs_time,secs,days) - call write_location(-1, base_obs_loc, charstring=base_loc_text) - - write(localization_unit,'(i12,1x,i5,1x,i8,1x,A,2(f14.5,1x,i12))') i, secs, days, & - trim(base_loc_text), cutoff_orig, total_num_close_obs, cutoff_rev, rev_num_close_obs - endif - endif - endif + endif - else if (output_localization_diagnostics) then - - ! if you aren't adapting but you still want to know how many obs are within the - ! localization radius, set the diag output. this could be large, use carefully. - - ! this does a cross-task sum, so all tasks must make this call. - total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & + if ( output_localization_diagnostics ) then + ! Warning, this can be costly and generate large output + ! This is referred to as revised in case adaptive localization was done + rev_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & close_obs_dist, cutoff_rev*2.0_r8) + ! Output diagnostic information about the number of close obs if (my_task_id() == 0) then - call get_obs_def(observation, obs_def) this_obs_time = get_obs_def_time(obs_def) call get_time(this_obs_time,secs,days) call write_location(-1, base_obs_loc, charstring=base_loc_text) - write(localization_unit,'(i12,1x,i5,1x,i8,1x,A,f14.5,1x,i12)') i, secs, days, & - trim(base_loc_text), cutoff_rev, total_num_close_obs + ! If adaptive localization did something, output info about what it did + ! Probably would be more consistent to just output for all observations + if(adaptive_localization_threshold > 0 .and. & + total_num_close_obs > adaptive_localization_threshold) then + write(localization_unit,'(i12,1x,i5,1x,i8,1x,A,2(f14.5,1x,i12))') i, secs, days, & + trim(base_loc_text), cutoff_orig, total_num_close_obs, cutoff_rev, rev_num_close_obs + else + ! Full diagnostics + write(localization_unit,'(i12,1x,i5,1x,i8,1x,A,f14.5,1x,i12)') i, secs, days, & + trim(base_loc_text), cutoff_rev, rev_num_close_obs + endif endif endif From 964d1a798e200f3cf5ded7a5835a7ae42bc2002d Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Mon, 1 Nov 2021 12:02:38 -0600 Subject: [PATCH 06/18] Moved the adaptive localization and localization output to a subroutine and removed the previous code. Believe to bitwise reproduce. --- .../modules/assimilation/assim_tools_mod.f90 | 145 ++++++++++-------- 1 file changed, 81 insertions(+), 64 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index eba2c876a2..f341ab6689 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -344,13 +344,12 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & integer :: obs_mean_index, obs_var_index integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group integer :: num_close_obs, obs_index, num_close_states -integer :: total_num_close_obs, last_num_close_obs, last_num_close_states +integer :: last_num_close_obs, last_num_close_states integer :: base_obs_kind, base_obs_type, nth_obs integer :: num_close_obs_cached, num_close_states_cached integer :: num_close_obs_calls_made, num_close_states_calls_made -integer :: localization_unit, secs, days, rev_num_close_obs integer :: whichvert_obs_in_localization_coord -integer :: istatus +integer :: istatus, localization_unit integer, allocatable :: close_obs_ind(:) integer, allocatable :: close_state_ind(:) integer, allocatable :: last_close_obs_ind(:) @@ -360,8 +359,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & integer, allocatable :: my_state_kind(:) integer, allocatable :: vstatus(:) -character(len = 200) :: base_loc_text ! longer than longest location formatting possible - type(location_type) :: base_obs_loc, last_base_obs_loc, last_base_states_loc type(location_type) :: dummyloc type(location_type), allocatable :: my_obs_loc(:) @@ -370,7 +367,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & type(get_close_type) :: gc_obs, gc_state type(obs_type) :: observation type(obs_def_type) :: obs_def -type(time_type) :: obs_time, this_obs_time +type(time_type) :: obs_time logical :: allow_missing_in_state logical :: local_single_ss_inflate @@ -435,10 +432,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & endif endif -!GSR open the dignostics file -if(output_localization_diagnostics .and. my_task_id() == 0) then +! Open the localization diagnostics file +if(output_localization_diagnostics .and. my_task_id() == 0) & localization_unit = open_file(localization_diagnostics_file, action = 'append') -endif ! For performance, make local copies of these settings which ! are really in the inflate derived type. @@ -693,15 +689,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & end do endif - ! If we are doing adaptive localization then we need to know the number of - ! other observations that are within the localization radius. We may need - ! to shrink it, and so we need to know this before doing get_close() for the - ! state space (even though the state space increments will be computed and - ! applied first). - - !****************************************** - - if (.not. close_obs_caching) then call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & my_obs_loc, my_obs_kind, my_obs_type, & @@ -738,48 +725,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & cutoff_orig = cutoff endif - cutoff_rev = cutoff_orig - - ! For adaptive localization, need number of other obs close to the chosen observation - if(adaptive_localization_threshold > 0) then - ! this does a cross-task sum, so all tasks must make this call. - total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & - close_obs_dist, cutoff_rev*2.0_r8) - - ! Want expected number of close observations to be reduced to some threshold; - ! accomplish this by cutting the size of the cutoff distance. - if(total_num_close_obs > adaptive_localization_threshold) then - cutoff_rev = revised_distance(cutoff_rev*2.0_r8, adaptive_localization_threshold, & - total_num_close_obs, base_obs_loc, & - adaptive_cutoff_floor*2.0_r8) / 2.0_r8 - endif - endif - - if ( output_localization_diagnostics ) then - ! Warning, this can be costly and generate large output - ! This is referred to as revised in case adaptive localization was done - rev_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & - close_obs_dist, cutoff_rev*2.0_r8) - - ! Output diagnostic information about the number of close obs - if (my_task_id() == 0) then - this_obs_time = get_obs_def_time(obs_def) - call get_time(this_obs_time,secs,days) - call write_location(-1, base_obs_loc, charstring=base_loc_text) - - ! If adaptive localization did something, output info about what it did - ! Probably would be more consistent to just output for all observations - if(adaptive_localization_threshold > 0 .and. & - total_num_close_obs > adaptive_localization_threshold) then - write(localization_unit,'(i12,1x,i5,1x,i8,1x,A,2(f14.5,1x,i12))') i, secs, days, & - trim(base_loc_text), cutoff_orig, total_num_close_obs, cutoff_rev, rev_num_close_obs - else - ! Full diagnostics - write(localization_unit,'(i12,1x,i5,1x,i8,1x,A,f14.5,1x,i12)') i, secs, days, & - trim(base_loc_text), cutoff_rev, rev_num_close_obs - endif - endif - endif + call adaptive_localization_and_diags(cutoff_orig, cutoff_rev, adaptive_localization_threshold, & + adaptive_cutoff_floor, num_close_obs, close_obs_ind, close_obs_dist, my_obs_type, & + i, base_obs_loc, obs_def, localization_unit) ! Now everybody updates their close states ! Find state variables on my process that are close to observation being assimilated @@ -1007,10 +955,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & !call test_state_copies(ens_handle, 'end') -!GSR close the localization diagnostics file -if(output_localization_diagnostics .and. my_task_id() == 0) then - call close_file(localization_unit) -end if +! Close the localization diagnostics file +if(output_localization_diagnostics .and. my_task_id() == 0) call close_file(localization_unit) ! get rid of mpi window call free_mean_window() @@ -2572,6 +2518,77 @@ function count_close(num_close, index_list, my_types, dist, maxdist) end function count_close +!---------------------------------------------------------------------- +! Revise the cutoff for this observation if adaptive localization is required +! Output diagnostics for localization if requested + +subroutine adaptive_localization_and_diags(cutoff_orig, cutoff_rev, adaptive_localization_threshold, & + adaptive_cutoff_floor, num_close_obs, close_obs_ind, close_obs_dist, my_obs_type, & + base_obs_index, base_obs_loc, obs_def, out_unit) + +real(r8), intent(in) :: cutoff_orig +real(r8), intent(out) :: cutoff_rev +integer, intent(in) :: adaptive_localization_threshold +real(r8), intent(in) :: adaptive_cutoff_floor +integer, intent(in) :: num_close_obs +integer, intent(in) :: close_obs_ind(:) +real(r8), intent(in) :: close_obs_dist(:) +integer, intent(in) :: my_obs_type(:) +integer, intent(in) :: base_obs_index +type(location_type), intent(in) :: base_obs_loc +type(obs_def_type), intent(in) :: obs_def +integer, intent(in) :: out_unit + +integer :: total_num_close_obs, rev_num_close_obs, secs, days +type(time_type) :: this_obs_time +character(len = 200) :: base_loc_text ! longer than longest location formatting possible + +! Default is that cutoff is not revised +cutoff_rev = cutoff_orig + +! For adaptive localization, need number of other obs close to the chosen observation +if(adaptive_localization_threshold > 0) then + ! this does a cross-task sum, so all tasks must make this call. + total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & + close_obs_dist, cutoff_rev*2.0_r8) + + ! Want expected number of close observations to be reduced to some threshold; + ! accomplish this by cutting the size of the cutoff distance. + if(total_num_close_obs > adaptive_localization_threshold) then + cutoff_rev = revised_distance(cutoff_rev*2.0_r8, adaptive_localization_threshold, & + total_num_close_obs, base_obs_loc, & + adaptive_cutoff_floor*2.0_r8) / 2.0_r8 + endif +endif + +if ( output_localization_diagnostics ) then + ! Warning, this can be costly and generate large output + ! This is referred to as revised in case adaptive localization was done + rev_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, & + close_obs_dist, cutoff_rev*2.0_r8) + + ! Output diagnostic information about the number of close obs + if (my_task_id() == 0) then + this_obs_time = get_obs_def_time(obs_def) + call get_time(this_obs_time,secs,days) + call write_location(-1, base_obs_loc, charstring=base_loc_text) + + ! If adaptive localization did something, output info about what it did + ! Probably would be more consistent to just output for all observations + if(adaptive_localization_threshold > 0 .and. & + total_num_close_obs > adaptive_localization_threshold) then + write(out_unit,'(i12,1x,i5,1x,i8,1x,A,2(f14.5,1x,i12))') base_obs_index, & + secs, days, trim(base_loc_text), cutoff_orig, total_num_close_obs, cutoff_rev, & + rev_num_close_obs + else + write(out_unit,'(i12,1x,i5,1x,i8,1x,A,f14.5,1x,i12)') base_obs_index, & + secs, days, trim(base_loc_text), cutoff_rev, rev_num_close_obs + endif + endif +endif + +end subroutine adaptive_localization_and_diags + !---------------------------------------------------------------------- !> gets the location of of all my observations subroutine get_my_obs_loc(obs_ens_handle, obs_seq, keys, my_obs_loc, my_obs_kind, my_obs_type, my_obs_time) From 5a286faa47392b3694b3b8875f268c830de8adec Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Tue, 2 Nov 2021 16:18:01 -0600 Subject: [PATCH 07/18] Moved caching code for obs and state get_close into subroutines to reduce main loop clutter. Removed some timeing routines that are no longer used. --- .../modules/assimilation/assim_tools_mod.f90 | 248 +++++++----------- 1 file changed, 102 insertions(+), 146 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index f341ab6689..d9ece542e1 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -77,7 +77,7 @@ module assim_tools_mod public :: filter_assim, & set_assim_tools_trace, & test_state_copies, & - update_ens_from_weights ! Jeff thinks this routine is in the wild. + update_ens_from_weights ! Indicates if module initialization subroutine has been called yet logical :: module_initialized = .false. @@ -688,31 +688,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & orig_obs_prior_var(group), obs(1), obs_err_var, grp_size, inflate_only) end do endif - - if (.not. close_obs_caching) then - call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & - my_obs_loc, my_obs_kind, my_obs_type, & - num_close_obs, close_obs_ind, close_obs_dist, ens_handle) - else - - if (base_obs_loc == last_base_obs_loc) then - num_close_obs = last_num_close_obs - close_obs_ind(:) = last_close_obs_ind(:) - close_obs_dist(:) = last_close_obs_dist(:) - num_close_obs_cached = num_close_obs_cached + 1 - else - call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & - my_obs_loc, my_obs_kind, my_obs_type, & - num_close_obs, close_obs_ind, close_obs_dist, ens_handle) - - last_base_obs_loc = base_obs_loc - last_num_close_obs = num_close_obs - last_close_obs_ind(:) = close_obs_ind(:) - last_close_obs_dist(:) = close_obs_dist(:) - num_close_obs_calls_made = num_close_obs_calls_made +1 - endif - endif - + + ! Adaptive localization needs number of other observations within localization radius. + ! Do get_close_obs first, even though state space increments are computed before obs increments. + ! JLA: ens_handle doesn't ever appear to be used. Get rid of it. Should be obs_ens_handle anyway? + call get_close_obs_cached(close_obs_caching, gc_obs, base_obs_loc, base_obs_type, & + my_obs_loc, my_obs_kind, my_obs_type, num_close_obs, close_obs_ind, close_obs_dist, & + ens_handle, last_base_obs_loc, last_num_close_obs, last_close_obs_ind, & + last_close_obs_dist, num_close_obs_cached, num_close_obs_calls_made) n_close_obs_items(i) = num_close_obs ! set the cutoff default, keep a copy of the original value, and avoid @@ -725,35 +708,16 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & cutoff_orig = cutoff endif + ! JLA, could also cache for adaptive_localization which may be expensive? call adaptive_localization_and_diags(cutoff_orig, cutoff_rev, adaptive_localization_threshold, & adaptive_cutoff_floor, num_close_obs, close_obs_ind, close_obs_dist, my_obs_type, & i, base_obs_loc, obs_def, localization_unit) - ! Now everybody updates their close states ! Find state variables on my process that are close to observation being assimilated - if (.not. close_obs_caching) then - call get_close_state(gc_state, base_obs_loc, base_obs_type, & - my_state_loc, my_state_kind, my_state_indx, & - num_close_states, close_state_ind, close_state_dist, ens_handle) - else - if (base_obs_loc == last_base_states_loc) then - num_close_states = last_num_close_states - close_state_ind(:) = last_close_state_ind(:) - close_state_dist(:) = last_close_state_dist(:) - num_close_states_cached = num_close_states_cached + 1 - else - call get_close_state(gc_state, base_obs_loc, base_obs_type, & - my_state_loc, my_state_kind, my_state_indx, & - num_close_states, close_state_ind, close_state_dist, ens_handle) - - last_base_states_loc = base_obs_loc - last_num_close_states = num_close_states - last_close_state_ind(:) = close_state_ind(:) - last_close_state_dist(:) = close_state_dist(:) - num_close_states_calls_made = num_close_states_calls_made + 1 - endif - endif - + call get_close_state_cached(close_obs_caching, gc_state, base_obs_loc, base_obs_type, & + my_state_loc, my_state_kind, my_state_indx, num_close_states, close_state_ind, close_state_dist, & + ens_handle, last_base_states_loc, last_num_close_states, last_close_state_ind, & + last_close_state_dist, num_close_states_cached, num_close_states_calls_made) n_close_state_items(i) = num_close_states !print*, 'num close state', num_close_states !call test_close_obs_dist(close_state_dist, num_close_states, i) @@ -2626,109 +2590,101 @@ subroutine get_my_obs_loc(obs_ens_handle, obs_seq, keys, my_obs_loc, my_obs_kind end subroutine get_my_obs_loc !-------------------------------------------------------------------- -!> wrappers for timers -!> -!> t_space is where we store the time information -!> itemcount is the running count of how many times we've been called -!> maxitems is a limit on the number of times we want this timer to print. -!> do_sync overrides the default for whether we want to do a task sync -!> or not. right now this code defaults to yes, sync before getting -!> the time. for very large processor counts this increases overhead. - -subroutine start_timer(t_space, itemcount, maxitems, do_sync) - real(digits12), intent(out) :: t_space - integer(i8), intent(inout), optional :: itemcount - integer(i8), intent(in), optional :: maxitems - logical, intent(in), optional :: do_sync - -logical :: sync_me - -if (present(itemcount) .and. present(maxitems)) then - itemcount = itemcount + 1 - if (itemcount > maxitems) then - ! if called enough, this can roll over the integer limit. - ! set itemcount to maxitems+1 here to avoid this. - ! also, go ahead and set the time because there is - ! an option to print large time values even if over the - ! number of calls limit in read_timer() - - itemcount = maxitems + 1 - call start_mpi_timer(t_space) - return - endif -endif - -sync_me = .true. -if (present(do_sync)) sync_me = do_sync - -if (sync_me) call task_sync() -call start_mpi_timer(t_space) - -end subroutine start_timer +!> Get close obs from cache if appropriate. Cache new get_close_obs info +!> if requested. + +subroutine get_close_obs_cached(close_obs_caching, gc_obs, base_obs_loc, base_obs_type, & + my_obs_loc, my_obs_kind, my_obs_type, num_close_obs, close_obs_ind, close_obs_dist, & + ens_handle, last_base_obs_loc, last_num_close_obs, last_close_obs_ind, & + last_close_obs_dist, num_close_obs_cached, num_close_obs_calls_made) + +logical, intent(in) :: close_obs_caching +type(get_close_type), intent(in) :: gc_obs +type(location_type), intent(inout) :: base_obs_loc, my_obs_loc(:) +integer, intent(in) :: base_obs_type, my_obs_kind(:), my_obs_type(:) +integer, intent(out) :: num_close_obs, close_obs_ind(:) +real(r8), intent(out) :: close_obs_dist(:) +type(ensemble_type), intent(in) :: ens_handle +type(location_type), intent(inout) :: last_base_obs_loc +integer, intent(inout) :: last_num_close_obs +integer, intent(inout) :: last_close_obs_ind(:) +real(r8), intent(inout) :: last_close_obs_dist(:) +integer, intent(inout) :: num_close_obs_cached, num_close_obs_calls_made + +! This logic could be arranged to make code less redundant +if (.not. close_obs_caching) then + call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & + my_obs_loc, my_obs_kind, my_obs_type, & + num_close_obs, close_obs_ind, close_obs_dist, ens_handle) +else + if (base_obs_loc == last_base_obs_loc) then + num_close_obs = last_num_close_obs + close_obs_ind(:) = last_close_obs_ind(:) + close_obs_dist(:) = last_close_obs_dist(:) + num_close_obs_cached = num_close_obs_cached + 1 + else + call get_close_obs(gc_obs, base_obs_loc, base_obs_type, & + my_obs_loc, my_obs_kind, my_obs_type, & + num_close_obs, close_obs_ind, close_obs_dist, ens_handle) -!-------------------------------------------------------------------- -!> -!> t_space is where we store the time information -!> label is the string to print out with the time. limited to ~60 chars. -!> itemcount is the running count of how many times we've been called -!> maxitems is a limit on the number of times we want this timer to print. -!> do_sync overrides the default for whether we want to do a task sync -!> or not. right now this code defaults to yes, sync before getting -!> the time. for very large processor counts this increases overhead. -!> elapsed is an optional return of the value instead of only printing here - -subroutine read_timer(t_space, label, itemcount, maxitems, do_sync, elapsed) - real(digits12), intent(in) :: t_space - character(len=*), intent(in) :: label - integer(i8), intent(inout), optional :: itemcount - integer(i8), intent(in), optional :: maxitems - logical, intent(in), optional :: do_sync - real(digits12), intent(out), optional :: elapsed - -real(digits12) :: interval -logical :: sync_me -character(len=132) :: buffer - -! if interval time (in seconds) is > this, go ahead and -! print even if item count is over limit. -integer(i8), parameter :: T_ALWAYS_PRINT = 1.0_r8 - -! get the time first and then figure out what we're doing - -interval = read_mpi_timer(t_space) - -sync_me = .true. -if (present(do_sync)) sync_me = do_sync - -! if there's a limit on number of prints don't allow sync -! because if we return early we could hang everyone else. -if (present(maxitems)) sync_me = .false. - -! print out large values no matter what -! (large is defined above locally in this routine) -if (present(itemcount) .and. present(maxitems)) then - if (interval < T_ALWAYS_PRINT .and. itemcount > maxitems) return + last_base_obs_loc = base_obs_loc + last_num_close_obs = num_close_obs + last_close_obs_ind(:) = close_obs_ind(:) + last_close_obs_dist(:) = close_obs_dist(:) + num_close_obs_calls_made = num_close_obs_calls_made +1 + endif endif -! if syncing, wait and read the timer again -if (sync_me) then - call task_sync() - interval = read_mpi_timer(t_space) -endif +end subroutine get_close_obs_cached -if (sync_me) then - if (my_task_id() == 0) then - write(buffer,'(A75,F15.8)') "timer: "//trim(label)//" time ", interval - write(*,*) buffer - endif +!-------------------------------------------------------------------- +!> Get close state from cache if appropriate. Cache new get_close_state info +!> if requested. + +subroutine get_close_state_cached(close_obs_caching, gc_state, base_obs_loc, base_obs_type, & + my_state_loc, my_state_kind, my_state_indx, num_close_states, close_state_ind, close_state_dist, & + ens_handle, last_base_states_loc, last_num_close_states, last_close_state_ind, & + last_close_state_dist, num_close_states_cached, num_close_states_calls_made) + +logical, intent(in) :: close_obs_caching +type(get_close_type), intent(in) :: gc_state +type(location_type), intent(inout) :: base_obs_loc, my_state_loc(:) +integer, intent(in) :: base_obs_type, my_state_kind(:) +integer(i8), intent(in) :: my_state_indx(:) +integer, intent(out) :: num_close_states, close_state_ind(:) +real(r8), intent(out) :: close_state_dist(:) +type(ensemble_type), intent(in) :: ens_handle +type(location_type), intent(inout) :: last_base_states_loc +integer, intent(inout) :: last_num_close_states +integer, intent(inout) :: last_close_state_ind(:) +real(r8), intent(inout) :: last_close_state_dist(:) +integer, intent(inout) :: num_close_states_cached, num_close_states_calls_made + +! This logic could be arranged to make code less redundant +if (.not. close_obs_caching) then + call get_close_state(gc_state, base_obs_loc, base_obs_type, & + my_state_loc, my_state_kind, my_state_indx, & + num_close_states, close_state_ind, close_state_dist, ens_handle) else - write(buffer,'(A75,F15.8,A6,I7)') "timer: "//trim(label)//" time ", interval, " rank ", my_task_id() - write(*,*) buffer + if (base_obs_loc == last_base_states_loc) then + num_close_states = last_num_close_states + close_state_ind(:) = last_close_state_ind(:) + close_state_dist(:) = last_close_state_dist(:) + num_close_states_cached = num_close_states_cached + 1 + else + call get_close_state(gc_state, base_obs_loc, base_obs_type, & + my_state_loc, my_state_kind, my_state_indx, & + num_close_states, close_state_ind, close_state_dist, ens_handle) + + last_base_states_loc = base_obs_loc + last_num_close_states = num_close_states + last_close_state_ind(:) = close_state_ind(:) + last_close_state_dist(:) = close_state_dist(:) + num_close_states_calls_made = num_close_states_calls_made +1 + endif endif -if (present(elapsed)) elapsed = interval - -end subroutine read_timer +end subroutine get_close_state_cached !-------------------------------------------------------------------- !> log what the user has selected via the namelist choices From eb4aa8f9dfab71a84d40bb0bb974bcee30dfa259 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Fri, 5 Nov 2021 12:03:48 -0600 Subject: [PATCH 08/18] Cleaned up the computation of the various localization methods in the state and obs update loops. Tested with obs_impact in addition to all other variants. --- .../modules/assimilation/assim_tools_mod.f90 | 117 +++++++++--------- 1 file changed, 61 insertions(+), 56 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index d9ece542e1..8401e700ef 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -324,7 +324,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! changed the ensemble sized things here to allocatable real(r8) :: obs_prior(ens_size), obs_inc(ens_size), increment(ens_size) -real(r8) :: reg_factor, impact_factor +real(r8) :: impact_factor, reg_factor, final_factor real(r8) :: net_a(num_groups), reg_coef(num_groups), correl(num_groups) real(r8) :: cov_factor, obs(1), obs_err_var, my_inflate, my_inflate_sd real(r8) :: obs_qc, cutoff_rev, cutoff_orig @@ -536,8 +536,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & end do endif -! The computations in the two get_close_maxdist_init are redundant - ! Initialize the method for getting state variables close to a given ob on my process if (has_special_cutoffs) then call get_close_init(gc_state, my_num_state, 2.0_r8*cutoff, my_state_loc, 2.0_r8*cutoff_list) @@ -727,28 +725,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & STATE_UPDATE: do j = 1, num_close_states state_index = close_state_ind(j) - ! 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. if ( allow_missing_in_state ) then - ! Some models can take evasive action if one or more of the ensembles have - ! a missing value. Generally means 'do nothing' (as opposed to DIE) + ! Don't allow update of state ensemble with any missing values if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE 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)) - - ! if external impact factors supplied, factor them in here - ! FIXME: this would execute faster for 0.0 impact factors if - ! we check for that before calling comp_cov_factor. but it makes - ! the logic more complicated - this is simpler if we do it after. - if (adjust_obs_impact) then - impact_factor = obs_impact_table(base_obs_type, my_state_kind(state_index)) - cov_factor = cov_factor * impact_factor - endif - - ! If no weight is indicated, no more to do with this state variable + ! Compute the covariance localization and adjust_obs_impact factors + cov_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_state_loc(state_index), & + my_state_kind(state_index), close_state_dist(j), cutoff_rev, adjust_obs_impact, obs_impact_table) if(cov_factor <= 0.0_r8) cycle STATE_UPDATE ! Loop through groups to update the state variable ensemble members @@ -762,23 +746,18 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) end do - ! Compute an information factor for impact of this observation on this state - if(num_groups == 1) then - reg_factor = 1.0_r8 - else - ! Pass the time along with the index for possible diagnostic output - ! Compute regression factor for this obs-state pair + if(num_groups <= 1) then + final_factor = cov_factor + else reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index)) + final_factor = min(cov_factor, reg_factor) endif - ! The final factor is the minimum of group regression factor and localization cov_factor - reg_factor = min(reg_factor, cov_factor) - !PAR NEED TO TURN STUFF OFF MORE EFFICEINTLY ! If doing full assimilation, update the state variable ensemble with weighted increments if(.not. inflate_only) then ens_handle%copies(1:ens_size, state_index) = & - ens_handle%copies(1:ens_size, state_index) + reg_factor * increment + ens_handle%copies(1:ens_size, state_index) + final_factor * increment endif ! Compute spatially-varying state space inflation @@ -789,7 +768,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ens_handle%copies(ENS_INF_SD_COPY, state_index), & ens_handle%copies(ENS_SD_COPY, state_index), & orig_obs_prior_mean(group), orig_obs_prior_var(group), obs(1), & - obs_err_var, grp_size, reg_factor, correl(group), inflate_only) + obs_err_var, grp_size, final_factor, correl(group), inflate_only) end do endif @@ -810,19 +789,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! update the unassimilated observations if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE - ! Compute the distance and the covar_factor - cov_factor = comp_cov_factor(close_obs_dist(j), cutoff_rev, & - base_obs_loc, base_obs_type, my_obs_loc(obs_index), my_obs_kind(obs_index)) - - ! if external impact factors supplied, factor them in here - ! FIXME: this would execute faster for 0.0 impact factors if - ! we check for that before calling comp_cov_factor. but it makes - ! the logic more complicated - this is simpler if we do it after. - if (adjust_obs_impact) then - impact_factor = obs_impact_table(base_obs_type, my_obs_kind(obs_index)) - cov_factor = cov_factor * impact_factor - endif - + ! Compute the covariance localization and adjust_obs_impact factors + cov_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & + my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev, adjust_obs_impact, obs_impact_table) if(cov_factor <= 0.0_r8) cycle OBS_UPDATE ! Loop through and update ensemble members in each group @@ -838,22 +807,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! FIXME: could we move the if test for inflate only to here? ! Compute an information factor for impact of this observation on this state - if(num_groups == 1) then - reg_factor = 1.0_r8 - else - ! Pass the time along with the index for possible diagnostic output - ! Compute regression factor for this obs-state pair - ! Negative indicates that this is an observation index - reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, -1*my_obs_indx(obs_index)) + if(num_groups <= 1) then + final_factor = cov_factor + else + reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index)) + final_factor = min(cov_factor, reg_factor) endif - ! Final weight is min of group and localization factors - reg_factor = min(reg_factor, cov_factor) - ! Only update state if indicated (otherwise just getting inflation) if(.not. inflate_only) then obs_ens_handle%copies(1:ens_size, obs_index) = & - obs_ens_handle%copies(1:ens_size, obs_index) + reg_factor * increment + obs_ens_handle%copies(1:ens_size, obs_index) + final_factor * increment endif endif end do OBS_UPDATE @@ -2220,6 +2184,47 @@ subroutine update_ens_from_weights(ens, ens_size, rel_weight, ens_inc) end subroutine update_ens_from_weights +!------------------------------------------------------------- + +function cov_and_impact_factors(base_obs_loc, base_obs_type, state_loc, state_kind, & +dist, cutoff_rev, adjust_obs_impact, obs_impact_table) + +! Computes the cov_factor and multiplies by obs_impact_factor if selected + +real(r8) :: cov_and_impact_factors +type(location_type), intent(in) :: base_obs_loc +integer, intent(in) :: base_obs_type +type(location_type), intent(in) :: state_loc +integer, intent(in) :: state_kind +real(r8), intent(in) :: dist +real(r8), intent(in) :: cutoff_rev +logical, intent(in) :: adjust_obs_impact +real(r8), intent(in) :: obs_impact_table(:, 0:) + +real(r8) :: impact_factor, cov_factor + +! Get external impact factors, cycle if impact of this ob on this state is zerio +if (adjust_obs_impact) then + ! Get the impact factor from the table if requested + impact_factor = obs_impact_table(base_obs_type, state_kind) + if(impact_factor <= 0.0_r8) then + ! Avoid the cost of computing cov_factor if impact is 0 + cov_and_impact_factors = 0.0_r8 + return + endif +else + impact_factor = 1.0_r8 +endif + +! Compute the covariance factor +cov_factor = comp_cov_factor(dist, cutoff_rev, & + base_obs_loc, base_obs_type, state_loc, state_kind) + +! Combine the impact_factor and the cov_factor +cov_and_impact_factors = cov_factor * impact_factor + +end function cov_and_impact_factors + !------------------------------------------------------------------------ From 7a2d1f25aa43a8a21bd9b6548236833fa00cba8d Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Fri, 5 Nov 2021 14:08:24 -0600 Subject: [PATCH 09/18] Entire observation update loop can be skipped for posterior inflate update_only case. Additional cleanup to compact the update loops. --- .../modules/assimilation/assim_tools_mod.f90 | 92 ++++++++----------- 1 file changed, 37 insertions(+), 55 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 8401e700ef..726805bd73 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -664,8 +664,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Compute observation space increments for each group do group = 1, num_groups - grp_bot = grp_beg(group) - grp_top = grp_end(group) + grp_bot = grp_beg(group); grp_top = grp_end(group) 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)) @@ -717,9 +716,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ens_handle, last_base_states_loc, last_num_close_states, last_close_state_ind, & last_close_state_dist, num_close_states_cached, num_close_states_calls_made) n_close_state_items(i) = num_close_states - !print*, 'num close state', num_close_states !call test_close_obs_dist(close_state_dist, num_close_states, i) - !call test_state_copies(ens_handle, 'beforeupdates') ! Loop through to update each of my state variables that is potentially close STATE_UPDATE: do j = 1, num_close_states @@ -737,8 +734,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Loop through groups to update the state variable ensemble members do group = 1, num_groups - grp_bot = grp_beg(group) - grp_top = grp_end(group) + grp_bot = grp_beg(group); grp_top = grp_end(group) ! Do update of state, correl only needed for varying ss inflate but compute for all call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & obs_prior_var(group), obs_inc(grp_bot:grp_top), & @@ -753,7 +749,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & final_factor = min(cov_factor, reg_factor) endif -!PAR NEED TO TURN STUFF OFF MORE EFFICEINTLY ! If doing full assimilation, update the state variable ensemble with weighted increments if(.not. inflate_only) then ens_handle%copies(1:ens_size, state_index) = & @@ -771,59 +766,46 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & obs_err_var, grp_size, final_factor, correl(group), inflate_only) end do endif - end do STATE_UPDATE - !call test_state_copies(ens_handle, 'after_state_updates') - - !------------------------------------------------------ - - ! Now everybody updates their obs priors (only ones after this one) - 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 - - ! If the forward observation operator failed, no need to - ! update the unassimilated observations - if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE - - ! Compute the covariance localization and adjust_obs_impact factors - cov_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & - my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev, adjust_obs_impact, obs_impact_table) - if(cov_factor <= 0.0_r8) cycle OBS_UPDATE - - ! Loop through and update ensemble members in each group - do group = 1, num_groups - grp_bot = grp_beg(group) - grp_top = grp_end(group) - call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & - 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 - - ! FIXME: could we move the if test for inflate only to here? - - ! Compute an information factor for impact of this observation on this state - if(num_groups <= 1) then - final_factor = cov_factor - else - reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index)) - final_factor = min(cov_factor, reg_factor) - endif + if(.not. inflate_only) then + ! Now everybody updates their obs priors (only ones after this one) + 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 + + ! If forward observation operator failed, no need to update unassimilated observations + if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE + + ! Compute the covariance localization and adjust_obs_impact factors + cov_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & + my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev, adjust_obs_impact, obs_impact_table) + if(cov_factor <= 0.0_r8) cycle OBS_UPDATE + + ! Loop through and update ensemble members in each group + do group = 1, num_groups + grp_bot = grp_beg(group); grp_top = grp_end(group) + call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & + 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 + + ! Compute an information factor for impact of this observation on this state + if(num_groups <= 1) then + final_factor = cov_factor + else + reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index)) + final_factor = min(cov_factor, reg_factor) + endif - ! Only update state if indicated (otherwise just getting inflation) - if(.not. inflate_only) then obs_ens_handle%copies(1:ens_size, obs_index) = & - obs_ens_handle%copies(1:ens_size, obs_index) + final_factor * increment + obs_ens_handle%copies(1:ens_size, obs_index) + final_factor * increment endif - endif - end do OBS_UPDATE - - !call test_state_copies(ens_handle, 'after_obs_updates') - + end do OBS_UPDATE + endif end do SEQUENTIAL_OBS ! Every pe needs to get the current my_inflate and my_inflate_sd back From b5d9aa18ef657de708c3ff4e558912a548ca67bc Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Tue, 9 Nov 2021 15:26:33 -0700 Subject: [PATCH 10/18] Removed the redundant outlier threshold check for posterior since posterior forward operator was already checking for outlier threshold. Note that the deleted code checked the outlier threshold for each group and so could result in different answers for posterior inflation with a group filter. This code bitwise reproduces the main branch once this check has also been removed there. --- .../assimilation/adaptive_inflate_mod.f90 | 33 +++---------------- 1 file changed, 4 insertions(+), 29 deletions(-) diff --git a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 index e20f3c1c3e..d7d81cf6ed 100644 --- a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 +++ b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 @@ -705,13 +705,13 @@ end subroutine update_single_state_space_inflation !------------------------------------------------------------------------------- !> Computes updatea inflation mean and inflation sd for varying state space inflation -subroutine update_varying_state_space_inflation(inflate, inflate_mean_inout, inflate_sd_inout, & +subroutine update_varying_state_space_inflation(inflate, inflate_mean, inflate_sd, & ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, & ens_size, reg_factor, correl, inflate_only) type(adaptive_inflate_type), intent(in) :: inflate -real(r8), intent(inout) :: inflate_mean_inout -real(r8), intent(inout) :: inflate_sd_inout +real(r8), intent(inout) :: inflate_mean +real(r8), intent(inout) :: inflate_sd real(r8), intent(in) :: ss_inflate_base real(r8), intent(in) :: orig_obs_prior_mean real(r8), intent(in) :: orig_obs_prior_var @@ -722,14 +722,10 @@ subroutine update_varying_state_space_inflation(inflate, inflate_mean_inout, inf real(r8), intent(in) :: correl logical, intent(in) :: inflate_only -real(r8) :: inflate_mean, inflate_sd, gamma, ens_var_deflate, r_var, r_mean +real(r8) :: gamma, ens_var_deflate, r_var, r_mean real(r8) :: diff_sd, outlier_ratio logical :: do_adapt_inf_update -! Copy the inflate_mean and inflate_sd which can be changed by update_inflation -inflate_mean = inflate_mean_inout -inflate_sd = inflate_sd_inout - if(inflate_mean > 0.0_r8 .and. inflate_sd > 0.0_r8) then ! Gamma is less than 1 for varying ss, see adaptive inflate module gamma = reg_factor * abs(correl) @@ -758,27 +754,6 @@ subroutine update_varying_state_space_inflation(inflate, inflate_mean_inout, inf ! Update the inflation values call update_inflation(inflate, inflate_mean, inflate_sd, & r_mean, r_var, ens_size, obs, obs_err_var, gamma) -else - ! if we don't go into the previous if block, make sure these - ! have good values going out for the block below - r_mean = orig_obs_prior_mean - r_var = orig_obs_prior_var -endif - -! Update adaptive values if posterior outlier_ratio test doesn't fail. -! Match code in obs_space_diags() in filter.f90 -do_adapt_inf_update = .true. -if (inflate_only) then - diff_sd = sqrt(obs_err_var + r_var) - if (diff_sd > 0.0_r8) then - outlier_ratio = abs(obs - r_mean) / diff_sd -! JLA: ********** Magic number of 3 sd outlier here needs to be addressed ********* - do_adapt_inf_update = (outlier_ratio <= 3.0_r8) - endif -endif -if (do_adapt_inf_update) then - inflate_mean_inout = inflate_mean - inflate_sd_inout = inflate_sd endif end subroutine update_varying_state_space_inflation From cb2711a100eac8f13d4e51a57d097fecce5a8337 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Wed, 10 Nov 2021 11:47:03 -0700 Subject: [PATCH 11/18] Moved the update of either state or obs ensembles by obs to a single subroutine. Removed net_a from the broadcast since it is computed by all processes under the revised communication that was implemented earlier. Still appears to bitwise reproduce. --- .../modules/assimilation/assim_tools_mod.f90 | 140 +++++++++++------- 1 file changed, 87 insertions(+), 53 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 726805bd73..050b38a4dc 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -323,8 +323,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! changed the ensemble sized things here to allocatable -real(r8) :: obs_prior(ens_size), obs_inc(ens_size), increment(ens_size) -real(r8) :: impact_factor, reg_factor, final_factor +real(r8) :: obs_prior(ens_size), obs_inc(ens_size), increment(ens_size), updated_ens(ens_size) +real(r8) :: reg_factor, final_factor real(r8) :: net_a(num_groups), reg_coef(num_groups), correl(num_groups) real(r8) :: cov_factor, obs(1), obs_err_var, my_inflate, my_inflate_sd real(r8) :: obs_qc, cutoff_rev, cutoff_orig @@ -638,7 +638,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! coordinate transformation whichvert_real = real(whichvert_obs_in_localization_coord, r8) call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, & - net_a, orig_obs_prior_mean, orig_obs_prior_var, & + orig_obs_prior_mean, orig_obs_prior_var, & scalar1=obs_qc, scalar2=vertvalue_obs_in_localization_coord, & scalar3=whichvert_real, scalar4=my_inflate, scalar5=my_inflate_sd) @@ -646,7 +646,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & !----------------------------------------------------------------------- else call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, & - net_a, orig_obs_prior_mean, orig_obs_prior_var, & + orig_obs_prior_mean, orig_obs_prior_var, & scalar1=obs_qc, scalar2=vertvalue_obs_in_localization_coord, & scalar3=whichvert_real, scalar4=my_inflate, scalar5=my_inflate_sd) whichvert_obs_in_localization_coord = nint(whichvert_real) @@ -727,36 +727,17 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE endif - ! Compute the covariance localization and adjust_obs_impact factors - cov_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_state_loc(state_index), & - my_state_kind(state_index), close_state_dist(j), cutoff_rev, adjust_obs_impact, obs_impact_table) - if(cov_factor <= 0.0_r8) cycle STATE_UPDATE - - ! Loop through groups to update the state variable ensemble members - do group = 1, num_groups - grp_bot = grp_beg(group); grp_top = grp_end(group) - ! Do update of state, correl only needed for varying ss inflate but compute for all - call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & - obs_prior_var(group), obs_inc(grp_bot:grp_top), & - ens_handle%copies(grp_bot:grp_top, state_index), grp_size, & - increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) - end do - - if(num_groups <= 1) then - final_factor = cov_factor - else - reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index)) - final_factor = min(cov_factor, reg_factor) - endif + call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & + updated_ens, my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & + obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & + close_state_dist(j), cutoff_rev, net_a, adjust_obs_impact, obs_impact_table, & + grp_size, grp_beg, grp_end, i, my_state_indx(state_index), final_factor, correl) ! If doing full assimilation, update the state variable ensemble with weighted increments - if(.not. inflate_only) then - ens_handle%copies(1:ens_size, state_index) = & - ens_handle%copies(1:ens_size, state_index) + final_factor * increment - endif + if(.not. inflate_only) ens_handle%copies(1:ens_size, state_index) = updated_ens ! Compute spatially-varying state space inflation - if(local_varying_ss_inflate) then + if(local_varying_ss_inflate .and. final_factor > 0.0_r8) then do group = 1, num_groups call update_varying_state_space_inflation(inflate, & ens_handle%copies(ENS_INF_COPY, state_index), & @@ -779,30 +760,13 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! If forward observation operator failed, no need to update unassimilated observations if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE - ! Compute the covariance localization and adjust_obs_impact factors - cov_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & - my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev, adjust_obs_impact, obs_impact_table) - if(cov_factor <= 0.0_r8) cycle OBS_UPDATE - - ! Loop through and update ensemble members in each group - do group = 1, num_groups - grp_bot = grp_beg(group); grp_top = grp_end(group) - call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & - 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 - - ! Compute an information factor for impact of this observation on this state - if(num_groups <= 1) then - final_factor = cov_factor - else - reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index)) - final_factor = min(cov_factor, reg_factor) - endif + call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & + updated_ens, my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & + obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & + close_obs_dist(j), cutoff_rev, net_a, adjust_obs_impact, obs_impact_table, & + grp_size, grp_beg, grp_end, i, -1*my_obs_indx(obs_index), final_factor, correl) - obs_ens_handle%copies(1:ens_size, obs_index) = & - obs_ens_handle%copies(1:ens_size, obs_index) + final_factor * increment + obs_ens_handle%copies(1:ens_size, obs_index) = updated_ens endif end do OBS_UPDATE endif @@ -2165,6 +2129,76 @@ subroutine update_ens_from_weights(ens, ens_size, rel_weight, ens_inc) end do end subroutine update_ens_from_weights +!--------------------------------------------------------------- + +subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_kind, & + obs_prior, obs_inc, obs_prior_mean, obs_prior_var, obs_loc, obs_type, obs_time, & + dist, cutoff_rev, net_a, adjust_obs_impact, obs_impact_table, & + grp_size, grp_beg, grp_end, reg_factor_obs_index, reg_factor_ens_index, & + final_factor, correl) + +integer, intent(in) :: ens_size +integer, intent(in) :: num_groups +real(r8), intent(in) :: ens(ens_size) +real(r8), intent(out) :: updated_ens(ens_size) +type(location_type), intent(in) :: ens_loc +integer, intent(in) :: ens_kind +real(r8), intent(in) :: obs_prior(ens_size) +real(r8), intent(in) :: obs_inc(ens_size) +real(r8), intent(in) :: obs_prior_mean(num_groups) +real(r8), intent(in) :: obs_prior_var(num_groups) +type(location_type), intent(in) :: obs_loc +integer, intent(in) :: obs_type +type(time_type), intent(in) :: obs_time +real(r8), intent(in) :: dist +real(r8), intent(in) :: cutoff_rev +real(r8), intent(inout) :: net_a(num_groups) +logical, intent(in) :: adjust_obs_impact +real(r8), intent(in) :: obs_impact_table(:, :) +integer, intent(in) :: grp_size +integer, intent(in) :: grp_beg(num_groups) +integer, intent(in) :: grp_end(num_groups) +integer, intent(in) :: reg_factor_obs_index +integer(i8), intent(in) :: reg_factor_ens_index +real(r8), intent(out) :: final_factor +real(r8), intent(out) :: correl(num_groups) + +real(r8) :: reg_coef(num_groups), increment(ens_size) +real(r8) :: cov_factor, reg_factor +integer :: group, grp_bot, grp_top + +! Compute the covariance localization and adjust_obs_impact factors +cov_factor = cov_and_impact_factors(obs_loc, obs_type, ens_loc, & + ens_kind, dist, cutoff_rev, adjust_obs_impact, obs_impact_table) + +! If no impact, don't do anything else +if(cov_factor <= 0.0_r8) then + final_factor = cov_factor + updated_ens = ens + return +endif + +! Loop through groups to update the state variable ensemble members +do group = 1, num_groups + grp_bot = grp_beg(group); grp_top = grp_end(group) + ! Do update of state, correl only needed for varying ss inflate but compute for all + call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & + obs_prior_var(group), obs_inc(grp_bot:grp_top), ens(grp_bot:grp_top), grp_size, & + increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) +end do + +if(num_groups <= 1) then + final_factor = cov_factor +else + reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, & + reg_factor_obs_index, reg_factor_ens_index) + final_factor = min(cov_factor, reg_factor) +endif + +! Get the updated ensemble +updated_ens = ens + final_factor * increment + +end subroutine obs_updates_ens !------------------------------------------------------------- From 660aabfed913fcb578cfd51c4e9cf1cb0413efc1 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Thu, 11 Nov 2021 08:13:56 -0700 Subject: [PATCH 12/18] Made net_a into an intent in argument for state updates to avoid confusion. Removed unused declarations. --- .../modules/assimilation/assim_tools_mod.f90 | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 050b38a4dc..e515995ccc 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -323,10 +323,10 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! changed the ensemble sized things here to allocatable -real(r8) :: obs_prior(ens_size), obs_inc(ens_size), increment(ens_size), updated_ens(ens_size) -real(r8) :: reg_factor, final_factor -real(r8) :: net_a(num_groups), reg_coef(num_groups), correl(num_groups) -real(r8) :: cov_factor, obs(1), obs_err_var, my_inflate, my_inflate_sd +real(r8) :: obs_prior(ens_size), obs_inc(ens_size), updated_ens(ens_size) +real(r8) :: final_factor +real(r8) :: net_a(num_groups), correl(num_groups) +real(r8) :: obs(1), obs_err_var, my_inflate, my_inflate_sd real(r8) :: obs_qc, cutoff_rev, cutoff_orig real(r8) :: orig_obs_prior_mean(num_groups), orig_obs_prior_var(num_groups) real(r8) :: obs_prior_mean(num_groups), obs_prior_var(num_groups) @@ -1374,7 +1374,7 @@ end subroutine obs_increment_kernel subroutine update_from_obs_inc(obs, obs_prior_mean, obs_prior_var, obs_inc, & - state, ens_size, state_inc, reg_coef, net_a, correl_out) + state, ens_size, state_inc, reg_coef, net_a_in, correl_out) !======================================================================== ! Does linear regression of a state variable onto an observation and @@ -1385,12 +1385,12 @@ subroutine update_from_obs_inc(obs, obs_prior_mean, obs_prior_var, obs_inc, & real(r8), intent(in) :: obs_prior_mean, obs_prior_var real(r8), intent(in) :: state(ens_size) real(r8), intent(out) :: state_inc(ens_size), reg_coef -real(r8), intent(inout) :: net_a +real(r8), intent(in) :: net_a_in real(r8), optional, intent(inout) :: correl_out real(r8) :: obs_state_cov, intermed real(r8) :: restoration_inc(ens_size), state_mean, state_var, correl -real(r8) :: factor, exp_true_correl, mean_factor +real(r8) :: factor, exp_true_correl, mean_factor, net_a ! For efficiency, just compute regression coefficient here unless correl is needed @@ -1465,7 +1465,7 @@ subroutine update_from_obs_inc(obs, obs_prior_mean, obs_prior_var, obs_inc, & ! Spread restoration algorithm option. if(spread_restoration) then ! Don't use this to reduce spread at present (should revisit this line) - if(net_a > 1.0_r8) net_a = 1.0_r8 + net_a = min(net_a_in, 1.0_r8) ! Default restoration increment is 0.0 restoration_inc = 0.0_r8 @@ -2152,14 +2152,14 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ type(time_type), intent(in) :: obs_time real(r8), intent(in) :: dist real(r8), intent(in) :: cutoff_rev -real(r8), intent(inout) :: net_a(num_groups) +real(r8), intent(in) :: net_a(num_groups) logical, intent(in) :: adjust_obs_impact real(r8), intent(in) :: obs_impact_table(:, :) integer, intent(in) :: grp_size integer, intent(in) :: grp_beg(num_groups) integer, intent(in) :: grp_end(num_groups) integer, intent(in) :: reg_factor_obs_index -integer(i8), intent(in) :: reg_factor_ens_index +integer(i8), intent(in) :: reg_factor_ens_index real(r8), intent(out) :: final_factor real(r8), intent(out) :: correl(num_groups) From 2d71cc8d24d1e83d83773c3493c9588514490557 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Thu, 18 Nov 2021 13:35:01 -0700 Subject: [PATCH 13/18] Removed unused code. Removed module storage variables from callling hierarchy. This fixed a problem with passing unallocated allocatable arrays to subroutines that was inconsistent with language spec. --- .../modules/assimilation/assim_tools_mod.f90 | 61 ++++--------------- 1 file changed, 11 insertions(+), 50 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index e515995ccc..c338de8cfc 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -48,7 +48,7 @@ module assim_tools_mod use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars, & compute_copy_mean_var, get_var_owner_index, & - prepare_to_update_copies, map_pe_to_task + map_pe_to_task use mpi_utilities_mod, only : my_task_id, broadcast_send, broadcast_recv, & sum_across_tasks, task_count, start_mpi_timer, & @@ -93,7 +93,6 @@ module assim_tools_mod real(r8), allocatable :: cutoff_list(:) logical :: has_special_cutoffs logical :: close_obs_caching = .true. -real(r8), parameter :: small = epsilon(1.0_r8) ! threshold for avoiding NaNs/Inf ! true if we have multiple vert choices and we're doing vertical localization ! (make it a local variable so we don't keep making subroutine calls) @@ -375,15 +374,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & logical :: local_ss_inflate logical :: local_obs_inflate -integer, allocatable :: n_close_state_items(:), n_close_obs_items(:) - -! Just to make sure multiple tasks are running for tests -write(*, *) 'my_task_id ', my_task_id() - -! how about this? look for imbalances in the tasks -allocate(n_close_state_items(obs_ens_handle%num_vars), & - n_close_obs_items( obs_ens_handle%num_vars)) - ! 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), & @@ -404,10 +394,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & my_state_loc( ens_handle%my_num_vars)) ! end alloc -! we are going to read/write the copies array -call prepare_to_update_copies(ens_handle) -call prepare_to_update_copies(obs_ens_handle) - ! Initialize assim_tools_module if needed if (.not. module_initialized) call assim_tools_init() @@ -689,11 +675,10 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Adaptive localization needs number of other observations within localization radius. ! Do get_close_obs first, even though state space increments are computed before obs increments. ! JLA: ens_handle doesn't ever appear to be used. Get rid of it. Should be obs_ens_handle anyway? - call get_close_obs_cached(close_obs_caching, gc_obs, base_obs_loc, base_obs_type, & + call get_close_obs_cached(gc_obs, base_obs_loc, base_obs_type, & my_obs_loc, my_obs_kind, my_obs_type, num_close_obs, close_obs_ind, close_obs_dist, & ens_handle, last_base_obs_loc, last_num_close_obs, last_close_obs_ind, & last_close_obs_dist, num_close_obs_cached, num_close_obs_calls_made) - n_close_obs_items(i) = num_close_obs ! set the cutoff default, keep a copy of the original value, and avoid ! looking up the cutoff in a list if the incoming obs is an identity ob @@ -711,11 +696,10 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & i, base_obs_loc, obs_def, localization_unit) ! Find state variables on my process that are close to observation being assimilated - call get_close_state_cached(close_obs_caching, gc_state, base_obs_loc, base_obs_type, & + call get_close_state_cached(gc_state, base_obs_loc, base_obs_type, & my_state_loc, my_state_kind, my_state_indx, num_close_states, close_state_ind, close_state_dist, & ens_handle, last_base_states_loc, last_num_close_states, last_close_state_ind, & last_close_state_dist, num_close_states_cached, num_close_states_calls_made) - n_close_state_items(i) = num_close_states !call test_close_obs_dist(close_state_dist, num_close_states, i) ! Loop through to update each of my state variables that is potentially close @@ -730,7 +714,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & updated_ens, my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - close_state_dist(j), cutoff_rev, net_a, adjust_obs_impact, obs_impact_table, & + close_state_dist(j), cutoff_rev, net_a, & grp_size, grp_beg, grp_end, i, my_state_indx(state_index), final_factor, correl) ! If doing full assimilation, update the state variable ensemble with weighted increments @@ -763,7 +747,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & updated_ens, my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - close_obs_dist(j), cutoff_rev, net_a, adjust_obs_impact, obs_impact_table, & + close_obs_dist(j), cutoff_rev, net_a, & grp_size, grp_beg, grp_end, i, -1*my_obs_indx(obs_index), final_factor, correl) obs_ens_handle%copies(1:ens_size, obs_index) = updated_ens @@ -787,21 +771,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! this is going to be per-task. so only print if something interesting ! shows up in the stats? maybe it would be worth a reduce() call here? -!>@todo FIXME: -! we have n_close_obs_items and n_close_state_items for each assimilated -! observation. what we really want to know is across the tasks is there -! a big difference in counts? so that means communication. maybe just -! the largest value? and the number of 0 values? and if the largest val -! is way off compared to the other tasks, warn the user? -! we don't have space or time to do all the obs * tasks but could we -! send enough info to make a histogram? compute N bin counts and then -! reduce that across all the tasks and have task 0 print out? -! still thinking on this idea. -! write(msgstring, *) 'max state items per observation: ', maxval(n_close_state_items) -! call error_handler(E_MSG, 'filter_assim:', msgstring) -! if i come up with something i like, can we use the same idea -! for the threed_sphere locations boxes? - ! Assure user we have done something if (print_trace_details >= 0) then write(msgstring, '(A,I8,A)') & @@ -854,8 +823,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & my_state_kind, & my_state_loc) -deallocate(n_close_state_items, & - n_close_obs_items) ! end dealloc end subroutine filter_assim @@ -2133,7 +2100,7 @@ end subroutine update_ens_from_weights subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_kind, & obs_prior, obs_inc, obs_prior_mean, obs_prior_var, obs_loc, obs_type, obs_time, & - dist, cutoff_rev, net_a, adjust_obs_impact, obs_impact_table, & + dist, cutoff_rev, net_a, & grp_size, grp_beg, grp_end, reg_factor_obs_index, reg_factor_ens_index, & final_factor, correl) @@ -2153,8 +2120,6 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ real(r8), intent(in) :: dist real(r8), intent(in) :: cutoff_rev real(r8), intent(in) :: net_a(num_groups) -logical, intent(in) :: adjust_obs_impact -real(r8), intent(in) :: obs_impact_table(:, :) integer, intent(in) :: grp_size integer, intent(in) :: grp_beg(num_groups) integer, intent(in) :: grp_end(num_groups) @@ -2167,9 +2132,9 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ real(r8) :: cov_factor, reg_factor integer :: group, grp_bot, grp_top -! Compute the covariance localization and adjust_obs_impact factors +! Compute the covariance localization and adjust_obs_impact factors (module storage) cov_factor = cov_and_impact_factors(obs_loc, obs_type, ens_loc, & - ens_kind, dist, cutoff_rev, adjust_obs_impact, obs_impact_table) + ens_kind, dist, cutoff_rev) ! If no impact, don't do anything else if(cov_factor <= 0.0_r8) then @@ -2203,7 +2168,7 @@ end subroutine obs_updates_ens !------------------------------------------------------------- function cov_and_impact_factors(base_obs_loc, base_obs_type, state_loc, state_kind, & -dist, cutoff_rev, adjust_obs_impact, obs_impact_table) +dist, cutoff_rev) ! Computes the cov_factor and multiplies by obs_impact_factor if selected @@ -2214,8 +2179,6 @@ function cov_and_impact_factors(base_obs_loc, base_obs_type, state_loc, state_ki integer, intent(in) :: state_kind real(r8), intent(in) :: dist real(r8), intent(in) :: cutoff_rev -logical, intent(in) :: adjust_obs_impact -real(r8), intent(in) :: obs_impact_table(:, 0:) real(r8) :: impact_factor, cov_factor @@ -2614,12 +2577,11 @@ end subroutine get_my_obs_loc !> Get close obs from cache if appropriate. Cache new get_close_obs info !> if requested. -subroutine get_close_obs_cached(close_obs_caching, gc_obs, base_obs_loc, base_obs_type, & +subroutine get_close_obs_cached(gc_obs, base_obs_loc, base_obs_type, & my_obs_loc, my_obs_kind, my_obs_type, num_close_obs, close_obs_ind, close_obs_dist, & ens_handle, last_base_obs_loc, last_num_close_obs, last_close_obs_ind, & last_close_obs_dist, num_close_obs_cached, num_close_obs_calls_made) -logical, intent(in) :: close_obs_caching type(get_close_type), intent(in) :: gc_obs type(location_type), intent(inout) :: base_obs_loc, my_obs_loc(:) integer, intent(in) :: base_obs_type, my_obs_kind(:), my_obs_type(:) @@ -2662,12 +2624,11 @@ end subroutine get_close_obs_cached !> Get close state from cache if appropriate. Cache new get_close_state info !> if requested. -subroutine get_close_state_cached(close_obs_caching, gc_state, base_obs_loc, base_obs_type, & +subroutine get_close_state_cached(gc_state, base_obs_loc, base_obs_type, & my_state_loc, my_state_kind, my_state_indx, num_close_states, close_state_ind, close_state_dist, & ens_handle, last_base_states_loc, last_num_close_states, last_close_state_ind, & last_close_state_dist, num_close_states_cached, num_close_states_calls_made) -logical, intent(in) :: close_obs_caching type(get_close_type), intent(in) :: gc_state type(location_type), intent(inout) :: base_obs_loc, my_state_loc(:) integer, intent(in) :: base_obs_type, my_state_kind(:) From 09858bde74bd11f2f215dac18c66813a3beb4d30 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Fri, 19 Nov 2021 08:05:31 -0700 Subject: [PATCH 14/18] Eliminated the unused computation of correlation between obs and state ensembles. This tested bitwise across 'all' relevant namelist selections. Had no measurable impact on run-time in L96. --- .../modules/assimilation/assim_tools_mod.f90 | 30 +++++++++++-------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index c338de8cfc..78970f05cf 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -714,8 +714,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & updated_ens, my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - close_state_dist(j), cutoff_rev, net_a, & - grp_size, grp_beg, grp_end, i, my_state_indx(state_index), final_factor, correl) + close_state_dist(j), cutoff_rev, net_a, grp_size, grp_beg, grp_end, i, & + my_state_indx(state_index), final_factor, correl, local_varying_ss_inflate) ! If doing full assimilation, update the state variable ensemble with weighted increments if(.not. inflate_only) ens_handle%copies(1:ens_size, state_index) = updated_ens @@ -747,8 +747,8 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & updated_ens, my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - close_obs_dist(j), cutoff_rev, net_a, & - grp_size, grp_beg, grp_end, i, -1*my_obs_indx(obs_index), final_factor, correl) + close_obs_dist(j), cutoff_rev, net_a, grp_size, grp_beg, grp_end, i, & + -1*my_obs_indx(obs_index), final_factor, correl, .false.) obs_ens_handle%copies(1:ens_size, obs_index) = updated_ens endif @@ -2099,10 +2099,9 @@ end subroutine update_ens_from_weights !--------------------------------------------------------------- subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_kind, & - obs_prior, obs_inc, obs_prior_mean, obs_prior_var, obs_loc, obs_type, obs_time, & - dist, cutoff_rev, net_a, & - grp_size, grp_beg, grp_end, reg_factor_obs_index, reg_factor_ens_index, & - final_factor, correl) + obs_prior, obs_inc, obs_prior_mean, obs_prior_var, obs_loc, obs_type, obs_time, & + dist, cutoff_rev, net_a, grp_size, grp_beg, grp_end, reg_factor_obs_index, & + reg_factor_ens_index, final_factor, correl, correl_needed) integer, intent(in) :: ens_size integer, intent(in) :: num_groups @@ -2127,6 +2126,7 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ integer(i8), intent(in) :: reg_factor_ens_index real(r8), intent(out) :: final_factor real(r8), intent(out) :: correl(num_groups) +logical, intent(in) :: correl_needed real(r8) :: reg_coef(num_groups), increment(ens_size) real(r8) :: cov_factor, reg_factor @@ -2146,10 +2146,16 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ ! Loop through groups to update the state variable ensemble members do group = 1, num_groups grp_bot = grp_beg(group); grp_top = grp_end(group) - ! Do update of state, correl only needed for varying ss inflate but compute for all - call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & - obs_prior_var(group), obs_inc(grp_bot:grp_top), ens(grp_bot:grp_top), grp_size, & - increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) + ! Do update of state, correl only needed for varying ss inflate + if(correl_needed) then + call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & + obs_prior_var(group), obs_inc(grp_bot:grp_top), ens(grp_bot:grp_top), grp_size, & + increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group)) + else + call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), & + obs_prior_var(group), obs_inc(grp_bot:grp_top), ens(grp_bot:grp_top), grp_size, & + increment(grp_bot:grp_top), reg_coef(group), net_a(group)) + endif end do if(num_groups <= 1) then From 09fb3420415f2c83961bab3fb381f0ab4faed320 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Wed, 24 Nov 2021 11:35:12 -0700 Subject: [PATCH 15/18] Increased speed a bit by being smarter about checking for an impact before doing a subroutine call and computation. --- .../modules/assimilation/assim_tools_mod.f90 | 66 ++++++++----------- 1 file changed, 28 insertions(+), 38 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 78970f05cf..8376ac70b5 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -711,14 +711,16 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE endif - call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & - updated_ens, my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & - obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - close_state_dist(j), cutoff_rev, net_a, grp_size, grp_beg, grp_end, i, & - my_state_indx(state_index), final_factor, correl, local_varying_ss_inflate) + ! Compute the covariance localization and adjust_obs_impact factors (module storage) + final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_state_loc(state_index), & + my_state_kind(state_index), close_state_dist(j), cutoff_rev) - ! If doing full assimilation, update the state variable ensemble with weighted increments - if(.not. inflate_only) ens_handle%copies(1:ens_size, state_index) = updated_ens + if(final_factor > 0.0_r8) & + call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & + my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & + obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & + net_a, grp_size, grp_beg, grp_end, i, & + my_state_indx(state_index), final_factor, correl, local_varying_ss_inflate, inflate_only) ! Compute spatially-varying state space inflation if(local_varying_ss_inflate .and. final_factor > 0.0_r8) then @@ -744,13 +746,16 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! If forward observation operator failed, no need to update unassimilated observations if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE - call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & - updated_ens, my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & - obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - close_obs_dist(j), cutoff_rev, net_a, grp_size, grp_beg, grp_end, i, & - -1*my_obs_indx(obs_index), final_factor, correl, .false.) + ! Compute the covariance localization and adjust_obs_impact factors (module storage) + final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & + my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev) - obs_ens_handle%copies(1:ens_size, obs_index) = updated_ens + if(final_factor > 0.0_r8) & + call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & + my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & + obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & + net_a, grp_size, grp_beg, grp_end, i, & + -1*my_obs_indx(obs_index), final_factor, correl, .false., inflate_only) endif end do OBS_UPDATE endif @@ -2098,15 +2103,14 @@ subroutine update_ens_from_weights(ens, ens_size, rel_weight, ens_inc) end subroutine update_ens_from_weights !--------------------------------------------------------------- -subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_kind, & +subroutine obs_updates_ens(ens_size, num_groups, ens, ens_loc, ens_kind, & obs_prior, obs_inc, obs_prior_mean, obs_prior_var, obs_loc, obs_type, obs_time, & - dist, cutoff_rev, net_a, grp_size, grp_beg, grp_end, reg_factor_obs_index, & - reg_factor_ens_index, final_factor, correl, correl_needed) + net_a, grp_size, grp_beg, grp_end, reg_factor_obs_index, & + reg_factor_ens_index, final_factor, correl, correl_needed, inflate_only) integer, intent(in) :: ens_size integer, intent(in) :: num_groups -real(r8), intent(in) :: ens(ens_size) -real(r8), intent(out) :: updated_ens(ens_size) +real(r8), intent(inout) :: ens(ens_size) type(location_type), intent(in) :: ens_loc integer, intent(in) :: ens_kind real(r8), intent(in) :: obs_prior(ens_size) @@ -2116,33 +2120,21 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ type(location_type), intent(in) :: obs_loc integer, intent(in) :: obs_type type(time_type), intent(in) :: obs_time -real(r8), intent(in) :: dist -real(r8), intent(in) :: cutoff_rev real(r8), intent(in) :: net_a(num_groups) integer, intent(in) :: grp_size integer, intent(in) :: grp_beg(num_groups) integer, intent(in) :: grp_end(num_groups) integer, intent(in) :: reg_factor_obs_index integer(i8), intent(in) :: reg_factor_ens_index -real(r8), intent(out) :: final_factor +real(r8), intent(inout) :: final_factor real(r8), intent(out) :: correl(num_groups) logical, intent(in) :: correl_needed +logical, intent(in) :: inflate_only real(r8) :: reg_coef(num_groups), increment(ens_size) -real(r8) :: cov_factor, reg_factor +real(r8) :: reg_factor integer :: group, grp_bot, grp_top -! Compute the covariance localization and adjust_obs_impact factors (module storage) -cov_factor = cov_and_impact_factors(obs_loc, obs_type, ens_loc, & - ens_kind, dist, cutoff_rev) - -! If no impact, don't do anything else -if(cov_factor <= 0.0_r8) then - final_factor = cov_factor - updated_ens = ens - return -endif - ! Loop through groups to update the state variable ensemble members do group = 1, num_groups grp_bot = grp_beg(group); grp_top = grp_end(group) @@ -2158,16 +2150,14 @@ subroutine obs_updates_ens(ens_size, num_groups, ens, updated_ens, ens_loc, ens_ endif end do -if(num_groups <= 1) then - final_factor = cov_factor -else +if(num_groups > 1) then reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, & reg_factor_obs_index, reg_factor_ens_index) - final_factor = min(cov_factor, reg_factor) + final_factor = min(final_factor, reg_factor) endif ! Get the updated ensemble -updated_ens = ens + final_factor * increment +if(.not. inflate_only) ens = ens + final_factor * increment end subroutine obs_updates_ens From 2635cabf81fc2f410eb9f0982683738f49718fbc Mon Sep 17 00:00:00 2001 From: jlaucar <54950434+jlaucar@users.noreply.github.com> Date: Tue, 7 Dec 2021 14:41:24 -0700 Subject: [PATCH 16/18] Apply suggestions from code review Co-authored-by: nancy collins Co-authored-by: hkershaw-brown <20047007+hkershaw-brown@users.noreply.github.com> --- .../modules/assimilation/adaptive_inflate_mod.f90 | 4 ++-- assimilation_code/modules/assimilation/assim_tools_mod.f90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 index d7d81cf6ed..b362b7b9f5 100644 --- a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 +++ b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 @@ -651,7 +651,7 @@ subroutine update_inflation(inflate_handle, inflate, inflate_sd, prior_mean, pri end subroutine update_inflation !------------------------------------------------------------------------------- -!> Computes updatea inflation mean and inflation sd for single state space inflation +!> Computes updated inflation mean and inflation sd for single state space inflation subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd, & ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, & @@ -703,7 +703,7 @@ subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd end subroutine update_single_state_space_inflation !------------------------------------------------------------------------------- -!> Computes updatea inflation mean and inflation sd for varying state space inflation +!> Computes updated inflation mean and inflation sd for varying state space inflation subroutine update_varying_state_space_inflation(inflate, inflate_mean, inflate_sd, & ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, & diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 8376ac70b5..2b57123f78 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -747,7 +747,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE ! Compute the covariance localization and adjust_obs_impact factors (module storage) - final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & + final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev) if(final_factor > 0.0_r8) & @@ -2178,7 +2178,7 @@ function cov_and_impact_factors(base_obs_loc, base_obs_type, state_loc, state_ki real(r8) :: impact_factor, cov_factor -! Get external impact factors, cycle if impact of this ob on this state is zerio +! Get external impact factors, cycle if impact of this ob on this state is zero if (adjust_obs_impact) then ! Get the impact factor from the table if requested impact_factor = obs_impact_table(base_obs_type, state_kind) From 757add213103e6cae7cd3ff4aa87695269bad303 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Tue, 7 Dec 2021 15:46:04 -0700 Subject: [PATCH 17/18] Modified some if statements to exit loops or return from functions rather that conditionally executing large blocks of code. Cleaned up erroneous comment. Reformatted indentation on several lines. --- .../assimilation/adaptive_inflate_mod.f90 | 104 +++++++++--------- .../modules/assimilation/assim_tools_mod.f90 | 32 +++--- 2 files changed, 68 insertions(+), 68 deletions(-) diff --git a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 index b362b7b9f5..151f0cd91e 100644 --- a/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 +++ b/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90 @@ -671,35 +671,35 @@ subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd real(r8) :: gamma, ens_var_deflate, r_var, r_mean ! If either inflation or sd is not positive, not really doing inflation -if(inflate_mean > 0.0_r8 .and. inflate_sd > 0.0_r8) then - ! For case with single spatial inflation, use gamma = 1.0_r8 - gamma = 1.0_r8 - ! Deflate the inflated variance; required for efficient single pass - ! This is one of many places that assumes linear state/obs relation - ! over range of ensemble; Essentially, we are removing the inflation - ! which has already been applied in filter to see what inflation should - ! have been needed. - ens_var_deflate = orig_obs_prior_var / & - (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 - - ! If this is inflate_only (i.e. posterior) remove impact of this obs. - ! This is simulating independent observation by removing its impact. - if(inflate_only .and. & - ens_var_deflate > small .and. & - obs_err_var > small .and. & - obs_err_var - ens_var_deflate > small ) then - r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) - r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var) - else - r_var = ens_var_deflate - r_mean = orig_obs_prior_mean - endif - - ! Update the inflation mean value and standard deviation - call update_inflation(inflate, inflate_mean, inflate_sd, & - r_mean, r_var, ens_size, obs, obs_err_var, gamma) +if(inflate_mean <= 0.0_r8 .or. inflate_sd <= 0.0_r8) return + +! For case with single spatial inflation, use gamma = 1.0_r8 +gamma = 1.0_r8 +! Deflate the inflated variance; required for efficient single pass +! This is one of many places that assumes linear state/obs relation +! over range of ensemble; Essentially, we are removing the inflation +! which has already been applied in filter to see what inflation should +! have been needed. +ens_var_deflate = orig_obs_prior_var / & + (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 + +! If this is inflate_only (i.e. posterior) remove impact of this obs. +! This is simulating independent observation by removing its impact. +if(inflate_only .and. & + ens_var_deflate > small .and. & + obs_err_var > small .and. & + obs_err_var - ens_var_deflate > small ) then + r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) + r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var) +else + r_var = ens_var_deflate + r_mean = orig_obs_prior_mean endif +! Update the inflation mean value and standard deviation +call update_inflation(inflate, inflate_mean, inflate_sd, & + r_mean, r_var, ens_size, obs, obs_err_var, gamma) + end subroutine update_single_state_space_inflation !------------------------------------------------------------------------------- @@ -726,36 +726,36 @@ subroutine update_varying_state_space_inflation(inflate, inflate_mean, inflate_s real(r8) :: diff_sd, outlier_ratio logical :: do_adapt_inf_update -if(inflate_mean > 0.0_r8 .and. inflate_sd > 0.0_r8) then - ! Gamma is less than 1 for varying ss, see adaptive inflate module - gamma = reg_factor * abs(correl) +if(inflate_mean <= 0.0_r8 .or. inflate_sd <= 0.0_r8) return - ! Remove the impact of inflation to allow efficient single pass with assim. - if ( abs(gamma) > small ) then - ens_var_deflate = orig_obs_prior_var / & - (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 - else - ens_var_deflate = orig_obs_prior_var - endif +! Gamma is less than 1 for varying ss, see adaptive inflate module +gamma = reg_factor * abs(correl) - ! If this is inflate only (i.e. posterior) remove impact of this obs. - if(inflate_only .and. & - ens_var_deflate > small .and. & - obs_err_var > small .and. & - obs_err_var - ens_var_deflate > small ) then - r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) - r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var) - else - r_var = ens_var_deflate - r_mean = orig_obs_prior_mean - endif +! Remove the impact of inflation to allow efficient single pass with assim. +if ( abs(gamma) > small ) then + ens_var_deflate = orig_obs_prior_var / & + (1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2 +else + ens_var_deflate = orig_obs_prior_var +endif - ! IS A TABLE LOOKUP POSSIBLE TO ACCELERATE THIS? - ! Update the inflation values - call update_inflation(inflate, inflate_mean, inflate_sd, & - r_mean, r_var, ens_size, obs, obs_err_var, gamma) +! If this is inflate only (i.e. posterior) remove impact of this obs. +if(inflate_only .and. & + ens_var_deflate > small .and. & + obs_err_var > small .and. & + obs_err_var - ens_var_deflate > small ) then + r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var) + r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var) +else + r_var = ens_var_deflate + r_mean = orig_obs_prior_mean endif +! IS A TABLE LOOKUP POSSIBLE TO ACCELERATE THIS? +! Update the inflation values +call update_inflation(inflate, inflate_mean, inflate_sd, & + r_mean, r_var, ens_size, obs, obs_err_var, gamma) + end subroutine update_varying_state_space_inflation diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 2b57123f78..36c3fe012f 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -674,7 +674,6 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Adaptive localization needs number of other observations within localization radius. ! Do get_close_obs first, even though state space increments are computed before obs increments. - ! JLA: ens_handle doesn't ever appear to be used. Get rid of it. Should be obs_ens_handle anyway? call get_close_obs_cached(gc_obs, base_obs_loc, base_obs_type, & my_obs_loc, my_obs_kind, my_obs_type, num_close_obs, close_obs_ind, close_obs_dist, & ens_handle, last_base_obs_loc, last_num_close_obs, last_close_obs_ind, & @@ -715,15 +714,16 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_state_loc(state_index), & my_state_kind(state_index), close_state_dist(j), cutoff_rev) - if(final_factor > 0.0_r8) & - call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & - my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & - obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - net_a, grp_size, grp_beg, grp_end, i, & - my_state_indx(state_index), final_factor, correl, local_varying_ss_inflate, inflate_only) + if(final_factor <= 0.0_r8) cycle STATE_UPDATE + + call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), & + my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, & + obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & + net_a, grp_size, grp_beg, grp_end, i, & + my_state_indx(state_index), final_factor, correl, local_varying_ss_inflate, inflate_only) ! Compute spatially-varying state space inflation - if(local_varying_ss_inflate .and. final_factor > 0.0_r8) then + if(local_varying_ss_inflate) then do group = 1, num_groups call update_varying_state_space_inflation(inflate, & ens_handle%copies(ENS_INF_COPY, state_index), & @@ -750,12 +750,13 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), & my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev) - if(final_factor > 0.0_r8) & - call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & - my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & - obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & - net_a, grp_size, grp_beg, grp_end, i, & - -1*my_obs_indx(obs_index), final_factor, correl, .false., inflate_only) + if(final_factor <= 0.0_r8) cycle OBS_UPDATE + + call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), & + my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, & + obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, & + net_a, grp_size, grp_beg, grp_end, i, & + -1*my_obs_indx(obs_index), final_factor, correl, .false., inflate_only) endif end do OBS_UPDATE endif @@ -778,8 +779,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Assure user we have done something if (print_trace_details >= 0) then -write(msgstring, '(A,I8,A)') & - 'Processed', obs_ens_handle%num_vars, ' total observations' + write(msgstring, '(A,I8,A)') 'Processed', obs_ens_handle%num_vars, ' total observations' call error_handler(E_MSG,'filter_assim:',msgstring) endif From 5ec034523eea8ab4248966bcaf48f84cbd4b1abb Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Tue, 7 Dec 2021 17:58:25 -0700 Subject: [PATCH 18/18] bump version and CHANGELOG note, I bumped the wrong number on the last release --- CHANGELOG.rst | 7 ++++++- conf.py | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 4e02cd5c59..03da62bd9d 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,7 +22,12 @@ individual files. The changes are now listed with the most recent at the top. -**November 22 2021 :: Bug fix for groups with posterior spatially-varying adaptive inflation. Tag: v9.12.13** +**December 7 2021 :: Refactored filter_assim. Tag: v9.12.0** + +- Filter_assim refactored so each process calcuates increments +- Code readability changes + +**November 22 2021 :: Bug fix for groups with posterior spatially-varying adaptive inflation. Tag: v9.11.13** - Removed the additional outlier threshold test for each group when using posterior spatially-varying adaptive inflation. The outlier test is done for the entire ensemble diff --git a/conf.py b/conf.py index 97111ea469..8294d0277e 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '9.11.13' +release = '9.12.0' master_doc = 'README' # -- General configuration ---------------------------------------------------