Improve DA RANDOMCV capability

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"      ""  "" 

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.

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")
 

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

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, &

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")