[Dart-dev] DART/branches Revision: 12870

[dart at ucar.edu]

thoar at ucar.edu
2018-10-01 16:21:46 -0600 (Mon, 01 Oct 2018)
339
Better detection method for determining outliers.
checks the ratio of the variance of the original with the 
variance of the same with the max value replaced by the median.
If the ratio is > 100.0, the max is an outlier.
Some of the fields have _tiny_ values ... replacing the outlier
with a useful value is still an area of exploration.




Modified: DART/branches/cesm_clm/models/clm/clean_forcing.f90
===================================================================
--- DART/branches/cesm_clm/models/clm/clean_forcing.f90	2018-10-01 22:17:03 UTC (rev 12869)
+++ DART/branches/cesm_clm/models/clm/clean_forcing.f90	2018-10-01 22:21:46 UTC (rev 12870)
@@ -11,7 +11,7 @@
 !             CAM_DATM.cpl_0055.ha2x1dx6h.2008.nc
 !             CAM_DATM.cpl_0055.ha2x1dx6h.2009.nc
 !             CAM_DATM.cpl_0055.ha2x1dx6h.2010.nc
-!              a2x6h_Faxa_rainc     rainc
+!               a2x6h_Faxa_rainc     rainc
 !               a2x6h_Faxa_rainl     rainl
 !               a2x6h_Faxa_snowc     snowc
 !               a2x6h_Faxa_snowl     snowl
@@ -44,6 +44,10 @@
 
 use             sort_mod, only : index_sort
 
+use       random_seq_mod, only : random_seq_type, &
+                                 init_random_seq, &
+                                 random_gaussian
+
 use        utilities_mod, only : initialize_utilities, &
                                  finalize_utilities, &
                                  find_namelist_in_file, &
@@ -82,10 +86,10 @@
 
 character(len=16) :: varname
 character(len=16) :: variables(10) = (/'a2x6h_Faxa_lwdn ', &
+                                       'a2x6h_Faxa_swndf', &
+                                       'a2x6h_Faxa_swvdf', &
                                        'a2x6h_Faxa_swndr', &
                                        'a2x6h_Faxa_swvdr', &
-                                       'a2x6h_Faxa_swndf', &
-                                       'a2x6h_Faxa_swvdf', &
                                        'null            ', &
                                        'a2x6h_Faxa_rainc', &
                                        'a2x6h_Faxa_rainl', &
@@ -92,14 +96,17 @@
                                        'a2x6h_Faxa_snowc', &
                                        'a2x6h_Faxa_snowl'/)
 
-integer :: indices(ensemble_size)
-integer ::    ncid(ensemble_size)
+integer  :: indices(ensemble_size)
+integer  ::    ncid(ensemble_size)
+real(r8) ::  sorted(ensemble_size)
 character(len=256) :: input_file(ensemble_size)
 
 real(r8), allocatable :: tensor(:,:,:) ! nx,ny,ensemble_size
 real(r8) :: minvalue, maxvalue, q1, q2, q3, iqr
-real(r8) :: newhigh, newlow, original_variance, new_variance
-logical :: suspicious, newmin, newmax
+real(r8) :: mean, stddev, noise
+real(r8) :: newhigh, newlow, variance, new_variance
+logical  :: suspicious, newmin, newmax
+type(random_seq_type) :: r
 
 integer :: imember, itime, iy, ix, nT, ny, nx
 integer :: ncstart(3)
@@ -110,7 +117,7 @@
 integer            :: iunit, io
 character(len=512) :: string1, string2, string3
 
-integer(i8) :: icount
+integer(i8) :: num_outliers
 
 !-------------------------------------------------------------------------------
 ! namelist
@@ -117,9 +124,9 @@
 !-------------------------------------------------------------------------------
 
 integer :: year = 2008
-real(r8) :: iqr_multiplier = 20.0_r8
+real(r8) :: criterion = 100.0_r8
 
-namelist /clean_forcing_nml/ year, iqr_multiplier
+namelist /clean_forcing_nml/ year, criterion
 
 !======================================================================
 
@@ -134,6 +141,8 @@
 
 call initialize_utilities(progname=routine)
 
+call init_random_seq(r,1)
+
 100 format('/glade/p_old/image/thoar/CAM_DATM/4xdaily/CAM_DATM.cpl_',i4.4,'.ha2x1dx6h.',i4.4,'.nc')
 
 do imember = 1,ensemble_size
@@ -142,12 +151,12 @@
    ncid(imember) = nc_open_file_readwrite(input_file(imember),routine)
 enddo
 
-varname = variables(1)
+varname = variables(5)
 
 call nc_get_variable_num_dimensions(ncid(1), varname, numdims)
 call nc_get_variable_size(ncid(1), varname, dimlens)
 
-if ( .false. ) write(*,*)'dimlens are ',dimlens(1:numdims)
+if ( .true. ) write(*,*)'dimlens are ',dimlens(1:numdims)