[Dart-dev] [8248] DART/trunk/models/bgrid_solo/model_mod.f90: add an explicit perturb routine which perturbs the T

[nancy at ucar.edu]

Revision: 8248
Author:   nancy
Date:     2015-07-24 13:49:03 -0600 (Fri, 24 Jul 2015)
Log Message:
-----------
add an explicit perturb routine which perturbs the T
field only by adding gaussian noise.  by default it is
NOT enabled; comment out the early return code in
pert_model_state() to use it.

Modified Paths:
--------------
    DART/trunk/models/bgrid_solo/model_mod.f90

-------------- next part --------------
Modified: DART/trunk/models/bgrid_solo/model_mod.f90
===================================================================
--- DART/trunk/models/bgrid_solo/model_mod.f90	2015-07-24 19:41:58 UTC (rev 8247)
+++ DART/trunk/models/bgrid_solo/model_mod.f90	2015-07-24 19:49:03 UTC (rev 8248)
@@ -74,6 +74,8 @@
 use          obs_kind_mod, only: KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT, &
                                  KIND_SURFACE_PRESSURE, KIND_TEMPERATURE
 
+use mpi_utilities_mod, only : my_task_id
+
 !-----------------------------------------------------------------------
 
 implicit none
@@ -188,7 +190,8 @@
 !-----------------------------------------------------------------------
 
 ! Stuff to allow addition of random 'sub-grid scale' noise
-!!!!!!type(random_seq_type) :: randseed
+! and also used in the pert_model_state() call.
+type(random_seq_type) :: randtype, randnoise
 logical :: first_call = .true.
 
 
@@ -275,7 +278,7 @@
 ! Need to initialize random sequence on first call
 if(first_call) then
    !write(*, *) 'NOISE_SD is ', noise_sd
-   !!!!!!call init_random_seq(randseed)
+   !!!!!!call init_random_seq(randnoise)
    first_call = .false.
 endif
 
@@ -286,7 +289,7 @@
 !      do k = 1, size(Var_dt%t, 3)
 !         temp_t = Var_dt%t(i, j, k)
          !!!!!!Var_dt%t(i, j, k) = &
-            !!!!!!(1.0_r8 + random_gaussian(randseed, 0.0_r8, dble(noise_sd))) * Var_dt%t(i, j, k)
+            !!!!!!(1.0_r8 + random_gaussian(randnoise, 0.0_r8, dble(noise_sd))) * Var_dt%t(i, j, k)
 !      end do
 !   end do
 !end do
@@ -1254,7 +1257,6 @@
    val(1, 2) =  get_val(x, lon_below, lat_above, nint(level), itype)
    val(2, 1) =  get_val(x, lon_above, lat_below, nint(level), itype)
    val(2, 2) =  get_val(x, lon_above, lat_above, nint(level), itype)
-
 else
 ! Case of pressure specified in vertical
    val(1, 1) =  get_val_pressure(x, lon_below, lat_below, pressure, itype, istatus)
@@ -1264,7 +1266,6 @@
 endif
 
 ! Do the weighted average for interpolation
-!write(*, *) 'fracts ', lon_fract, lat_fract
 do i = 1, 2
    a(i) = lon_fract * val(2, i) + (1.0 - lon_fract) * val(1, i)
 end do
@@ -1836,7 +1837,7 @@
 !-------------------------------------------------------------------------------
 call check(nf90_sync(ncFileID),"atts sync")
 
-write (*,*)'nc_write_model_atts: netCDF file ',ncFileID,' is synched ...'
+!write (*,*)'nc_write_model_atts: netCDF file ',ncFileID,' is synched ...'
 
 contains
 
@@ -2054,17 +2055,43 @@
 real(r8), intent(out) :: pert_state(:)
 logical,  intent(out) :: interf_provided
 
+logical, save :: first_call = .true.
+integer :: i, j, k
+
 if ( .not. module_initialized ) call static_init_model
 
+! option 1: let filter do the perturbs
+! (comment this out to select other options below)
 interf_provided = .false.
+return
+ 
+! (debug) option 2: tell filter we are going to perturb, but don't.
+! interf_provided = .true.
+! pert_state = state
+! return
 
-! you *cannot* set this to junk. in at least one
-! location the caller is passing the same
-! array into the first arg as the second.  doing this
-! corrupts the state vector completely. set it to the
-! incoming state if you must do something.
-pert_state = state ! Just to satisfy INTENT(OUT)
+! option 3: really perturb, t only
 
+interf_provided = .true.
+
+if (first_call) then
+   call init_random_seq(randtype, my_task_id()+1)
+   first_call = .false.
+endif
+
+call vector_to_prog_var(state, get_model_size(), global_Var)
+
+do k = 1, size(global_Var%t, 3)
+   do j = 1, size(global_Var%t, 2)
+      do i = 1, size(global_Var%t, 1)
+         global_Var%t(i, j, k) = global_Var%t(i, j, k) + &
+             random_gaussian(randtype, 0.0_r8, 0.01_r8)
+      end do
+   end do
+end do
+
+call prog_var_to_vector(global_Var, pert_state, get_model_size())
+
 end subroutine pert_model_state