Merge pull request #363 from NCAR/mitgcm-perturb

Fix for  MITgcm_ocean pert_model_copies

CHANGELOG.rst

@@ -22,6 +22,16 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
+**June 24 2022 :: Bug-fixes for MITgcm_ocean and Var obs converter. Tag: v10.0.2**
+
+- MITgcm_ocean pert_model_copies routine fixed to use the correct variable clamping
+  value and indices for each element of the copies array. 
+- Var obs converter quicklbuild.sh fixed to correctly locate the required 
+  3DVAR_OBSPROC code.
+- Documentation for Var obs converter updated with information for where to 
+  get the latest WRF 3DVAR_OBSPROC code.
+
+
 **June 2 2022 :: Bug-fixes for ps_rand_local in the Bgrid Model. Tag: v10.0.1**
 
 - performs the missing call for initialize_utilities() 

conf.py

@@ -21,7 +21,7 @@
 author = 'Data Assimilation Research Section'
 
 # The full version, including alpha/beta/rc tags
-release = '10.0.1'
+release = '10.0.2'
 master_doc = 'README'
 
 # -- General configuration ---------------------------------------------------

models/MITgcm_ocean/model_mod.f90

@@ -97,8 +97,6 @@ module model_mod
 character(len=512) :: string1, string2, string3
 logical, save      :: module_initialized = .false.
 
-! Storage for a random sequence for perturbing a single initial state
-type(random_seq_type) :: random_seq
 
 !------------------------------------------------------------------
 !
@@ -1304,9 +1302,11 @@ subroutine pert_model_copies(state_ens_handle, ens_size, pert_amp, interf_provid
 real(r8),            intent(in)    :: pert_amp
 logical,             intent(out)   :: interf_provided
 
-integer  :: i, ivar
-integer  :: copy, start_ind, end_ind
+integer  :: ivar
+integer  :: copy
+integer  :: i
 real(r8) :: pertval, clamp_min_val 
+integer  :: iloc, jloc, kloc ! not used, but required for get_model_variable_indices
 
 type(random_seq_type) :: random_seq
 
@@ -1315,33 +1315,28 @@ subroutine pert_model_copies(state_ens_handle, ens_size, pert_amp, interf_provid
 interf_provided = .true.
 
 call init_random_seq(random_seq, my_task_id())
+
+INDICES : do i = 1, state_ens_handle%my_num_vars
 
-VARIABLES : do ivar = 1, get_num_variables(domain_id)
-
-   start_ind = get_index_start(domain_id, ivar)
-   end_ind   = get_index_end(  domain_id, ivar)
-
+   call get_model_variable_indices(state_ens_handle%my_vars(i), iloc, jloc, kloc, ivar)
    clamp_min_val = get_io_clamping_minval(domain_id, ivar)
 
-   INDICES : do i = start_ind, end_ind
-      MEMBERS : do copy = 1, ens_size
-
-         ! Only perturb the actual ocean cells;
-         ! Leave the land and ocean floor values alone.
-         if( state_ens_handle%copies(copy, i) == FVAL ) cycle MEMBERS
+   MEMBERS : do copy = 1, ens_size
 
-         pertval = random_gaussian(random_seq, state_ens_handle%copies(copy, i), &
-                                         model_perturbation_amplitude)
+      ! Only perturb the actual ocean cells;
+      ! Leave the land and ocean floor values alone.
+      if( state_ens_handle%copies(copy, i) == FVAL ) cycle MEMBERS
 
-         ! Clamping: Samples obtained from truncated Gaussian dist. 
-         if (.not. log_transform) pertval = max(clamp_min_val, pertval)
+      pertval = random_gaussian(random_seq, state_ens_handle%copies(copy, i), &
+                                      model_perturbation_amplitude)
 
-         state_ens_handle%copies(copy, i) = pertval
+      ! Clamping: Samples obtained from truncated Gaussian dist. 
+      if (.not. log_transform) pertval = max(clamp_min_val, pertval)
 
-      enddo MEMBERS
-   enddo INDICES
+      state_ens_handle%copies(copy, i) = pertval
 
-enddo VARIABLES
+   enddo MEMBERS
+enddo INDICES
 
 end subroutine pert_model_copies
 
@@ -1361,12 +1356,12 @@ function read_model_time(filename)
 
 read_model_time = model_time
 
-!if (do_output() .and. debug > 0 .and. present(last_time)) then
+if (do_output()) then
    call print_time(read_model_time, str='MITgcm_ocean time is ',iunit=logfileunit)
    call print_time(read_model_time, str='MITgcm_ocean time is ')
    call print_date(read_model_time, str='MITgcm_ocean date is ',iunit=logfileunit)
    call print_date(read_model_time, str='MITgcm_ocean date is ')
-!endif
+endif
 
 end function read_model_time