[Dart-dev] DART/branches Revision: 12901

[dart at ucar.edu]

thoar at ucar.edu
2018-10-10 15:05:32 -0600 (Wed, 10 Oct 2018)
834
This checks each gridcell for outliers and replaces them with a noisy estimate.

The methodology is that all the ensemble members for that gridcell and at a particular time
have their variance calculated. The largest value is replaced by the median and the variance
is recalculated. If the ratio of the variances exceeds some threshhold (in this case a factor of 100.0)
the largest value is considered to be an outlier and is replaced by the value at the 75th quantile 
that has a random_gaussian() noise from a distribution with the interquartile range as the stddev.
The process is iterated until there are no outliers in the gridcell. Preliminary tests for all of 
2008 with the downwelling longwave radiation indicates that _at_most_ there were 3 outliers 
in a gridcell. This is loosely based on the Modified Thompson Tau test.




Modified: DART/branches/cesm_clm/models/clm/clean_forcing.f90
===================================================================
--- DART/branches/cesm_clm/models/clm/clean_forcing.f90	2018-10-10 20:45:02 UTC (rev 12900)
+++ DART/branches/cesm_clm/models/clm/clean_forcing.f90	2018-10-10 21:05:32 UTC (rev 12901)
@@ -82,7 +82,7 @@
 character(len=*), parameter :: routine  = 'clean_forcing'
 
 integer, parameter :: ensemble_size = 80
-character(len=*), parameter :: directory = '/glade/p_old/image/thoar/CAM_DATM/4xdaily'
+character(len=*), parameter :: directory = '/glade/p/cisl/dares/thoar/CAM_DATM/4xdaily'
 
 character(len=16) :: varname
 character(len=16) :: variables(10) = (/'a2x6h_Faxa_lwdn ', &
@@ -128,7 +128,7 @@
 
 namelist /clean_forcing_nml/ year, criterion
 
-!======================================================================
+!===============================================================================
 
 ! Read the namelist entry
 call find_namelist_in_file('input.nml', 'clean_forcing_nml', iunit)
@@ -143,15 +143,15 @@
 
 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')
+100 format(A,'/CAM_DATM.cpl_',i4.4,'.ha2x1dx6h.',i4.4,'.nc')
 
 do imember = 1,ensemble_size
-   write(input_file(imember),100) imember, year
-
+   write(input_file(imember),100) directory, imember, year
    ncid(imember) = nc_open_file_readwrite(input_file(imember),routine)
+   write(*,*)'Opened '//trim(input_file(imember))
 enddo
 
-varname = variables(5)
+varname = variables(1)
 
 call nc_get_variable_num_dimensions(ncid(1), varname, numdims)
 call nc_get_variable_size(ncid(1), varname, dimlens)
@@ -165,7 +165,7 @@
 
 TIMESTEP : do itime = 1,nT
 
-   call fill_tensor(itime)
+   call read_tensor(itime)
 
    num_outliers = 0_i8
 
@@ -238,11 +238,10 @@
 
          num_outliers = num_outliers + 1_i8
 
-         mean   = sum(sorted(10:70))/61.0_r8
-         stddev = sqrt(sum(sorted(10:70) - mean)**2)/60.0_r8
+!        mean   = sum(sorted(10:70))/61.0_r8
+!        stddev = sqrt(sum(sorted(10:70) - mean)**2)/60.0_r8
+!        write(*,*)'ix,iy,stddev,iqr = ',ix,iy,stddev,iqr
 
-         write(*,*)'ix,iy,stddev,iqr = ',ix,iy,stddev,iqr
-
 !        noise = random_gaussian(r, 0.0_r8, 3.0_r8*stddev)
          noise = random_gaussian(r, 0.0_r8, iqr)
 
@@ -260,12 +259,13 @@
 
    write(*,*)'timestep ',itime,' had ',num_outliers, &
              ' values replaced when criterion is ',criterion
+   call write_tensor(itime)
 
 enddo TIMESTEP
 
-
 do imember = 1,ensemble_size
    call nc_close_file(ncid(imember),routine)
+   write(*,*)'Closed '//trim(input_file(imember))
 enddo
 
 call finalize_utilities(routine)
@@ -274,10 +274,10 @@
 contains
 !-------------------------------------------------------------------------------
 
-subroutine fill_tensor(itime)
+subroutine read_tensor(itime)
 integer, intent(in) :: itime
 
-character(len=*), parameter :: routine = 'fill_tensor'
+character(len=*), parameter :: routine = 'read_tensor'
 integer :: imember, io
 
 MEMBER : do imember = 1,ensemble_size
@@ -292,15 +292,32 @@
                       start=ncstart, count=nccount)
    call nc_check(io, routine, 'get_var '//trim(varname))
 
-!     minvalue = minval(tensor(imember,:,:))
-!     maxvalue = maxval(tensor(imember,:,:))
-!