diff --git a/Registry/registry.var b/Registry/registry.var index 9ecfe1fee0..44b4b38937 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -472,6 +472,7 @@ rconfig logical hybrid_dual_res namelist,wrfvar16 1 .false. - "hy rconfig integer dual_res_upscale_opt namelist,wrfvar16 1 3 - "dual_res_upscale_opt" "" "" rconfig logical use_4denvar namelist,wrfvar16 1 .false. - "4D-En-Var" "switch for activating 4D-Ensemble-Var" "" rconfig character analysis_type namelist,wrfvar17 1 "3D-VAR" - "analysis_type" "" "" +rconfig integer n_randomcv namelist,wrfvar17 1 1 - "n_randomcv" "number of realizations of random_cv" "" rconfig integer sensitivity_option namelist,wrfvar17 1 -1 - "sensitivity_option" "" "" rconfig logical adj_sens namelist,wrfvar17 1 .false. - "adj_sens" "" "" rconfig character analysis_date namelist,wrfvar18 1 "2002-08-03_00:00:00.0000" - "analysis_date" "" "" diff --git a/var/README.namelist b/var/README.namelist index 015b345ee6..07c395870c 100644 --- a/var/README.namelist +++ b/var/README.namelist @@ -510,6 +510,8 @@ Description of WRFDA namelist variables, defined in Registry/registry.var ; "VERIFY": verification mode. ; WRFDA resets check_max_iv=.false. and ntmax=0 ; "RANDOMCV": for creating ensemble perturbations + n_randomcv = 1 ; number of realizations of randomcv + ; set to > 1 to get multiple perturbed output in one execution adj_sens = .false. ; .true.: write out gradient of Jo for adjoint sensitivity / &wrfvar18 ; needs to set &wrfvar21 and &wrfvar22 as well if ob_format=1 and/or radiances are used. diff --git a/var/da/da_define_structures/da_initialize_cv.inc b/var/da/da_define_structures/da_initialize_cv.inc index 770d587e97..93f8494d3a 100644 --- a/var/da/da_define_structures/da_initialize_cv.inc +++ b/var/da/da_define_structures/da_initialize_cv.inc @@ -18,37 +18,7 @@ subroutine da_initialize_cv(cv_size, cv) ! [1.0] Initialize cv: !--------------------------------------------------------------------------- - if (anal_type_randomcv) then - - ! [2.1] Initialize random number generator and scalars: - - call da_random_seed - - if( use_rf )then - - ! [2.2] Calculate random numbers with Gaussian distribution: - - do i = 1, cv_size - call da_gauss_noise(z) - cv(i) = z - end do - else - write(unit=message(1),fmt='(a)')'Need to inject CV into wavelet space' - call da_error(__FILE__,__LINE__,message(1:1)) - endif - - mean_cv = sum(cv) / real(cv_size) - rms_cv = sqrt(sum(cv*cv) / real(cv_size)) - std_dev_cv = sqrt(rms_cv * rms_cv - mean_cv * mean_cv) - - write(unit=message(1),fmt='(a)')' Gaussian (Normal) noise statistics:' - write(unit=message(2),fmt='(a,f15.5)')' Mean = ',mean_cv - write(unit=message(3),fmt='(a,f15.5)')' RMS = ', rms_cv - write(unit=message(4),fmt='(a,f15.5)')' STD DEV = ', std_dev_cv - call da_message(message(1:4)) - else - cv = 0.0 - end if + cv = 0.0 if (trace_use) call da_trace_exit("da_initialize_cv") diff --git a/var/da/da_main/da_solve.inc b/var/da/da_main/da_solve.inc index 22a9655531..58cc5d6161 100644 --- a/var/da/da_main/da_solve.inc +++ b/var/da/da_main/da_solve.inc @@ -194,9 +194,9 @@ endif if (anal_type_randomcv) then - ntmax = 0 - write(unit=stdout,fmt='(a)') & - ' Resetting ntmax = 0 for analysis_type = randomcv' + write_gts_omb_oma = .false. + write_rej_obs_conv = .false. + write_unpert_obs = .false. end if if ( cv_options == 3 ) then @@ -596,6 +596,35 @@ ! allocate (full_eignvec(cv_size)) ! end if + !------------------------------------------------------ + ! set CV to random noise ("RANDOMCV") + !------------------------------------------------------ + if (anal_type_randomcv) then + it = 1 + ! Initialize random number generator and scalars + call da_random_seed + do i= 1, n_randomcv + write(ci,'(i3.3)') i + write(unit=message(1),fmt='(a,a)') & + ' Setting randomcv for wrfvar_output_randomcv.e', trim(ci) + call da_message(message(1:1)) + call da_set_randomcv (cv_size, cvt) + call da_transform_vtox (grid,cv_size,xbx,be,grid%ep,cvt,grid%vv,grid%vp) + call da_transform_xtoxa (grid) + call da_transfer_xatoanalysis (it, xbx, grid, config_flags) + call da_update_firstguess(grid,'wrfvar_output_randomcv.e'//trim(ci)) + if ( i < n_randomcv ) then + ! restore the original grid info for the next realization of randomcv + call da_med_initialdata_input( grid , config_flags, 'fg') + call da_setup_firstguess(xbx, grid, config_flags, .false.) + end if + end do + ! Done with randomcv. + ! Set the following to skip some code to get to the deallocation part. + max_ext_its = 0 + outer_loop_restart = .false. + end if !anal_type_randomcv + if ( outer_loop_restart ) then !call da_get_unit(cvt_unit) cvt_unit=600 @@ -847,11 +876,6 @@ end if !------------------------------------------------------------------------ - ! reset cv to random noise - if (anal_type_randomcv) then - call da_set_randomcv (cv_size, xhat) - end if - ! [8.5] Update latest analysis solution: if (.not. var4d) then diff --git a/var/da/da_main/da_wrfvar_top.f90 b/var/da/da_main/da_wrfvar_top.f90 index c794911105..39ddee57d4 100644 --- a/var/da/da_main/da_wrfvar_top.f90 +++ b/var/da/da_main/da_wrfvar_top.f90 @@ -47,7 +47,7 @@ module da_wrfvar_top use da_define_structures, only : y_type, j_type, iv_type, be_type, & xbx_type,da_deallocate_background_errors,da_initialize_cv, & da_zero_vp_type,da_allocate_y,da_deallocate_observations, & - da_deallocate_y, da_zero_x + da_deallocate_y, da_zero_x, da_random_seed use da_minimisation, only : da_get_innov_vector,da_minimise_cg, & da_minimise_lz, da_write_diagnostics, da_calculate_residual, & da_calculate_grady, da_sensitivity, da_lanczos_io, da_calculate_j, & diff --git a/var/da/da_tools/da_set_randomcv.inc b/var/da/da_tools/da_set_randomcv.inc index 65829574df..a3581c29ca 100644 --- a/var/da/da_tools/da_set_randomcv.inc +++ b/var/da/da_tools/da_set_randomcv.inc @@ -19,7 +19,8 @@ subroutine da_set_randomcv(cv_size, rcv) ! [1] Initialize random number generator and scalars: - call da_random_seed + !moved to da_solve.inc + !call da_random_seed sum_cv = 0.0 sum_cv2 = 0.0 @@ -39,8 +40,8 @@ subroutine da_set_randomcv(cv_size, rcv) write(unit=stdout,fmt=*) write(unit=stdout,fmt='(a)')' Gaussian (Normal) noise statistics:' - write(unit=stdout,fmt='(a,f15.5)')' Mean = ',mean_cv - write(unit=stdout,fmt='(a,f15.5)')' RMS = ', rms_cv + write(unit=stdout,fmt='(a,f15.5)')' Mean = ', mean_cv + write(unit=stdout,fmt='(a,f15.5)')' RMS = ', rms_cv write(unit=stdout,fmt='(a,f15.5)')' STD DEV = ', std_dev_cv if (trace_use) call da_trace_exit("da_set_randomcv")