Merge branch 'master' into sphinx-documentation

keeping sphinx-branch upto date with master

assimilation_code/modules/utilities/utilities_mod.f90

@@ -67,15 +67,20 @@ module utilities_mod
           find_first_occurrence, &
           array_dump, &
           dump_unit_attributes, &
-          ! lowest level routines
+          ! lowest level routines follow.
+          ! these need to be cautious about logging, error handlers, etc
           initialize_utilities, &
           finalize_utilities, &
           register_module, &
           find_namelist_in_file, &
           check_namelist_read, &
+          ! this should follow string_to_logical
+          get_value_from_string, &
           do_nml_file, &
           do_nml_term, &
           log_it, &
+          ! these two routines should move up to after get_value_from_string
+          ! they shouldn't be grouped with these other low level routines.
           interactive_r, &
           interactive_i
 
@@ -480,6 +485,120 @@ subroutine check_namelist_read(iunit, iostat_in, nml_name)
 
 end subroutine check_namelist_read
 
+!-----------------------------------------------------------------------
+! TODO: the next 2 routines belong right after the string_to_logical function.
+!> convert integers or strings describing options into integer. 
+!> 
+!> Some input items have historically been integers which do not help describe the
+!> option they are selecting.  They are being converted to descriptive strings.
+!> During a transition period before the integers are deprecated this routine will return 
+!> an integer value if the input string is either an integer or a string.  
+!> For strings, the corresponding integer values are input as a 1-to-1 array with the
+!> valid strings.   On error this routine prints an error message and does not return.
+!> If needed, this routine could be changed to take an optional return code, which if 
+!> present would return to the calling code without printing or calling the error 
+!> handler so the caller can take an alternative code path.
+
+function get_value_from_string(input,integer_options,string_options,context)
+
+character(len=*), intent(in) :: input                 ! value from namelist, eg
+integer,          intent(in) :: integer_options(:)    ! possible integer values
+character(len=*), intent(in) :: string_options(:)     ! matching string values
+character(len=*), intent(in), optional :: context 
+integer                      :: get_value_from_string
+
+character(len=len_trim(input)) :: uppercase
+character(len=len(string_options)) :: possibility
+integer :: ios, iopt, candidate
+
+if ( .not. module_initialized ) call initialize_utilities
+
+get_value_from_string = MISSING_I
+
+! Try to read the input as an ASCII-coded integer ( i.e. '3') 
+
+read(input,*,iostat=ios) candidate
+if (ios == 0) then
+
+   get_value_from_string = candidate
+
+   ! Check for valid value
+   if (any(integer_options == get_value_from_string)) then
+      return
+   else
+      call get_value_error(input, integer_options, string_options, context)
+   endif
+endif
+
+! numeric conversion not possible, try to interpret as a string
+
+if (size(integer_options) /= size(string_options)) then
+   write(msgstring2,*)'size(integer_options) =',size(integer_options), &
+                  ' /= size( string_options) =',size( string_options)
+   call error_handler(E_ERR,'get_value_from_string',msgstring2,source,text2=context)
+endif
+
+uppercase = trim(input)
+call to_upper(uppercase)
+
+LOOP : do iopt = 1,size(integer_options)
+
+   possibility = string_options(iopt)
+   call to_upper(possibility)
+
+   if (uppercase == possibility) then
+      get_value_from_string = integer_options(iopt)
+      exit LOOP
+   endif
+
+enddo LOOP
+
+if (get_value_from_string == MISSING_I) &
+   call get_value_error(input, integer_options, string_options, context)
+
+end function get_value_from_string
+
+
+!-----------------------------------------------------------------------
+!> report any errors from the get_value_from_string routine
+
+subroutine get_value_error(input, integer_options, string_options, whofrom)
+
+character(len=*),           intent(in) :: input
+integer,                    intent(in) :: integer_options(:)
+character(len=*),           intent(in) :: string_options(:)
+character(len=*), optional, intent(in) :: whofrom
+
+character(len=256) :: string1
+integer :: iopt
+
+if (present(whofrom)) then
+   msgstring1 = trim(whofrom)//' no valid option found.' 
+else
+   msgstring1 = 'No valid option found.' 
+endif
+write(msgstring2,*)'input is "'//trim(input)//'"'
+write(msgstring3,*)'valid values are the following integers or matching character strings:'
+
+call error_handler(E_MSG, 'get_value_from_string', msgstring1, source, &
+           text2=msgstring2, text3=msgstring3)
+
+do iopt = 1,size(integer_options)
+   write(string1,*) integer_options(iopt), ' = "'//trim(string_options(iopt))//'"' 
+   call error_handler(E_MSG,'get_value_from_string',string1)
+enddo
+
+! Repeat the leading message (as an E_ERR) to help delineate the problem.
+call error_handler(E_ERR, 'get_value_from_string', msgstring1, source)
+
+end subroutine get_value_error
+
+!----------------------------------------------------------------------
+! TODO: next pull request, put interactive_r and interactive_i after
+! these routines, just before array_1d_dump.
+! they don't belong in the debug section at the end of this file.
+!----------------------------------------------------------------------
+
 !-----------------------------------------------------------------------
 !> if trying to write an unformatted string, like "write(*,*)"
 !> to both standard output and the logfile, call this routine instead.
@@ -2015,6 +2134,8 @@ function string_to_logical(inputstring, string_to_match)
 
 end function string_to_logical
 
+
+
 !-----------------------------------------------------------------------
 !> dump the contents of a 1d array with a max of N items per line.
 !> optional arguments allow the caller to restrict the output to 
@@ -2789,6 +2910,9 @@ subroutine output_unit_attribs(ios, label, cvalue, ivalue, lvalue)
 
 end subroutine output_unit_attribs
 
+! TODO: these next 2 routines should move to after the get_value_from_string
+! routine.  they aren't debug routines which is what routines in this section are.
+
 !----------------------------------------------------------------------
 !> prompt for a real value, optionally setting min and/or max limits
 !> loops until valid value input.
@@ -2877,7 +3001,6 @@ function interactive_i(str1,minvalue,maxvalue)
 
 end function interactive_i
 
-
 !=======================================================================
 ! End of utilities_mod
 !=======================================================================

assimilation_code/programs/obs_utils/obs_info.f90

@@ -66,9 +66,9 @@ program obs_info
    type(time_type) :: last_time
 end type
 
-! in spite of the name, this is the number of specific types.
+! an array to hold counts for each obs type.
 ! also one for all obs types.
-type(obs_info_type) :: oinfo(0:max_defined_types_of_obs)
+type(obs_info_type) :: oinfo(max_defined_types_of_obs)
 type(obs_info_type) :: all_obs
 
 type(location_type) :: location
@@ -140,7 +140,7 @@ program obs_info
 do fnum = 1, num_input_files
 
    ! initialize the bookkeeping structures
-   do i=0, max_defined_types_of_obs
+   do i=1, max_defined_types_of_obs
       call initialize(oinfo(i))
    enddo
    call initialize(all_obs)
@@ -213,7 +213,7 @@ program obs_info
    endif
 
    ! print out the results
-   ALLTYPES: do i=0, max_defined_types_of_obs
+   ALLTYPES: do i=1, max_defined_types_of_obs
       if (oinfo(i)%count == 0) cycle ALLTYPES
       if (counts_only) then
          call compute_times(oinfo(i)%first_time, oinfo(i)%last_time, avg_string=mid_string)
@@ -321,7 +321,7 @@ subroutine print_obs_seq(seq_in, filename)
 
 
 ! Initialize input obs_types
-do i = 0, max_defined_types_of_obs
+do i = 1, max_defined_types_of_obs
    type_count(i) = 0
 enddo
 identity_count = 0
@@ -395,7 +395,7 @@ subroutine print_obs_seq(seq_in, filename)
 call error_handler(E_MSG, '', msgstring)
 write(msgstring, *) '---------------------------------------------------------'
 call error_handler(E_MSG, '', msgstring)
-do i = 0, max_defined_types_of_obs
+do i = 1, max_defined_types_of_obs
    if (type_count(i) > 0) then 
       write(msgstring, '(a32,i8,a)') trim(get_name_for_type_of_obs(i)), &
                                      type_count(i), ' obs'

guide/DART_LAB/matlab/oned_model_inf.m

@@ -1386,10 +1386,15 @@ function step_ahead(~, ~)
                     handles.inflation = str2double(get(handles.ui_edit_fixed_inflation,'String'));
 
                 case 'Adaptive Inflation'
+                    % Here, only stick to one algorithm (Gaussian; flavor 2). Users may manually 
+                    % switch between 'Gaussian' [Anderson 2009] and 'Gamma'
+                    % [El Gharamti 2018]. The GUI option is only available
+                    % in the Lorenz'96 section. 
                     [lambda, handles.adap_inf_Std] = ...
                         update_inflate(mean(ens), var(ens), observation, obs_error_sd^2, inf_prior, ...
-                        handles.inflation, handles.adap_inf_Min, handles.adap_inf_Max, ...
-                        1, handles.adap_inf_Std, handles.adap_inf_Std_Min);
+                        handles.inflation, handles.adap_inf_Std, handles.adap_inf_Min, handles.adap_inf_Max, ...
+                        1, handles.adap_inf_Std_Min, handles.ens_size, 'Gaussian');
+
                     % Damping is placed unusually here to obtain a less messy code
                     % It won't matter because it's a single variable case!
                     handles.inflation = 1.0 + handles.adap_inf_Damp * ( lambda - 1.0 );

guide/DART_LAB/matlab/private/change_GA_IG.m

@@ -0,0 +1,37 @@
+function beta = change_GA_IG(mode, var)
+
+% Routine to change the Gaussian prior into an inverse gamma (IG).
+% The Gaussian prior is represented by a mode (:= mean) and a variance; var
+
+
+%% DART software - Copyright UCAR. This open source software is provided
+% by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+
+% Computation savers 
+var_p    = zeros(1, 3); 
+var_p(1) = var;
+for i = 2:3
+   var_p(i) = var_p(i-1)*var;
+end
+
+mode_p    = zeros(1, 9);
+mode_p(1) = mode;
+for i = 2:9
+  mode_p(i) = mode_p(i-1)*mode;
+end
+
+% Calculate the rate parameter for IG distribution.
+% It's a function of both the prior mean and variannce, 
+% obtained as a "real" solution to a cubic polynomial.
+AA = mode_p(4) * sqrt((var_p(2) + 47*var*mode_p(2) + 3*mode_p(4)) / var_p(3));
+BB = 75*var_p(2)*mode_p(5);
+CC = 21*var*mode_p(7);
+DD = var_p(3)*mode_p(3);
+EE = (CC + BB + DD + mode_p(9) + 6*sqrt(3)*AA*var_p(3)) / var_p(3);
+
+beta = (7*var*mode + mode_p(3))/(3*var)                               + ...
+       EE^(1/3)/3 + mode_p(2)*(var_p(2) + 14*var*mode_p(2) + ...
+       mode_p(4)) / (3*var_p(2)*EE^(1/3));
+

guide/DART_LAB/matlab/private/enh_compute_new_density.m

@@ -0,0 +1,22 @@
+function enh_density = enh_compute_new_density(dist_2, sigma_p_2, sigma_o_2, alpha, beta, gamma_corr, lambda, ens_size)
+
+%% DART software - Copyright UCAR. This open source software is provided
+% by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+
+% Compute probability of this lambda being correct
+exp_prior = - beta / lambda;
+
+% Compute probability that observation would have been observed given this lambda
+fac1 = (1 + gamma_corr * (sqrt(lambda) - 1.0))^2;
+fac2 = -1 / ens_size;
+if fac1 < abs(fac2), fac2 = 0; end
+
+theta    = sqrt( (fac1+fac2) * sigma_p_2 + sigma_o_2 );
+exp_like = - 0.5 * dist_2 / theta^2;
+
+% Compute the updated probability density for lambda
+enh_density = beta^alpha / gamma(alpha) * lambda^(- alpha - 1) / ... 
+              (sqrt(2.0 * pi) * theta) * exp(exp_like + exp_prior);
+